summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/dvips/pst-ode/pst-ode.pro
blob: cc9f66cb4274fd245e7b536ba6c631103a393d27 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PostScript prologue for pst-ode.tex.
% Version 0.2, 2012/09/14
%
% Alexander Grahn (C) 2012
%
% This program can be redistributed and/or modified under the terms
% of the LaTeX Project Public License Distributed from CTAN archives
% in directory macros/latex/base/lppl.txt.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/tx@odeDict 1 dict def
tx@odeDict begin
%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]
  aload pop xlength1 -1 roll {xlength1 -1 roll add} forall
  xlength array astore
} bind def
/subvect { % [1 2 3] [4 5 6] subvect => [-3 -3 -3]
  aload pop xlength1 -1 roll {xlength1 -1 roll sub} forall
  xlength array astore
} bind def
/mulvect { % [1 2 3] 4 mulvect => [4 8 12]
  /mul cvx 2 array astore cvx forall xlength array astore
} bind def
/edivvect { % [1 2 3] [4 5 6] edivvect => [0.25 0.4 0.5]
  aload pop xlength1 -1 roll {xlength1 -1 roll div} forall
  xlength array astore
} bind def
/eabsvect { % [-1 2 -3] eabsvect => [1 2 3]
  {abs} forall xlength array astore
} 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
%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 {
  /t tcur def
  dup ODESET ddt mulvect /k1 exch def
  /t tcur ddt 4 div add def
  dup k1 b21 mulvect addvect ODESET ddt mulvect /k2 exch def
  /t tcur ddt 3 mul 8 div add def
  dup k1 b31 mulvect addvect k2 b32 mulvect addvect ODESET ddt mulvect
    /k3 exch def
  /t tcur ddt 12 mul 13 div add def
  dup k1 b41 mulvect addvect k2 b42 mulvect addvect k3 b43 mulvect addvect
    ODESET ddt mulvect /k4 exch def
  /t tcur ddt add def
  dup k1 b51 mulvect addvect k2 b52 mulvect addvect k3 b53 mulvect addvect
    k4 b54 mulvect addvect ODESET ddt mulvect /k5 exch def
  /t tcur ddt 2 div add def
  dup k1 b61 mulvect addvect k2 b62 mulvect addvect k3 b63 mulvect addvect
    k4 b64 mulvect addvect k5 b65 mulvect addvect ODESET 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]
  %maximum rel. error
  /errmax 0 def
  edivvect {abs dup errmax gt {/errmax exch def}{pop} ifelse} forall errmax
  % => [x(t)] [x(t+ddt) by RKF4] errmax
} 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
  tcur ddt add tout sub tcur ddt add tout sub dt add mul 0 gt
  {/ddt tout tcur sub def} if
  RKF45
  ode@tol div dup 1 gt {
    %failed step -> reduce step size
    exch pop pshrink exp 0.1 max sfty mul ddt mul /ddt exch def
    (-) odeprint tcur ddt add tcur eq {
      (error: step size underflow in ODEINT) print quit
    } if
    ODEINT %repeat step with new ddt
  }{
    %success
    3 -1 roll pop % remove previous state vector
    /tcur tcur ddt add def
    pgrow exp 5 min sfty mul ddt mul /ddt exch def
    %output step completed?
    tcur tout sub dt mul 0 ge {
      (o) odeprint /tcur tout def /tout tout dt add def
    }{
      (+) odeprint ODEINT %continue integration
    } ifelse
  }ifelse
} bind def
/writeresult { %write output vector to file
  {
    20 string cvs statefile exch writestring
    statefile ( ) writestring
  } forall
  statefile (\n) writestring
} bind def
/assembleresult { %assembles state vector for building table of results
  {
    dup (t) eq {
      pop tcur
    }{
      ode@laststate exch get
    } ifelse
  } forall
} bind def
end