Romberg Integration

Romberg integration applies Richardson extrapolation to the trapezoidal rule, building a table that converges rapidly for smooth integrands.

Extrapolating the trapezoidal rule

The composite trapezoidal rule has an error expansion in even powers of the step size h. Romberg integration exploits this by computing trapezoidal estimates at h, h/2, h/4, ... and applying Richardson extrapolation to cancel successive error terms, producing a triangular table whose entries climb quickly in accuracy.

The table

The first column holds trapezoidal estimates with halved step sizes. Each next column applies R(i,j) = R(i,j-1) + (R(i,j-1) - R(i-1,j-1))/(4^j - 1), cancelling the next h^2 term. The diagonal entries converge fastest; a few rows often reach machine precision for smooth integrands.

def romberg(f, a, b, k):
    R = [[0.0]*(k) for _ in range(k)]
    h = b-a
    R[0][0] = 0.5*h*(f(a)+f(b))
    for i in range(1, k):
        h *= 0.5
        s = sum(f(a+(2*j-1)*h) for j in range(1, 2**(i-1)+1))
        R[i][0] = 0.5*R[i-1][0] + h*s
        for j in range(1, i+1):
            R[i][j] = R[i][j-1] + (R[i][j-1]-R[i-1][j-1])/(4**j - 1)
    return R[k-1][k-1]

Reusing evaluations

A key efficiency is that halving the step reuses all previous function values; only the new midpoints are evaluated. This makes each new row cheap. The method is ideal when the integrand is smooth and evaluations are moderately expensive.

Limits

Romberg's rapid convergence depends on the smooth error expansion, which requires an integrand with many continuous derivatives. Singularities, discontinuities, or endpoint blowups break the expansion, and adaptive Gauss-Kronrod methods handle those cases better. For clean, smooth integrands, Romberg remains one of the most efficient one-dimensional methods.

Romberg integration is a convenient high-accuracy default for smooth profile integrals in the breeder Hyperion analyses.