%% syntaxe : a f x0 y0 h basemilne --> depose les points, calcules par
%% la methode de Milne, de la courbe sur
%% [x0, a] de la fonction h solution de l'equa diff y' = f (x, y)
%% verifiant h (x0) = y0
%% Plus precisement :
%% predicteur : tyn+1 = yn-3 + (4h/3) (2y'n - y'n-1 + 2y'n-2)
%% correcteur : yn+1 = yn-1 + (h/3) (ty'n+1 + 4 y'n + y'n-1)
%% ou y'n = f (xn, yn) et ty'n+1 = f (xn+1, tyn+1)
%% attention : on peut avoir a < x0, mais dans ce cas h doit etre negatif
%% autre syntaxe : a f x0 y0 n basemilne ou n est le nombre de points
%% calcules. Il y aura alors n+1 points deposes sur la pile : (x0, y0)
%% plus les n suivants
/basemilne {
14 dict begin
/v@r exch def
/y0 exch def
/x0 exch def
/@@@f exch def %% y'= @@@f (x, y)
/a exch def
/tyn+1 0 def %% valeur temporaire de yn+1
/yn+1 0 def
v@r isinteger { %% si v@r est entier
/n v@r def %% c'est que c'est le nombre de points
/h a x0 sub n div def
} { %% sinon, v@r n'est pas entier
/h v@r def %% et c'est le pas
/n a x0 sub h div truncate cvi def
} ifelse
%% sur [x0, a]
x0 3 h mul add {@@@f} x0 y0 3 baserungekutta
/yn exch def
/xn exch def
/yn-1 exch def pop
/yn-2 exch def pop
/yn-3 exch def pop
x0 y0
x0 h add yn-2
x0 2 h mul add yn-1
x0 3 h mul add yn
n 3 sub {
/tyn+1
yn-3
4 h mul 3 div
xn yn @@@f 2 mul
xn h sub yn-1 @@@f sub
xn 2 h mul sub yn-2 @@@f 2 mul add
mul add
store
/yn+1
yn-1 h 3 div
xn h add tyn+1 @@@f
xn yn @@@f 4 mul add
xn h sub yn-1 @@@f add
mul add
store
%% ok, maintenant on translate
/yn-3 yn-2 store
/yn-2 yn-1 store
/yn-1 yn store
/yn yn+1 store
/xn xn h add store
xn yn+1
} repeat
end
} def
%% syntaxe : a b f x0 y0 h Milne --> depose les points, calcules par
%% la methode de Milne, de la courbe sur
%% [a, b] de la fonction h solution de l'equa diff y' = f (x, y)
%% verifiant h (x0) = y0
%% Plus precisement :
%% predicteur : tyn+1 = yn-3 + (4h/3) (2y'n - y'n-1 + 2y'n-2)
%% correcteur : yn+1 = yn-1 + (h/3) (ty'n+1 + 4 y'n + y'n-1)
%% ou y'n = f (xn, yn) et ty'n+1 = f (xn+1, tyn+1)
/Milne {
6 dict begin
/@h exch def %% le pas, ou le nb de points
/y@ exch def
/x@ exch def
/@@f@ exch def %% y'= @f (x, y)
/@b exch def
/@a exch def
%% sur [a, x0]
@a x@ lt {
[
@a {@@f@} x@ y@ @h dup isinteger not {neg} if basemilne
]
reversep aload pop
x@ @b lt {
pop pop
} if
} if
%% sur [x0, b]
x@ @b lt {
@b {@@f@} x@ y@ @h basemilne
} if
end
} def
%% syntaxe : f x0 y0 h milne --> depose les points, calcules par
%% la methode de Milne, de la courbe sur [xmin, xmax] de la fonction h solution
%% de l'equa diff y' = f (x, y) verifiant h (x0) = y0
/milne {
xmax 5 1 roll
xmin 6 1 roll
Milne
} def
%% %% syntaxe : a f x0 y0 h basemilne --> depose les points, calcules par
%% %% la methode de Milne, de la courbe sur
%% %% [x0, a] de la fonction h solution de l'equa diff y' = f (x, y)
%% %% verifiant h (x0) = y0
%% %% Plus precisement :
%% %% predicteur : tyn+1 = yn-3 + (4h/3) (2y'n - y'n-1 + 2y'n-2)
%% %% correcteur : yn+1 = yn-1 + (h/3) (ty'n+1 + 4 y'n + y'n-1)
%% %% ou y'n = f (xn, yn) et ty'n+1 = f (xn+1, tyn+1)
%% %% attention : on peut avoir a < x0, mais dans ce cas h doit etre negatif
%% /basemilne {
%% 6 dict begin
%% /h exch def %% le pas
%% /y0 exch def
%% /x0 exch def
%% /@@@f exch def %% y'= @@@f (x, y)
%% /a exch def
%% /tyn+1 0 def %% valeur temporaire de yn+1
%% /yn+1 0 def
%%
%% %% sur [x0, a]
%% x0 3 h mul add {@@@f} x0 y0 3 baserungekutta
%% /yn exch def
%% /xn exch def
%% /yn-1 exch def pop
%% /yn-2 exch def pop
%% /yn-3 exch def pop
%%
%% x0 y0
%%
%% x0 h add yn-2
%% x0 2 h mul add yn-1
%% x0 3 h mul add yn
%%
%% a x0 3 h mul add sub h div truncate cvi {
%% /tyn+1
%% yn-3
%% 4 h mul 3 div
%% xn yn @@@f 2 mul
%% xn h sub yn-1 @@@f sub
%% xn 2 h mul sub yn-2 @@@f 2 mul add
%% mul add
%% store
%% /yn+1
%% yn-1 h 3 div
%% xn h add tyn+1 @@@f
%% xn yn @@@f 4 mul add
%% xn h sub yn-1 @@@f add
%% mul add
%% store
%%
%% %% ok, maintenant on translate
%% /yn-3 yn-2 store
%% /yn-2 yn-1 store
%% /yn-1 yn store
%% /yn yn+1 store
%% /xn xn h add store
%% xn yn+1
%% } repeat
%%
%% end
%% } def
%%
%% %% syntaxe : a b f x0 y0 h Milne --> depose les points, calcules par
%% %% la methode de Milne, de la courbe sur
%% %% [a, b] de la fonction h solution de l'equa diff y' = f (x, y)
%% %% verifiant h (x0) = y0
%% %% Plus precisement :
%% %% predicteur : tyn+1 = yn-3 + (4h/3) (2y'n - y'n-1 + 2y'n-2)
%% %% correcteur : yn+1 = yn-1 + (h/3) (ty'n+1 + 4 y'n + y'n-1)
%% %% ou y'n = f (xn, yn) et ty'n+1 = f (xn+1, tyn+1)
%% /Milne {
%% 6 dict begin
%% /@h exch def %% le pas
%% /y@ exch def
%% /x@ exch def
%% /@@f@ exch def %% y'= @f (x, y)
%% /@b exch def
%% /@a exch def
%%
%% %% sur [a, x0]
%% @a x@ lt {
%% [
%% @a {@@f@} x@ y@ @h neg basemilne
%% ]
%% reversep aload pop
%% x@ @b lt {
%% pop pop
%% } if
%% } if
%%
%% %% sur [x0, b]
%% x@ @b lt {
%% @b {@@f@} x@ y@ @h basemilne
%% } if
%%
%% end
%% } def
%%
%% %% syntaxe : f x0 y0 h milne --> depose les points, calcules par
%% %% la methode de Milne, de la courbe sur [xmin, xmax] de la fonction h solution
%% %% de l'equa diff y' = f (x, y) verifiant h (x0) = y0
%% /milne {
%% xmax 5 1 roll
%% xmin 6 1 roll
%% Milne
%% } def
|