Romberg Integration.

We have now seen how to develop the Romberg integration method-- we use the trapezoid rule with (2^n) subdivisions as the initial estimator, and then apply Richardson extrapolation to get the best performance possible from this estimator.

Coding Romberg integration

Romberg integration can be defined recursively as follows:

RombergIntegration[f_,{a_,b_},k_] := 
Block[{Romberg,j,n},
Romberg[0,0] = ((b-a)/2)(f[a] + f[b]);
Romberg[n_,0] := Romberg[n,0] = (1/2)Romberg[n-1,0]+((b-a)/(2^(n+1)))((∑_(j=1)^(2^(n-1)))2f[a+(2j-1)(((b-a)/(2^n)))]);
Romberg[n_,m_] :=Romberg[n,m] =  ((4^m)/((4^m)-1))Romberg[n,m-1] - (1/((4^m)-1))Romberg[n-1,m-1];
Romberg[k,k-1]
]

Coding notes

You may have noticed that we used a weird-looking construction above, in the form

f[x_] := f[x] = (something)

The effect of this construction is to produce the same result as (something), but also to permanently define f[x] as that value as a side effect.  The reason to do this is that it caches the result of the computation for the rest of the run (or the rest of the enclosing Block or Module). After the first call, repeating the calculation is essentially instantaneous and doesn’t rerun the original code.

Here’s a simple example:

f[x_] := Sum[Sin[1.0 n^2 + 3x],{n,1,1000000}]

Running this with x = 3.0 takes my laptop a little more than one second:

Timing[f[3.0]]

and that’s about the same if I run it again (given that my laptop is also backing itself up to an external disk, downloading my email, and playing an old Arcade Fire album as I write this):

Timing[f[3.0]]

On the other hand, if I define

Clear[g]

g[x_] := g[x] = Sum[Sin[1.0 n^2 + 3x],{n,1,1000000}]

Then the first call takes about the same amount of time

Timing[g[3.0]]

But if I run it again with the same input, it’s basically instant:

Timing[g[3.0]]

Coding the trapezoid rule

We now recover our previous trapezoid rule code, but since we’re comparing with Romberg, 
I’ll alter it to make the number of subdivisions a power of 2:

Trap[f_,{a_,b_},k_] := ((b-a)/(2^(k+1)))(f[a] + (∑_(i=1)^((2^k)-1))2f[a+i(((b-a)/(2^k)))] + f[b])

Comparing Trap and Romberg

We’re now going to compare the performance of the trapezoid rule and Romberg integration on some of our favorite integrands:

TrapezoidError[f_,{a_,b_},k_]:= 
Block[{exact,x,h,$MaxExtraPrecision=20000},
exact = Integrate[f[x],{x,a,b}];
h = ((b-a)/(2^k)); 
{N@Log10[h],N@Log10[Abs[exact-Trap[f,{N[a,1000],N[b,1000]},k]]]}
]

RombergError[f_,{a_,b_},k_]:= 
Block[{exact,x,h,$MaxExtraPrecision=20000},
exact = Integrate[f[x],{x,a,b}];
h = ((b-a)/(2^k)); 
{N@Log10[h],N@Log10[Abs[exact-RombergIntegration[f,{N[a,1000],N[b,1000]},k]]]}
]

We’ll start with the integral

(∫_0^π)Sin[x]\[DifferentialD]x

and gather data from (2^2) intervals to (2^7) intervals:

TrapezoidErrorData1 = Table[TrapezoidError[Sin,{0,π},k],{k,2,7}]

RombergErrorData1 = Table[RombergError[Sin,{0,π},k],{k,2,7}]

We know that the Trapezoid rule error is O((h^(2 ))), so the log-log data should fit to a line with slope 2:

lm = LinearModelFit[TrapezoidErrorData1,lgh,lgh][lgh]

However, the error order of the Romberg data depends on the number of subdivisions, and improves as the number of subdivisions increases. It turns out that the log-log data fits very well to a quadratic function:

Rm = NonlinearModelFit[RombergErrorData1,a + b lgh + c lgh^2,{a,b,c},lgh][lgh]

Now we can plot the two sets of error data together:

Show[{Plot[{lm,Rm},{lgh,-2,0}],
ListPlot[{TrapezoidErrorData1,RombergErrorData1},]},]

This is amazing when you think about it: Romberg integration with only (2^7)=128 intervals has already surpassed full-precision in double arithmetic!

The total picture: error including roundoff.

We have deliberately suppressed roundoff error in the calculation above so that we could see only the approximation error. We know that the computations involved in the trapezoid rule are numerically stable (certainly for this few intervals!), but we don’t know the same about Romberg. Perhaps, we might worry, the subtractions involved in Richardson extrapolation will cause us numerical problems.

We now take a look at these effects, recomputing everything without additional precision to see the total error (roundoff error plus approximation error):

MachineTrapezoidError[f_,{a_,b_},k_]:= 
Block[{exact,x,h},
exact = 1.0Integrate[f[x],{x,a,b}];
h = (((1.0b)-(1.0a))/(2.0^k)); 
{1.0Log10[h],1.0Log10[Abs[exact-Trap[f,{1.0a,1.0b},k]]]}
]

MachineRombergError[f_,{a_,b_},k_]:= 
Block[{exact,x,h},
exact = 1.0Integrate[f[x],{x,a,b}];
h = (((1.0b)-(1.0a))/(2.0^k)); 
{1.0Log10[h],1.0Log10[Abs[exact-RombergIntegration[f,{1.0a,1.0b},k]]]}
]

and gather data from (2^2) intervals to (2^8) intervals:

MachineTrapezoidErrorData1 = Table[MachineTrapezoidError[Sin,{0,π},k],{k,2,8}]

MachineRombergErrorData1 = Table[MachineRombergError[Sin,{0,π},k],{k,2,8}]

We won’t try to fit error to theoretical functions for the machine errors, we’ll just connect the dots in the plot.

Show[{Plot[-16,{lgh,-2,0},PlotLegends->{"Min error (double precision)"}],
ListLinePlot[{TrapezoidErrorData1,RombergErrorData1,MachineTrapezoidErrorData1,MachineRombergErrorData1},PlotMarkers->{Automatic,Small},PlotLegends->{"Trapezoid approx. error","Romberg approx. error","Trapezoid total error","Romberg total error"}]},]

This should reassure us that Romberg isn’t losing anything to roundoff error-- it really is that much better than the trapezoid rule!

Computation speed.

For this example, the number of computations involved is so small that the computation is basically instant. However, we can see that Romberg integration takes only a little more time than the Trapezoid rule itself. (This is largely due to the way we coded it.)

Timing[Trap[Sin,{0.0,1.0π},7]]

Timing[RombergIntegration[Sin,{0.0,1.0π},7]]

For many cases of numerical integration, it’s evaluating the original function which is the time-consuming part, either because you are actually reading values from a sensor or because the function is the result of running some other complicated code. So we can also measure time in terms of function evaluations.

So it’s worthwhile to count function evaluations. Now you might think that we should do this by rewriting the Trap and RombergIntegration functions and incrementing a counter whenever f is evaluated, or perhaps by redefining f to increment a global variable every time it is called.

Both of these approaches have drawbacks: in the first, it’s true here that we wrote the Trap and RombergIntegration functions, so we can alter their code. But that might not be the case, especially if we’re using a built-in function like NIntegrate. In the second, globals are inherently awkward in programming. We either risk overwriting a variable from somewhere else in our notebook or have to limit the scope of the variable by wrapping everything in a Block or Module. Mathematica provides an easier way to do this with Sow and Reap.

Consider the code:

Table[Sin[x],{x,1,10}]

which evaluates Sin 10 times. If, every time we evaluate Sin, we Sow a “1”, the results are the same,

Table[(Sow[1]; Sin[#])&[x],{x,1,10}]

but if we enclose the whole thing in a Reap, we get the results together with a List of everything we Sow’d during the enclosed code:

Reap[Table[(Sow[1]; Sin[#])&[x],{x,1,10}]]

The last step to counting is simply to take the length of the list in the second argument

Length[Reap[Table[(Sow[1]; Sin[#])&[x],{x,1,10}]][[2,1]]]

Now we do it for our code! And we see that in each case, the number of function evaluations is the same:

Length[Reap[N@Trap[(Sow[1];Sin[#])&,{0,π},7]][[2,1]]]

Length[Reap[N@RombergIntegration[(Sow[1];Sin[#])&,{0,π},7]][[2,1]]]

This means that when the time spent evaluating the integrand is dominant, the two methods are basically the same speed.