%% syntaxe : a b {f} epsilon romberg --> une approximation de l'integrale de
%% f (x) entre a et b, calculee avec la methode de Romberg
%% voir http://www.mat.ulaval.ca/anum/ch5/html/node12.html
/romberg {
5 dict begin
/epsil@n exch def
/@f exch def
/@b exch def
/@a exch def
%% initialisation
/i 0 def
/j 0 def
/hi-1 @b @a sub def
/hi 0 def
hi-1 0 ne {
/Ti-1 [ @b @f @a @f add 2 mul hi-1 div ] def
{
/i i 1 add store
/hi hi-1 2 div store
/Ti i 1 add array def
Ti 0
Ti-1 0 get
0
1 1 2 i 1 sub exp {
2 mul 1 sub hi-1 2 div mul
@a add @f
add
} for
hi-1 mul
add
.5 mul
put
{
/j j 1 add store
j i gt {
exit
} if
Ti j
4 j exp
Ti j 1 sub get mul
Ti-1 j 1 sub get sub
4 j exp 1 sub div
put
} loop
/j 0 store
Ti i get Ti-1 i 1 sub get sub abs epsil@n le {
exit
} if
currentdict /Ti-1 undef
/Ti-1 Ti def
currentdict /Ti undef
/hi-1 hi store
} loop
Ti i get
} {
0
} ifelse
end
} def
|