Prev
Unit 8 · Numerical Differentiation and Integration
Romberg and Adaptive Quadrature
55m100 XP
Next
∫Numerical Differentiation and Integration

Romberg and Adaptive Quadrature

Romberg integration, adaptive Simpson, error estimation and scipy.integrate.quad

QUADPACK, 1983
In the late 1970s Robert Piessens and Elise de Doncker at the University of Leuven, Christoph Überhuber in Vienna and David Kahaner at the US National Bureau of Standards set out to write a definitive library of one-dimensional integration routines in Fortran. Each routine had to take a function, an interval and a requested accuracy, and return the integral and an estimate of its own error — without asking the user for a grid, a step size or a number of points. The library, QUADPACK, appeared with a book in 1983. Its routines bisect intervals where the integrand is difficult, apply 21-point Gauss–Kronrod rules on each piece, estimate errors from the difference between two rules, and even extrapolate their way past endpoint singularities. Four decades later it is still running everywhere. When you type scipy.integrate.quad(f, 0, 1), SciPy calls QUADPACK's routine QAGS — Fortran from 1983 — and it is the reason quad usually just works.
💡Romberg integration
The composite trapezoid rule's error has an expansion in even powers of hh: T(h)=I+c2h2+c4h4+…T(h) = I + c_2h^2 + c_4h^4 + \dots (the Euler–Maclaurin formula). So Richardson extrapolation (lesson 8-1) applies repeatedly. Compute TT with 1,2,4,8,…1, 2, 4, 8, \dots panels — each reusing all previous function values — and build the table Rk,j=Rk,j−1+Rk,j−1−Rk−1,j−14j−1.R_{k,j} = R_{k,j-1} + \dfrac{R_{k,j-1} - R_{k-1,j-1}}{4^j - 1}. Column jj has error O(h2j+2)O(h^{2j+2}); the first extrapolated column is exactly Simpson's rule. Werner Romberg proposed the scheme in 1955. For ∫01e−x2dx\int_0^1 e^{-x^2}dx, the diagonal errors go 6×10−26 \times 10^{-2}, 4×10−44 \times 10^{-4}, 10−510^{-5}, 10−710^{-7}, 3×10−103 \times 10^{-10}, 2×10−132 \times 10^{-13} — thirteen digits from 33 function values.
Derivation
The Romberg recurrence is Richardson, column by column
Suppose the entries of column j−1j - 1 have leading error c h2jc\,h^{2j}, where hh is the panel width of row kk. The entry one row up used panels twice as wide, so its leading error is c (2h)2j=4j c h2jc\,(2h)^{2j} = 4^j\,c\,h^{2j}. Richardson's combination with 2p=4j2^p = 4^j (lesson 8-1) removes it:
Rk,j=4jRk,j−1−Rk−1,j−14j−1=Rk,j−1+Rk,j−1−Rk−1,j−14j−1.R_{k,j} = \frac{4^jR_{k,j-1} - R_{k-1,j-1}}{4^j - 1} = R_{k,j-1} + \frac{R_{k,j-1} - R_{k-1,j-1}}{4^j - 1}.
The first column is Simpson. On one double panel [0,2h][0, 2h], T(h)=h2(f0+2f1+f2)T(h) = \tfrac h2(f_0 + 2f_1 + f_2) and T(2h)=h(f0+f2)T(2h) = h(f_0 + f_2), so
4T(h)−T(2h)3=2h(f0+2f1+f2)−h(f0+f2)3=h3(f0+4f1+f2).\frac{4T(h) - T(2h)}{3} = \frac{2h(f_0 + 2f_1 + f_2) - h(f_0 + f_2)}{3} = \frac h3(f_0 + 4f_1 + f_2).
Richardson discovers Simpson's rule by itself. Each further column cancels the next power of h2h^2, which the Euler–Maclaurin formula guarantees is all the trapezoid rule's error contains for a smooth integrand.
Error estimates for free. Any two approximations of different accuracy give an estimate of the error of the worse one: their difference. Romberg's table, Simpson on one panel versus two half panels, a Gauss rule versus its Kronrod extension — all use the same idea. It is the foundation of every automatic integrator, and of step-size control in ODE solvers (Unit 9).
💡Adaptive Simpson
Integrate a panel with Simpson's rule (S1S_1) and again as two half panels (S2S_2). Because Simpson's error scales like h4h^4, S2S_2's error is about (S2−S1)/15(S_2 - S_1)/15. If ∣S2−S1∣<15 tol|S_2 - S_1| < 15\,\text{tol}, accept S2+(S2−S1)/15S_2 + (S_2 - S_1)/15 — the extrapolated value, at no extra cost. Otherwise apply the same test to each half with tolerance tol/2\text{tol}/2. The recursion spends evaluations only where the test fails. For ∫01x dx\int_0^1\sqrt{x}\,dx — infinite slope at 0, like the displacement near a crack tip — adaptive Simpson reaches an error of 3×10−143 \times 10^{-14} with 985 evaluations, nearly all of them crowded near x=0x = 0. Uniform Simpson with 65,537 evaluations is still off by 5×10−95 \times 10^{-9}.
Derivation
Where the 15 comes from
Let the panel have width HH. One Simpson rule on it has error I−S1≈K H5I - S_1 \approx K\,H^5, with KK proportional to f(4)f^{(4)} (lesson 8-2). The two half panels have width H/2H/2, so each has error K(H/2)5=KH5/32K(H/2)^5 = KH^5/32, and together
I−S2≈2×KH532=KH516=I−S116.I - S_2 \approx 2 \times \frac{KH^5}{32} = \frac{KH^5}{16} = \frac{I - S_1}{16}.
We do not know II, but we know the difference of the two results:
S2−S1=(I−S1)−(I−S2)≈16(I−S2)−(I−S2)=15 (I−S2).S_2 - S_1 = (I - S_1) - (I - S_2) \approx 16(I - S_2) - (I - S_2) = 15\,(I - S_2).
So the error of the better estimate is about (S2−S1)/15(S_2 - S_1)/15. Testing it against the tolerance is the acceptance test, and adding it to S2S_2 is a free Richardson step — the corrected value is exact for quintic polynomials. The argument assumed f(4)f^{(4)} roughly constant over the panel, which is true for small panels on smooth functions and false across a spike: that is how adaptive routines are fooled.
I wish to God these calculations had been executed by steam
– Charles Babbage (1821), recalled in Memoir of the Life and Labours of the Late Charles Babbage
⚠How adaptive integrators are fooled
An adaptive routine only knows the function at the points it sampled. A narrow spike that falls between all of its first samples is invisible: the two rules agree, the error estimate is small, and the answer is confidently wrong. Help it: pass known trouble spots with points=[...] (for a step at 0.3, quad needs 42 evaluations instead of 357), split the interval at discontinuities, transform away singularities, and use the special weights for oscillatory integrands (weight='cos', wvar=50 integrates e−xcos⁡50xe^{-x}\cos 50x over [0,∞)[0, \infty) to twelve digits). And always read the error estimate quad returns.
💡Where you will use this
Orbital Mechanics: averaged perturbation integrals and probability-of-collision integrals. Structures and Geotechnics: influence-line integrals, stress integrals around singular points, Boussinesq-type load integrals. Propulsion: integrating heat flux and pressure distributions over nozzle contours. Statistics and reliability: probabilities as integrals of densities with long tails over infinite intervals — quad handles ±∞\pm\infty by a change of variables.
📚Further reading
  • Robert Piessens, Elise de Doncker-Kapenga, Christoph Überhuber and David Kahaner, QUADPACK: A Subroutine Package for Automatic Integration (Springer, 1983) — the book behind scipy.integrate.quad.
  • Walter Gander and Walter Gautschi, 'Adaptive quadrature — revisited', BIT Numerical Mathematics 40(1), 84–101 (2000) — adaptive Simpson and Gauss–Lobatto, their error estimates and how to stop them.
  • Philip Davis and Philip Rabinowitz, Methods of Numerical Integration (2nd ed., Academic Press, 1984) — Romberg and the other extrapolation methods, and the theory of automatic integration.
↓But there is a problem
Romberg and adaptive Simpson are clever ways of using equally spaced points within each panel. Every rule so far fixed the nodes first and then chose the weights. In 1814 Gauss asked the obvious question nobody had asked: if the nodes are free too, where should they go? With nn points chosen freely, he found, you can integrate polynomials of degree 2n−12n - 1 exactly — twice what equally spaced points achieve.