Mathematical foundations of Composite Trapezoidal & Romberg Quadrature.
Numerical quadrature—the algorithmic approximation of definite integrals—is essential across scientific computing, from evaluating quantum expectation values to solving hydrodynamic conservation laws. The composite trapezoidal rule is the simplest and most robust Newton-Cotes formula, replacing the integrand with piecewise linear secants. While its basic second-order O(h²) convergence is modest, the Euler-Maclaurin summation formula reveals that its truncation error expansion proceeds exclusively in even powers of the mesh spacing h. Werner Romberg seized upon this mathematical structure to devise Romberg extrapolation: by taking successive linear combinations of trapezoidal sums with halved step sizes, one systematically cancels error terms, boosting precision to O(h⁴), O(h⁶), O(h⁸), and beyond. In this laboratory, we implement the composite trapezoidal rule, verify its O(h²) scaling, and construct the Romberg tableau.
Viva voce questions & analytical theory
How is the composite trapezoidal rule derived geometrically?
The interval [a, b] is divided into N sub-panels of width h = (b - a)/N. Over each sub-panel, the curve y = f(x) is approximated by a linear secant line segment, forming a trapezoid of area (h/2)[f(x_i) + f(x_{i+1})]. Summing these yields T(h) = h [f(a)/2 + ∑ f(x_i) + f(b)/2].
What is the Euler-Maclaurin summation formula and what does it reveal about trapezoidal error?
The Euler-Maclaurin formula states that T(h) - I = c₂ h² + c₄ h⁴ + c₆ h⁶ + ... where c_{2k} = [B_{2k} / (2k)!] · [f^{(2k-1)}(b) - f^{(2k-1)}(a)] and B_{2k} are Bernoulli numbers. This proves that the error series contains only even powers of h.
How does Romberg integration exploit the Euler-Maclaurin expansion?
Because the error contains only even powers of h (h², h⁴, h⁶, ...), Richardson extrapolation can eliminate the leading O(h²) error by taking R(1, 1) = [4 T(h/2) - T(h)] / 3, yielding Simpson's rule O(h⁴). Repeating this process eliminates O(h⁴), O(h⁶), etc.
Why does the trapezoidal rule exhibit miraculous exponential convergence for smooth periodic functions?
If f(x) is smooth and periodic on [a, b] with all derivatives matching at endpoints (f^{(k)}(a) = f^{(k)}(b) for all k), all Euler-Maclaurin coefficients c_{2k} vanish identically. The error then decays exponentially fast as O(e^{-c/h}) by the Poisson summation formula.
How does Romberg integration reuse function evaluations efficiently?
When halving the step size from h to h/2, all existing function evaluation points x_i are retained. The new trapezoidal sum is computed as T(h/2) = 0.5 T(h) + (h/2) ∑_{odd} f(x_j), requiring only N new function evaluations rather than 2N.
What is the algebraic precision of the composite trapezoidal rule?
The trapezoidal rule is exact for polynomials of degree ≤ 1 (affine linear functions). Its degree of algebraic precision is d = 1.
What happens if Romberg extrapolation is applied to an integrand with a boundary derivative singularity (e.g., √x)?
If f(x) lacks smooth derivatives at boundaries, the Euler-Maclaurin series does not hold with integer powers of h (fractional powers like h^{1.5} appear). Standard Romberg weights (4^m - 1) fail, degrading convergence unless singularity-adapted extrapolation or variable transformations are used.
Algorithm precautions and stability rules
- Reuse previously evaluated function values when halving panel spacing to avoid redundant floating point operations.
- Sum function evaluations in increasing order of magnitude or use Kahan compensated summation to minimize accumulated roundoff error.
- Ensure panel spacing h is computed as (b - a)/N rather than repeatedly adding h, preventing cumulative floating point drift.
- Verify that the integrand is sufficiently smooth (at least C^{2m} for m extrapolation levels) before relying on Romberg results.
- Monitor differences between successive Romberg diagonal elements |R(k, k) - R(k-1, k-1)| as a rigorous practical error estimator.
- Ensure floating point precision buffer remains 64-bit IEEE 754 to prevent precision collapse during Richardson division.
Theoretical limitations
- The composite trapezoidal rule is inefficient for oscillatory or peaked integrands unless adaptive mesh refinement is employed.
- Standard Romberg extrapolation fails to achieve high-order speedup if the integrand has endpoint singularities (e.g. x^{-1/2} ln(x)).
- High-order Romberg extrapolation can become sensitive to roundoff error when panel counts N exceed 2^{12} in double precision.