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);

Math

Math
Math
Math

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 

Math

Math

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 

Math

Math

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.