Legendre Polynomials and Gaussian Quadrature.

The Vector Space L^2(-1,1).

We are going to find this set of orthogonal polynomials with respect to the L^2 inner product on [-1,1]. We define this inner product in Mathematica:

IP[f_,g_] := Integrate[f*g,{x,-1,1}];

We try a few simple examples:

IP[x^2,x^3 + 4x + 2]

IP[x,1]

Gram-Schmidt Orthogonalization

We now need to define Gram-Schmidt orthogonalization given a set of polynomials. The big idea is that given a list of orthogonal polynomials 
q0, ...., qn and some p of degree n+1 we should transform p to 

             p - (<p,q0>/<q0,q0>) q0 - (<p,q1>/<q1,q1>) q1 - ...

GramSchmidtCoefficient[p_,q_] :=IP[p,q]/IP[q,q]; 
GramSchmidt[p_,Q_] := p - Sum[GramSchmidtCoefficient[p,Q[[i]]]*Q[[i]],{i,1,Length[Q]}];

We now try this out:

GramSchmidt[x,{1}]

GramSchmidt[x^2,{1,x}]

Putting everything together.

We now recursively apply GramSchmidt.

LegendrePolylist = {1};
Do[LegendrePolylist = Append[LegendrePolylist,GramSchmidt[x^n,LegendrePolylist]],{n,1,15}];
  
Simplify[Collect[LegendrePolylist,x]]

This gives us the list of Legendre polynomials.

Generating a Quadrature Rule from the Legendre Polynomials

In order to generate a (numerical) quadrature rule from the collection of Legendre polynomials, we need to find their roots. Luckily, Mathematica is pretty good at this.

RL = Table[Solve[LegendrePolylist[[i]] == 0, x],{i,2,15}];

As you can see, these are not rational numbers, but they are numbers that we can evaluate numerically relatively easily with Mathematica. To do this, we'll need to actually get these guys out of rule form and expressed just as numbers. Mathematica makes this easy by intelligently applying /. across the nested list which is RL.

RootList = x /. RL;

Evaluating numerically, we get

GaussQuadratureNodes = N[RootList,20];

Weights for Gaussian Quadrature.

We now have to figure out the appropriate (exact and numerical) weights for Gaussian quadrature. By construction, in the Gauss quadrature theorem, we know that 

A_i = Integral_a^b l_i(x). 

We could be more clever about the algebra if we were doing this by hand, but Mathematica is happy to construct and integrate the l_i for us. Suppose that x is the list of roots. Then we have

L_i = (∏_(j=1)^n)((x - (Subscript[x, j]))/((Subscript[x, i])-(Subscript[x, j]))),

where j=i is supposed to be missing from the sum, but I don't know how to write that!

L[i_,X_] :=
If[i>1,Product[(x - X[[j]])/(X[[i]] - X[[j]]),{j,1,i-1}],1]*
If[i<Length[X],Product[(x - X[[j]])/(X[[i]] - X[[j]]),{j,i+1,Length[X]}],1];

We can now try this out:

RootList[[3]]

L[1,RootList[[3]]]

L[2,RootList[[3]]]

L[3,RootList[[3]]]

We can now build a table of all the L_i with a couple of nested Table commands (you don't want to look at it):

LTable = Table[
Table[L[i,RootList[[j]]],{i,1,Length[RootList[[j]]]}],
{j,2,6}];

We now want to integrate the L_i to get the weights A_i for Gauss Quadrature. For instance,

Integrate[L[3,RootList[[3]]],{x,-1,1}]

In general, we should make a table of integrals, of course:

ATable = Integrate[LTable,{x,-1,1}];

Remember that the RootList table starts with the polynomial x, with single root 0, while the LTable and ATable start one degree higher. To keep our notation straight, we just prepend a fake weight to the table.

ATable = Prepend[ATable,0];

Testing the actual Gaussian Quadrature Rule (exactly)

We now have tables of nodes (RootList and GaussQuadratureNodes) and weights (ATable) for the Gauss quadrature rules with between 2 and 5 nodes. The n node formula should correctly integrate all polynomials up to degree 2n -1. 

Of course, it will be no fun whatsoever unless we actually test this hypothesis on some real examples. First, we code the actual Gauss quadrature rule:

GaussQuadrature[f_,n_] := Sum[ATable[[n,i]]*(f /. {x ->RootList[[n,i]]}),{i,1,Length[RootList[[n]]]}];

Now we try it out!

GaussQuadrature[x+4,3]

Simplify[GaussQuadrature[x+4,3]]

Integrate[x+4,{x,-1,1}]

So the three node GQ rule should integrate polynomials up to degree 5. Let's try it!

P = x^5 - 4*x^2 + 17*x^3 - 2*x + 1;
Simplify[GaussQuadrature[P,3]]
Integrate[P,{x,-1,1}]

Trying out Gauss Quadrature on a non-polynomial integrand

We now return to our favorite example: integrating trig functions. Let's try Cos[x].

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];

GaussQuadrature[Cos[x],3]

Ok, we'll have to N it:

N[GaussQuadrature[Cos[x],5],20]

N[Integrate[Cos[x],{x,-1,1}],20]

Of course, in order to make a serious test of these methods, we ought to count digits and make a graph, as we usually do:

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];
GQDigits = Table[Digits[N[GaussQuadrature[Cos[x],i],20],N[Integrate[Cos[x],{x,-1,1}],20]],{i,2,5}]
ListLinePlot[GQDigits,Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

We now compare these digits to the original NewtonCotes rules:

FitData[n_,a_, h_,f_] := Table[{a + i*h,f[a + i*h]},{i,0,n}];
Poly[n_,a_,h_,f_] := Collect[Simplify[InterpolatingPolynomial[FitData[n,a,h,f],x]],x];
NewtonCotes[n_,a_,h_,f_] := Integrate[Poly[n,a,h,f],{x,a,a+n*h}];
CompositeNewtonCotes[k_,n_,a_,h_,f_] := Sum[NewtonCotes[n,a+i*n*h,h,f],{i,0,k-1}];
NCDigits = Table[Digits[N[NewtonCotes[i,-1,2/i,Cos],20],N[Integrate[Cos[x],{x,-1,1}],20]],{i,2,5}]

ListLinePlot[{GQDigits,NCDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

We can see that the Gauss Quadrature rules (top) are substantially better than the equally spaced Newton-Cotes rules (bottom).

Integrands for which polynomial interpolation is particularly poorly behaved -- does Gauss Quadrature do better?

Recall that when we were first studying polynomial interpolation, we saw that some functions were unusually difficult to model with an interpolating polynomial based on equally spaced gridpoints. How do they perform under Gauss Quadrature?

ff[x_] := 1/(1 + (10*x)^2);

a = -1;
b =    1;
n  =   40;

EqualNodes := Table[{a + (b-a)i/n,ff[a+(b-a)i/n]},{i,0,n}]
EqP := InterpolatingPolynomial[EqualNodes,x];

Here's the equally-spaced nodes interpolating polynomial:

Show[Plot[{ff[x],EqP},{x,a,b}],ListPlot[EqualNodes,PlotStyle->{PointSize[Large]}],PlotRange->{{-1,1},{-2,2}}]

If we reduce the number of nodes to six, we get a better looking polynomial graph, but you can still see there's substantial error.

ff[x_] := 1/(1 + (10*x)^2);

a = -1;
b = 1;
n  =  6;

EqualNodes := Table[{a + (b-a)i/n,ff[a+(b-a)i/n]},{i,0,n}]
EqP := InterpolatingPolynomial[EqualNodes,x];

Show[Plot[{ff[x],EqP},{x,-1.1,1.1},Filling->{2}],ListPlot[EqualNodes,PlotStyle->{PointSize[Large]}]]

We need to integrate this over [-1,1] using our various rules, so let's reestablish the Newton-Cotes formulas:

It's immediately clear that the regular NC rules are a disaster:

N[Table[NewtonCotes[i,-1,2/i,ff],{i,2,13}]]

N[Integrate[ff[x],{x,-1,1}],20]

The composite rules fare better (with a lot more function evaluations):

N[Table[CompositeNewtonCotes[60/i,i,-1,2/60,ff],{i,2,5}],20]

We now compare the GaussQuadrature performance:

GQData = N[Table[GaussQuadrature[ff[x],i],{i,2,5}],20]

This is still a disaster! Clearly, we're going to have to come up with a new idea in order to handle this integrand correctly.

Numerical Gauss Quadrature Rules.

It's not hard to see that even Mathematica is going to have a hard time computing the weights for Gauss Quadrature exactly after a certain point (I tried to get the first 15 Gauss Quadrature rules exactly, which takes more than 1/2 hour on my desktop machine). 

But we can easily compute floating point approximations to the weights:

NLTable = Table[
Table[L[i,GaussQuadratureNodes[[j]]],{i,1,Length[GaussQuadratureNodes[[j]]]}],
{j,2,14}];

NATable = Integrate[NLTable,{x,-1,1}];

NATable = Prepend[NATable,0];

Comparing the numerical weight table and the ordinary weight table, we see:

NATable[[5]]
N[ATable,18][[5]]

GaussQuadratureNodes[[5]]
N[RootList[[5]],20]

NGaussQuadrature[f_,n_] := Sum[NATable[[n,i]]*(f /. {x ->GaussQuadratureNodes[[n,i]]}),{i,1,Length[GaussQuadratureNodes[[n]]]}];

N[Integrate[Cos[x],{x,-1,1}],20]

NGaussQuadrature[Cos[x],4]

N[GaussQuadrature[Cos[x],4],20]

NGQDigits = Table[Digits[NGaussQuadrature[Cos[x],i],N[Integrate[Cos[x],{x,-1,1}],40]],{i,2,8}]

ListLinePlot[{GQDigits,NCDigits,NGQDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

We see from this graph the expected result that Gauss Quadrature with numerically computed weights and nodes performs quite well until we reach 8 nodes. At that point, there seems to be a precipitous loss of accuracy!

Extra Credit: Explain this phenomenon. What's going wrong in our code?

Singular integrands and Gauss Quadrature.

You often have to deal with integrals such as 

∫(1/(Sqrt[x+1]))dx
which are certainly well-defined, but have a singularity. Does Gauss quadrature help? Not really, as the examples below show. The real problem is that Gauss quadrature is based on integration with result

Integrate[1/Sqrt[x+1],{x,-1,1}]

N[Integrate[1/Sqrt[x+1],{x,-1,1}],20]

GQData = N[Table[GaussQuadrature[1/Sqrt[x+1],i],{i,2,5}],20]

Gauss Quadrature with respect to a weight function.

Note: The following code doesn't quite work. Since the class is moving on, I'm posting it for historical value, but I won't debug it now.

Suppose we have a collection of integrals in the form 

polynomial(x) * exp(x) 

or 

polynomial(x) * 1/sqrt(1+x)

or 

polynomial(x) * sin(x)

As long as the nonpolynomial (weight) function has a single sign, we can construct specialized Gauss quadrature rules for functions in this form without any new ideas. We start with a new inner product

The weighted inner product of functions.

wIP[f_,g_,w_] := Integrate[f*g*w,{x,-1,1}];

We try a few simple examples:

wIP[x^2,x^3 + 4x + 2,1/Sqrt[1+x]]

wIP[x,1,1/Sqrt[1+x]]

Weighted Gram-Schmidt Orthogonalization

We now need to define weighted Gram-Schmidt orthogonalization given a set of polynomials. This is exactly the same as ordinary Gram-Schmidt orthogonalization, except that we use the weighted inner product wIP.

wGramSchmidtCoefficient[p_,q_,w_] :=wIP[p,q,w]/wIP[q,q,w]; 
wGramSchmidt[p_,Q_,w_] := p - Sum[wGramSchmidtCoefficient[p,Q[[i]],w]*Q[[i]],{i,1,Length[Q]}];

We now try this out:

wGramSchmidt[x,{1},1/Sqrt[1+x]]

wGramSchmidt[x^2,{1,x},1/Sqrt[1+x]]

Putting everything together.

We now recursively apply GramSchmidt.

wLegendrePolylist = {1};
Do[wLegendrePolylist = Append[wLegendrePolylist,wGramSchmidt[x^n,wLegendrePolylist,1/Sqrt[1+x]]],{n,1,5}];
  
Simplify[Collect[wLegendrePolylist,x]]

This gives us the list of Legendre polynomials.

Generating a Quadrature Rule from the Legendre Polynomials

In order to generate a (numerical) quadrature rule from the collection of weighted Legendre polynomials, we need to find their roots. Luckily, Mathematica is pretty good at this.

wRL = Table[Solve[wLegendrePolylist[[i]] == 0, x],{i,2,5}];

As you can see, these are not rational numbers, but they are numbers that we can evaluate numerically relatively easily with Mathematica. To do this, we'll need to actually get these guys out of rule form and expressed just as numbers. Mathematica makes this easy by intelligently applying /. across the nested list which is RL.

wRootList = x /. wRL;

Evaluating numerically, we get

wGaussQuadratureNodes = N[wRootList,20]

Weights for Gaussian Quadrature.

We now have to figure out the appropriate (exact and numerical) weights for Gaussian quadrature. By construction, in the Gauss quadrature theorem, we know that 

A_i = Integral_a^b l_i(x). 

We could be more clever about the algebra if we were doing this by hand, but Mathematica is happy to construct and integrate the l_i for us. Suppose that x is the list of roots. Then we have

L_i = (∏_(j=1)^n)((x - (Subscript[x, j]))/((Subscript[x, i])-(Subscript[x, j]))),

where j=i is supposed to be missing from the sum, but I don't know how to write that!

L[i_,X_] :=
If[i>1,Product[(x - X[[j]])/(X[[i]] - X[[j]]),{j,1,i-1}],1]*
If[i<Length[X],Product[(x - X[[j]])/(X[[i]] - X[[j]]),{j,i+1,Length[X]}],1];

We can now try this out:

wRootList[[3]]

L[1,wRootList[[3]]]

L[2,wRootList[[3]]]

L[3,wRootList[[3]]]

We can now build a table of all the L_i with a couple of nested Table commands (you don't want to look at it):

wLTable = Table[
Table[L[i,wRootList[[j]]],{i,1,Length[wRootList[[j]]]}],
{j,2,4}];

We now want to integrate the L_i to get the weights A_i for Gauss Quadrature. For instance,

N[Integrate[L[3,wRootList[[3]]],{x,-1,1}],20]

In general, we should make a table of integrals, of course:

wATable = Integrate[wLTable,{x,-1,1}];

Remember that the RootList table starts with the polynomial x, with single root 0, while the LTable and ATable start one degree higher. To keep our notation straight, we just prepend a fake weight to the table.

wATable = Prepend[wATable,0];
N[wATable,20]

Testing the Gaussian Quadrature Rule (exactly)

We now have tables of nodes (wRootList and wGaussQuadratureNodes) and weights (wATable) for the weighted Gauss quadrature rules with between 2 and 5 nodes. The n node formula should correctly integrate all polynomials up to degree 2n -1. 

Of course, it will be no fun whatsoever unless we actually test this hypothesis on some real examples. First, we code the actual Gauss quadrature rule:

N[wATable[[3]]]
N[wRootList[[3]]]

wGaussQuadrature[f_,n_] := Sum[wATable[[n,i]]*(f /. {x ->wRootList[[n,i]]}),{i,1,Length[wRootList[[n]]]}];

N[wGaussQuadrature[(x^2 + 5x + 2)/Sqrt[1+x],2],20]

N[Integrate[(x^2+5x+2)/Sqrt[1+x],{x,-1,1}],20]