Calculating an integral by the Romberg Algorithm

In my "Numerical Analysis" course, I learned the Romberg Algorithm to numerically calculate the integral.

The Romberg Algorithm as shown below:

$$T_{2n}(f)+\frac{1}{4^1-1}[T_{2n}(f)-T_{n}(f)]=S_n(f) \\ S_{2n}(f)+\frac{1}{4^2-1}[S_{2n}(f)-S_{n}(f)]=C_n(f) \\ C_{2n}(f)+\frac{1}{4^3-1}[C_{2n}(f)-C_{n}(f)]=R_n(f)$$

where $$T_n=\frac{h}{2}\left[ f(a)+ 2\sum_{i=1}^{n-1}f(x_i) +f(b)\right]$$ and $h=\frac{b-a}{n}$.

My code:

 trapezium[func_, n_, {a_, b_}] :=
With[{h = (b - a)/n},
1/2 h (func@a + 2 Sum[func[a + i h], {i, 1, n - 1}] + func@b)
]

rombergCalc[func_,iter_, {a_, b_}] :=
Module[{m = 1},
Nest[
MovingAverage[#, {-1,4^(m++)}] &,
Table[trapezium[func, 2^i, {a, b}], {i, 0., iter}], 3]
]


The calculation process comes from my textbook

Fixed one bug 1

Update

Fixed bug 2

Test

 rombergCalc[Exp, 5, {0, 1}]//InputForm

 {1.7182818287945305, 1.7182818284603885, 1.7182818284590502}


My Question update

• In function rombergCalc, I utilized the usage m++ that I believe is not suitable in Mathematica Programming. Is there any other method to replace m++ or implement Romberg Algorithm?

- Why Block[{$MinPrecision = precision,$MaxPrecision = precision}..] cannot give the result that contain significance digit that I gave(seeing the graphic of textbook)?

(Thanks for @xzczd's solution for dealing with precision problem)

N[{a, b}, precision]


and replace

trapezium[func, 2^i, {a, b}], {i, 0., iter}]


with

trapezium[func, 2^i, {a, b}], {i, 0, iter}]


- Except for SetOptions[SelectedNotebook[], PrintPrecision -> 16] or InputForm, is there other solutions to set precision conveniently?

• Have you read this ? reference.wolfram.com/language/tutorial/… – Dr. belisarius Nov 14 '14 at 4:10
• I @belisarius haven't read that page. My teacher suggested us implementing that algorithm by a computer language and improve our programming ability by this practice. – xyz Nov 14 '14 at 5:05
• I do not understand why you are dealing with Precision stuff in the implementation of this algorithm. I've implemented this also for a HW assignment but do not have time to post it now. I think you are making things much more complicated than needed with adding Precision in there. – Nasser Nov 14 '14 at 6:49
• @Nasser,This trick come from @RunnyKine's answer A nice trick to force Mathematica to use a given precision is to use Block and make $MinPrecision equal to $MaxPrecision. – xyz Nov 14 '14 at 7:01
• I understand what this is doing. My point is that, you should first work on the algorithm to make sure it works, but eliminating any un-needed extra fancy complexity. You can always fine tune things later if needed. I'll try to look up my HW and post my code. – Nasser Nov 14 '14 at 7:04

Nice to meet you, Mr. Shu.

Bug fix first. Your function doesn't work under desired precision because:

Table[trapezium[func, 2^i, {a, b}], {i, **0.**, iter}]


Changing it to

Table[trapezium[func, 2^i, {a, b}], {i, 0, iter}]


still doesn't fix the problem, because all the numbers taking part in the calculation have infinite precision. Adding an

N[…, precision]


somewhere still doesn't fix the problem, because your utilization of m++ is not only unsuitable, but also wrong. Try changing your m++ into (Print[i = m++]; i) and see what will happen.

Fixed code:

trapezium[func_, n_, {a_, b_}] :=
With[{h = (b - a)/n},
1/2 h (func@a + 2 Sum[func[a + i h], {i, 1, n - 1}] + func@b)]

rombergCalc[func_, iter_, {a_, b_}, precision_] :=
Block[{$MinPrecision = precision,$MaxPrecision = precision},
Module[{m = 1},
NestList[
With[{n = 4^(m++)},
Flatten@(MovingAverage[#, {-1, n}] & /@ Partition[#, 2, 1])] &,
Table[trapezium[func, 2^(i - 1), N[{a, b}, precision]], {i,
iter}], 3]]]


It's so ugly now that I'd like to turn to:

trapezium[func_, n_, {a_, b_}] :=
With[{h = (b - a)/n},
Module[{f = Function[x, func@x, Listable]}, h (Total@f@Range[a, b, h] - 1/2 (f@a + f@b))]]

rombergCalc[func_, iter_, {a_, b_}, precision_] :=
FoldList[Rest@# + 1/(4^#2 - 1) Differences@# &,
trapezium[func, 2^(# - 1), N[{a, b}, precision]] & /@ Range@iter, Range@3]


Finally "visualize" the result:

TableForm[Flatten[rombergCalc[Exp, 6, {0, 1}, 13], {{2}, {1}}],
{"\!$$\*SubscriptBox[\(T$$, $$n$$]\)",
"\!$$\*SubscriptBox[\(S$$, $$n$$]\)",
"\!$$\*SubscriptBox[\(C$$, $$n$$]\)",
"\!$$\*SubscriptBox[\(R$$, $$n$$]\)"}}, TableAlignments -> Center]


For completeness, here's a compiled version of the function above. Notice that it only speeds up Listable compilable internal functions or pure functions formed by Listable compilable internal function.

trapezium[f_, n_, {a_, b_}] :=
With[{h = (b - a)/n},
h (Total[f@Range[a, b, h]] - 1/2 (f@a + f@b))]

compiledrombergCalc[f_, {a_, b_}] :=
ReleaseHold[
Hold@Compile[{{i, _Integer}},
Fold[Rest@# + 1/(4^#2 - 1) Differences@# &,
trapezium[f, 2^(# - 1), {a, b}] & /@ Range@i, Range@3]] /.
DownValues@trapezium]

g = Exp[#] &
NumberForm[compiledrombergCalc[g, {0, 1}][18], 13] // AbsoluteTiming


It's a 23X speedup compared to the uncompiled version.

Also notice that inside the compiled function, the calculation is under MachinePrecision so the Precision of the result is also MachinePrecision, though its appearance is changed by NumberForm. However, I think this treatment may be closer to what your text book has done: as said in the comment below, significance arithmetic is the secret recipe of Mathematica anyway.

• +1, (Seeing mu Update)Good solution! One of my question is that I didn't know how to deal with $$B_n=\frac{4^m}{4^m-1}A_{2n}+\frac{-1}{4^m-1}A_n$$ where, $m=1,2,3$, so I introduced the intermediate variable m and the usage m++(I think it is not suitable for Function Programming) :-) – xyz Nov 14 '14 at 11:48
• My name is 汤书桃, @xzczd so it should be Mr.Tang:D – xyz Nov 14 '14 at 12:00
• @ShutaoTang Oh, I think that besides some minorities, Chinese always put family name ahead. For your new question, you should notice that InputForm and PrintPrecision and NumberForm etc. Only control the display of the number, it doesn't influence Precision nor Accuracy. Nevertheless, calculate with MachinePrecision and control the display with NumberForm etc. may be more close to what your text book has done, significant arithmetic is the secret recipe of Mathematica anyway. – xzczd Nov 14 '14 at 12:46
• OK, thanks a lot. I am alaways having confusion about precison. But I think I will know it gradually by practicing and learning:) – xyz Nov 14 '14 at 13:19
• BTW, you can replace DeleteCases[Transpose@PadRight@rombergCalc[Exp, 6, {0, 1}, 13], 0, 2] with Flatten[rombergCalc[Exp, 6, {0, 1}, 13], {{2}, {1}}].I know that Flatten[lst,{{2},{1}}] can transpose the two-dimensional data that owns different length in the 2nd level. :-) – xyz Nov 15 '14 at 2:00

Here is an old HW assignment. The code is not very functional at all. I even used Break[] in there (OMG!), but I get the same result as the table in the book shows.

I put them side-by-side with the magic of cut/paste:

Code:

nmaRomberg[c_] := Module[{len = Length[c], r, k = 1, f, j, i},
r = Table[0, {len}, {len}];
r[[All, 1]] = c;
Do[
k++;
f = 4^(k - 1);
Do[r[[j, k]] = (f*r[[j + 1, k - 1]] - r[[j, k - 1]])/(f - 1), {j, 1, i - 1}],
{i, len, 1, -1}
];
r
];

nmaTrapz[func_, from_, to_, nStrips_] := Module[{nPoints = nStrips + 1,
h, x, int = 0, f},
h = Abs[from - to]/nStrips;
x = Range[from, to, h];
Do[f = func[x[[i]]];
int += If[i == 1 || i == nPoints, f, 2 f],
{i, 1, nPoints}
];
int*h/2
]


To use:

ClearAll[i, c];
fun = Exp; from = 0; to = 1; nStrips = 32; c = Table[0, {nStrips + 1}];
Do[
c[[i]] = nmaTrapz[fun, from, to, 2^(i - 1)];
If[2^(i - 1) > nStrips, Break[]],
{i, 1, nStrips}];
SetOptions[SelectedNotebook[], PrintPrecision -> 16];
r = nmaRomberg[N[c]][[1 ;; 6, 1 ;; 3]]; (*pull out 3 columns from table*)
Insert[Transpose[r], Table[2^(i - 1), {i, Length[r]}], 1]; (*add n label*)
Grid[Transpose[%], Frame -> All, Spacings -> {1, 1}]

• Reading this code is like looking at a picture of young Nasser and his pacifier. Thanks for sharing ;) – Dr. belisarius Nov 14 '14 at 11:37
• It is first for me to see the usage SetOptions[SelectedNotebook[], PrintPrecision -> 16]. Thanks a lot!:-) – xyz Nov 14 '14 at 11:41
• @belisarius it is not functional because the original HW here and here was in Matlab actually, and I just did now a quick 1:1 translation on the fly. (did not want to change many things so not to break it :) – Nasser Nov 14 '14 at 11:44
• @Nasser Oh! What a pity. I thought you were sharing the recording of your first baby steps with us :) – Dr. belisarius Nov 14 '14 at 11:53
• @ShutaoTang one hr? No, it runs in Matlab just as fast as with Mathematica. I can't imagine why this would take one hr. May be your friend has a bug in the code. They can run the Matlab profiler on it to find where the slow down is. – Nasser Nov 14 '14 at 11:59

Here's a straightforward implementation of the formulas posted in the question.

Clear[t, s, c, r];
t[f_, {a_, b_}, n_] := (b - a)/(2 n) *
ReplacePart[ConstantArray[2, n + 1], {1 -> 1, n + 1 -> 1}] .
(f /@ Rescale[Range[0, n], {0, n}, {a,b}])
s[f_, {a_, b_}, n_] := (4^1 t[f, {a, b}, 2 n] - t[f, {a, b}, n])/(4^1 - 1);
c[f_, {a_, b_}, n_] := (4^2 s[f, {a, b}, 2 n] - s[f, {a, b}, n])/(4^2 - 1);
r[f_, {a_, b_}, n_] := (4^3 c[f, {a, b}, 2 n] - c[f, {a, b}, n])/(4^3 - 1);


The OP's code suggests rather that what is desired is an implementation of Romberg quadrature via Richardson extrapolation of the trapezoidal rule. It's pretty much the same thing, except that instead of using three functions s, c, r to perform the three extrapolation, one can do it iteratively in one function.

Here is a way to implement a Richardson extrapolation of the trapezoidal rule. It adds the new coefficients and abscissas at each step to the previous arrays, so that function evaluations are not repeated. It could be adapted to a recursive method that could stop when some convergence criteria were met. It takes some steps to be efficient. The coefficients of the function values rx[n, m] is a packed array of integers, and abscissas will be a packed array of reals, when using machine precision. Not much checking is included, but for an industrial-strength integrator, use NIntegrate. The test listableQ is an unnecessary enhancement: it assumes vectorization/Listable attributes of the integrand, if it passes a simple test; set listableQ = False &, if there's a problem, or cut it out of the code. It saves time on functions like Exp.

ClearAll[rx, romberg];

rx[0, 0] = ConstantArray[1, 2];
rx[n_, 0] := ReplacePart[ConstantArray[2, 2^n + 1], {1 -> 1, 2 -> 1}];
rx[n_, m_] := Module[{rx0 = 2^(2 m - 1) rx[n, m - 1]},
rx0[[;; 2^(n - 1) + 1]] -= rx[n - 1, m - 1];
rx[n, m] = rx0];

rx[n_, m_, {a_, b_}, "abscissas"] := Fold[Join[#, Rest[#] - (b - a)/2^#2] &,
DeveloperToPackedArray@{a, b}, Range[1, n]];

rx[n_, m_, {a_, b_}, "scale"] :=
Fold[#1*2/(4^#2 - 1) &, (b - a)/2^n/2, Range[m]];

listableQ[f_, x0_] := Quiet@MatchQ[f[{N@x0}], {_Real}];  (* weak test of Listable ? *)
romberg[f_, {a_, b_}, n_, m_, prec_: MachinePrecision] /; n >= m :=
Module[{fvalues},
If[listableQ[f, {x, a}],
fvalues = f[rx[n, m, N[{a, b}, prec], "abscissas"]],
fvalues = f /@ rx[n, m, N[{a, b}, prec], "abscissas"]];
N[rx[n, m, {a, b}, "scale"], prec] rx[n, m].fvalues
];


Comparison:

N[E - 1, 20]
r[Exp, {0.20, 1.20}, 2^(7 - 3)]
romberg[Exp, {0, 1}, 7, 3, 20]
(*
1.7182818284590452354
1.718281828459045235     <-- r loses a little more precision than romberg
1.7182818284590452354
*)


The following shows the equivalence between romberg and t, s, c, and r:

Block[{n = 3},
t[f, {a, b}, 2^(n - 0)] == romberg[f, {a, b}, n, 0, Infinity] // Simplify // Print;
s[f, {a, b}, 2^(n - 1)] == romberg[f, {a, b}, n, 1, Infinity] // Simplify // Print;
c[f, {a, b}, 2^(n - 2)] == romberg[f, {a, b}, n, 2, Infinity] // Simplify // Print;
r[f, {a, b}, 2^(n - 3)] == romberg[f, {a, b}, n, 3, Infinity] // Simplify // Print;
]
(*
True
True
True
True
*)


Here is a quick and dirty trick for obtaining only the top row of the Romberg triangle (which contains the (usually) most accurate estimates of the integral) by exploiting the connexion between Richardsonian extrapolation and the interpolating polynomial. First, we write a simple routine to generate the trapezoidal sums:

trapezoidal[f_, {x_, a_, b_}, n_Integer?Positive] :=
Module[{fn = Function[x, f], h = b - a, m = 1, s},
s = (h/2) (fn[a] + fn[b]);
Table[If[k > 1, h /= 2; m *= 2;
s = s/2 + h Sum[fn[a + j h], {j, 1, m - 1, 2}]];
s, {k, n}]]


Let's use $\int_0^1 \exp(x)\,\mathrm dx$ as the example, just like in the OP:

n = 5;
ts = trapezoidal[N[Exp[x], 20], {x, 0, 1}, n]
{1.8591409142295226177, 1.7539310924648253823, 1.7272219045575167293,
1.7205185921643018614, 1.7188411285799943937}


Now, here's the Richardsonian magic:

pts = Transpose[{4^Range[0, 1 - n, -1], ts}];
rom = Table[InterpolatingPolynomial[Take[pts, k], 0], {k, n}]
{1.8591409142295226177, 1.7188611518765929705, 1.7182826879247574588,
1.7182818287945304232, 1.7182818284590783227}


Compare the relative errors:

correct = E - 1;
Round[-Log10[(ts - correct)/correct]]
{1, 2, 2, 3, 3}

Round[-Log10[(rom - correct)/correct]]
{1, 3, 6, 10, 14}