Skip to main content

Numerical Analysis A look at rootfinding, quadrature and approximation using analytic and computing tools

Peter Staab

Section 5.4 Gaussian Quadrature

Recall that a quadrature rule approximates \(\int_{a}^{b} f(x) \, dx\) with a formula
\begin{equation*} \begin{aligned}I_{n} (f)\amp = \sum_{i=0}^{n} w_{i} f(x_{i}),\end{aligned} \end{equation*}
for \(x_{i}\) nodes (or quadrature points) and \(w_{i}\) are the weights. The concept of Gaussian Quadrature is the formula that determines the weights and nodes with the highest degree of precision.
🔗
Also, in the derivations below, we will find the integral with \(a=-1\) and \(b=1\) for simplicity. For general intervals, a linear conversion can be made.
🔗

Subsection 5.4.1 2-point Gaussian Quadrature Formula

In this case the quadrature formula will be
\begin{equation*} I_{2} (f) = w_{0} f(x_{0}) + w_{1} f(x_{1}) \end{equation*}
and we want to find \(x_{0}, x_{1}, w_{0}, w_{1}\text{,}\) which is in contrast to Newton-Cotes formulas where we knew \(x_{0}\) and \(x_{1}\text{.}\)
🔗
\begin{equation} \int_{-1}^{1} x^{n} \, dx = \begin{cases}\frac{2}{n+1} \amp \text{for $n$ even} \\ 0 \amp \text{for $n$ odd}\end{cases}\tag{5.4.1} \end{equation}
🔗
Applying this formula for \(n=0,1,2,3\) results in the following nonlinear system:
\begin{align} 2\amp = w_{0} + w_{1} \tag{5.4.2}\\ 0\amp = w_{0} x_{0} + w_{1} x_{1} \tag{5.4.3}\\ \frac{2}{3}\amp = w_{0} x_{0}^{2} + w_{1} x_{1}^{2} \tag{5.4.4}\\ 0\amp = w_{0} x_{0}^{3} + w_{1} x_{1}^{3}\tag{5.4.5} \end{align}
🔗
Next, we solve this nonlinear system. First, find \(w_{0}\) in (5.4.3),
\begin{equation} w_{0} = -\frac{w_{1} x_{1}}{x_{0}}\tag{5.4.6} \end{equation}
and substitute this into (5.4.5) to arrive at
\begin{equation*} \begin{aligned}0\amp = - \frac{w_{1} x_{1}}{x_{0}}x_{0}^{3} + w_{1} x_{1}^{3} \\\amp = w_{1} x_{1} ( x_{0}^{2}-x_{1}^{2}).\end{aligned} \end{equation*}
🔗
This results in either \(w_{1}=0\text{,}\) \(x_{1}=0\text{,}\) \(x_{0}=x_{1}\) or \(x_{0}=-x_{1}\text{.}\)
🔗
If (5.4.6) is substituted into (5.4.4), one gets
\begin{equation*} \begin{aligned}\frac{2}{3}\amp = - w_{1} \frac{x_{1}}{x_{0}}x_{0}^{2} + w_{1} x_{1}^{2}, \\ \frac{2}{3}\amp = w_{1} x_{1} (x_{1}-x_{0}).\end{aligned} \end{equation*}
🔗
If \(w_{1}=0\text{,}\) \(x_{1}=0\) or \(x_{1}=x_{0}\text{,}\) there is no solution to this, therefore, \(x_{0}=-x_{1}\text{.}\) Using this result in (5.4.3) means that \(w_{0}=w_{1}\) and from (5.4.2), both weights must be 1.
🔗
Lastly, substituting \(x_{0}=-x_{1}\text{,}\) \(w_{0}=w_{1}=1\) into (5.4.4), we get:
\begin{equation} \begin{aligned}\frac{2}{3}\amp = x_{1}^{2} + x_{1}^{2} \\ \text{or} \\ x_{1}^{2} - \frac{1}{3}\amp = 0 \end{aligned}\tag{5.4.7} \end{equation}
🔗
so \(x_{1} = \sqrt{1/3}\) and \(x_{0}=-\sqrt{1/3}\text{.}\)
🔗
The 2-point Gaussian Quadrature formula for \(f(x)\) on \([-1,1]\) is given by
\begin{equation} I_{2} (f) = f\biggl( \sqrt{\frac{1}{3}}\biggr) + f\biggl( -\sqrt{\frac{1}{3}}\biggr)\tag{5.4.8} \end{equation}
🔗
We will see in Chapter 6 that the quadratic in (5.4.7) is special. It is called a Legendre polynomial and the roots of it are the nodes for Gaussian Quadrature. The next example shows an application of Gaussian quadrature.
🔗

Example 5.4.1.

Use the two-point Gaussian Quadrature formula in (5.4.8) to approximate \(\int_{-1}^{1} \cos (\pi x/2) \, dx\text{.}\)
🔗
In this case,
\begin{equation*} \begin{aligned}I_{2} (f)\amp = \cos \biggl( \frac{\pi}{2}\sqrt{\frac{1}{3}}\biggr) + \cos \biggl( -\frac{\pi}{2}\sqrt{\frac{1}{3}}\biggr) \approx 1.23238101695912\end{aligned} \end{equation*}
🔗
Note that the actual value of the integral is \(4/\pi \approx 1.27323954473516\text{.}\)
🔗
🔗
The above example shows the usefulness of Gaussian quadrature. Since \(f\) evaluated at each endpoint is 0, the trapezoid rule applied to this is 0 and Simpson’s rule is 4/3. By selecting other nodes, Gaussian quadrature is often much more accurate than the Newton-Cotes formulas.
🔗
🔗

Subsection 5.4.2 Gaussian Quadrature on \([a,b]\)

In this case, we translate the interval \([-1,1]\) to \([a,b]\) as we have seen before. If \(x \in [-1,1]\) and \(z \in [a,b]\) then
\begin{equation} z = \frac{a+b}{2}+ x \frac{b-a}{2}.\tag{5.4.9} \end{equation}
🔗
Therefore the 2-point Gaussian quadrature formula in (5.4.8) can be written as
\begin{align} \int_{a}^{b} f(z) \, dz\amp = \frac{b-a}{2}\int_{-1}^{1} f\biggl(\frac{a+b}{2}+ x \frac{b-a}{2}\biggr) \, dx, \notag\\ \amp \approx \frac{b-a}{2}\biggl( f\biggl(\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}\biggr) + f\biggl(\frac{a+b}{2}+ \frac{b-a}{2\sqrt{3}}\biggr) \tag{5.4.10} \end{align}
🔗

Example 5.4.2.

Use (5.4.10) to approximate \(\int_{0}^{1} e^{-x^2}\, dx\text{.}\)
🔗
Solution.
In this case \(a=0\text{,}\) \(b=1\) and we get
\begin{equation*} \begin{aligned}\int_{0}^{1} e^{-x^2}\, dx\amp \approx \frac{1}{2}\biggl( f\biggl(\frac{1}{2}- \frac{1}{2\sqrt{3}}\biggr) + f\biggl(\frac{1}{2}+ \frac{1}{2\sqrt{3}}\biggr) \biggr) \\\amp = 0.7465946882\end{aligned} \end{equation*}
which has error of 2.29 \(\times 10^{-4}\text{.}\) Note that this error is better than the composite trapezoid rule with \(n=10\) seen in the previous section in Example 5.3.1.
🔗
🔗
🔗
🔗

Subsection 5.4.3 3-point Gaussian Quadrature Formula

The next more-accurate Gaussian Quadrature Formula is using 3 nodes. That is, if we let
\begin{equation*} \begin{aligned}I_{3}(f)\amp = w_{0} f(x_{0}) + w_{1} f(x_{1}) + w_{2} f(x_{2})\end{aligned} \end{equation*}
to have the highest precision on the interval \([-1,1]\text{.}\) If we apply (5.4.1) for \(n=0,1,\ldots,5\text{,}\) the result is:
\begin{equation*} \begin{aligned}2\amp = w_{0} +w_{1}+w_{2}\amp 0\amp = w_{0} x_{0} + w_{1} x_{1} + w_{2} x_{2} \\ \frac{2}{3}\amp = w_{0} x_{0}^{2} + w_{1} x_{1}^{2} + w_{2} x_{2}^{2}\amp 0\amp = w_{0} x_{0}^{3} + w_{1} x_{1}^{3} + w_{2} x_{2}^{3} \\ \frac{2}{5}\amp = w_{0} x_{0}^{4} + w_{1} x_{1}^{4} + w_{2} x_{2}^{4}\amp 0\amp = w_{0} x_{0}^{5} + w_{1} x_{1}^{5} + w_{2} x_{2}^{5} \\\end{aligned} \end{equation*}
🔗
Although this may be a nightmare to try to solve by hand, a CAS (like Maple) can generate the solution as:
\begin{equation*} \begin{aligned}w_{0}\amp = \frac{5}{9}\amp w_{1}\amp = \frac{8}{9}\amp w_{2}\amp = \frac{5}{9}\\ x_{0}\amp = -\sqrt{\frac{3}{5}}\amp x_{1}\amp = 0\amp x_{2}\amp = \sqrt{\frac{3}{5}}.\end{aligned} \end{equation*}
The 3-point Gaussian Quadrature formula for approximating \(\int_{-1}^{1} f(x) \, dx\) is:
\begin{equation} \int_{-1}^{1}f(x) \, dx \approx \frac{5}{9}f\biggl( -\sqrt{\frac{3}{5}}\biggr) + \frac{8}{9}f(0) + \frac{5}{9}f\biggl( \sqrt{\frac{3}{5}}\biggr)\tag{5.4.11} \end{equation}
🔗
Note, the solutions for \(x_{i}\) above are the roots of the polynomial
\begin{equation*} P(x) = x^{3} - \frac{3}{5}x \end{equation*}
which is also a Legendre polynomial. We will discuss this is more generality in Chapter 6.
🔗

Example 5.4.3.

Use the 3-point Gaussian Quadrature formula in (5.4.11) as well as the interval conversion in (5.4.9) to find an approximation to
\begin{equation*} \int_{0}^{1} e^{-x^2}\, dx \end{equation*}
🔗
Solution.
Using (5.4.9), the three points \(z_{0},z_{1},z_{2}\) become:
\begin{equation*} \begin{aligned}z_{0}\amp = \frac{1}{2}+ \frac{1}{2}x_{0} = \frac{1}{2}- \frac{1}{2}\sqrt{\frac{3}{5}}\\ z_{1}\amp = \frac{1}{2}+ \frac{1}{2}x_{1} = \frac{1}{2}\\ z_{2}\amp = \frac{1}{2}+ \frac{1}{2}x_{2} = \frac{1}{2}+ \frac{1}{2}\sqrt{\frac{3}{5}}\\\end{aligned} \end{equation*}
and using \(f(x)=e^{-x^2}\text{,}\)
\begin{equation*} \begin{aligned}\int_{0}^{1} e^{-x^2}\amp = \frac{1}{2}\biggl(\frac{5}{9}f\biggl(\frac{1}{2}- \frac{1}{2}\sqrt{\frac{3}{5}}\biggr) + \frac{8}{9}f\biggl(\frac{1}{2}\biggr) + \frac{5}{9}f\biggl(\frac{1}{2}+ \frac{1}{2}\sqrt{\frac{3}{5}}\biggr) \biggr) \\\amp \approx 0.7468145842\end{aligned} \end{equation*}
🔗
Note that the \(1/2\) that multiplies in front is the term \((b-a)/2\) due to the change in length of the interval. Also, the quadrature result is within \(10^{-5}\) of the correct answer.
🔗
🔗
🔗
🔗

Subsection 5.4.4 Error Formula of Gaussian Quadrature

The error associated with Gaussian quadrature can be found in much the same manner that we found for the Newton-Cotes formulas in Section 5.2. In short, the error associated with the interpolated polynomial is used and upon integration an error formula can be found. We won’t go into the details here mainly because the irregular spacing of the nodes (roots of Legendre polynomials) make the calculation more difficult than we saw for Newton-Cotes.
🔗
🔗

Subsection 5.4.5 Higher-Order Gaussian Quadrature

In Subsection 6.5.2, we will show why using roots of Legendre Polynomials for the quadrature notes results in the highest degree of precision possible. For this section, let’s assume that the nodes \(x_k\) are the roots of \(P_n(x)\text{,}\) found using (6.5.1). In that section, we also derive the formula for the weights or
\begin{equation} w_k = \frac{2}{(1-x_k^2)(P_n'(x_k))^2}\tag{5.4.12} \end{equation}
🔗

Example 5.4.4.

Find the nodes from the roots of the Legendre polynomial and the weights using (5.4.12) for both the \(n=2\) and \(n=3\) cases for Gaussian Quadrature formula.
🔗
Solution.
First, for the \(n=2\) case, we need \(P_2\text{,}\) which can be derived from (6.5.1)
\begin{equation*} P_2(x) = \frac{2(2)-1}{2} x P_1(x) - \frac{2-1}{2} P_0(x) = \frac{3}{2} x (x) - \frac{1}{2} (1) = \frac{3}{2}x^2-\frac{1}{2} \end{equation*}
🔗
The nodes are the root of this or
\begin{equation*} \begin{aligned} x^2 = \frac{1}{3} \qquad \Longrightarrow \quad x = \pm \frac{1}{\sqrt{3}} \end{aligned} \end{equation*}
which is the same as the nodes found above. Next, we’ll find the weights.
🔗
First, we need the derivative so \(P'_2(x) = 3x\text{,}\) then
\begin{equation*} w_0 = \frac{2}{(1-(1/3))(3(1/\sqrt{3}))^2} = 1 \end{equation*}
and similarly (mainly due to symmetry), \(w_1=1\text{.}\)
🔗
For the \(n=3\) case, we start with the 3rd Legendre Polynomial which can be derived as
\begin{equation*} \begin{aligned} P_3(x) \amp = \frac{2(3)-1}{3} x P_2(x) - \frac{3-1}{3} P_1(x) \\ \amp = \frac{5}{3} x \left(\frac{3}{2}x^2-\frac{1}{2}\right) - \frac{2}{3} x \\ \amp = \frac{5}{2} x^3 -\frac{3}{2} x = \frac{1}{2} \left(5x^3-3x^2 \right) \end{aligned} \end{equation*}
🔗
The quadrature nodes are the roots of this polynomial or
\begin{equation*} x = 0, \pm \sqrt{\frac{3}{5}}. \end{equation*}
🔗
To find the weights, we first need the derivative or
\begin{equation*} P'_3 (x) = \frac{1}{2} \left(15 x^2 - 6\right) \end{equation*}
and then
\begin{equation*} \begin{aligned} w_0 \amp = \frac{2}{(1-(\sqrt{3/5})^2)((15(\sqrt{3/5})^2 -6)/2)^2} \\ \amp = \frac{2}{(2/5)(3)^2} = \frac{5}{9} \\ w_1 \amp = \frac{2}{(1-0^2)((-3/2)^2)} = \frac{8}{9} \\ w_2 \amp = w_0 = \frac{5}{9} \end{aligned} \end{equation*}
🔗
🔗
🔗
Let’s look at a higher-order Gauss-Quadrature formula
🔗

Example 5.4.5.

Find the Gauss-Quadrature formula for \(n=7\) and use it to find an estimate for \(\int_0^1 e^{-x^2} \, dx\text{.}\)
🔗
Solution.
First, we need, \(P_7\text{,}\) the seventh Legendre Polynomial. Again, using (6.5.1), we can generate this polynomial. Skipping the details, we get
\begin{equation*} P_7(x) = \frac{1}{16} \left( 429x^7 - 693 x^5 + 315 x^3-35x\right) \end{equation*}
and we seek the roots. Note that one of the roots is \(x=0\text{,}\) which is always the fact for odd Legendre polynomials. Use \(z=x^2\text{,}\) we can reduce the rest to a cubic polynomial and generate the exact solutions, however they are quite complicated. Instead, we will generate the roots numerically with Newton’s method to high precision.
🔗
It’s helpful to have the graph of \(P_7\) to estimate the roots.
🔗
and as expected, the roots are symmetric about \(x=0\text{.}\) If we use Newton’s method as described in Section 3.4, we can find all 7 roots. Note though that they are symmetric, so we only need to find 3 of them. Also the weights are found from (5.4.12) and again there is symmetry to these. The following table shows 20 decimals of each of these.
🔗
Table 5.4.6. Nodes and Weights for \(n=7\) Gaussian Quadrature
🔗
\(i\) \(x_i\) \(w_i\)
1 \(-0.94910791234275852453\) \(0.12948496616886969327\)
2 \(-0.74153118559939443986\) \(0.27970539148927666790\)
3 \(-0.40584515137739716691\) \(0.38183005050511894495\)
4 \(0.00000000000000000000\) \(0.41795918367346938776\)
5 \(0.40584515137739716691\) \(0.38183005050511894495\)
6 \(0.74153118559939443986\) \(0.27970539148927666790\)
7 \(0.94910791234275852453\) \(0.12948496616886969327\)
Since the interval of integration is not \([-1,1]\text{,}\) we perform the transformation in (5.4.9)
\begin{equation*} z = \frac{0+1}{2} + x \frac{1-0}{2} = \frac{1+x}{2} \end{equation*}
Lastly, we use the formula
\begin{equation*} I_7(f) = \sum_{i=1}^7 w_i f((1+x_i)/2) \end{equation*}
with \(f(x) = e^{-x^2}\) to get the estimation: \(I_7(f) = 0.74682413281163830879066365\text{.}\) This is within \(10^{-12}\) of the true answer.
🔗
🔗
🔗
🔗
🔗
PrevTopNext