How to compute Gaussian Quadrature weights?

[ektrules]

My numerical analysis book doesn't explain it. It just tells you to use precomputed tables, and directs you to an out of print book from the 80's that I can't find anywhere.

After searching, I found http://en.wikipedia.org/wiki/Gaussian_quadrature#Computation_of_Gaussian_quadrature_rules" in the "Gaussian Quadrature" Wikipedia article. But I don't really understand it, and don't have access to the referenced book either. Specifically, in the "Jacobi" matrix, what are An and Bn supposed to be?

 

[lpetrich]

Here's what orthogonal polynomials are:

A set of polynomials p_i(x) with degree i having values 0, 1, 2, 3, ... that satisfy

<p_i(x)*p_j(x)> = integral over x of w(x)*p_i(x)*p_j(x) = 0 for j != i

and some weight function w(x).

They can be found from a recurrence relation; this one assumes monic polynomials:

p_i+1(x) + (B_i - x)*p_i(x) + A_i*p_i-1 = 0

with p₀(x) = 1 and p₁(x) = x - B₀

That's where the A's and B's come from. One can calculate them using orthogonality:

B_i = <x*(p_i(x))²>/<(p_i(x))²>
A_i = <x*p_i(x)*p_i-1(x)>/<(p_i-1(x))²>

To derive them, multiply the recurrence relation by p_i(x) and p_i-1(x) and integrate.

Once you get to there, the Wikipedia article explains what next.


Some orthogonal polynomials satisfy second-order differential equations. There are several families of them, and some mathematicians have found formulas for Gaussian-quadrature weights for them.

Jacobi polynomials - (1-x)^a*(1+x)^b over (-1,1)

Laguerre polynomials - exp(-x) over (0,infinity)

Associated Laguerre polynomials - x^a*exp(x) over (0,infinity)

Hermite polynomials - exp(-x²) over (-infinity,infinity)