26 Numerical Quadrature
Approximating \(I = \int_a^b f(x)\,dx\) via weighted sums: \(\sum_{i=1}^n w_i f(x_i)\).
A quadrature rule approximates an integral by function values: \[ \begin{align} \int_a^b f(x)\,dx \approx \sum_{i=1}^p w_i f(x_i). \end{align} \] The nodes \(x_i\) choose where \(f\) is sampled, and the weights \(w_i\) choose how those samples are averaged.
A quadrature rule has degree of exactness \(m\) if it integrates every polynomial of degree at most \(m\) exactly, and fails for at least one polynomial of degree \(m+1\).
(Interpolation view) Newton-Cotes rules interpolate \(f\) at selected nodes and integrate the interpolating polynomial. The quadrature error is therefore controlled by interpolation error and smoothness of \(f\).
26.1 Newton-Cotes Rules
Uses equally-spaced nodes. Higher-order rules are prone to Runge’s phenomenon.
On \(n\) panels of width \(h = (b-a)/n\): \[ \begin{align} T_n = h \left[ \frac{f(a) + f(b)}{2} + \sum_{i=1}^{n-1} f(a + ih) \right]. \end{align} \]
The one-panel trapezoidal rule is exact for polynomials of degree at most \(1\).
On \([a,b]\), the trapezoidal rule integrates the line through \((a,f(a))\) and \((b,f(b))\). If \(f\) itself is constant or linear, this interpolant equals \(f\), so the integral is exact.
If \(f\) has two continuous derivatives on \([a,b]\), then the composite trapezoidal rule satisfies \[ \begin{align} I-T_n = O(h^2). \end{align} \] Thus doubling \(n\) reduces the leading error by a factor of about \(4\).
Derive the one-panel trapezoidal rule by integrating the linear interpolant through \((a,f(a))\) and \((b,f(b))\).
Use Taylor expansion on one panel to show the local error is \(O(h^3)\).
Sum over \(O(1/h)\) panels to prove the global \(O(h^2)\) statement in the result above.
Requires \(n\) to be even. Uses piecewise quadratic interpolation: \[ \begin{align} S_n = \frac{h}{3} \left[ f(a) + f(b) + 4\sum_{\text{odd}\;i} f(x_i) + 2\sum_{\text{even}\;i} f(x_i) \right]. \end{align} \]
Simpson’s rule is exact for polynomials of degree at most \(3\). For sufficiently smooth \(f\), the composite Simpson rule has global error \(O(h^4)\).
(Degree effect) Simpson’s rule uses a quadratic interpolant, but symmetry cancels the cubic error term. This is why its exactness degree is \(3\), not merely \(2\).
Integrate the quadratic interpolant through \(x_0,x_1,x_2\) to derive the Simpson weights \(1,4,1\).
Verify exactness for \(1,x,x^2,x^3\) on a symmetric panel.
Use the exactness degree to explain why the leading error involves the fourth derivative.
Compare the observed error reduction after halving \(h\) with the result above.
26.2 Gaussian Quadrature
Optimizes both nodes and weights to achieve the highest possible degree of exactness.
Integrating over \([-1, 1]\) with \(p\) points: \[ \begin{align} \int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^p w_i f(x_i). \end{align} \] The nodes are the roots of the degree-\(p\) Legendre polynomial \(P_p(x)\), and the rule is exact for all polynomials of degree at most \(2p-1\).
(Gauss efficiency) Gaussian quadrature spends nodes where orthogonality says they are most useful. For smooth, non-periodic integrands, a small number of Gauss points often beats many equally spaced Newton-Cotes points.
For \(p=2\), use the roots of \(P_2(x)=(3x^2-1)/2\) to find the nodes.
Solve for weights that integrate \(1\) and \(x^2\) exactly on \([-1,1]\).
Verify that the resulting rule also integrates \(x\) and \(x^3\) exactly by symmetry.
Explain why exactness through degree \(2p-1\) is the maximum possible with \(2p\) free parameters.
26.3 Richardson Extrapolation and Romberg Integration
Given an \(O(h^2)\) estimate \(A(h)\), build an \(O(h^4)\) estimate by: \[ \begin{align} A_{new} = \frac{4A(h/2) - A(h)}{3}. \end{align} \]
Let \(A(h)\) be an \(O(h^2)\) approximation of the true value \(I\). Using a Taylor expansion in \(h\): \[ \begin{align} A(h) = I + C_1 h^2 + C_2 h^4 + O(h^6). \end{align} \] Halving the step size gives: \[ \begin{align} A(h/2) &= I + C_1 (h/2)^2 + C_2 (h/2)^4 + O(h^6) \\ &= I + \frac{1}{4} C_1 h^2 + \frac{1}{16} C_2 h^4 + O(h^6). \end{align} \] To eliminate the \(O(h^2)\) error term, multiply \(A(h/2)\) by 4 and subtract \(A(h)\): \[ \begin{align} 4A(h/2) - A(h) &= (4I + C_1 h^2 + \frac{1}{4} C_2 h^4) - (I + C_1 h^2 + C_2 h^4) + O(h^6) \\ &= 3I - \frac{3}{4} C_2 h^4 + O(h^6). \end{align} \] Dividing by 3 yields the \(O(h^4)\) estimate: \[ \begin{align} I = \frac{4A(h/2) - A(h)}{3} + O(h^4). \end{align} \]
Iterative application of Richardson extrapolation to the trapezoidal rule.
(Romberg idea) Romberg integration repeatedly cancels the leading even powers in the trapezoidal error expansion. It is effective when \(f\) is smooth enough for that expansion to be accurate.
26.4 Adaptive Quadrature
(Adaptive quadrature) Adaptive quadrature focuses function evaluations where \(f\) is difficult and uses large panels where \(f\) is smooth. The usual mechanism compares one large panel with two smaller panels; if the difference is above tolerance \(\varepsilon\), the method recurses.
(Nonsmooth integrands) Algebraic error rates assume enough derivatives. Corners, jumps, and endpoint singularities reduce the observed order. Adaptive quadrature helps by isolating nonsmooth regions, but it cannot restore smoothness.
(Monte Carlo contrast) Deterministic quadrature is efficient in low dimension for smooth functions. Monte Carlo convergence is only \(O(N^{-1/2})\), but its rate is largely dimension-independent, which makes it useful for high-dimensional integrals.
26.5 Periodic Functions and Spectral Accuracy
For smooth periodic functions on \([0,2\pi]\), the trapezoidal rule can converge faster than any algebraic power of \(h\). For analytic periodic functions, the convergence is exponential.
Apply the trapezoidal rule to \(e^{ikx}\) on \([0,2\pi]\) and determine when the discrete sum is exactly zero.
Use a Fourier series argument to explain why smooth periodic functions are integrated unusually well.
Compare errors for \(\sin(x)\), \(\exp(\cos x)\), and \(|x-\pi|\) as \(n\) increases.
26.6 Exercises
Doubling nodes in Trapezoidal rule reduces error by \(4\times\). Use the result above to predict Simpson’s reduction factor.
Derive the \(p=2\) Gauss points for \([-1, 1]\) using the result above.
Use Richardson extrapolation from the result above to improve \(T(h)\) and \(T(h/2)\) to get Simpson’s rule.
Integration of \(\sin(x)\) on \([0, 2\pi]\) with \(n=4\) nodes. Is the error \(0\)? Why?
Compare composite Simpson on \(f(x)=|x-1/2|\) with a smooth quartic. Estimate the observed order.