%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PostScript prologue for pst-ode.tex.
% Version 0.19, 2024/01/04
%
% Alexander Grahn (C) 2012--today
%
% This program can be redistributed and/or modified under the terms
% of the LaTeX Project Public License
%
% http://www.latex-project.org/lppl.txt
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/tx@odeDict 1 dict def
tx@odeDict begin
/ode@@dict 1 dict def
/ode@dict {ode@@dict begin} bind def
ode@dict
 %some constants for step size calculation
 /sfty 0.9 def /pgrow -0.2 def /pshrink -0.25 def
 %helper functions
 /addvect { % [1 2 3] [4 5 6] addvect => [5 7 9]
   ode@dict
     aload pop xlength1 -1 roll {xlength1 -1 roll add} forall
     xlength array astore
   end
 } bind def
 /subvect { % [1 2 3] [4 5 6] subvect => [-3 -3 -3]
   ode@dict
     aload pop xlength1 -1 roll {xlength1 -1 roll sub} forall
     xlength array astore
   end
 } bind def
 /mulvect { % [1 2 3] 4 mulvect => [4 8 12]
   ode@dict /mul cvx 2 array astore cvx forall xlength array astore end
 } bind def
 /edivvect { % [1 2 3] [4 5 6] edivvect => [0.25 0.4 0.5]
   ode@dict
     aload pop xlength1 -1 roll {xlength1 -1 roll div} forall
     xlength array astore
   end
 } bind def
 /eabsvect { % [-1 2 -3] eabsvect => [1 2 3]
   ode@dict {abs} forall xlength array astore end
 } bind def
 %/revstack { % (a) (b) (c) (d) 3 revstack => (a) (d) (c) (b)
 %  -1 2 {dup 1 sub neg roll} for
 %} bind def
 /min { 2 copy gt { exch } if pop } bind def
 /max { 2 copy lt { exch } if pop } bind def
 %coefficient table (Butcher table) of RKF45
 /c2 1 4 div def /c3 3 8 div def /c4 12 13 div def /c6 1 2 div def
 /a21 1 4 div def /a31 3 32 div def /a32 9 32 div def
 /a41 1932 2197 div def /a42 7200 2197 div neg def /a43 7296 2197 div def
 /a51 439 216 div def /a52 8 neg def /a53 3680 513 div def
 /a54 845 4104 div neg def /a61 8 27 div neg def /a62 2 def
 /a63 3544 2565 div neg def /a64 1859 4104 div def /a65 11 40 div neg def
 /b1 16 135 div def /b3 6656 12825 div def
 /b4 28561 56430 div def /b5 9 50 div neg def
 /b6 2 55 div def
 /b1* 25 216 div def /b3* 1408 2565 div def
 /b4* 2197 4104 div def /b5* 1 5 div neg def
end
%Runge-Kutta-Fehlberg (RKF45) method
%performs one integration step over tentative step size ddt
%[state vector x(t)] RKF45 => [x(t)] [x(t+ddt) by RKF4] errmax
/RKF45 {
 dup
 /t ode@dict tcur end def
 ODESET
 ode@dict
   ddt mulvect /k1 exch def
   dup k1 a21 mulvect addvect
   tcur ddt c2 mul add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k2 exch def
   dup k1 a31 mulvect addvect k2 a32 mulvect addvect
   tcur ddt c3 mul add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k3 exch def
   dup k1 a41 mulvect addvect k2 a42 mulvect addvect k3 a43 mulvect addvect
   tcur ddt c4 mul add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k4 exch def
   dup k1 a51 mulvect addvect k2 a52 mulvect addvect k3 a53 mulvect addvect
   k4 a54 mulvect addvect
   tcur ddt add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k5 exch def
   dup k1 a61 mulvect addvect k2 a62 mulvect addvect k3 a63 mulvect addvect
   k4 a64 mulvect addvect k5 a65 mulvect addvect
   tcur ddt c6 mul add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k6 exch def % => [x(t)]
   %fourth order solution (increment dx)
   dup dup k1 b1* mulvect k3 b3* mulvect addvect k4 b4* mulvect addvect
     k5 b5* mulvect addvect dup
     % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [dx by RKF4]
   %fifth order solution (abs. error)
   k1 b1 mulvect k3 b3 mulvect addvect k4 b4 mulvect addvect
     k5 b5 mulvect addvect k6 b6 mulvect addvect subvect
     % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [err]
   5 1 roll addvect 4 -2 roll % => [x(t)] [x(t+ddt) by RKF4] [err] [x(t)]
   %scaling vector for step size adjustment (Numerical Recipies)
   eabsvect k1 eabsvect addvect {1e-30 add} forall xlength array astore
   % => [x(t)] [x(t+ddt) by RKF4] [err] [xscale]
   edivvect eabsvect 1e-30 exch {max} forall %maximum rel. error
   % => [x(t)] [x(t+ddt) by RKF4] errmax
 end
} bind def
%classical Runge-Kutta (RK4) method
%one integration step over step size ddt; no error estimate
%[state vector x(t)] RK4 => [x(t+ddt) by RK4]
/RK4 {
 dup
 ode@dict tcur end /t exch def
 ODESET
 ode@dict
   ddt mulvect /k1 exch def
   dup k1 0.5 mulvect addvect
   tcur ddt 2 div add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k2 exch def
   dup k2 0.5 mulvect addvect
 end
 ODESET
 ode@dict
   ddt mulvect /k3 exch def
   dup k3 addvect
   tcur ddt add
 end
 /t exch def
 ODESET
 ode@dict
   ddt mulvect /k4 exch def % => [x(t)]
   k1 k2 k3 addvect 2 mulvect addvect k4 addvect 1 6 div mulvect addvect
   % => [x(t+ddt) by RK4]
 end
} bind def
/ODEINT { % performs integration over output step [t,t+dt]
         % [ state vector x(t) ] ODEINT => [ state vector x(t+dt) ]
 rk4 {
   RK4 (o) odeprint
   ode@dict
     /tcur tout def
     /outStepCnt outStepCnt 1 add def
     /tout tStart dt outStepCnt mul add def
   end
 }{
   %decrease overshooting step size
   ode@dict
     tout tcur sub dup abs ddt abs lt {/ddt exch def}{pop} ifelse
   end
   RKF45
   ode@tol div dup 1 gt {
     %failed step -> reduce step size
     ode@dict
       exch pop pshrink exp 0.1 max sfty mul ddt mul /ddt exch def
       tcur ddt add tcur eq {
         % error: step size underflow in ODEINT
         (! t=) odeprint tcur odeprint
         % print & remove previous state vector
         (, x=[) odeprint /ode@spc () def
         {ode@spc odeprint /ode@spc ( ) def odeprint} forall (]) odeprint
         true
       }{
         (-) odeprint
         false
       } ifelse
     end
     % on step size underflow, leave loop over output steps (pst-ode.tex)
     {exit} if
     ODEINT %repeat step with new ddt
   }{
     %success
     3 -1 roll pop % remove previous state vector
     ode@dict
       /tcur tcur ddt add def
       dup 0 ne {pgrow exp 5 min sfty mul ddt mul /ddt exch def}{pop} ifelse
       %output step completed?
       tcur tout sub
     end
     0 eq {
       (o) odeprint
       ode@dict
         /tcur tout def
         /outStepCnt outStepCnt 1 add def
         /tout tStart dt outStepCnt mul add def
       end
     }{
       (+) odeprint ODEINT %continue integration
     } ifelse
   } ifelse
 } ifelse
} bind def
/writeresult { %write output vector to file
 /loopproc load forall
 statefile (\n) writestring
} bind def
/loopproc {
 0 cvs statefile exch writestring
 statefile ( ) writestring
} bind def
/loopproc load 0 256 string put
/assembleresult { %assembles state vector for building table of results
 {
   dup (t) eq {
     pop ode@dict tcur end
   }{
     ode@laststate exch get
   } ifelse
 } forall
} bind def
end

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.