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}


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


 \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.


the Romberg formula

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


the Romberg integral
a Romberg integral
					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