Gaussian Quadrature

and formula generation with Sympy

A key point of Simpson's rule is that it is exact for cubic polynomials. While this is not at all hard to show by hand, this verification is one example of somthing that can be done with Sympy. Furthermore, there are related formulae that can be derived by solving systems exactly and Sympy provides a convenient way to do this exactly.

We start by importing the whole SymPy name space and declaring a variable, as we've done before:

In [1]:
from sympy import *
x = symbols('x')

Checking Simpson's rule

Again, one way to think of Simpson's rule is that is exact for cubic polynomials. Let's check this for the general cubic $f(x)=ax^3+bx^2+cx+d$ over the interval $[-1,1]$, i.e. $$\int_{-1}^1 f(x) \, dx = \frac{1}{3}f(-1) + \frac{4}{3}f(0) + \frac{1}{3}f(1).$$

In [2]:
a,b,c,d = symbols('a,b,c,d')
def f(x): return a*x**3 + b*x**2 + c*x + d
print(integrate(f(x),(x,-1,1)))
print(expand(f(-1)/3+4*f(0)/3+f(1)/3))

2*b/3 + 2*d
2*b/3 + 2*d

Now suppose we turn this around and ask, suppose we want $$\int_{-1}^1 f(x) \, dx \approx w_1\,f(-1) + w_2\,f(0) + w_3\,f(1).$$ What choice of weights $w_1$, $w_2$, and $w_3$ work best? Well, if it's exact for polynomials of high degree, that seems good and, the higher the better. Thus, might ask "what values of $w_1$, $w_2$, and $w_3$ make $$\int_{-1}^1 f(x) \, dx = w_1\,f(-1) + w_2\,f(0) + w_3\,f(1)$$ for $f(x) = ax^3+bx^2+cx+d$? Well, let's compute the right hand side:

In [3]:
w1,w2,w3 = symbols('w1,w2,w3')
collect(expand(w1*f(-1)+w2*f(0)+w3*f(1)), (a,b,c,d))
Out[3]:
a*(-w1 + w3) + b*(w1 + w3) + c*(-w1 + w3) + d*(w1 + w2 + w3)

This suggests that we solve the following:

In [4]:
solve([Eq(w3-w1,0), Eq(w1+w3,2/S(3)), Eq(w3-w1,0), Eq(w1+w2+w3,2)])
Out[4]:
[{w2: 4/3, w3: 1/3, w1: 1/3}]

As expected!

Gaussian quadrature

Gauss came up with the following idea: rather than choosing just the weights $w_1$, $w_2$, and $w_3$, let's also choose the nodes $x_1$, $x_2$, and $x_3$ at which the function is evaluated. That leads us to the following equations:

In [5]:
x1,x2,x3 = symbols('x1,x2,x3')
w1,w2,w3 = symbols('w1,w2,w3')
eqs = [Eq(integrate(x**k,(x,-1,1)),w1*x1**k + w2*x2**k + w3*x3**k) for k in range(6)]
eqs
Out[5]:
[Eq(2, w1 + w2 + w3),
 Eq(0, w1*x1 + w2*x2 + w3*x3),
 Eq(2/3, w1*x1**2 + w2*x2**2 + w3*x3**2),
 Eq(0, w1*x1**3 + w2*x2**3 + w3*x3**3),
 Eq(2/5, w1*x1**4 + w2*x2**4 + w3*x3**4),
 Eq(0, w1*x1**5 + w2*x2**5 + w3*x3**5)]

OK, let's solve 'em.

In [6]:
sols = solve(eqs)
sols
Out[6]:
[{w1: 5/9, w2: 5/9, w3: 8/9, x1: -sqrt(15)/5, x2: sqrt(15)/5, x3: 0},
 {w1: 5/9, w2: 5/9, w3: 8/9, x1: sqrt(15)/5, x2: -sqrt(15)/5, x3: 0},
 {w1: 5/9, w2: 8/9, w3: 5/9, x1: -sqrt(15)/5, x2: 0, x3: sqrt(15)/5},
 {w1: 5/9, w2: 8/9, w3: 5/9, x1: sqrt(15)/5, x2: 0, x3: -sqrt(15)/5},
 {w1: 8/9, w2: 5/9, w3: 5/9, x1: 0, x2: -sqrt(15)/5, x3: sqrt(15)/5},
 {w1: 8/9, w2: 5/9, w3: 5/9, x1: 0, x2: sqrt(15)/5, x3: -sqrt(15)/5}]

Looks like a lot but, in reality, it's just the same solution with the $x_i$s permuted. Let's grab the one in the implied order.

In [7]:
the_sol = [sol for sol in sols if sol[x1]<sol[x2]<sol[x3]]
the_sol
Out[7]:
[{w1: 5/9, w2: 8/9, w3: 5/9, x1: -sqrt(15)/5, x2: 0, x3: sqrt(15)/5}]

Checking Gauss's rule

OK, let's use it.

In [8]:
gauss_nodes = [the_sol[0][x] for x in (x1,x2,x3)]
gauss_weights = [the_sol[0][w] for w in (w1,w2,w3)]
def gauss_rule(f):
    x1,x2,x3 = gauss_nodes
    w1,w2,w3 = gauss_weights
    return w1*f(x1)+w2*f(x2)+w3*f(x3)

First, let's see if it works exactly for 5th degree polynomials. Here's a fifth degree polynomial. The degree is just high enough that we'll represent the coefficient with a function $c_k$.

In [9]:
c = symbols('c', cls=Function)
def f(x): return sum([c(k)*x**k for k in range(6)])
f(x)
Out[9]:
x**5*c(5) + x**4*c(4) + x**3*c(3) + x**2*c(2) + x*c(1) + c(0)

Let's compare the gauss_rule and the actual integral.

In [10]:
print(gauss_rule(f))
print(integrate(f(x), (x,-1,1)))

2*c(0) + 2*c(2)/3 + 2*c(4)/5
2*c(0) + 2*c(2)/3 + 2*c(4)/5

Looks good, let's do a numerical example.

In [11]:
def f(x): return exp(-x**2)
N(gauss_rule(f))
Out[11]:
1.49867959566003

I wonder if that's right? Let's compare to Numpy.

In [12]:
import numpy as np
from scipy.integrate import quad, simps
def f(x): return np.exp(-x**2)
quad(f, -1,1)
Out[12]:
(1.493648265624854, 1.6582826951881447e-14)

Pretty close for a three point evaluation! I wonder how it compares to Simpson's rule?

In [13]:
simps([f(-1), f(0), f(1)])
Out[13]:
1.5785862941142952