Supplement: Romberg Quadrature re-implemented
This section is primarily intended to exercise our coding with Maple. It has no direct connection to our exploration of Euler's 'Pi-equation'.
TangentTrapez := proc (ig::procedure, ab::range,
dig::integer)
global a, b, h, f, tangentSum, reqDig;
f := ig; a := op(1,ab);
b := op(2,ab); h := b-a;
reqDig := dig;
tangentSum := proc (subdiv::integer)
local k, p, H;
p := 2^subdiv;
H := h/p;
evalf(H*sum(f(a+H*(k+1/2)),k=0..p-1),reqDig) end;
tangentSum end:
tangentSin := TangentTrapez(x -> sin(x), 0..Pi, 16):
seq(tangentSin(i),i=0..4);
RombergTangent(x->sin(x),0..Pi,4,16,2);
i T[i]
0 +3.1415926536
1 +2.2214414691 +1.9147244076
2 +2.0523443060 +1.9959785849 +2.0013955301
3 +2.0129090856 +1.9997640121 +2.0000163740 +1.9999944826
4 +2.0032163782 +1.9999854757 +2.0000002399 +1.9999999838 +2.0000000054
RombergTangent(x -> 2/(x^2+1),-1..1,5,16,2): 'Err'=evalf(Pi-%,16);
This differs from Richardson(piTangent) only in the first term, because piTangent suppresses the first approximation '4'.
i T[i]
0 +4.0000000000
1 +3.2000000000 +2.9333333333
2 +3.1623529412 +3.1498039216 +3.1642352941
3 +3.1468005184 +3.1416163775 +3.1410705412 +3.1407028467
4 +3.1428947296 +3.1415928000 +3.1415912282 +3.1415994930 +3.1416030093
5 +3.1419181743 +3.1415926559 +3.1415926463 +3.1415926688 +3.1415926420 +3.1415926319
I can't help, but it is much more easy for me to think about Romberg integration as applying two operators one after the other, first an integral-operator, then an extrapolation-operator. But what is more natural than to code the procedures in correspondence to this mental picture? Is it only for some irrelevant small performance gain that most implementations mix these two operations together, making a reuse in a new combination impossible? (To be faire to the above implementation of W. Werner, we add that our implementation uses a feature of Maple, which was not available prior to Release 5).
And the last section shows, that this approach does pay off. At some point in our development we could switch from the tangent approximation 'piTangent' to a simplified form of this function without even thinking about, whether we are still performing an integration, and yet kept getting meaningful results. In fact all we need is to input a sequence, which has an asymptotic behavior, the Richardson procedure is designed to cope with.