Burden & Faires, 9th Ed. - §4.8–4.9
Chapter 10: Adaptive & Gaussian Quadrature
Romberg integration · Adaptive Simpson's rule · Gaussian quadrature · Optimal node placement
Chapter 9's rules all fix the nodes first (equally spaced) and then integrate whatever interpolant those nodes produce. This chapter asks two sharper questions: can we squeeze near-machine-precision accuracy out of the trapezoidal rule without ever using a dangerously small $h$ (Romberg)? Can we spend more effort exactly where the function actually needs it instead of everywhere uniformly (adaptive quadrature)? And, most radically, what if the nodes themselves aren't required to be equally spaced at all - what's the best possible place to put them (Gaussian quadrature)?
Section 4.8
Romberg Integration
Richardson Extrapolation Applied to the Trapezoidal Rule
The composite trapezoidal rule's error is a series in even powers of $h$: $T(h) = I + K_1h^2+K_2h^4+\cdots$ - exactly the structure Chapter 8's Richardson extrapolation needs. Define $R_{k,1}=$ composite trapezoidal with $2^{k-1}$ subintervals, then recursively
$$R_{k,j} = R_{k,j-1} + \frac{R_{k,j-1}-R_{k-1,j-1}}{4^{j-1}-1}$$
Column $j$ is accurate to $O(h^{2j})$ - column 2 matches Simpson's rule exactly, column 3 matches Boole's rule, and so on, all generated from nothing but repeated trapezoidal evaluations and this one recursion.
📘 Example 10.1 - Romberg Table for $\int_0^\pi \sin x\,dx$
Build the Romberg table using $R_{1,1},\ldots,R_{5,1}$ from composite trapezoidal with $1,2,4,8,16$ subintervals.
Solution
$R_{1,1}=0$ is a single (useless) trapezoid; by $R_{5,5}$, the answer is correct to 9 decimal places - reached using only trapezoidal-rule evaluations that individually range from poor to mediocre, with all the real accuracy manufactured by the extrapolation itself.
| $k$ | $R_{k,1}$ | $R_{k,2}$ | $R_{k,3}$ | $R_{k,4}$ | $R_{k,5}$ |
|---|---|---|---|---|---|
| 1 | 0.0000000 | ||||
| 2 | 1.5707963 | 2.0943951 | |||
| 3 | 1.8961189 | 2.0045597 | 1.9985707 | ||
| 4 | 1.9742316 | 2.0002692 | 1.9999831 | 2.0000056 | |
| 5 | 1.9935703 | 2.0000166 | 1.9999998 | 2.0000000 | 1.9999999999 |
Section 4.8 continued
Adaptive Quadrature
The Idea - Spend Effort Where the Function Needs It
Composite rules with a uniform $h$ are wasteful when $f$ is nearly flat over most of $[a,b]$ but has a sharp feature somewhere - a uniform grid fine enough for the sharp region wastes huge effort on the flat region. Adaptive quadrature instead estimates the local error on each subinterval and only subdivides further where that estimate exceeds tolerance.
Adaptive Simpson's Algorithm
On interval $[a,b]$ with midpoint $c$: compute Simpson's estimate $S[a,b]$ over the whole interval, and $S[a,c]+S[c,b]$ using two half-width panels. Since Simpson's error is $O(h^5)$ per panel, halving $h$ should shrink the error by a factor of $\approx16$, so
$$|S[a,c]+S[c,b] - S[a,b]| \approx 15\times(\text{true error of } S[a,c]+S[c,b])$$
gives a cheap, computable error estimate. If it exceeds tolerance, recurse on $[a,c]$ and $[c,b]$ separately (each with half the original tolerance); otherwise accept $S[a,c]+S[c,b]$ as the answer on this panel.
$f(x) = e^{-100(x-0.5)^2}$ - a sharp Gaussian spike centered at $x=0.5$, nearly zero elsewhere on $[0,1]$. Adaptive Simpson's Rule automatically clusters small panels (bars) tightly around the spike and leaves the flat tails covered with just one or two wide panels - no uniform grid could do this without wasting enormous effort on the flat regions.
Section 4.9
Gaussian Quadrature
A Different Question Entirely
Every method so far fixed the nodes (equally spaced) and only chose weights. Gaussian quadrature chooses both the $n$ nodes $x_1,\ldots,x_n$ and the $n$ weights $c_1,\ldots,c_n$ to maximize the degree of polynomial the rule integrates exactly. With $2n$ free parameters, the best possible degree of precision is $2n-1$ - and it is achieved exactly when the nodes are the roots of the $n$th-degree Legendre polynomial on $[-1,1]$.
Gauss-Legendre Quadrature on $[-1,1]$
$$\int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n c_i f(x_i)$$
| $n$ | Nodes $x_i$ | Weights $c_i$ | Exact for degree $\le$ |
|---|---|---|---|
| 2 | $\pm1/\sqrt3 \approx \pm0.57735$ | $1,\,1$ | 3 |
| 3 | $0,\;\pm\sqrt{3/5}\approx\pm0.77460$ | $8/9,\;5/9,\;5/9$ | 5 |
📘 Example 10.2 - Gaussian Quadrature Beats Simpson at Equal Cost
Approximate $\int_{-1}^1 e^x\,dx$ (true value $e-e^{-1}=2.350402$) using 2-point and 3-point Gauss-Legendre, and compare to Simpson's Rule using the same number of function evaluations.
Solution
$$\text{2-point:} \quad e^{-1/\sqrt3}+e^{1/\sqrt3} = 2.342696, \quad \text{error}=7.706\times10^{-3}$$
$$\text{3-point:} \quad \tfrac89 e^0 + \tfrac59\left(e^{-\sqrt{3/5}}+e^{\sqrt{3/5}}\right) = 2.350337, \quad \text{error}=6.546\times10^{-5}$$
3-point Gauss-Legendre (3 evaluations) beats even 5-point Simpson (5 evaluations) by more than a factor of 10 - because Gauss also optimizes where to evaluate $f$, not just how to weight fixed, equally spaced points.
| Method (# evaluations) | Error |
|---|---|
| Simpson, $n=2$ (3 points) | $1.165\times10^{-2}$ |
| 3-point Gauss-Legendre | $6.546\times10^{-5}$ |
| Simpson, $n=4$ (5 points) | $7.924\times10^{-4}$ |
The Trade Gaussian Quadrature Makes
Gauss-Legendre's nodes are irrational and not reusable between different $n$ (unlike Romberg's nested trapezoidal grids) and are harder to look up or hand-compute - the accuracy gain comes at the cost of needing precomputed tables (or a root-finder from Chapter 4, applied to Legendre polynomials) rather than simple equally spaced arithmetic. It is the method of choice when $f$ is smooth and evaluations are expensive; adaptive Simpson is usually preferred when $f$'s behavior is unknown or irregular in advance.
Practice
Practice Problem
Try this problem yourself before expanding the solution below.
✏️ Practice Problem 10.1 - Romberg Table for $\int_0^1 e^{-x^2}\,dx$
Build the Romberg table $R_{1,1}$ through $R_{3,3}$ for $\displaystyle\int_0^1 e^{-x^2}\,dx$ (no elementary antiderivative exists, but the true value is $0.746824$ to six decimals). Use trapezoidal estimates with $n=1,2,4$ as the first column.
Solution
First column (composite trapezoidal, $n=1,2,4$):
$$R_{1,1}=0.683940, \qquad R_{2,1}=0.731370, \qquad R_{3,1}=0.742984.$$
Applying $R_{i,j} = R_{i,j-1} + \dfrac{R_{i,j-1}-R_{i-1,j-1}}{4^{j-1}-1}$ down the diagonals:
$R_{3,3}=0.746834$, within $10^{-5}$ of the true value $0.746824$, achieved from trapezoidal estimates that individually had errors as large as $0.063$ - Romberg's Richardson-extrapolation cascade is doing essentially all of the work.
| $n$ | $j=1$ | $j=2$ | $j=3$ |
|---|---|---|---|
| 1 | 0.683940 | ||
| 2 | 0.731370 | 0.747180 | |
| 3 | 0.742984 | 0.746855 | 0.746834 |
Connections Across the Course
- Romberg integration is literally Chapter 8's Richardson extrapolation, applied to a different base sequence.
- Adaptive quadrature's error-doubling-check echoes the a posteriori error estimates of Chapter 3's fixed-point bounds - estimate error from two successive approximations rather than a worst-case bound.
- Finding Gauss-Legendre nodes as roots of Legendre polynomials is itself a root-finding problem, solvable with Newton's Method from Chapter 4.