Romberg's Method

Romberg's Method (also called Romberg integration) is an extrapolation-based numerical integration algorithm in computation and algorithms and numerical analysis that approximates the definite integral I = ∫_a^b f(x) dx of a smooth function by applying Neville's algorithm (recursive Richardson extrapolation) to sequence refinements of the composite trapezoidal rule. Formulated by German-Norwegian physicist Werner Romberg in 1955 building on Lewis Fry Richardson's deferred approach to the limit (1911, 1927), the method organizes approximations into a lower-triangular array R(k, m) where each column eliminates the next leading even power h^2m in the asymptotic Euler–Maclaurin error expansion.

By evaluating the integrand on nested dyadic grids with step sizes h_k = (b − a) / 2^k, Romberg's method reuses 100% of previously computed function evaluations when halving the mesh spacing from h_k−1 to h_k. The first extrapolation column R(k, 1) reproduces composite Simpson's rule of order O(h⁴), the second column R(k, 2) reproduces Boole's rule of order O(h⁶), and the diagonal entries R(k, k) achieve superlinearly accelerating algebraic precision O(h_k^2k+2) without the negative weights that destabilize high-order equispaced Newton–Cotes formulas—complementing optimal non-equispaced schemes such as Gauss–Legendre quadrature and underpinning Bulirsch–Stoer extrapolation for Runge–Kutta methods.

The Asymptotic Even-Power Error Series of the Trapezoidal Rule

Let f ∈ C^2M+2([a, b]) be a sufficiently smooth integrand on a finite interval [a, b]. Partitioning [a, b] into N_k = 2^k uniform subintervals of step size h_k = (b − a) / 2^k (k = 0, 1, 2, ...), the composite trapezoidal rule T(h_k) = R(k, 0) is defined as:

By the Euler–Maclaurin summation formula, the difference between the true integral I = ∫_a^b f(x) dx and the trapezoidal approximation T(h) expands in even powers of the step size h:

where each coefficient c_m = −(B_2*m / (2m)!) [ f^(2m−1)(b) − f^(2m−1)(a) ] depends on the Bernoulli numbers B_2m (B₂ = 1/6, B₄ = −1/30, B₆ = 1/42) and the odd endpoint derivatives of f(x), but is independent of h.

Non-Redundant Nested Function Evaluation

Starting from the single-panel trapezoid R(0, 0) = ((b − a) / 2)[f(a) + f(b)] with h₀ = b − a, each subsequent refinement R(k, 0) at half the step size h_k = h_k−1 / 2 shares all 2^k−1 + 1 grid nodes of R(k − 1, 0) and requires evaluating f(x) only at the 2^k−1 new odd-indexed midpoints:

Thus, computing the entire first column R(0, 0), R(1, 0), ..., R(K, 0) requires 2^K + 1 total function evaluations rather than 2^K+1 + K.

Elimination of Leading Error Terms and the 4^m Recurrence

Because I = R(k − 1, 0) + c₁h_k−1² + c₂h_k−1⁴ + ··· and h_k−1 = 2*h_k*, comparing two consecutive trapezoidal levels R(k − 1, 0) and R(k, 0) reveals that the h² error term shrinks by a factor of 2² = 4 when h is halved:

Subtracting R(k − 1, 0) from 4R(k, 0) and dividing by 3 cancels the c₁h_k² error term, yielding R(k, 1) = R(k, 0) + (R(k, 0) − R(k − 1, 0)) / 3 = I + O(h_k⁴). By mathematical induction across column index m = 1, 2, ..., k, eliminating the O(h_k^2m) term between R(k − 1, m − 1) and R(k, m − 1) produces the universal Romberg extrapolation recurrence:

Each entry R(k, m) has truncation error O(h_k^2m+2) and integrates polynomials of degree up to 2m + 1 with zero truncation error:

Positive Linear Weights and Absence of Runge Instability

A remarkable analytical property proved by Bauer, Rutishauser, and Stiefel (1963) is that whereas closed Newton–Cotes formulas of degree ≥ 8 possess large alternating positive and negative weights whose absolute sum ∑|w_i| diverges as n → ∞ (amplifying IEEE 754 floating-point rounding errors via the Runge phenomenon), every entry R(k, m) in the Romberg tableau is a linear combination of function values f(x_j) with strictly positive weights summing to b − a. Consequently, Romberg's diagonal sequence R(k, k) converges to ∫_a^b f(x) dx as k → ∞ for every continuous function f ∈ C([a, b]), and achieves superlinear convergence whenever f ∈ C^∞([a, b]).

Bulirsch–Stoer Step Sequences and Adaptive Stopping Criteria

In practice, Romberg integration terminates when successive diagonal approximations satisfy the relative error criterion |R(k, k) − R(k − 1, k − 1)| ≤ ε_tol |R(k, k)|. One limitation of pure dyadic halving (N_k = 2^k = 1, 2, 4, 8, 16, 32, 64, ...) is that the evaluation count doubles at every row. Roland Bulirsch and Josef Stoer (1964) introduced the gentler step-count sequence N_k = {1, 2, 3, 4, 6, 8, 12, 16, 24, 32, ...} (alternating powers of 2 and 3 · 2^j), which replaces the factor 4^m in Neville's extrapolation formula with (h_k−m / h_k)² to restrain function evaluations in both numerical quadrature and Bulirsch–Stoer ordinary differential equation solvers alongside Runge–Kutta methods and spectral fast Fourier transform Clenshaw–Curtis rules.

Why does Romberg's method use 4^m − 1 instead of 2^m − 1 in its extrapolation formula?

Because the composite trapezoidal rule is symmetric under reversing the interval orientation (h → −h), all odd powers of h (h¹, h³, h⁵, ...) vanish identically in the Euler–Maclaurin expansion. Each extrapolation step therefore eliminates an even power h^2m rather than h^m, so halving h reduces the leading error term by (2²)^m = 4^m at each column m.

When does Romberg's method fail to accelerate convergence?

Romberg extrapolation relies on the validity of the Euler–Maclaurin expansion, which requires the integrand f(x) to possess continuous higher-order derivatives at the endpoints a and b. If f(x) has an integrable endpoint singularity—such as f(x) = √x on [0, 1], whose error expands in fractional powers h^3/2, h², h^5/2—standard 4^m extrapolation does not cancel the leading h^3/2 error term unless modified powers or a smoothing variable transformation are applied.

How does Romberg's method compare in efficiency with Gauss–Legendre quadrature?

To achieve polynomial degree 2k + 1, the Romberg diagonal entry R(k, k) requires 2^k + 1 equispaced function evaluations, whereas an (n = k + 1)-point Gauss–Legendre quadrature rule achieves the same polynomial degree 2k + 1 using only k + 1 evaluations (orders of magnitude fewer points for high k). However, Romberg's method requires no precomputed irrational nodes or weights, reuses 100% of function values across all refinement levels, and provides built-in error estimates along every row and diagonal.

1. Romberg, W. (1955). "Vereinfachte numerische Integration." Det Kongelige Norske Videnskabers Selskabs Forhandlinger, 28(7), 30–36. https://mathscinet.ams.org/mathscinet-getitem?mr=0073884
2. Bauer, F. L., Rutishauser, H., & Stiefel, E. (1963). "New aspects in numerical quadrature." Proceedings of Symposia in Applied Mathematics, 15, 199–218. https://doi.org/10.1090/psapm/015/0158539
3. Bulirsch, R., & Stoer, J. (1964). "Fehlerabschätzungen und Extrapolation mit rationalen Funktionen bei Verfahren vom Richardson-Typus." Numerische Mathematik, 6(1), 413–427. https://doi.org/10.1007/BF01386092
4. Stoer, J., & Bulirsch, R. (2002). Introduction to Numerical Analysis (3rd ed., Section 3.4: Integration by Extrapolation — Romberg's Method, pp. 160–167). Springer. https://doi.org/10.1007/978-0-387-22750-4
5. Davis, P. J., & Rabinowitz, P. (1984). Methods of Numerical Integration (2nd ed., Chapter 6.3: Romberg Integration). Academic Press. https://doi.org/10.1016/C2013-0-10566-1