MATH/APPM 4650
Class notes, March 12

Gaussian quadrature

(Review): The basic idea of Newton-Cotes formulas is to use some combination of function values to get integrals of low-degree polynomials exactly.

Recall for example to get Simpson's rule, we set h=(b-a)/2 and started with a linear combination
∫ab f(x) dx ∼ p f(a) + q f(a+h) + r f(a+2h),
for some unknown coefficients p, q, and r.

We then demanded that the left and right sides be exactly equal when f(x) = 1, f(x) = x-a, or f(x) = (x-a)2. (Three conditions for three unknowns.)

This gave us three equations: and we solved them to find p=h/3, q=4h/3, and r=h/3.

Thus we obtained Simpson's rule
∫ab f(x) dx ∼ h/3 [f(a) + 4 f(a+h) + f(a+2h)]


(New):
Now we get "meta," which is how all really interesting math arises.

Did we really have to use points a, a+h, and a+2h?

What happens if we let those points be unknowns as well?
Then we'd have six unknowns, so we could demand six conditions on them.
The logical thing to do is to demand that the rule be exact on f(x) = 1, f(x) = x, ..., f(x) = x5.
It's not totally obvious that this will work (after all, now the equations will be nonlinear, so there may be no solution).
But let's try and see what happens.

To make things simple, let's assume that a=-1 and b=1. (Otherwise it's a mess.)
So we have three unknown points u, v, and w, all lying in [-1,1], and three coefficients p, q, r.

∫-11 f(x) dx ∼ p f(u) + q f(v) + r f(w)

OMG WHAT DO WE DO NOW????

All right, calm down.

Remember, the equations are linear in p, q, and r.
Take the odd equations (the ones with zeroes on the left side) and write them as a matrix.

/ u v w \ / p \ / 0 \
| u3 v3 w3 | | q | = | 0 |
\ u5 v5 w5 / \ r / \ 0 /


That means the determinant of this matrix must be zero.

So we get the equation
uvw(v-u)(u-w)(v-w)(u+v)(u+w)(v+w) = 0.

This implies that either one of the points is zero, or that two of the points are negatives of each other.

In fact both conditions turn out to be required. (We'll skip the details.)
So let's say w=-u and v=0. Then the six equations become Since u cannot be zero, the second, fourth, and sixth equations all give p=r.
We now have two equations for p and u:
So u=-sqrt(3/5) and p=5/9.
That means w=sqrt(3/5) and r=5/9 and q=8/9.

So the Gaussian 3-point quadrature rule is
∫-11 f(x) ∼ 1/9 [ 5 f(-sqrt(3/5)) + 8 f(0) + 5 f(sqrt(3/5)) ].


If we wanted to use this on the interval x ∈ [a, b], we'd have to change variables to rescale.
u-substitution: u = 2x/(b-a) - (b+a)/(b-a).

We expect the error to be roughly the integral of x6 (since there is no error from any power below that), which is roughly O(h7) for a one-interval application.
This means the error is O(h6) for a composite version.

So for about the same amount of work as Simpson's method, we should have a much better error estimate.

Obviously nobody actually figures these values out by hand.
They have all been tabulated, and there are techniques involving orthogonal polynomials which would allow you to get more if you wanted.


Gaussian quadrature is a little safer with high-order formulas, because you're not really using interpolating polynomials and Taylor series.
Instead you're using orthogonal polynomials, which act a lot more like Fourier series. It's generally easier to understand how these work.


(Not in book):
The problem with Gaussian quadrature is that using a two-point Gaussian formula doesn't help you at all if you decide to use a three-point formula.
The nodes are totally different in every order.

So unlike the trapezoid rule, where it's not much additional work to double the number of steps (because half the calculations have already been done), it's a lot of trouble to increase the number of steps in a Gaussian method.

But it's also important to be able to do this, because how else would you practically estimate the error?

This is a lot like the problem for interpolating polynomials.

Recall from February 18, when we studied the optimal placement of points to interpolate.
Evenly-spaced points did badly near the ends, while Chebyshev points did much better throughout.
But Chebyshev points are totally different for every step, so you can't really add in an extra point to see how you're doing.


"Gauss-Kronrod" quadrature is a technique due to Kronrod, which takes the Gaussian points and adds in extra points to get a higher-order estimate.

Suppose we had already figured out the evaluation points (-sqrt(3/5), 0, sqrt(3/5)) and the weights (5/9, 8/9, 5/9).
Let's add four more points: -y, -z, z, and y, and figure out what new weights s and t we'd need.
New rule:
∫-11 f(x) dx ∼ 1/9 [ 5 f(-sqrt(3/5)) + 8 f(0) + 5 f(sqrt(3/5)) ]
+ s f(-y) + s f(y) + t f(-z) + t f(z).


(Since we know how this is going to work, choosing the coefficients and points symmetrically will obviously make all the odd integrals zero, so we'll just have to worry about the even integrals.)
So now with four more unknowns to work with, we'll want the integral to be correct on x6, x8, x10, and x12.
Because it was already correct on x0, x2, and x4.
Right?




Wrong!
Adding these new points and coefficients probably breaks the formula on the old polynomials.
So we'll have to correct for that too.

We can therefore only impose four constraints total, which means we can only make it correct on x0, x2, x4, and x6. (We'll get x7 for free, of course.)

This hardly seems worth it.

Instead, compromise and change the weights on the old function values (but keep those old function values).

∫-11 f(x) dx ∼ p f(-sqrt(3/5)) + q f(0) + p f(sqrt(3/5)) ]
+ s f(-y) + s f(y) + t f(-z) + t f(z).

Now we have six unknowns to determine, so we can impose six conditions.

By symmetry, all odd powers will be correct for free, so we just need to worry about even powers.
We can get x0, x2, ..., x10 correct, which would make the first error term the integral of x12 (which gives an error of h13 on the one-interval method, or h12 on the composite method).


We can then use the Kronrod formula to check the accuracy of the Gaussian formula. Hopefully the error is small!

Of course, these equations are kind of a pain to solve, which is why people just use tables of pre-calculated values.

General rules: