Contents

20. Definite Integrals, Part 4: Romberg Integration

Updated on Thursday, April 1 adding a note about an error estimate that can be used in an algorithm that computes to a specified error tolerance.

References:

  • Section 5.3 Romberg Integration of Sauer

  • Section 4.5 Romberg Integration of Burden&Faires

20.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} \]

With the “Romberg notation” \(R_{i,0} = T_{2^i}\) so that \(n=2^i\), this becomes

\[R_{i+1,0} = \frac{R_{i,0} + M_{2^i}}{2}\]

or more conveniently for below,

\[R_{i,0} = \frac{R_{i-1,0} + M_{2^{i-1}}}{2}\]

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}\]

where

\[E_{i,j-1} = \frac{R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}\]

is an error estimate. This can be used as a stopping condition in an algorithm that computes to a specified error tolerance.

Aside: the first level of extrapolation gives the Composite Simpson’s Rule:

\[ R_{i,1} = S_{2n} = \frac{4 T_{2 n} - T_{n}}{4 - 1}, \; n=2^{i-1} \]

20.2. A first 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

19. Definite Integrals, Part 3: The (Composite) Simpson’s Rule and Richardson Extrapolation 21. Initial Value Problems for Ordinary Differential Equations, Part 1: Basic Concepts and Euler’s Method