Gauss-Laguerre Quadrature Evaluation Points and Weights

This example shows how to solve polynomial equations and systems of equations, and work with the results using Symbolic Math Toolbox™.

Gaussian quadrature rules approximate an integral by sums abf(t)w(t)dti=1nf(xi)αi. Here, the xi and αi are parameters of the method, depending on n but not on f. They follow from the choice of the weight function w(t), as follows. Associated to the weight function is a family of orthogonal polynomials. The polynomials' roots are the evaluation points xi. Finally, the weights αi are determined by the condition that the method be correct for polynomials of small degree. Consider the weight function w(t)=exp(-t) on the interval [0,]. This case is known as Gauss-Laguerre quadrature.

syms t
n = 4;
w(t) = exp(-t);

Assume you know the first n members of the family of orthogonal polynomials. In case of the quadrature rule considered here, they turn out to be the Laguerre polynomials.

F = laguerreL(0:n-1, t)
F = 

(11-tt22-2t+1-t36+3t22-3t+1)

Let L be the n+1 st polynomial, the coefficients of which are still to be determined.

X = sym('X', [1, n+1])
X = (X1X2X3X4X5)
L = poly2sym(X, t)
L = X1t4+X2t3+X3t2+X4t+X5

Represent the orthogonality relations between the Laguerre polynomials F and L in a system of equations sys.

sys = [int(F.*L.*w(t), t, 0, inf) == 0]
sys = (24X1+6X2+2X3+X4+X5=0-96X1-18X2-4X3-X4=0144X1+18X2+2X3=0-96X1-6X2=0)

Add the condition that the polynomial have norm 1.

sys = [sys, int(L^2.*w(t), 0, inf) == 1]
sys = (24X1+6X2+2X3+X4+X5=0-96X1-18X2-4X3-X4=0144X1+18X2+2X3=0-96X1-6X2=040320X12+10080X1X2+1440X1X3+240X1X4+48X1X5+720X22+240X2X3+48X2X4+12X2X5+24X32+12X3X4+4X3X5+2X42+2X4X5+X52=1)

Solve for the coefficients of L.

S = solve(sys, X)
S = struct with fields:
    X1: [2x1 sym]
    X2: [2x1 sym]
    X3: [2x1 sym]
    X4: [2x1 sym]
    X5: [2x1 sym]

solve returns the two solutions in a structure array. Display the solutions.

structfun(@display, S)
ans = 

(-124124)

ans = 

(23-23)

ans = 

(-33)

ans = 

(4-4)

ans = 

(-11)

Make the solution unique by imposing an extra condition that the first coefficient be positive:

sys = [sys, X(1)>0];
S = solve(sys, X)
S = struct with fields:
    X1: [1x1 sym]
    X2: [1x1 sym]
    X3: [1x1 sym]
    X4: [1x1 sym]
    X5: [1x1 sym]

Substitute the solution into L.

L = subs(L, S)
L = 

t424-2t33+3t2-4t+1

As expected, this polynomial is the |n|th Laguerre polynomial:

laguerreL(n, t)
ans = 

t424-2t33+3t2-4t+1

The evaluation points xi are the roots of the polynomial L. Solve L for the evaluation points. The roots are expressed in terms of the root function.

x = solve(L)
x = 

(root(σ1,z,1)root(σ1,z,2)root(σ1,z,3)root(σ1,z,4))where  σ1=z4-16z3+72z2-96z+24

The form of the solutions might suggest that nothing has been achieved, but various operations are available on them. Compute floating-point approximations using vpa:

vpa(x)
ans = 

(0.322547689619392311800361459104371.74576110115834657568681671251794.53662029692112798327928538495719.3950709123011331292335364434205)

Some spurious imaginary parts might occur. Prove symbolically that the roots are real numbers:

isAlways(in(x, 'real'))
ans = 4x1 logical array

   1
   1
   1
   1

For polynomials of degree less than or equal to 4, you can use MaxDegree to obtain the solutions in terms of nested radicals instead in terms of root. However, subsequent operations on results of this form would be slow.

xradical = solve(L, 'MaxDegree', 4)
xradical = 

(4-σ1-σ34+σ1-σ3σ3-σ2+4σ3+σ2+4)where  σ1=96σ6σ4-3σ5σ4-288σ4-512366+36i2768+12836i1/6144σ6+9σ5+8641/4  σ2=96σ6σ4-3σ5σ4-288σ4+512366+36i2768+12836i1/6144σ6+9σ5+8641/4  σ3=σ42768+12836i1/6  σ4=16σ6+σ5+96  σ5=768+12836i2/3  σ6=768+12836i1/3

The weights αi are given by the condition that for polynomials of degree less than n, the quadrature rule must produce exact results. It is sufficient if this holds for a basis of the vector space of these polynomials. This condition results in a system of four equations in four variables.

y = sym('y', [n, 1]);
sys = sym(zeros(n));
for k=0:n-1 
    sys(k+1) = sum(y.*(x.^k)) == int(t^k * w(t), t, 0, inf); 
end
sys
sys = 

(y1+y2+y3+y4=1000y1root(σ1,z,1)+y2root(σ1,z,2)+y3root(σ1,z,3)+y4root(σ1,z,4)=1000y1root(σ1,z,1)2+y2root(σ1,z,2)2+y3root(σ1,z,3)2+y4root(σ1,z,4)2=2000y1root(σ1,z,1)3+y2root(σ1,z,2)3+y3root(σ1,z,3)3+y4root(σ1,z,4)3=6000)where  σ1=z4-16z3+72z2-96z+24

Solve the system both numerically and symbolically. The solution is the desired vector of weights αi.

[a1, a2, a3, a4] = vpasolve(sys, y)
a1 = 0.60315410434163360163596602381808
a2 = 0.35741869243779968664149201745809
a3 = 0.03888790851500538427243816815621
a4 = 0.00053929470556132745010379056762059
[alpha1, alpha2, alpha3, alpha4] = solve(sys, y)
alpha1 = 

-σ3σ2+σ3σ1+σ2σ1-σ3σ2σ1-2σ3-2σ2-2σ1+6σ4-σ3σ4σ2+σ4σ1-σ2σ1-σ42where  σ1=root(z4-16z3+72z2-96z+24,z,4)  σ2=root(z4-16z3+72z2-96z+24,z,3)  σ3=root(z4-16z3+72z2-96z+24,z,2)  σ4=root(z4-16z3+72z2-96z+24,z,1)

alpha2 = 

root(σ1,z,1)root(σ1,z,3)+root(σ1,z,1)root(σ1,z,4)+root(σ1,z,3)root(σ1,z,4)-root(σ1,z,1)root(σ1,z,3)root(σ1,z,4)-2root(σ1,z,1)-2root(σ1,z,3)-2root(σ1,z,4)+6root(σ1,z,2)-root(σ1,z,1)root(σ1,z,2)-root(σ1,z,3)root(σ1,z,2)-root(σ1,z,4)where  σ1=z4-16z3+72z2-96z+24

alpha3 = 

σ3σ2+σ3σ1+σ2σ1-σ3σ2σ1-2σ3-2σ2-2σ1+6σ4-σ1σ3σ2-σ3σ4-σ2σ4+σ42where  σ1=root(z4-16z3+72z2-96z+24,z,4)  σ2=root(z4-16z3+72z2-96z+24,z,2)  σ3=root(z4-16z3+72z2-96z+24,z,1)  σ4=root(z4-16z3+72z2-96z+24,z,3)

alpha4 = 

-σ3σ2+σ3σ1+σ2σ1-σ3σ2σ1-2σ3-2σ2-2σ1+6σ42σ3+σ42σ2+σ42σ1-σ43+σ3σ2σ1-σ3σ2σ4-σ3σ1σ4-σ2σ1σ4where  σ1=root(z4-16z3+72z2-96z+24,z,3)  σ2=root(z4-16z3+72z2-96z+24,z,2)  σ3=root(z4-16z3+72z2-96z+24,z,1)  σ4=root(z4-16z3+72z2-96z+24,z,4)

Alternatively, you can also obtain the solution as a structure by giving only one output argument.

S = solve(sys, y)
S = struct with fields:
    y1: [1x1 sym]
    y2: [1x1 sym]
    y3: [1x1 sym]
    y4: [1x1 sym]

structfun(@double, S)
ans = 4×1

    0.6032
    0.3574
    0.0389
    0.0005

Convert the structure S to a symbolic array:

Scell = struct2cell(S);
alpha = transpose([Scell{:}])
alpha = 

(-σ3+σ15+σ2-σ6-σ12-σ11-σ10+6root(σ16,z,1)-root(σ16,z,2)σ1+σ4-σ2-root(σ16,z,1)2σ1+σ4+σ2-σ7-σ13-σ11-σ10+6root(σ16,z,2)-root(σ16,z,1)root(σ16,z,2)-root(σ16,z,3)root(σ16,z,2)-root(σ16,z,4)σ5+σ4+σ15-σ8-σ13-σ12-σ10+6root(σ16,z,3)-root(σ16,z,4)σ5-σ1-σ3+root(σ16,z,3)2-σ5+σ1+σ3-σ9-σ13-σ12-σ11+6σ14root(σ16,z,1)+σ14root(σ16,z,2)+σ14root(σ16,z,3)-root(σ16,z,4)3+σ9-σ8-σ7-σ6)where  σ1=root(σ16,z,1)root(σ16,z,3)  σ2=root(σ16,z,3)root(σ16,z,4)  σ3=root(σ16,z,2)root(σ16,z,3)  σ4=root(σ16,z,1)root(σ16,z,4)  σ5=root(σ16,z,1)root(σ16,z,2)  σ6=root(σ16,z,2)root(σ16,z,3)root(σ16,z,4)  σ7=root(σ16,z,1)root(σ16,z,3)root(σ16,z,4)  σ8=root(σ16,z,1)root(σ16,z,2)root(σ16,z,4)  σ9=root(σ16,z,1)root(σ16,z,2)root(σ16,z,3)  σ10=2root(σ16,z,4)  σ11=2root(σ16,z,3)  σ12=2root(σ16,z,2)  σ13=2root(σ16,z,1)  σ14=root(σ16,z,4)2  σ15=root(σ16,z,2)root(σ16,z,4)  σ16=z4-16z3+72z2-96z+24

The symbolic solution looks complicated. Simplify it, and convert it into a floating point vector:

alpha = simplify(alpha)
alpha = 

(root(σ1,z,1)272-29root(σ1,z,1)144+23root(σ1,z,2)272-29root(σ1,z,2)144+23root(σ1,z,3)272-29root(σ1,z,3)144+23root(σ1,z,4)272-29root(σ1,z,4)144+23)where  σ1=z4-16z3+72z2-96z+24

vpa(alpha)
ans = 

(0.603154104341633601635966023818080.357418692437799686641492017458090.038887908515005384272438168156210.00053929470556132745010379056762059)

Increase the readability by replacing the occurrences of the roots x in alpha by abbreviations:

subs(alpha, x, sym('R', [4, 1]))
ans = 

(R1272-29R1144+23R2272-29R2144+23R3272-29R3144+23R4272-29R4144+23)

Sum the weights to show that their sum equals 1:

simplify(sum(alpha))
ans = 1

A different method to obtain the weights of a quadrature rule is to compute them using the formula αi=abw(t)jit-xjxi-xjdt. Do this for i=1. It leads to the same result as the other method:

int(w(t) * prod((t - x(2:4)) ./ (x(1) - x(2:4))), t, 0, inf)
ans = 

root(z4-16z3+72z2-96z+24,z,1)272-29root(z4-16z3+72z2-96z+24,z,1)144+23

The quadrature rule produces exact results even for all polynomials of degree less than or equal to 2n-1, but not for t2n.

simplify(sum(alpha.*(x.^(2*n-1))) -int(t^(2*n-1)*w(t), t, 0, inf))
ans = 0
simplify(sum(alpha.*(x.^(2*n))) -int(t^(2*n)*w(t), t, 0, inf))
ans = -576

Apply the quadrature rule to the cosine, and compare with the exact result:

vpa(sum(alpha.*(cos(x))))
ans = 0.50249370546067059229918484198931
int(cos(t)*w(t), t, 0, inf)
ans = 

12

For powers of the cosine, the error oscillates between odd and even powers:

errors = zeros(1, 20);
for k=1:20 
    errors(k) = double(sum(alpha.*(cos(x).^k)) -int(cos(t)^k*w(t), t, 0, inf));
end
plot(real(errors))