%%% rk4 ------------------------------------------------------------------------
% y n h t rk4 --> y_{0} y_{1} ... y_{n-1}
% y : tableau
% n : taille de y
% h : pas de temps
% t : temps
%%%
% La procédure SM (second membre) doit être définie
/rk4Dict 9 dict def
/rk4 {
rk4Dict begin
/t exch def % La variable
/h exch def % Le pas
/n exch def % Le nombre de fonctions
/y exch def % Le tableau de fonctions
/k1 y aload pop % k1 = h * SM(t, y)
t SM n array astore { h mul } forall n array astore def
/t t h 2 div add def
/k2 0 1 n 1 sub % k2 = h * SM(t + h / 2, y + k1 / 2)
{ dup y exch get exch k1 exch get 2 div add } for
t SM n array astore { h mul } forall n array astore def
/k3 0 1 n 1 sub % k3 = h * SM(t + h / 2, y + k2 / 2)
{ dup y exch get exch k2 exch get 2 div add } for
t SM n array astore { h mul } forall n array astore def
/t t h 2 div add def
/k4 0 1 n 1 sub % k4 = h * SM(t + h, y + k3)
{ dup y exch get exch k3 exch get add } for
t SM n array astore { h mul } forall n array astore def
/i 0 def
y { % RETOUR : y + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6
k1 i get 6 div add
k2 i get 3 div add
k3 i get 3 div add
k4 i get 6 div add
/i i 1 add def
} forall
end
} def
%%% rk4ap ----------------------------------------------------------------------
%%%
% La variable EPSILON (le seuil de décision) doit être définie
% Ce codage est largement inspiré des Numerical Recipes en fortran.
/rk4apDict 10 dict def
rk4apDict begin
/CONS 4 0.9 div -5 exp def % La consistance
end
/rk4ap {
rk4apDict begin
/t exch def
/h exch def
/n exch def
/y exch def
/hh h 2 div def
% Approximation de y(t + h) en passant par celle de y(t + h / 2)
/yA y n hh t rk4 n array astore n hh t hh add rk4 n array astore def
% Approximation directe de y(t + h)
/yB y n h t rk4 n array astore def
% Mesure de la distance entre les deux approximations (norme inf)
/i 0 def
0 yA { yB i get sub abs
2 copy gt {pop} {exch pop} ifelse
/i i 1 add def
} forall
% On fait le rapport avec le seuil
EPSILON div dup
1 gt {
% L'écart est trop important entre les deux, on recommence avec un pas
% estimé suffisant
/hn exch -0.25 exp h mul 0.9 mul def
y n hn t rk4ap
} {
dup
CONS gt {
% Diminution du pas pour tomber sous la consistance - accord avec
% le fait que la méthode est d'ordre 5 , l'erreur est équivalente
% à alpha x h^5, (1/5 = 0.2)
/hn exch -0.20 exp h mul 0.9 mul def
} {
% On est sous la consistance, on peux multiplier le pas par 4
pop
/hn h 4 mul def
} ifelse
/i 0 def
% On profite du fait que deux valeurs sont présentes pour gagner un
% ordre en les combinants pour éliminer l'ordre 5
yA { 16 mul yB i get sub 15 div /i i 1 add def } forall
n hn t h add
} ifelse
end
} def
%%% rk4ap ----------------------------------------------------------------------
|