Gauss-Laguerre Quadrature Evaluation Points and Weights

Open Live Script

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 $\int_{a}^{b} f (t) w (t) d t \approx \sum_{i = 1}^{n} f (x_{i}) α_{i}$ . Here, the $x_{i}$ 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 $x_{i}$ . 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, \infty]$ . This case is known as the 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 = 
 $(\begin{array}{c} 1 & 1 - t & \frac{t^{2}}{2} - 2 t + 1 & - \frac{t^{3}}{6} + \frac{3 t^{2}}{2} - 3 t + 1 \end{array})$

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

X = sym('X', [1, n+1])

X =  $(\begin{array}{c} X_{1} & X_{2} & X_{3} & X_{4} & X_{5} \end{array})$

L = poly2sym(X, t)

L =  $X_{1} t^{4} + X_{2} t^{3} + X_{3} t^{2} + X_{4} t + X_{5}$

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 =  $(\begin{array}{c} 24 X_{1} + 6 X_{2} + 2 X_{3} + X_{4} + X_{5} = 0 & - 96 X_{1} - 18 X_{2} - 4 X_{3} - X_{4} = 0 & 144 X_{1} + 18 X_{2} + 2 X_{3} = 0 & - 96 X_{1} - 6 X_{2} = 0 \end{array})$

Add the condition that the polynomial have norm 1.

sys = [sys, int(L^2.*w(t), 0, inf) == 1]

sys =  $(\begin{array}{c} 24 X_{1} + 6 X_{2} + 2 X_{3} + X_{4} + X_{5} = 0 & - 96 X_{1} - 18 X_{2} - 4 X_{3} - X_{4} = 0 & 144 X_{1} + 18 X_{2} + 2 X_{3} = 0 & - 96 X_{1} - 6 X_{2} = 0 & 40320 {X_{1}}^{2} + 10080 X_{1} X_{2} + 1440 X_{1} X_{3} + 240 X_{1} X_{4} + 48 X_{1} X_{5} + 720 {X_{2}}^{2} + 240 X_{2} X_{3} + 48 X_{2} X_{4} + 12 X_{2} X_{5} + 24 {X_{3}}^{2} + 12 X_{3} X_{4} + 4 X_{3} X_{5} + 2 {X_{4}}^{2} + 2 X_{4} X_{5} + {X_{5}}^{2} = 1 \end{array})$

Solve for the coefficients of L.

S = solve(sys, X)

S = struct with fields:
    X1: [2×1 sym]
    X2: [2×1 sym]
    X3: [2×1 sym]
    X4: [2×1 sym]
    X5: [2×1 sym]

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

structfun(@display, S)

ans = 
 $(\begin{array}{c} - \frac{1}{24} \\ \frac{1}{24} \end{array})$

ans = 
 $(\begin{array}{c} \frac{2}{3} \\ - \frac{2}{3} \end{array})$

ans = 
 $(\begin{array}{c} - 3 \\ 3 \end{array})$

ans = 
 $(\begin{array}{c} 4 \\ - 4 \end{array})$

ans = 
 $(\begin{array}{c} - 1 \\ 1 \end{array})$

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: 1/24
    X2: -2/3
    X3: 3
    X4: -4
    X5: 1

Substitute the solution into L.

L = subs(L, S)

L = 
 $\frac{t^{4}}{24} - \frac{2 t^{3}}{3} + 3 t^{2} - 4 t + 1$

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

laguerreL(n, t)

ans = 
 $\frac{t^{4}}{24} - \frac{2 t^{3}}{3} + 3 t^{2} - 4 t + 1$

The evaluation points $x_{i}$ 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 = 
 $\begin{array}{l} (\begin{array}{c} root (σ_{1}, z, 1) \\ root (σ_{1}, z, 2) \\ root (σ_{1}, z, 3) \\ root (σ_{1}, z, 4) \end{array}) \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

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 = 
 $(\begin{array}{c} 0.32254768961939231180036145910437 \\ 1.7457611011583465756868167125179 \\ 4.5366202969211279832792853849571 \\ 9.3950709123011331292335364434205 \end{array})$

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

isAlways(in(x, 'real'))

ans = 4×1 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 = 
 $\begin{array}{l} (\begin{array}{c} 4 - σ_{1} - σ_{3} \\ 4 + σ_{1} - σ_{3} \\ σ_{3} - σ_{2} + 4 \\ σ_{3} + σ_{2} + 4 \end{array}) \\ where \\ σ_{1} = \frac{\sqrt{96 σ_{6} σ_{4} - 3 σ_{5} σ_{4} - 288 σ_{4} - 512 \sqrt{3} \sqrt{6} \sqrt{6 + \sqrt{3} \sqrt{6} i}}}{2 {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 6} {(144 σ_{6} + 9 σ_{5} + 864)}^{1 / 4}} \\ σ_{2} = \frac{\sqrt{96 σ_{6} σ_{4} - 3 σ_{5} σ_{4} - 288 σ_{4} + 512 \sqrt{3} \sqrt{6} \sqrt{6 + \sqrt{3} \sqrt{6} i}}}{2 {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 6} {(144 σ_{6} + 9 σ_{5} + 864)}^{1 / 4}} \\ σ_{3} = \frac{σ_{4}}{2 {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 6}} \\ σ_{4} = \sqrt{16 σ_{6} + σ_{5} + 96} \\ σ_{5} = {(768 + 128 \sqrt{3} \sqrt{6} i)}^{2 / 3} \\ σ_{6} = {(768 + 128 \sqrt{3} \sqrt{6} i)}^{1 / 3} \end{array}$

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 = 
 $\begin{array}{l} (\begin{array}{c} y_{1} + y_{2} + y_{3} + y_{4} = 1 & 0 & 0 & 0 \\ y_{1} root (σ_{1}, z, 1) + y_{2} root (σ_{1}, z, 2) + y_{3} root (σ_{1}, z, 3) + y_{4} root (σ_{1}, z, 4) = 1 & 0 & 0 & 0 \\ y_{1} {root (σ_{1}, z, 1)}^{2} + y_{2} {root (σ_{1}, z, 2)}^{2} + y_{3} {root (σ_{1}, z, 3)}^{2} + y_{4} {root (σ_{1}, z, 4)}^{2} = 2 & 0 & 0 & 0 \\ y_{1} {root (σ_{1}, z, 1)}^{3} + y_{2} {root (σ_{1}, z, 2)}^{3} + y_{3} {root (σ_{1}, z, 3)}^{3} + y_{4} {root (σ_{1}, z, 4)}^{3} = 6 & 0 & 0 & 0 \end{array}) \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

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 = 
 $\begin{array}{l} - \frac{σ_{3} σ_{2} + σ_{3} σ_{1} + σ_{2} σ_{1} - σ_{3} σ_{2} σ_{1} - 2 σ_{3} - 2 σ_{2} - 2 σ_{1} + 6}{(σ_{4} - σ_{3}) (σ_{4} σ_{2} + σ_{4} σ_{1} - σ_{2} σ_{1} - {σ_{4}}^{2})} \\ where \\ σ_{1} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 4) \\ σ_{2} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 3) \\ σ_{3} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 2) \\ σ_{4} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1) \end{array}$

alpha2 = 
 $\begin{array}{l} \frac{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) - 2 root (σ_{1}, z, 1) - 2 root (σ_{1}, z, 3) - 2 root (σ_{1}, z, 4) + 6}{(root (σ_{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} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

alpha3 = 
 $\begin{array}{l} \frac{σ_{3} σ_{2} + σ_{3} σ_{1} + σ_{2} σ_{1} - σ_{3} σ_{2} σ_{1} - 2 σ_{3} - 2 σ_{2} - 2 σ_{1} + 6}{(σ_{4} - σ_{1}) (σ_{3} σ_{2} - σ_{3} σ_{4} - σ_{2} σ_{4} + {σ_{4}}^{2})} \\ where \\ σ_{1} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 4) \\ σ_{2} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 2) \\ σ_{3} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1) \\ σ_{4} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 3) \end{array}$

alpha4 = 
 $\begin{array}{l} - \frac{σ_{3} σ_{2} + σ_{3} σ_{1} + σ_{2} σ_{1} - σ_{3} σ_{2} σ_{1} - 2 σ_{3} - 2 σ_{2} - 2 σ_{1} + 6}{{σ_{4}}^{2} σ_{3} + {σ_{4}}^{2} σ_{2} + {σ_{4}}^{2} σ_{1} - {σ_{4}}^{3} + σ_{3} σ_{2} σ_{1} - σ_{3} σ_{2} σ_{4} - σ_{3} σ_{1} σ_{4} - σ_{2} σ_{1} σ_{4}} \\ where \\ σ_{1} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 3) \\ σ_{2} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 2) \\ σ_{3} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1) \\ σ_{4} = root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 4) \end{array}$

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: -(root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 2)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 3) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 2)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 4) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 3)*root(z^4…
    y2: (root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 1)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 3) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 1)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 4) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 3)*root(z^4 …
    y3: (root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 1)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 2) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 1)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 4) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 2)*root(z^4 …
    y4: -(root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 1)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 2) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 1)*root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 3) + root(z^4 - 16*z^3 + 72*z^2 - 96*z + 24, z, 2)*root(z^4…

structfun(@double, S)

Convert the structure S to a symbolic array:

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

alpha = 
 $\begin{array}{l} (\begin{array}{c} - \frac{σ_{3} + σ_{15} + σ_{2} - σ_{6} - σ_{12} - σ_{11} - σ_{10} + 6}{(root (σ_{16}, z, 1) - root (σ_{16}, z, 2)) (σ_{1} + σ_{4} - σ_{2} - {root (σ_{16}, z, 1)}^{2})} \\ \frac{σ_{1} + σ_{4} + σ_{2} - σ_{7} - σ_{13} - σ_{11} - σ_{10} + 6}{(root (σ_{16}, z, 2) - root (σ_{16}, z, 1)) (root (σ_{16}, z, 2) - root (σ_{16}, z, 3)) (root (σ_{16}, z, 2) - root (σ_{16}, z, 4))} \\ \frac{σ_{5} + σ_{4} + σ_{15} - σ_{8} - σ_{13} - σ_{12} - σ_{10} + 6}{(root (σ_{16}, z, 3) - root (σ_{16}, z, 4)) (σ_{5} - σ_{1} - σ_{3} + {root (σ_{16}, z, 3)}^{2})} \\ - \frac{σ_{5} + σ_{1} + σ_{3} - σ_{9} - σ_{13} - σ_{12} - σ_{11} + 6}{σ_{14} root (σ_{16}, z, 1) + σ_{14} root (σ_{16}, z, 2) + σ_{14} root (σ_{16}, z, 3) - {root (σ_{16}, z, 4)}^{3} + σ_{9} - σ_{8} - σ_{7} - σ_{6}} \end{array}) \\ 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} = 2 root (σ_{16}, z, 4) \\ σ_{11} = 2 root (σ_{16}, z, 3) \\ σ_{12} = 2 root (σ_{16}, z, 2) \\ σ_{13} = 2 root (σ_{16}, z, 1) \\ σ_{14} = {root (σ_{16}, z, 4)}^{2} \\ σ_{15} = root (σ_{16}, z, 2) root (σ_{16}, z, 4) \\ σ_{16} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

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

alpha = simplify(alpha)

alpha = 
 $\begin{array}{l} (\begin{array}{c} \frac{{root (σ_{1}, z, 1)}^{2}}{72} - \frac{29 root (σ_{1}, z, 1)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z, 2)}^{2}}{72} - \frac{29 root (σ_{1}, z, 2)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z, 3)}^{2}}{72} - \frac{29 root (σ_{1}, z, 3)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z, 4)}^{2}}{72} - \frac{29 root (σ_{1}, z, 4)}{144} + \frac{2}{3} \end{array}) \\ where \\ σ_{1} = z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24 \end{array}$

vpa(alpha)

ans = 
 $(\begin{array}{c} 0.60315410434163360163596602381808 \\ 0.35741869243779968664149201745809 \\ 0.03888790851500538427243816815621 \\ 0.00053929470556132745010379056762059 \end{array})$

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

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

ans = 
 $\begin{array}{l} (\begin{array}{c} \frac{{R_{1}}^{2}}{72} - \frac{29 R_{1}}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z_{1}, 2)}^{2}}{72} - \frac{29 root (σ_{1}, z_{1}, 2)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z_{1}, 3)}^{2}}{72} - \frac{29 root (σ_{1}, z_{1}, 3)}{144} + \frac{2}{3} \\ \frac{{root (σ_{1}, z_{1}, 4)}^{2}}{72} - \frac{29 root (σ_{1}, z_{1}, 4)}{144} + \frac{2}{3} \end{array}) \\ where \\ σ_{1} = {z_{1}}^{4} - 16 {z_{1}}^{3} + 72 {z_{1}}^{2} - 96 z_{1} + 24 \end{array}$

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} = \int_{a}^{b} w (t) \prod_{j \neq i} \frac{t - x_{j}}{x_{i} - x_{j}} d t$ . 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 = 
 $\frac{{root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1)}^{2}}{72} - \frac{29 root (z^{4} - 16 z^{3} + 72 z^{2} - 96 z + 24, z, 1)}{144} + \frac{2}{3}$

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

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 = 
 $\frac{1}{2}$

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))

Figure contains an axes object. The axes object contains an object of type line.