The Romberg quadrature formula is also called the successive half-acceleration method. It is based on the relationship between the trapezoidal formula, the Simpson’s formula and the higher-order Newton-Cotes formula to construct a method to accelerate the calculation of integral. The Romberg algorithim is an extrapolation method of combining the previous approximations to generate more accurate apprximations in the process of successively dividing the integration interval into half. It improves the accuracy of the integral without increasing the amount of calculation.
According to the error estimation of the trapezoidal rule, it can be seen that the truncation error of the integral or tyhe approximation value \( T_n \) is roughly proportional to \( h^2 \). So, when the step size is divided by two sub-intervals which are double than that of the numerator, the truncation error \( 1-T_n \) will them be reduced to \( \frac{1}{4} \) of the original error. Mathematically,
\( \displaystyle \frac{1-T_{2n}}{1-T_n}\approx\frac{1}{4} \)
Rearranging the terms, we get
\( \displaystyle 1-T_{2n}\approx\frac{1}{3}(T_{2n}-T_n) \)
It can be seen that as long as the two successive integral values \( T_n \) and \( T_{2n} \), before and after the further division are close enough, the error of \( T_{2n} \) will alos be small. The error of \( T_{2n} \) is roughly equal to \( \frac{1}{3}(T_{2n}T_n) \). We can use this error to compensate for \( T_{2n} \) to get a better approximation, that is a more accurate approximation
\( \displaystyle \Bar{T}=T_{2n}+\frac{1}{3}(T_{2n}-T_n)=\frac{4}{3}4T_{2n}-\frac{1}{3}T_n \)
In other words, we can say that we linearly combine two integral values \( T_n \) and \( T_{2n} \) computed using the trapezoidal rule to get more accurate approximations. We can go further using the similar approach.
We can generate a Simpson sequence, \( \{S_{2k}\} \), by linearly ombining the values from the trapezoidal sequence \( \{T_{2k}\} \). The Simpson sequence has a faster convergence rate than that of the trapezoidal sequence.
\( \displaystyle \{S_{2k}\}:S_1,S_2,S_4\hdots \)
\( \displaystyle S_n=\Bar{T}=\frac{4}{3}T_{2n}-\frac{1}{3}T_n=\frac{4T_{2n}-T_n}{4-1} \)
Similary,
\( \displaystyle I\approx S_{2n}+\frac{1}{15}(S_{2n}-S_n) \)
We can go further by combining the values from the Simpson sequence to produce a Newton-Cotes sequence \( \{C_{2k}\} \) with a faster onvergence rate. It can be shown thart
\( \displaystyle C_{2n}=\frac{16}{15}S_{2n}-\frac{1}{15}S_n=\frac{4^2S_{2n}-S_n}{4^2-1} \)
Similarly,
\( \displaystyle I\approx C_{2n}+\frac{1}{63}(C_{2n}-C_n) \)
Similarly, like or the above three sequences, we can go further from the Newton-Cotes sequence to produce a Romberg sequence \( \{R_{2k}\} \) with an even faster convergence rate.
\( \displaystyle \{R_{2k}\}:R-1,R-2,R-4\hdots \)
\( \displaystyle R_n=\frac{64}{63}C_{2n}-\frac{1}{63}C_n=\frac{4^3C_{2n}-C_n}{4^3-1} \)
By using the acceleration formula in the process of variable step size, we gradually process the rough trapezoidal values \( T_n \) into the fine Simpson values \( S_n \), then into the even finer Newton-Cotes sequence \( C_ n \) and then further into the Romberg sequence \( R_n \) which has the highest precision. The following diagram illustrates the following sequences.
Now consider the example of calculating \( \pi \) using the integral given,
\( \displaystyle I=\int_{0}^{1}ydx=\int_{0}^{1}\frac{4}{1+x^2}dx=\pi \)
Below are the sequence of Trapezoidal values, Simpson alues, Newton-Cotes values and the Romberg values produced.
\( \displaystyle f(x)=\frac{4}{1+x^2} \)
\( \displaystyle T_1=\frac{1}{2}[f(0)+f(1)]=3 \)
\( \displaystyle T_2=\frac{1}{2}T_1+\frac{1}{2}f\Big(\frac{1}{2}\Big)=3.1 \)
\( \displaystyle S_1=\frac{4}{3}T_2-\frac{1}{3}T_1=3.13333 \)
\( \displaystyle T_4=\frac{1}{2}T_2+\frac{1}{4}\Big[f\Big(\frac{1}{4}\Big)+f\Big(\frac{3}{4}\Big)\Big]=3.131177 \)
\( \displaystyle S_2=\frac{4}{3}T_4-\frac{1}{3}T_2=3.141569 \)
\( \displaystyle T_8=\frac{1}{2}T_4+\frac{1}{8}\Big[f\Big(\frac{1}{8}\Big)+f\Big(\frac{3}{8}\Big)+f\Big(\frac{5}{8}\Big)+f\Big(\frac{7}{8}\Big)\Big]=3.138989 \)
\( \displaystyle S_4=\frac{8}{3}T_8-\frac{1}{3}T_4=3.141593 \)
\( \displaystyle C_1=\frac{16}{15}S_2-\frac{1}{15}S_1=3.142118\)
\( \displaystyle C_2=\frac{16}{15}S_4-\frac{1}{15}S_2=3.1415946 \)
\( \displaystyle R_1=\frac{64}{63}C_2-\frac{1}{63}C_1=3.141586292 \)
val f: final UnivariateRealFunction f = new AbstractUnivariateRealFunction() {
override public double evaluate(double x) {
return exp(2 * x) - 4 * x - 7;
}
};
IterativeIntegrator integrator1 = new Trapezoidal(1e-8, 20); // using the trapezoidal rule
Integrator integrator1 = new Simpson(1e-8, 20); // using the Simpson rule
Integrator integrator3 = new Romberg(integrator1);
double a = 0, b = 1; //the limit
// the integrations
double I1 = integrator1.integrate(f, a, b);
Sysytem.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %.16f, using the trapezoidal rule", a, b, I1));
double I2 = integrator2.integrate(f, a, b);
Sysytem.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %.16f, using the Simpson rule", a, b, I2));
double I3 = integrator3.integrate(f, a, b);
Sysytem.out.println(String.format("S_[%.0f,%.0f] f(x) dx = %.16f, using the Romberg formula", a, b, I3));
The output is:
S_[0,1] f(x) dx = -5.8054719346672840, sing the trapezoidal rule
S_[0,1] f(x) dx = -5.8054719494768790, sing the Simpson rule
S_[0,1] f(x) dx = -5.8054719505327520, sing the Romberg formula