8.5 Gaussian quadrature

Newton–Cotes formulas fix the nodes in advance. If the nodes are chosen as well as the weights, an \(n\)-point formula \(\sum _{i=1}^nc_if(x_i)\) has \(2n\) free parameters, and can be made exact for all polynomials of degree up to \(2n-1\). On \([-1,1]\) the optimal nodes are the roots of the Legendre polynomial \(P_n\).

Theorem 8.5.1. Let \(x_1,\dots ,x_n\) be the roots of the \(n\)th Legendre polynomial \(P_n\) and \(c_i = \int _{-1}^1\prod _{j\neq i}\frac {x-x_j}{x_i-x_j}\,dx\). Then \(\int _{-1}^1P(x)\,dx = \sum _{i=1}^nc_iP(x_i)\) for every polynomial \(P\) of degree less than \(2n\).

Proof. For degree less than \(n\), the formula integrates the interpolating polynomial at the \(n\) nodes, which is \(P\) itself. For degree less than \(2n\), divide by \(P_n\): \(P = QP_n+R\) with \(Q\) and \(R\) of degree less than \(n\). The Legendre polynomial \(P_n\) is orthogonal on \([-1,1]\) to every polynomial of lower degree, so \(\int _{-1}^1QP_n = 0\) and \(\int P = \int R\). At the nodes \(P_n(x_i) = 0\), so \(P(x_i) = R(x_i)\) and \(\sum c_iP(x_i) = \sum c_iR(x_i) = \int R\), by the first case. □

For two and three points the nodes and weights are

\(n\) nodes \(x_i\) weights \(c_i\)
\(2\) \(\pm \frac {1}{\sqrt 3} = \pm 0.5773502692\) \(1,\ 1\)
\(3\) \(0,\ \pm \sqrt {0.6} = \pm 0.7745966692\) \(\frac 89,\ \frac 59,\ \frac 59\)

An integral over \([a,b]\) is transformed to \([-1,1]\) by \(x = \frac {(b-a)t+(a+b)}{2}\): \[\int _a^bf(x)\,dx = \frac {b-a}{2}\int _{-1}^1f\!\left (\frac {(b-a)t+(a+b)}{2}\right )dt .\]

Example 8.5.2. Approximate \(\int _1^{1.5}e^{-x^2}\,dx = 0.1093643\) by Gaussian quadrature with \(n = 2\) and \(n = 3\), and compare with Simpson’s rule.

Solution. Here \(x = \frac {t+5}{4}\) and \(dx = \frac {dt}{4}\), so \(\int _1^{1.5}e^{-x^2}dx = \frac 14\int _{-1}^1e^{-(t+5)^2/16}\,dt\). With \(n = 2\), \[\frac 14\left [e^{-(5-0.5773503)^2/16}+e^{-(5+0.5773503)^2/16}\right ] = 0.1094003,\] error \(3.6\times 10^{-5}\). With \(n = 3\), using the nodes \(0,\pm 0.7745967\) and weights \(\frac 89,\frac 59,\frac 59\), the result is \(0.1093642\), error \(6.5\times 10^{-8}\). Simpson’s rule, with the same three function values’ worth of work at \(1, 1.25, 1.5\), gives \(0.1093104\), error \(5.4\times 10^{-5}\). For smooth integrands Gaussian quadrature is the most accurate rule for a given number of evaluations.

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.