blob: f9d8fdcf50baa88a907d112c77e57d07f6fd24f7 (
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
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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
|
