5.3 Gaussian quadrature#

Newton-Cotes formulas were obtained by integrating interpolating polynomials with equally-spaced nodes. The equal spacing is convenient in deriving simple expressions for the composite rules. However, this placement of nodes is not necessarily the optimal placement. For example, the trapezoidal rule approximates the integral by integrating a linear function that joins the endpoints of the function. The fact that this is not the optimal choice can be seen by sketching a simple parabola.

The idea of Gaussian quadrature is the following: in the numerical quadrature rule

\[\begin{equation*} \int_{a}^{b}f(x)dx\approx\sum_{i=1}^{n}w_{i}f(x_{i}) \end{equation*}\]

choose \(x_{i}\) and \(w_{i}\) in such a way that the quadrature rule has the highest possible accuracy. Note that unlike Newton-Cotes formulas where we started labeling the nodes with \(x_0\), in the Gaussian quadrature the first node is \(x_1\). This difference in notation is common in the literature, and each choice makes the subsequent equations in the corresponding theory easier to read.

Example 46

Let \((a,b)=(-1,1),\) and \(n=2.\) Find the ``best’’ \(x_{i}\) and \(w_{i}\).

There are four parameters to determine: \(x_{1},x_{2},w_{1},w_{2}.\) We need four constraints. Let’s require the rule to integrate the following functions exactly: \(f(x)=1,f(x)=x,f(x)=x^{2}\), and \(f(x)=x^{3}\).

If the rule integrates \(f(x)=1\) exactly, then \(\int_{-1}^{1}dx=\sum_{i=1}^{2}w_{i}\), i.e., \(w_{1}+w_{2}=2.\) If the rule integrates \(f(x)=x\) exactly, then \(\int_{-1}^{1}xdx=\sum_{i=1}^{2}w_{i}x_{i}\), i.e., \(w_{1}x_{1}+w_{2}x_{2}=0.\) Continuing this for \(f(x)=x^{2}\), and \(f(x)=x^{3},\) we obtain the following equations:

\[\begin{align*} w_{1}+w_{2} & =2\\ w_{1}x_{1}+w_{2}x_{2} & =0\\ w_{1}x_{1}^{2}+w_{2}x_{2}^{2} & =\frac{2}{3}\\ w_{1}x_{1}^{3}+w_{2}x_{2}^{3} & =0. \end{align*}\]

Solving the equations gives: \(w_{1}=w_{2}=1,x_{1}=\frac{-\sqrt{3}}{3},x_{2}=\frac{\sqrt{3}}{3}\). Therefore the quadrature rule is:

\[\begin{equation*} \int_{-1}^{1}f(x)dx\approx f\left(\frac{-\sqrt{3}}{3}\right)+f\left(\frac{\sqrt{3}}{3}\right). \end{equation*}\]

Observe that:

  • The two nodes are not equally spaced on \((-1,1)\).

  • The accuracy of the rule is three, and it uses only two nodes. Recall that the accuracy of Simpson’s rule is also three but it uses three nodes. In general Gaussian quadrature gives a degree of accuracy of \(2n-1\) using only \(n\) nodes.

We were able to solve for the nodes and weights in the simple example above; however, as the number of nodes increases, the resulting non-linear system of equations will be very difficult to solve. There is an alternative approach using the theory of orthogonal polynomials, a topic we will discuss in detail later. Here, we will use a particular set of orthogonal polynomials \(\{L_{0}(x),L_{1}(x),...,L_{n}(x),...\}\) called Legendre polynomials. We will give a definition of these polynomials in the next chapter. For this discussion, we just need the following properties of these polynomials:

  • \(L_{n}(x)\) is a monic polynomial of degree \(n\) for each \(n\).

  • \(\int_{-1}^{1}P(x)L_{n}(x)dx=0\) for any polynomial \(P(x)\) of degree less than \(n\).

The first few Legendre polynomials are

\[\begin{align*} L_{0}(x) =1,L_{1}(x)=x,L_{2}(x)=x^{2}-\frac{1}{3}, L_{3}(x) =x^{3}-\frac{3}{5}x,L_{4}(x)=x^{4}-\frac{6}{7}x^{2}+\frac{3}{35}. \end{align*}\]

How do these polynomials help with finding the nodes and the weights of the Gaussian quadrature rule? The answer is short: the roots of the Legendre polynomials are the nodes of the quadrature rule!

To summarize, the Gauss-Legendre quadrature rule for the integral of \(f\) over \((-1,1)\) is

\[\begin{equation*} \int_{-1}^1 f(x)dx = \sum_{i=1}^n w_i f(x_i) \end{equation*}\]

where \(x_{1},x_{2},...,x_{n}\) are the roots of the \(n\)th Legendre polynomial, and the weights are computed using the following theorem.

Theorem 31

Suppose that \(x_{1},x_{2},...,x_{n}\) are the roots of the \(n\)th Legendre polynomial \(L_{n}(x)\) and the weights are given by

\[\begin{equation*} w_{i}=\int_{-1}^{1}\prod_{j=1,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}dx. \end{equation*}\]

Then the Gauss-Legendre quadrature rule has degree of accuracy \(2n-1\). In other words, if \(P(x)\) is any polynomial of degree less than or equal to \(2n-1\), then

\[\begin{equation*} \int_{-1}^{1}P(x)dx=\sum_{i=1}^{n}w_{i}P(x_{i}). \end{equation*}\]

Proof. Let’s start with a polynomial \(P(x)\) with degree less than \(n\). Construct the Lagrange interpolant for \(P(x),\) using the nodes as \(x_{1},...,x_{n}\)

\[\begin{equation*} P(x)=\sum_{i=1}^{n}P(x_{i})l_{i}(x)=\sum_{i=1}^{n}\prod_{j=1,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}P(x_{i}). \end{equation*}\]

There is no error term above because the error term depends on the \(n\)th derivative of \(P(x)\), but \(P(x)\) is a polynomial of degree less than \(n\), so that derivative is zero. Integrate both sides to get

\[\begin{align*} \int_{-1}^{1}P(x)dx =\int_{-1}^{1}\left[\sum_{i=1}^{n}\prod_{j=1,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}P(x_{i})\right]dx & =\int_{-1}^{1}\left[\sum_{i=1}^{n}P(x_{i})\prod_{j=1,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}\right]dx\\ & =\sum_{i=1}^{n}\left[\int_{-1}^{1}P(x_{i})\prod_{j=1,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}dx\right]\\ & =\sum_{i=1}^{n}\left[\int_{-1}^{1}\prod_{j=1,j\neq i}^{n}\frac{x-x_{j}}{x_{i}-x_{j}}dx\right]P(x_{i})\\ & =\sum_{i=1}^{n}w_{i}P(x_{i}). \end{align*}\]

Therefore the theorem is correct for polynomials of degree less than \(n\). Now let \(P(x)\) be a polynomial of degree greater than or equal to \(n\), but less than or equal to \(2n-1\). Divide \(P(x)\) by the Legendre polynomial \(L_{n}(x)\) to get

\[\begin{equation*} P(x)=Q(x)L_{n}(x)+R(x). \end{equation*}\]

Note that

\[\begin{equation*} P(x_{i})=Q(x_{i})L_{n}(x_{i})+R(x_{i})=R(x_{i}) \end{equation*}\]

since \(L_{n}(x_{i})=0\) , for \(i=1,2,...,n.\) Also observe the following facts:

  1. \(\int_{-1}^{1}Q(x)L_{n}(x)dx=0,\) since \(Q(x)\) is a polynomial of degree less than \(n,\) and from the second property of Legendre polynomials.

  2. \(\int_{-1}^{1}R(x)dx=\sum_{i=1}^{n}w_{i}R(x_{i}),\) since \(R(x)\) is a polynomial of degree less than \(n,\) and from the first part of the proof.

Putting these facts together we get:

\[\begin{equation*} \int_{-1}^{1}P(x)dx=\int_{-1}^{1}\left[Q(x)L_{n}(x)+R(x)\right]dx=\int_{-1}^{1}R(x)dx=\sum_{i=1}^{n}w_{i}R(x_{i})=\sum_{i=1}^{n}w_{i}P(x_{i}). \end{equation*}\]

Table 1 displays the roots of the Legendre polynomials \(L_2,L_3,L_4,L_5\) and the corresponding weights.

Table 1 Roots of Legendre polynomials \(L_2\) through \(L_5\)#

\(n\)

Roots

Weights

2

\(\frac{1}{\sqrt{3}}=0.5773502692\)

\(1\)

\(-\frac{1}{\sqrt{3}}=-0.5773502692\)

\(1\)

\(3\)

\(-(\frac{3}{5})^{1/2}=-0.7745966692\)

\(\frac{5}{9}=0.5555555556\)

\(0.0\)

\(\frac{8}{9}=0.8888888889\)

\((\frac{3}{5})^{1/2}=0.7745966692\)

\(\frac{5}{9}=0.5555555556\)

\(4\)

\(0.8611363116\)

\(0.3478548451\)

\(0.3399810436\)

\(0.6521451549\)

\(-0.3399810436\)

\(0.6521451549\)

\(-0.8611363116\)

\(0.3478548451\)

\(5\)

\(0.9061798459\)

\(0.2369268850\)

\(0.5384693101\)

\(0.4786286705\)

\(0.0\)

\(0.5688888889\)

\(-0.5384693101\)

\(0.4786286705\)

\(-0.9061798459\)

\(0.2369268850\)

Example 47

Approximate \(\int_{-1}^1 \cos x dx\) using Gauss-Legendre quadrature with \(n=3\) nodes.

Solution.

From Table 1, and using two-digit rounding, we have

\[\begin{equation*} \int_{-1}^1 \cos x dx \approx 0.56 \cos (-0.77) + 0.89 \cos 0 + 0.56 \cos (0.77)=1.69 \end{equation*}\]

and the true solution is \(\sin(1)-\sin(-1)=1.68\).

So far we discussed integrating functions over the interval \((-1,1).\) What if we have a different integration domain? The answer is simple: change of variables! To compute \(\int_{a}^{b}f(x)dx\) for any \(a<b\), we use the following change of variables:

\[\begin{equation*} t=\frac{2x-a-b}{b-a}\Leftrightarrow x=\frac{1}{2}\left[(b-a)t+a+b\right] \end{equation*}\]

With this substitution, we have

\[\begin{equation*} \int_a^b f(x) dx = \frac{b-a}{2} \int_{-1}^1 f \left( \frac{1}{2}\left[(b-a)t+a+b\right] \right)dt. \end{equation*}\]

Now we can approximate the integral on the right-hand side as before.

Example 48

Approximate \(\int_{0.5}^{1}x^x dx\) using Gauss-Legendre quadrature with \(n=2\) nodes.

Solution.

Transform the integral using \(x=\frac{1}{2}\left(0.5t+1.5\right)=\frac{1}{2}(\frac{t}{2}+\frac{3}{2})=\frac{t+3}{4},dx=\frac{dt}{4}\) to get:

\[\begin{equation*} \int_{0.5}^{1}x^x dx=\frac{1}{4}\int_{-1}^{1} \left(\frac{t+3}{4}\right)^{\frac{t+3}{4}}dt. \end{equation*}\]

For \(n=2\)

\[\begin{align*} \frac{1}{4}\int_{-1}^{1} \left(\frac{t+3}{4}\right)^{\frac{t+3}{4}}dt\approx \frac{1}{4} \left[ \left( \frac{1}{4\sqrt{3}}+\frac{3}{4}\right)^{\left( \frac{1}{4\sqrt{3}}+\frac{3}{4}\right)}+\left(- \frac{1}{4\sqrt{3}}+\frac{3}{4}\right)^{\left( -\frac{1}{4\sqrt{3}}+\frac{3}{4}\right)} \right]=0.410759, \end{align*}\]

using six significant digits. We will next use Python for a five-node computation of the integral.

Python code for Gauss-Legendre rule with five nodes#

The following code computes the Gauss-Legendre rule for \(\int_{-1}^1 f(x)dx\) using \(n=5\) nodes. The nodes and weights are from Table 1.

def gauss(f):
    return 0.2369268851*f(-0.9061798459) + 0.2369268851*f(0.9061798459) + \
           0.5688888889*f(0) + 0.4786286705*f(0.5384693101) + \
           0.4786286705*f(-0.5384693101)

Now we compute \(\frac{1}{4}\int_{-1}^{1} \left(\frac{t+3}{4}\right)^{\frac{t+3}{4}}dt\) using the code:

0.25*gauss(lambda t: (t/4+3/4)**(t/4+3/4))

0.41081564812239885

The next theorem is about the error of the Gauss-Legendre rule. Its proof can be found in Atkinson [1989]. The theorem shows, in particular, that the degree of accuracy of the quadrature rule, using \(n\) nodes, is \(2n-1\).

Theorem 32

Let \(f\in C^{2n}[-1,1].\) The error of Gauss-Legendre rule satisfies

\[\begin{equation*} \int_{-1}^{1}f(x)dx-\sum_{i=1}^{n}w_{i}f(x_{i})=\textcolor{blue}{\frac{2^{2n+1}(n!)^{4}}{(2n+1)\left[(2n)!\right]^{2}}}\frac{f^{(2n)}(\xi)}{(2n)!} \end{equation*}\]

for some \(\xi\in(-1,1)\).

Using Stirling’s formula \(n! \sim e^{-n} n^n (2\pi n)^{1/2}\), where the symbol \(\sim\) means the ratio of the two sides converges to 1 as \(n\rightarrow \infty\), it can be shown that

\[\begin{equation*} \textcolor{blue}{\frac{2^{2n+1}(n!)^{4}}{(2n+1)\left[(2n)!\right]^{2}}}\sim \frac{\pi}{4^{n}}. \end{equation*}\]

This means the error of Gauss-Legendre rule decays at an exponential rate of \(1/4^{n}\) as opposed to, for example, the polynomial rate of \(1/n^{4}\) for composite Simpson’s rule.

Exercise 5.3-1

Prove that the sum of the weights in Gauss-Legendre quadrature is 2, for any \(n\).

Exercise 5.3-2

Take \(f(x)=1\) in \(\int_{-1}^1 f(x)dx=\sum_{i=1}^n w_i f(x_i)\) (the error term is zero because the derivative of \(f\) is zero) to get \(2=\sum_{i=1}^n w_i\).

Exercise 5.3-3

Approximate \(\int_{1}^{1.5}x^{2}\log xdx\) using Gauss-Legendre rule with \(n=2\) and \(n=3\). Compare the approximations to the exact value of the integral.

Atk89

K.E. Atkinson. An Introduction to Numerical Analysis. John Wiley & Sons, 2nd edition, 1989.