Romberg Integration.

This demonstration shows some examples of computing numerical integrals via the trapezoid rule. We see some interesting phenomena.

We first reproduce our Trapezoid rule code from a while back, then follow it with a Romberg algorithm. Our Romberg code is not super-efficient, but it avoids needless recomputing of previous values, at least. 

(One could write a short and sweet recursive Romberg algorithm, but at the cost of re-evaluating the Trapezoid rule terms far too many times.)

Trap[{f_,a_,b_,n_}] := Module[ {h},

h = (b - a)/n; 
h*((1/2)*(f[b] + f[a]) + Sum[f[a + i*h],{i,1,n-1}])

];

RombergArray[{f_,a_,b_,NN_,digits_}] := Module[{A,RA,m,i,n},

(* Code to generate an nxn array of R(n,m) terms. *)

RA = ConstantArray[0.0,{NN,NN}];

(* Mathematica arrays are indexed by 1--which is annoying from our point of view. We solve this by computing the 0th column of the array containing the R(n,0) terms obtained from the Trapezoid rule separately *)

R0 = Table[N[Trap[{f,a,b,2^(i-1)}],digits],{i,1,NN}];

(* Now we fill in the first column of the array by hand *)
	
RA[[2;;NN,1]] = 
Table[
(4/3)*R0[[n]] - (1/3)*R0[[n-1]],
{n,2,NN}];

(* The funny [[ 2 ;; N, 1]] means "the submatrix from row 2 to row N and column 1" *)

(* Now we can fill the remaining columns with a loop *)

Do[

RA [[m+1;;NN,m]] = 

Table[

(4^m/(4^m -1))*RA[[n,m-1]] - (1/(4^m - 1))*RA[[n-1,m-1]],

{n,m+1,NN}] ,

{m,2,NN}];

(* We now want to shift the columns of the array. After all,    the 0th column is stored in R0 and the Nth column is empty. *)

Do[ RA[[All,m]]= RA[[All,m-1]],{m,NN,2,-1}];
RA[[All,1]] = R0;

RA

];

Romberg[{f_,a_,b_,NN_,digits_}]  := RombergArray[{f,a,b,NN,digits}][[NN,NN]];

f[x_] := Sin[x];
Grid[RombergArray[{f,0,Pi,5,10}] ]

Romberg[{f,0,Pi,5,10}]

We now compare performance of the Romberg algorithm with varying numbers of steps, just to enjoy the superior performance of this method!

Digits[a_,b_] := N[Abs[Log[10,Abs[a -b]]]];

Clear[f];
Clear[Ans];

f[x_]:= Sin[x];
Ans = N[Integrate[f[x],{x,0,Pi}],100]

RomSinDigits = Table[Digits[Romberg[{f,0,Pi,i,100}],Ans],{i,1,10}]
TrapSinDigits =Table[Digits[N[Trap[{f,0,Pi,2^(i-1)}]],Ans],{i,1,10}]

ListLinePlot[{RomSinDigits,TrapSinDigits},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0}, PlotMarkers->Automatic]

How many function evaluations would it take to get to full machine precision with the Romberg and Trapezoid rules? Well, we try

Digits[Trap[{f,0,Pi,2^(15-1)}],2]

It's clear that we're not going to want to wait that long. How about we use PredictedDigits to guesstimate?

PredictedDigits[{M_,a_,b_,n_}] := N[Abs[Log[10,Abs[(1/12)*(b-a)^3/n^2]]]];

PredictedDigits[{1,0,Pi,2^28}]

There you have it: 2^28 function evaluations to compute the integral to full machine precision with the trapezoid rule vs 2^7 function evaluations with the Romberg method. That's the awesome power of Richardson extrapolation at work!