pst-ode.pro

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PostScript prologue for pst-ode.tex.
% Version 0.11, 2017/08/15
%
% Alexander Grahn (C) 2012--today
%
% This program can be redistributed and/or modified under the terms
% of the LaTeX Project Public License Distributed from CTAN archives:
% http://mirrors.ctan.org/tex-archive/macros/latex/base/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
  /a41 25 216 div def /a43 1408 2565 div def
  /a44 2197 4104 div def /a45 1 5 div neg def
  /a51 16 135 div def /a53 6656 12825 div def
  /a54 28561 56430 div def /a55 9 50 div neg def
  /a56 2 55 div def
  /b21 1 4 div def /b31 3 32 div def /b32 9 32 div def
  /b41 1932 2197 div def /b42 7200 2197 div neg def /b43 7296 2197 div def
  /b51 439 216 div def /b52 8 neg def /b53 3680 513 div def
  /b54 845 4104 div neg def /b61 8 27 div neg def /b62 2 def
  /b63 3544 2565 div neg def /b64 1859 4104 div def /b65 11 40 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
  ode@dict tcur end /t exch def
  ODESET
  ode@dict
    ddt mulvect /k1 exch def
    dup k1 b21 mulvect addvect
    tcur ddt 4 div add
  end
  /t exch def
  ODESET
  ode@dict
    ddt mulvect /k2 exch def
    dup k1 b31 mulvect addvect k2 b32 mulvect addvect
    tcur ddt 3 mul 8 div add
  end
  /t exch def
  ODESET
  ode@dict
    ddt mulvect /k3 exch def
    dup k1 b41 mulvect addvect k2 b42 mulvect addvect k3 b43 mulvect addvect
    tcur ddt 12 mul 13 div add
  end
  /t exch def
  ODESET
  ode@dict
    ddt mulvect /k4 exch def
    dup k1 b51 mulvect addvect k2 b52 mulvect addvect k3 b53 mulvect addvect
    k4 b54 mulvect addvect
    tcur ddt add
  end
  /t exch def
  ODESET
  ode@dict
    ddt mulvect /k5 exch def
    dup k1 b61 mulvect addvect k2 b62 mulvect addvect k3 b63 mulvect addvect
    k4 b64 mulvect addvect k5 b65 mulvect addvect
    tcur ddt 2 div add
  end
  /t exch def
  ODESET
  ode@dict
    ddt mulvect /k6 exch def % => [x(t)]
    %fourth order solution (increment dx)
    dup dup k1 a41 mulvect k3 a43 mulvect addvect k4 a44 mulvect addvect
      k5 a45 mulvect addvect dup
      % => [x(t)] [x(t)] [x(t)] [dx by RKF4] [dx by RKF4]
    %fifth order solution (abs. error)
    k1 a51 mulvect k3 a53 mulvect addvect k4 a54 mulvect addvect
      k5 a55 mulvect addvect k6 a56 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
/ODEINT { % performs integration over output step [t,t+dt]
          % [ state vector x(t) ] ODEINT => [ state vector x(t+dt) ]
  %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
      ode@dict tcur ddt add tcur end 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
      tcur tout sub %output step completed?
    end
    0 ge {
      (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
} 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