25. Definite Integrals, Part 4: Romberg Integration#

References:

  • Section 5.3 Romberg Integration of Sauer

  • Section 4.5 Romberg Integration of Burden&Faires

25.1. Introduction#

Romberg Integration is based on repeated Richardson extrapolalation from the composite trapezoidal rule, starting with one interval and repeatedly doubling. Our notation starts with

\[R_{i,0} = T_{2^i}, \quad i=0, 1, 2, \dots\]

where

\[ T_n = \left(\frac{f(a)}{2} + \sum_{k=1}^{n-1} f(a + kh) + \frac{f(b)}{2}\right) h, \quad h = \frac{b-a}{n} \]

and the second index will indicate the number of extapolation steps done (none so far!)

Actually we only need this \(T_n\) formula for the single trapezoidal rule, to get

\[ R_{0,0} = T_1 = \frac{f(a) + f(b)}{2}(b-a), \]

because the most efficient way to get the other values is recursively, with

\[ T_{2n} = \frac{T_n + M_n}{2} \]

where \(M_n\) is the composite midpoint rule,

\[ M_n = h \sum_{k=1}^n f(a + (k-1/2)h), \quad h = \frac{b-a}{n} \]

Extrapolation is then done with the formula

\[ R_{i,j} = \frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}, \quad j = 1, 2, \dots, i \]

which can also be expressed as

\[ R_{i,j} = R_{i,j-1} + E_{i,j-1}, \;\text{where}\; E_{i,j-1} = \frac{R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}\quad\text{is an error estimate.} \]
\[ R_{i,1} = S_{2n} = \frac{4 T_{2 n} - T_{n}}{4 - 1}, \; n=2^{i-1} \]

25.2. An algorithm, in pseudocode#

The above can now be arranged into a basic algorithm. It does a fixed number \(M\) of levels of extrapolation so using \(2^M\) intervals; a refinement would be to use the above error estimate \(E_{i,j-1}\) as the basis for a stopping condition.

\(n \leftarrow 1\)
\(h \leftarrow b-a\)
\(\displaystyle R_{0,0} = \frac{f(a) + f(b)}{2} h\)
for i from 1 to M:
\(\quad\) \(R_{i,0} = \left( R_{i-1,0} + h \sum_{k=1}^n f(a + (i - 1/2)h) \right)/2\)
\(\quad\) for j from 1 to i:
\(\quad\quad \displaystyle R_{i,j} = \frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}\)
\(\quad\) end for
\(\quad n \leftarrow 2 n\)
\(\quad h \leftarrow h/2\)
end for

This work is licensed under Creative Commons Attribution-ShareAlike 4.0 International