Romberg rule

Richardson extrapolation

Richardson extrapolation is an algorithm that conceptually uses successive applications of the trapezoid rule with differing step size to achieve superior results with less effort.

The exact integral can always be expressed: $\( I = I'(h) + E(h)\)\( where \)I’(h)\( is the approximation with step \)h\( and the associated error \)E(h)\(. We know that \)E(h) \propto h^2\( for the trapezoid rule. In fact, it is \)E(h) \propto f’’ h^2$.

Let’s sample the interval twice with step sizes \(h_1\) and \(h_2\). If we assume \(f''\) doesn’t change much we can say,

\[ \begin{align} \frac{E(h_1)}{E(h_2)} = \frac{h_1^2}{h_2^2}. \end{align}\]

Now, since the exact integral is the same in both cases, $\(\begin{align} I'(h_1) + E(h_1) &= I'(h_2) + E(h_2) \\ I'(h_1) + E(h_2) \bigg(\frac{h_1}{h_2}\bigg)^2 &= I'(h_2) + E(h_2) \\ E(h_2) &=\frac{I'(h_1) -I'(h_2)}{1-(\frac{h_1}{h_2})^2} \\ \end{align} \)$

and inserted into the formula for \(h_2\),

\[\begin{split} \begin{align} I&\approx I(h_2)+ \frac{I'(h_1) -I'(h_2)}{1-(\frac{h_1}{h_2})^2} \\ \end{align} \end{split}\]

which can be shown to be accurate to \(O(h^4)\)!

For the special case where \(h_1 = 2 h_2\) (which has advantages for overlapping point evaluations)

\[ I \approx \frac{4}{3} I(h_2) -\frac{1}{3} I(h_1) \]

This is an interesting result! Effectively what we have done is use a second estimate to estimate the next power in our expansion, leading to a higher order estimate!

Romberg Integration Algorithm

In fact we can repeat this proceedure arbitrarily! Above we combined two order \(O(h^2)\) to make \(O(h^4)\). We can take this results, combine it with another sampling at \(h_3 < h_2\) and combine to get an \(O(h^6)\) estimate and so on! If we successively halve the step size, we can get:

\[\begin{split}\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}\end{split}\]

where the last line is the Romberg Integration Algorithm. The structure lends itself to redundant programming and parallelizes nicely!

This algorithm is able to integrate to an arbitrary accuracy and does so remarkably efficiently compared to the alternatives.

Example

Approximate \(\int_0^\pi \sin(x) dx\) using the Rhomberg rule and compare with Simpson’s 1/3 rule

import numpy as np
import scipy as sp
def f(x):
  return np.sin(x)

tolerance = 1e-6

for n in [4,8,16,32,64, 128]:
  x = np.linspace(0, np.pi, n+1)
  f_x = f(x)
  rhomberg = sp.integrate.romb(f_x, dx = np.pi/n, show=True)
  print(f"Romberg with {n} intervals: {abs(rhomberg-2)}")
  simpson = sp.integrate.simpson(f_x, x=x)
  print(f"Simpson with {n} intervals: {abs(simpson-2)}")

Richardson Extrapolation Table for Romberg Integration
======================================================
 0.00000 
 1.57080  2.09440 
 1.89612  2.00456  1.99857 
======================================================
Romberg with 4 intervals: 0.001429268176164289
Simpson with 4 intervals: 0.0045597549844207386
Richardson Extrapolation Table for Romberg Integration
======================================================
 0.00000 
 1.57080  2.09440 
 1.89612  2.00456  1.99857 
 1.97423  2.00027  1.99998  2.00001 
======================================================
Romberg with 8 intervals: 5.549979670949057e-06
Simpson with 8 intervals: 0.00026916994838765973
Richardson Extrapolation Table for Romberg Integration
======================================================
 0.00000 
 1.57080  2.09440 
 1.89612  2.00456  1.99857 
 1.97423  2.00027  1.99998  2.00001 
 1.99357  2.00002  2.00000  2.00000  2.00000 
======================================================
Romberg with 16 intervals: 5.412709835894702e-09
Simpson with 16 intervals: 1.6591047935499148e-05
Richardson Extrapolation Table for Romberg Integration
======================================================
 0.00000 
 1.57080  2.09440 
 1.89612  2.00456  1.99857 
 1.97423  2.00027  1.99998  2.00001 
 1.99357  2.00002  2.00000  2.00000  2.00000 
 1.99839  2.00000  2.00000  2.00000  2.00000  2.00000 
======================================================
Romberg with 32 intervals: 1.3216094885137863e-12
Simpson with 32 intervals: 1.0333694131503535e-06
Richardson Extrapolation Table for Romberg Integration
======================================================
 0.00000 
 1.57080  2.09440 
 1.89612  2.00456  1.99857 
 1.97423  2.00027  1.99998  2.00001 
 1.99357  2.00002  2.00000  2.00000  2.00000 
 1.99839  2.00000  2.00000  2.00000  2.00000  2.00000 
 1.99960  2.00000  2.00000  2.00000  2.00000  2.00000  2.00000 
======================================================
Romberg with 64 intervals: 4.440892098500626e-16
Simpson with 64 intervals: 6.453000178652246e-08
Richardson Extrapolation Table for Romberg Integration
======================================================
 0.00000 
 1.57080  2.09440 
 1.89612  2.00456  1.99857 
 1.97423  2.00027  1.99998  2.00001 
 1.99357  2.00002  2.00000  2.00000  2.00000 
 1.99839  2.00000  2.00000  2.00000  2.00000  2.00000 
 1.99960  2.00000  2.00000  2.00000  2.00000  2.00000  2.00000 
 1.99990  2.00000  2.00000  2.00000  2.00000  2.00000  2.00000  2.00000 
======================================================
Romberg with 128 intervals: 0.0
Simpson with 128 intervals: 4.032257194808153e-09

Note the error became zero! What does that mean?