Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Oct 2016
Lecture topic: Gaussian Quadrature

Constructing Legendre Polynomials: Gram-Schmidt Process

Given an interval and a weight function, we are going to construct an orthogonal polynomial using the Gram-Schmidt process. This gives us another way of obtaining orthogonal polynomial that does not depend on the series solution to a differential equation!

Here, the weight function is w(x)=1and the interval is [-1,1]. We already know the orthogonal polynomials associated with this weight and interval are the Legendre polynomials, so that's what we should find.

Clear[ϕ,w,B,A,a,b,num,k,n]

a=-1;
b=1;
w[x_]=1;

(Subscript[ϕ, 0])[x_]=1;
num=4;

(Subscript[B, k_]):=(Subscript[B, k])=(((∫_a^b)x*w[x]*(((Subscript[ϕ, k-1])[x])^2)\[DifferentialD]x)/((∫_a^b)w[x]*(((Subscript[ϕ, k-1])[x])^2)\[DifferentialD]x))

(Subscript[A, k_]):=(Subscript[A, k])=(((∫_a^b)x*w[x]*(Subscript[ϕ, k-1])[x]*(Subscript[ϕ, k-2])[x]\[DifferentialD]x)/((∫_a^b)w[x]*(((Subscript[ϕ, k-2])[x])^2)\[DifferentialD]x))

(Subscript[ϕ, 1])[x_]=x-(Subscript[B, 1]);

Do[(Subscript[ϕ, k])[x_]=Expand[(x-(Subscript[B, k]))(Subscript[ϕ, k-1])[x]-(Subscript[A, k])(Subscript[ϕ, k-2])[x]],{k,2,num}]
TableForm[Table[{n,(Subscript[ϕ, n])[x]},{n,0,num}]]

TableForm[Table[{k,LegendreP[k,x]},{k,0,num}]]

We are missing the standardization (Subscript[P, n])(1)=1, but other than that, what we have found is the Legendre Polynomials!

Solving the System for Gaussian Quadrature: w(x)=1, n=2, and n=3

Simplest Case of getting weights and nodes: w(x)=1, (a,b)=(-1,1) and n=3

Here is the result for n=3:

Clear[eq,k,c,x,i,eqlist,varlist,n];
n=3;
Do[eq[k]=(∑_(i=1)^n)(Subscript[c, i])((((Subscript[x, i])))^k)==(∫_(-1)^1)(t^k)\[DifferentialD]t,{k,0,2n-1}];
eqlist={};
Do[eqlist=Append[eqlist,eq[i]],{i,0,2n-1}];

varlist={};
Do[varlist=Append[varlist,{(Subscript[c, i]),(Subscript[x, i])}],{i,1,n}];
varlist=Flatten[varlist];

eqlist//TableForm
varlist

Solve[eqlist,varlist]
N[%]

Since the equations we nonlinear we get multiple solutions, but notice all the solutions have the same (Subscript[x, i]),(Subscript[c, i]) pairs, and we call (Subscript[c, i]) the weight associated with node (Subscript[x, i]).

The important fact is that the nodes are always the roots of the orthogonal polynomial (Subscript[ϕ, n])(x) for interval [a,b] and weight function w(x). In this example, we had [a,b]=[-1,1] and w(x)=1 so we would want the roots of the Legendre polynomial.

Solve[LegendreP[n,x]==0,x]
N[%]

Examples of integrals

Here are the weights and nodes from above:

(Subscript[c, 1])=5/9;
(Subscript[c, 2])=8/9;
(Subscript[c, 3])=5/9;
(Subscript[x, 1])=-(Sqrt[(3/5)]);
(Subscript[x, 2])=0;
(Subscript[x, 3])=(Sqrt[(3/5)]);

Here is an integral, and our approximation:

f[t_]=123(t^4)-12(t^3)+13(t^2)-t+67
Integrate[f[t],{t,-1,1}];
N[%,24]
SetPrecision[(∑_(i=1)^3)(Subscript[c, i])*f[(Subscript[x, i])],24]

This is exact since the integrand was a polynomial of degree less than 2n-1=5 and we used the exact values for the weights and nodes.

Here is another example:

f[t_]=Sqrt[1-t^2]
Integrate[f[t],{t,-1,1}];
N[%,8]
SetPrecision[(∑_(i=1)^3)(Subscript[c, i])*f[(Subscript[x, i])],8]

Determining the weights via the formula in Theorem w(x)=1, (a,b)=(-1,1)

Clear[weight,xp]

w[x_]=1
weight[n_,k_]:=(∫_(-1)^1)w[x]((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))\[DifferentialD]x

Here are the weights. Notice it still depends on the xp[i].

weight[3,1]
weight[3,2]
weight[3,3]

The xp[i] are the roots of the Legendre polynomial.

NSolve[LegendreP[3,x]==0,x]

xp[1]=(Subscript[x, 1])=-0.774597;
xp[2]=(Subscript[x, 2])=0.0;
xp[3]=(Subscript[x, 3])=0.774597;

Now we can get the weights since we know the nodes.

weight[3,1]
weight[3,2]
weight[3,3]

Example using Laguerre orthogonal polynomial 
w(x)=(e^(-x)), (a,b)=(0,∞)

Clear[f]

f[x_]=Cos[x]
w[x_]=Exp[-x]
Integrate[w[x]f[x],{x,0,Infinity}]

Let's do this with n=5 and the Laguerre polynomial (I am cheating and using Mathematica to find the roots) :)

NSolve[LaguerreL[5,x]==0,x]

xp[1]=(Subscript[x, 1])=0.2635603197181409;
xp[2]=(Subscript[x, 2])=1.4134030591065165;
xp[3]=(Subscript[x, 3])=3.5964257710407197;
xp[4]=(Subscript[x, 4])=7.085810005858851;
xp[5]=(Subscript[x, 5])=12.64080084427577;

Now we calculate the weights (again, I will cheat and use Mathematica to do the integrals).

weight[n_,k_]:=(∫_0^∞)((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))w[x]\[DifferentialD]x

(Subscript[w, 1])=weight[5,1]
(Subscript[w, 2])=weight[5,2]
(Subscript[w, 3])=weight[5,3]
(Subscript[w, 4])=weight[5,4]
(Subscript[w, 5])=weight[5,5]

Notice that you have the sum of the weights being one. This is true for any of of n for the Laguerre polynomials. This isn't true for all orthogonal polynomials.

(∑_(i=1)^5)(Subscript[w, i])

And here is our estimate of the integral:

(∑_(i=1)^5)(Subscript[w, i])*f[(Subscript[x, i])]

Once this is set up you can do more complicated integrals using your formula.

f[x_]=1/(x^2+3)
(∑_(i=1)^5)(Subscript[w, i])*f[(Subscript[x, i])]
NIntegrate[w[x]f[x],{x,0,Infinity}]

Making n larger will improve the result. Typically you might have n as high as 16.

Example using Hermite orthogonal polynomial w(x)=(e^(-(x^2))), (a,b)=(-∞,∞)

Example 1

Clear[f]

f[x_]=1
w[x_]=Exp[-(x^2)]
Integrate[w[x]f[x],{x,-Infinity,Infinity}]

Let's do this with n=2 and the Hermite polynomial:

NSolve[HermiteH[2,x]==0,x]

xp[1]=(Subscript[x, 1])=-0.7071067811865475;
xp[2]=(Subscript[x, 2])=0.7071067811865475;

Now we calculate the weights:

weight[n_,k_]:=(∫_(-∞)^∞)((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))w[x]\[DifferentialD]x

(Subscript[w, 1])=weight[2,1]
(Subscript[w, 2])=weight[2,2]

And here is our estimate of the integral:

SetPrecision[(∑_(i=1)^2)(Subscript[w, i])*f[(Subscript[x, i])],16]

SetPrecision[(Sqrt[π]),15]

In fact, we will have the following for any value of n:

(∑_(i=1)^n)(Subscript[w, i])=(Sqrt[π])

Example 2: Maxwell Velocity Distribution

I have used this in the past when doing molecular calculations, and there is an integral over the velocity distribution of the gas you are studying. Things have to be modified slightly, but it works! And remember, the important use of this technique is that in application you don't always know the form of the function f(v) in your integrand. It is calculated once you know the velocities you are interested in.

Clear[dis,f,n,m,k,T]

The Hermite polynomials can be used for the Maxwell-Boltzmann Velocity Distribution in quantum chemsitry.

dis[v_]:=4π n((((m/(2π k T))))^(3/2))Exp[-((m (v^2))/(2k T))]*(v^2)

where

The mean speed of an ideal gas is given by the following integral:

Integrate[((dis[v]*v)/n),{v,0,∞},Assumptions->Re[(m/(k T))]>0]

For Helium gas at 70 K (units m[kg], k[J/K],T[K], v[m/s]), we have

m=4.0026*1.6605402*(10^(-27))
k=1.3806503*(10^(-23 ))
T=120

Here is the mean velocity (in meters/second):

result=SetPrecision[((2 (Sqrt[(2/π)]))/(Sqrt[(m/(k T))])),16]

Let's verify this integration numerically, using Hermite polynomials. First, what is the integrand?

((dis[v]*v)/n)

Let's split it into the weight part that depends on the exponential of velocity squared, which is the reason we approached this through the Hermite polynomials in the first place:

f[v_]=6.411062819705699*(10^(-9)) (v^3)
w[v_]= (E^(-2.0058416809455652*(10^(-6)) (v^2)))

Let's do this with n=6 and the Hermite polynomial. We need the zeros of the Hermite Polynomial:

n=6;
NSolve[HermiteH[n,x]==0,x]

xp[1]=(Subscript[x, 1])=-2.3506049736744923`;
xp[2]=(Subscript[x, 2])=-1.3358490740136968`;
xp[3]=(Subscript[x, 3])=-0.4360774119276165`;
xp[4]=(Subscript[x, 4])=0.4360774119276165`;
xp[5]=(Subscript[x, 5])=1.3358490740136968`;
xp[6]=(Subscript[x, 6])=2.3506049736744923`;

Now we calculate the weights. Notice that the integral has been modified to go from (0,∞).

weight[n_,k_]:=(∫_0^∞)((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))w[x]\[DifferentialD]x

(Subscript[w, 1])=weight[n,1]
(Subscript[w, 2])=weight[n,2]
(Subscript[w, 3])=weight[n,3]
(Subscript[w, 4])=weight[n,4]
(Subscript[w, 5])=weight[n,5]
(Subscript[w, 6])=weight[n,6]

And here is our estimate of the integral:

result
(Subscript[result, 1])=SetPrecision[(∑_(i=1)^n)(Subscript[w, i])*f[(Subscript[x, i])],16]

Surprisingly good results! You can see that you have to be careful, however, as the weights are large in absolute value and alternate in sign, so round off error can be a big factor. This results from the exponential function

w[v]

and so a change of variable might help fix this potentially damaging round off error problem. Also, since this weight function is not (e^(-(v^2))), we cannot use the weights other researchers have found! A change of variable will help with both of these issues.

Change of Variable Improves  Maxwell Velocity Distribution

((dis[v]*v)/n)

Let

-2.0058416809455652`*^-6 (v^2)==-(y^2)
Solve[%,v]

((dis[v]*v)/n)/.v-> 706.076364846701* y

Let's split it into the weight part that depends on the exponential of velocity squared, which is the reason we approached this through the Hermite polynomials in the first place:

f[y_]=2.256758334191026`  (y^3)* 706.076364846701
w[y_]=(E^(- (y^2)))

We can see that we have removed some small numbers. The important change is the weight function is (e^(- (y^2))). Note that extra 706 comes in from the change of the differential in the integral.

Let's do this with n=6 and the Hermite polynomial. This part is exactly the same as it was before we made out transformation!

n=6;
NSolve[HermiteH[n,x]==0,x]

xp[1]=(Subscript[x, 1])=-2.3506049736744923`;
xp[2]=(Subscript[x, 2])=-1.3358490740136968`;
xp[3]=(Subscript[x, 3])=-0.4360774119276165`;
xp[4]=(Subscript[x, 4])=0.4360774119276165`;
xp[5]=(Subscript[x, 5])=1.3358490740136968`;
xp[6]=(Subscript[x, 6])=2.3506049736744923`;

Now we calculate the weights. Notice that the integral has been modified to go from (0,∞), and the weight function is different from before:

weight[n_,k_]:=(∫_0^∞)((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))w[x]\[DifferentialD]x

(Subscript[w, 1])=weight[n,1]
(Subscript[w, 2])=weight[n,2]
(Subscript[w, 3])=weight[n,3]
(Subscript[w, 4])=weight[n,4]
(Subscript[w, 5])=weight[n,5]
(Subscript[w, 6])=weight[n,6]

And here is our estimate of the integral:

(Subscript[result, 2])=SetPrecision[(∑_(i=1)^n)(Subscript[w, i])*f[(Subscript[x, i])],16]

result
(Subscript[result, 1])
(Subscript[result, 2])

Same as before, but now the weights are not huge numbers, and we get a significantly better result. This is why you always want to change the vaiables to work with exactly the same weight function as the orthogonal polynomial is defined for.

The Mathematica Package for nodes and weights

We still need to load a package to get the weights.

<<"NumericalDifferentialEquationAnalysis`"

This returns the nodes, weights for a given n  in the interval a to b only for Legendre Gaussian quadrature.

GaussianQuadratureWeights[3,-1,1]

You can always get the nodes by using the Orthogonal polynomials

Clear[xp]
n=10;
roots=NSolve[ChebyshevU[n,x]==0]
Do[xp[i]=roots[[i,1,2]],{i,1,n}]

And from the nodes work out the weights:

Clear[w,weight]
w[x_]=(Sqrt[1-(x^2)]);
weight[n_,k_]:=weight[n,k]=(∫_(-1)^1)((∏_(i=1)^(k-1))((x-xp[i])/(xp[k]-xp[i])) )( (∏_(i=k+1)^n)((x-xp[i])/(xp[k]-xp[i])))w[x]\[DifferentialD]x
Do[weight[n,i],{i,1,n}]

Now do your integrals, once you have the nodes/weights all the hard work is done!

f[x_]=Cos[x];
SetPrecision[Sum[weight[n,i]*f[xp[i]],{i,1,n}],16]
SetPrecision[Integrate[w[x]f[x],{x,-1,1}],16]
Integrate[w[x]f[x],{x,-1,1}]

f[x_]=Cos[(x^2)];
SetPrecision[Sum[weight[n,i]*f[xp[i]],{i,1,n}],16]
SetPrecision[Integrate[w[x]f[x],{x,-1,1}],16]
Integrate[w[x]f[x],{x,-1,1}]

f[x_]=Cos[(x^2)]Sin[x];
SetPrecision[Sum[weight[n,i]*f[xp[i]],{i,1,n}],16]
SetPrecision[Integrate[w[x]f[x],{x,-1,1}],16]
Integrate[w[x]f[x],{x,-1,1}]

f[x_]=Exp[x]Tan[(x^2)];
SetPrecision[Sum[weight[n,i]*f[xp[i]],{i,1,n}],16]
SetPrecision[Integrate[w[x]f[x],{x,-1,1}],16]
Integrate[w[x]f[x],{x,-1,1}]