# 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 ${\int }_{a}^{b}f\left(t\right)w\left(t\right)dt\approx \sum _{i=1}^{n}f\left({x}_{i}\right){\alpha }_{i}$. Here, the ${x}_{i}$ and ${\alpha }_{i}$ are parameters of the method, depending on $n$ but not on $f$. They follow from the choice of the weight function $w\left(t\right)$, 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 ${\alpha }_{i}$ are determined by the condition that the method be correct for polynomials of small degree. Consider the weight function $w\left(t\right)=\mathrm{exp}\left(-t\right)$ on the interval $\left[0,\infty \right]$. 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 =  $\left(\begin{array}{cccc}1& 1-t& \frac{{t}^{2}}{2}-2 t+1& -\frac{{t}^{3}}{6}+\frac{3 {t}^{2}}{2}-3 t+1\end{array}\right)$```

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

`X = sym('X', [1, n+1])`
`X = $\left(\begin{array}{ccccc}{X}_{1}& {X}_{2}& {X}_{3}& {X}_{4}& {X}_{5}\end{array}\right)$`
`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 = $\left(\begin{array}{cccc}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}\right)$`

Add the condition that the polynomial have norm 1.

`sys = [sys, int(L^2.*w(t), 0, inf) == 1]`
`sys = $\left(\begin{array}{ccccc}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}\right)$`

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 =  $\left(\begin{array}{c}-\frac{1}{24}\\ \frac{1}{24}\end{array}\right)$```
```ans =  $\left(\begin{array}{c}\frac{2}{3}\\ -\frac{2}{3}\end{array}\right)$```
```ans =  $\left(\begin{array}{c}-3\\ 3\end{array}\right)$```
```ans =  $\left(\begin{array}{c}4\\ -4\end{array}\right)$```
```ans =  $\left(\begin{array}{c}-1\\ 1\end{array}\right)$```

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 =  $\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 =  ```

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

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 =  ```

The weights ${\alpha }_{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 =  ```

Solve the system both numerically and symbolically. The solution is the desired vector of weights ${\alpha }_{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 =  ```
```alpha2 =  ```
```alpha3 =  ```
```alpha4 =  ```

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 =  ```

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

`alpha = simplify(alpha)`
```alpha =  ```
`vpa(alpha)`
```ans =  $\left(\begin{array}{c}0.60315410434163360163596602381808\\ 0.35741869243779968664149201745809\\ 0.03888790851500538427243816815621\\ 0.00053929470556132745010379056762059\end{array}\right)$```

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

`subs(alpha, x, sym('R', [4, 1]))`
```ans =  $\left(\begin{array}{c}\frac{{{R}_{1}}^{2}}{72}-\frac{29 {R}_{1}}{144}+\frac{2}{3}\\ \frac{{{R}_{2}}^{2}}{72}-\frac{29 {R}_{2}}{144}+\frac{2}{3}\\ \frac{{{R}_{3}}^{2}}{72}-\frac{29 {R}_{3}}{144}+\frac{2}{3}\\ \frac{{{R}_{4}}^{2}}{72}-\frac{29 {R}_{4}}{144}+\frac{2}{3}\end{array}\right)$```

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 ${\alpha }_{i}={\int }_{a}^{b}w\left(t\right)\prod _{j\ne i}\frac{t-{x}_{j}}{{x}_{i}-{x}_{j}}dt$. 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{{\mathrm{root}\left({z}^{4}-16 {z}^{3}+72 {z}^{2}-96 z+24,z,1\right)}^{2}}{72}-\frac{29 \mathrm{root}\left({z}^{4}-16 {z}^{3}+72 {z}^{2}-96 z+24,z,1\right)}{144}+\frac{2}{3}$```

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

`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))``` ## Support

#### Mathematical Modeling with Symbolic Math Toolbox

Get examples and videos