Romberg rule - Executable Engineering: Interactive Numerical Methods in Python

import numpy as np
import scipy as sp
import plotly.graph_objects as go

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

Richardson Extrapolation¶

An alternative to the Newton-Cotes formulae uses successive sampling step sizes to estimate and cancel error. We begin with Richardson extrapolation and extend this to Romberg integration.

Richardson extrapolation is an algorithm that uses successive applications of the trapezoid rule with differing step sizes to achieve superior accuracy without computing higher derivatives.

The exact integral can always be expressed as:

I = I(h) + E(h)

(1)

where $I(h)$ is the approximation with step $h$ and $E(h)$ is the associated truncation error. Writing this for step $h_1$ and $h_2 = h_1 / 2$ :

\begin{aligned} I &= I(h_1) + c h_1^2 + O(h_1^4) \\ I &= I(h_2) + c h_2^2 + O(h_2^4) = I(h_2) + c \frac{h_1^2}{4} + O(h_1^4) \end{aligned}

(2)

NB $c$ here, which depends on $f''$ , is assumed to be the same between steps. This subtly assumtion is only a good one for well-behaved functions.

Multiplying the second equation by 4 and subtracting the first eliminates the leading $h_1^2$ error term:

I = \frac{4 I(h_2) - I(h_1)}{3} + O(h_1^4)

(3)

This is an interesting result! Effectively what we have done is use a second estimate to cancel the next power in our expansion, leading to a higher-order result without evaluating new derivatives.

The Romberg Integration Algorithm¶

In fact, we can repeat this procedure arbitrarily! Above we combined two $O(h^2)$ estimates to make an $O(h^4)$ estimate. We can take this result, combine it with another, and make an $O(h^6)$ estimate, and so on.

This recursive extrapolation generates the Romberg tableau:

\begin{align} I^4 &\approx \frac{4}{3} I(h_2) -\frac{1}{3} I(h_1) \\ I^6 &\approx \frac{16}{15} I(h_3) - \frac{1}{15} I^4 \\ I^8 &\approx \frac{64}{63} I(h_4) - \frac{1}{63} I^6 \\ \vdots \\ I_{j,k} &\approx \frac{4^{k-1} I_{j+1,k-1} - I_{j,k-1}}{4^{k-1}-1} \end{align}

(4)

where the last line is the Romberg Integration Algorithm which has some very nice properties:

lends itself to recursive programming (elegant)
parallelizes nicely since each interval is independent.
can integrate to arbitrary accuracy (if a callable function is provided)
Can reuse function evaluations if the step size is chosen well.

However, it requires either:

the ability to sample the function at will
a structured sampling grid

Remember however that we assumed the second derivative doesn’t change much between levels. If there is high frequency noise and a low frequency trend, the cancellation of the errors isn’t so robust. For these reasons, heuristics based on Newton-Cotes (Decision trees of applications of Simpson’s rules and trapezoid depending on sampling structure) remains a common robust choice.

Example: Romberg Integration vs. Simpson’'s Rule¶

Approximate $\int_0^\pi \sin(x) \,\mathrm{d}x = 2$ using the Romberg rule and compare with Simpson’'s 1/3 rule:

def f(x):
    return np.sin(x)

for n in [4, 8, 16, 32, 64, 128]:
    x = np.linspace(0, np.pi, n + 1)
    f_x = f(x)
    romberg_val = sp.integrate.romb(f_x, dx=np.pi / n, show=True)
    simpson_val = sp.integrate.simpson(f_x, x=x)
    print(f"Intervals n={n:3d} | Romberg error: {abs(romberg_val - 2):.2e} | Simpson error: {abs(simpson_val - 2):.2e}")

Convergence to Machine Precision¶

Note the error became zero (or reached double-precision floating point limits $\approx 10^{-16}$ )! What does that mean?

By systematically canceling higher-order polynomial truncation terms, Romberg integration converges exponentially for smooth functions with very few function evaluations.