Subject : The following demonstrates the convergence of the Romberg method of estimating integrals of continuous functions.
Authors : Nathan Collier, Autar Kaw, Loubna Guennoun
Date : 27 October 2005

Introduction

Romberg integration is based on the trapezoidal rule, where we use two estimates of an integral to compute a value that is more accurate than the previous estimates. This is called Richardson's extrapolation. Thus,
	I = I(h) + E(h)
	h = (b-a) / n
where I is the exact value of the integral, I(h) is the approximate integral using the trapezoidal rule with n segments, and E(h) is the truncation error. A general form of Romberg integration is:
	\(TraditionalForm\`I\_\(j, k\)\  = \(\(\ \)\(4\^\(k - 1\)\ I\_\(j + 1, k - 1\)\  - \ I\_\(j, k - 1\)\)\)\/\(\(\ \)\(4\^\(k - 1\)\  - \ 1\)\)\)
where the index j is the order of the estimate integral, and k is the level of integration. [click here for textbook notes] [click here for power point presentation].

Inputs

The following simulation will illustrate Romberg integration. This section is the only section where the user may interacts with the program. The user may enter any function f(x) and the lower and upper limit for the function. By entering this data, the program will calculate the exact value of the integral, followed by the results using the trapezoidal rule with n = 1, 2, 4, 8 segments, and the Romberg integration for each segments. The program will also display the approximate error, the absolute relative approximate percentage error, the least correct significant digits, and the least number of significant digits correct in the approximation.

Integrand f[x] = 0

f[x_] := 300.0*x/(1.0 + Exp[x]);

Lower and upper limit of the integral, a and b respectively.

a = 0.0;
b = 10.0;

Maximum number of iterations. This number must be even.

(Subscript[n, maximum]) = 8;

curve = 
        Plot[f[x], {x, a, b}, 
          PlotLabel -> "Entered function on given interval", 
          TextStyle -> {FontSize -> 
                11}];
General::spelll  Off;

True Solution

This is the solution found by Mathematica

Actual = Integrate[f[x], {x, a, b}]

Value of integral as a function of iterations

Trapezoidal Rule

Array[TR, 2^((Subscript[n, maximum]) - 1)];

For[i = 1,  i ≤ 2^((Subscript[n, maximum]) - 1), i++, 
      If[i ≤ 1, TR[i] = (f[a] + f[b])/2*(b - a), 
        h = (b - a)/i; 
        TR[i] = h/
              2*(f[a] + 2*Sum[f[a + j*h], {j, 1, i - 1}] + 
                f[b])]]

Romberg Method

Array[AV, (Subscript[n, maximum])];

For[i = 1, 
      i ≤  (Subscript[n, maximum]), i++, 
For[j = 2, 
        j ≤ i, j++, 
nn = 
          2^(j - 1); 
(Subscript[I, i, 1]) = 
          TR[nn]; 
For[k = 2, 
          k ≤ j, k++, 
l = 
            1 + j - k; 
(Subscript[I, l, k]) = (4^(k - 1)*
                    (Subscript[I, l + 1, k - 1]) - 
                  (Subscript[I, l, k - 1]))/(4^(k - 1) - 
                  1)]]; 
(Subscript[I, 1, 1]) = 
        TR[1]; 
AV[i] = (Subscript[I, 1, i])]

Absolute true error

Array[Et, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      Et[i] = Abs[Actual - AV[i]]]

Absolute relative true error

Array[ϵt, (Subscript[n, maximum])];

For[i = 1, 
      i <= (Subscript[n, maximum]), i++, ϵt[i] = 
        Abs[Et[i]/Actual*100]]

Absolute approximate error

Array[Ea, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      If[i <= 1,  Ea[i] =  0, Ea[i] = Abs[AV[i] - AV[i - 1]]]]

Absolute relative approximate error

Array[ϵa, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      If[i <= 1, ϵa[i] = 0, ϵa[i] = 
          Abs[Ea[i]/AV[i]*100]]]

Significant digits at least correct

Array[sigdig, (Subscript[n, maximum])];

For[i = 1, i <= (Subscript[n, maximum]), i++, 
      If[(ϵa[i] ≥ 5)  ||  (i <= 1), 
        sigdig[i] = 0, 
        sigdig[i] = 
          Floor[(2 - Log[10, Abs[ϵa[i]/0.5]])]]]

Graphs

xrplot = Table[AV[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[xrplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[xrplot]}, PlotLabel -> "Estimated root as a function of number \
of iterations"];

Etplot = Table[Et[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[Etplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[Etplot]}, PlotLabel -> "Absolute true error as a function of \
number of iterations"];

ϵtplot = 
        Table[ϵt[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[ϵtplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[ϵtplot]}, PlotLabel -> "Absolute relative true error as a \
function of number of iterations"];

Eaplot = Table[Ea[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[Eaplot,PlotJoined -> True, PlotRange -> All,AxesOrigin -> {1, Min[Eaplot]}, PlotLabel -> "Absolute approximate error as a function \
of number of iterations"];

ϵaplot = 
        Table[ϵa[i], {i, 1, (Subscript[n, maximum])}];

ListPlot[ϵaplot,PlotJoined -> True, PlotRange -> All, AxesOrigin -> {1, Min[ϵaplot]}, PlotLabel -> "Absolute relative approximate error as a \
function of number of iterations"];

sigdigplot = Table[sigdig[i], {i, 1, (Subscript[n, maximum])}];

<< Graphics`Graphics`

BarChart[sigdigplot];