blob: 28d0550dda2a136544703f91b974d2703d91ad63 (
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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PostScript prologue for pst-ode.tex.
% Version 0.10, 2017/06/16
%
% 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
% in directory 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]
/errmax 1e-30 def %maximum rel. error
edivvect {abs dup errmax gt {/errmax exch def}{pop} ifelse} forall errmax
% => [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
tcur ddt add tout sub tcur ddt add tout sub dt add mul 0 gt
{/ddt tout tcur sub def} if
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 {odeprint ( ) 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 dt mul
end
0 ge {
(o) odeprint ode@dict /tcur tout def /tout tout dt 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
|
