Romberg’s method

Motivation

Romberg’s method combines the composite trapezoidal rule with Richardson extrapolation to systematically cancel the leading error terms from the Euler-Maclaurin expansion. Each extrapolation step eliminates one error term, so with \(n\) columns the final estimate has accuracy \(O(h^{2n})\) — matching a Gauss rule of order \(2n\) but using only \(2^{n-1}+1\) function evaluations. It is an excellent choice when high precision is required for smooth integrands on finite intervals.

Example

use integrate::romberg::romberg_method;

fn square(x: f64) -> f64 {
    x.powi(2)
}

let a = 0.0;
let b = 1.0;

let num_steps: usize = 10;

let integral = romberg_method(square, a, b, num_steps);
println!("{}", integral);

Understanding Romberg’s method

Euler-Maclaurin expansion

The Euler-Maclaurin summation formula expresses the composite trapezoidal estimate \(T_h(f)\) with step size \(h = (b-a)/n\) as

\[ \begin{align} T_h(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^2}{12} \left[ f^\prime(b) - f^\prime(a) \right] - \frac{h^4}{720} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &+ \cdots + K h^{2p-2} \left[ f^{(2p-3)}(b) - f^{(2p-3)}(a) \right] + O(h^{2p}) \end{align} \]

Two features make this expansion ideal for Richardson extrapolation:

1. The error is a pure even-power series in \(h\): terms \(h^2, h^4, h^6, \ldots\) with coefficients that depend on the endpoints but not on \(h\).
2. When the step size is halved to \(h/2\), each coefficient is divided by the corresponding power of 4 — which allows exact cancellation.

Richardson extrapolation

Halving \(h\) gives a new estimate

\[ \begin{split} T_{h/2}(f) &= \int_{a}^{b} f(x)\,dx + \frac{h^2}{4 \cdot 12} \left[ f^\prime(b) - f^\prime(a) \right] \\ &- \frac{h^4}{16 \cdot 720} \left[ f^{(3)}(b) - f^{(3)}(a) \right] + \cdots \end{split} \]

Forming the linear combination \(\frac{4 T_{h/2}(f) - T_h(f)}{3}\) cancels the \(h^2\) term exactly:

\[ \begin{split} \frac{4 T_{h/2}(f) - T_h(f)}{3} &= \int_{a}^{b} f(x)\,dx + \frac{h^4}{2880} \left[ f^{(3)}(b) - f^{(3)}(a) \right] \\ &- \frac{h^6}{96768} \left[ f^{(5)}(b) - f^{(5)}(a) \right] + \cdots \end{split} \]

This extrapolation step improves the order from \(O(h^2)\) to \(O(h^4)\). Applying it again to two such \(O(h^4)\) estimates cancels the \(h^4\) term and yields \(O(h^6)\), and so on.

The Romberg table

The full process is organised into a triangular table \(R[i, j]\), where \(i\) is the refinement level (step size \(h_i = (b-a)/2^i\)) and \(j\) is the extrapolation column:

\[ R[i, j] = \frac{4^j , R[i, j-1] - R[i-1, j-1]}{4^j - 1}, \qquad 1 \leq j \leq i \]

The first column \(R[i, 0]\) is the composite trapezoidal estimate at step size \(h_i\). Column \(j\) contains estimates that are \(O(h^{2(j+1)})\) accurate. The divisors \(4^j - 1\) take the values \(3, 15, 63, 255, \ldots\) as \(j\) increases.

For \(n\) columns the table has \(n(n+1)/2\) entries and the bottom-right entry \(R[n-1, n-1]\) is the final estimate, with error \(O(h^{2n})\). For example, with \(n = 10\) columns the error behaves as \(O(h^{20})\), which for smooth integrands gives double-precision accuracy with far fewer evaluations than any fixed Newton-Cotes rule of that order.

Computation

Midpoint reuse. Each first-column entry \(R[i, 0]\) is obtained from the previous one by adding only the \(2^{i-1}\) new midpoints — the \(2^{i-1}+1\) points already evaluated at level \(i-1\) lie exactly on the level-\(i\) grid:

\[ \begin{split} R[i, 0] &= \frac{1}{2}, R[i-1, 0] \\ &+ h_i \sum_{k=0}^{2^{i-1}-1} f\left(a + (2k+1),h_i\right) \end{split} \]

This means the total number of function evaluations across all \(n\) levels is \(2^{n-1}+1\), not \(O(n \cdot 2^{n-1})\) as would be required if each level recomputed all its points independently.

Rolling buffer. The Richardson extrapolation depends only on the previous row and the current row of the table. The implementation therefore keeps only two rows of length \(n\) in memory at any time, rather than the full \(n \times n\) table, keeping the memory footprint at \(O(n)\).

Summary of cost for \(n\) columns:

Limitations

Romberg’s method requires the integrand to be infinitely differentiable on \([a, b]\) for the Euler-Maclaurin expansion to hold and the cancellation to be effective. It is not appropriate for functions with singularities, discontinuities, or even kinks (non-differentiable points) in the integration interval — at such points the error expansion breaks down and the extrapolation accelerates to the wrong limit.

The number of columns \(n\) controls both the accuracy and the number of evaluations. For double-precision arithmetic, \(n \approx 10\) is typically sufficient for smooth integrands; using larger \(n\) beyond the point where \(R[n-1, n-1]\) has converged to machine precision adds evaluations without improving the result. For functions with near-singularities or rapid oscillation, prefer the adaptive Simpson method.