%% http://homeomath.imingo.net/msimpson.htm
%% \int_a^b f(x) dx =
%% (b-a)/(3n) * (f (x0) + f (xn) + 4 \sum_{i=1,3,...,n-1} f(xi) + 2\sum_{i = 2,4,...,n-2 } f(xi))
%% voir egalement http://www.mat.ulaval.ca/anum/ch5/html/node10b.html
%% syntaxe : a b {f} n simpson --> une approximation de l'integrale de
%% f (x) entre a et b, calculee avec la methode de Simpson pour n+1
%% points (n est pair)
/simpson {
5 dict begin
/@n exch def
/@f exch def
/@b exch def
/@a exch def
/pas @b @a sub @n div def
0
0 1 @n 1 sub 2 div truncate cvi {
2 mul 1 add %% 2i + 1
pas mul @a add %% x_(2i+1)
@f add %% f (x_(2i+1))
} for
4 mul
0
1 1 @n 1 sub 2 div truncate cvi {
2 mul %% 2i
pas mul @a add %% x_(2i)
@f add %% f (x_(2i))
} for
2 mul
add
@a @f add @b @f add
@b @a sub @n 3 mul div mul
end
} def
|