8.4 Romberg integration
The composite trapezoidal rule \(R_{k,1}\) with \(2^{k-1}\) subintervals has an error expansion in even powers of \(h\) (the Euler–Maclaurin formula), so Richardson extrapolation applies: \[R_{k,j} = R_{k,j-1}+\frac {R_{k,j-1}-R_{k-1,j-1}}{4^{j-1}-1},\qquad j = 2,\dots ,k .\] This is Romberg integration. The trapezoidal values themselves are cheap to compute recursively, because each new one reuses all the previous function values: with \(h_k = \frac {b-a}{2^{k-1}}\), \[R_{k,1} = \frac 12\left [R_{k-1,1}+h_{k-1}\sum _{i=1}^{2^{k-2}}f\bigl (a+(2i-1)h_k\bigr )\right ] .\]
Example 8.4.1. Apply Romberg integration to \(\int _0^\pi \sin x\,dx\), using \(1, 2, 4, \dots , 32\) subintervals.
Solution. The Romberg table, with \(R_{k,1}\) in the first column, is
| \(0\) | |||||
| \(1.570796\) | \(2.094395\) | ||||
| \(1.896119\) | \(2.004560\) | \(1.998571\) | |||
| \(1.974232\) | \(2.000269\) | \(1.999983\) | \(2.000006\) | ||
| \(1.993570\) | \(2.000017\) | \(2.000000\) | \(2.000000\) | \(2.000000\) | |
| \(1.998393\) | \(2.000001\) | \(2.000000\) | \(2.000000\) | \(2.000000\) | \(2.000000\) |
The first column is the composite trapezoidal rule, which after \(33\) function evaluations is still in error by \(1.6\times 10^{-3}\). The extrapolated values are correct to eight decimals. (The second column is exactly the composite Simpson’s rule.)
Questions on this section
Stuck on something here? Ask below and it stays attached to this topic.