en.wikipedia.org/wiki/Romberg%27s_method
young.physics.ucsc.edu/115/romberg.pdf
The trapezium rule is a simple numerical approximation to the integral \begin{equation}I=\int_{a}^{b} f(x) d x\end{equation}For reasons that will become clear later we will write the trapezium rule for $n$ intervals as $I_{n}^{(0)}$, i.e.
\begin{equation}I_{n}^{(0)} \equiv T_{n}=h\left[\frac{1}{2} f_{0}+f_{1}+f_{2}+\cdots+f_{n-1}+\frac{1}{2} f_{n}\right]\end{equation}where\begin{equation}h=\frac{b-a}{n}\end{equation}is the width of one interval, $f_{i} \equiv f\left(x_{i}\right)$, with $x_{i}=x_{0}+i h$ and $x_{0}=a, x_{n}=b$. We have seen in class that, for small $h$, the error is proportional to $h^2$. It turns out that if one writes an expression
for the error as a power series in $h$ then only even powers of $h$ appear (assuming that the integrand doesn’t have any singularities). In other words,
\begin{equation}I=I_{n}^{(0)}+A h^{2}+B h^{4}+C h^{6}+\cdots\end{equation}where $A, B, C, · · ·$ are numbers which are related to derivatives of $f (x)$ at the endpoints (which of course we don’t know in the numerics). We already showed that the leading error in the trapezium rule is of order $h^2$, i.e. there is no term linear in $h$. Equation (4), with explicit expressions for $A, B, C, · · ·$ , is called the Euler-Maclaurin formula.
To obtain a more accurate estimate for $I$ we will eliminate the leading contribution to the error in Eq. (4), the term of order $h^2$, by taking $n$ to be even and determining the trapezium rule for $n/2$ intervals as well as for $n$ intervals. Since the width of one interval is now $2h$ we have\begin{equation}I_{n / 2}^{(0)}=2 h\left[\frac{1}{2} f_{0}+f_{2}+f_{4} \cdots+f_{n-2}+\frac{1}{2} f_{n}\right]\end{equation}The error in this expression is of the same form as Eq. (4) but with $h$ replaced by $2h$, i.e.\begin{equation}I=I_{n / 2}^{(0)}+A(2 h)^{2}+B(2 h)^{4}+C(2 h)^{6}+\cdots\end{equation}We can eliminate the leading contribution to the error by subtracting 4 times Eq. (4) from Eq. (6).
Dividing by 3 we get\begin{equation}I=\frac{4 I_{n}^{(0)}-I_{n / 2}^{(0)}}{3}-4 B h^{4}-20 C h^{6}+\cdots\end{equation}Hence we obtain the estimate $I_{n}^{(1)}$, where\begin{equation}I_{n}^{(1)}=\frac{4 I_{n}^{(0)}-I_{n / 2}^{(0)}}{3}\end{equation}From Eq. (7) this is related to the exact value of the integral $I$ by\begin{equation}I=I_{n}^{(1)}+B^{\prime} h^{4}+C^{\prime} h^{6}+\cdots\end{equation}where $B'= −4B$ and $C' = −20C$.
From Eqs.(2), (5) and (8) we find that $I_{n}^{(1)}$ is given in terms of the function values by\begin{equation}I_{n}^{(1)}=\frac{h}{3}\left[f_{0}+4 f_{1}+2 f_{2}+4 f_{3}+2 f_{4}+\cdots 2 f_{n-2}+4 f_{n-1}+f_{n}\right]\end{equation}which is just Simpson’s rule for $n$ intervals, i.e. $I_{n}^{(1)} \equiv S_{n}$.
Hence the first level of error elimination in Romberg integration gives Simpson’s rule. However, the higher levels of error elimination that we discuss next will not produce familiar approximations
1.
If $n$ is a multiple of 4, we can clearly repeat the trick of eliminating the leading contribution to the error, which is now the $B'h^4$ term in Eq. (9), by forming the appropriate linear combination of $I^{(1)}_n$ and $I^{(1)}_{n/2}$. Hence the second level of error extrapolation gives $I^{(2)}_n$ where\begin{equation}I_{n}^{(2)}=\frac{16 I_{n}^{(1)}-I_{n / 2}^{(1)}}{15}\end{equation}which has an error of order $h^6$. In general, if $n$ is a multiple of $2k$ we can form iteratively a sequence of higher order approximations, $I^{(k)}_n$, where\begin{equation}I_{n}^{(k)}=\frac{4^{k} I_{n}^{(k-1)}-I_{n / 2}^{(k-1)}}{4^{k}-1}\end{equation}for $k = 1, 2, 3, · · ·$, which have an error of order $h^{2k+2}$.
Romberg integration is a convenient way of carrying out this procedure of successively eliminating terms in the expression for the error. One proceeds as follows. Start with 1 interval and compute $I^{(0)}_1$. Then do two intervals and compute $I^{(0)}_2$ from which $I^{(1)}_2$ can be determined using Eq. (8). Then double the number of intervals again to 4 and calculate $I^{(0)}_4$, from which $I^{(1)}_4$ can be obtained using Eq. (8), and hence $I^{(2)}_4$ obtained from Eq. (12). The sequence of quantities thus obtained can be conveniently represented as a table,
Recall that $n$ is the number of intervals used in the trapezium rule, and $k$ is the number of levels of error elimination. The left hand column of estimates for the integral, labeled 0, contains results for the trapezium rule determined by evaluating the function at appropriate points. Each subsequent column is obtained, not from any more function evaluations, but rather by manipulating the data
in the column to the left of it. These columns contain results for a higher order approximation so, for example, the $k = 1$ column is Simpson’s rule. Note that the function values obtained for the trapezium rule with $n$ intervals also occur in the trapezium rule for $2n$ intervals, and so would not be recomputed, in a well written code when evaluating $T_{2 n}\left(\equiv I_{2 n}^{(0)}\right)$
