summaryrefslogtreecommitdiff
path: root/graphics/pstricks/contrib/pst-ode/pst-ode.pro
diff options
context:
space:
mode:
Diffstat (limited to 'graphics/pstricks/contrib/pst-ode/pst-ode.pro')
-rw-r--r--graphics/pstricks/contrib/pst-ode/pst-ode.pro164
1 files changed, 103 insertions, 61 deletions
diff --git a/graphics/pstricks/contrib/pst-ode/pst-ode.pro b/graphics/pstricks/contrib/pst-ode/pst-ode.pro
index 7358e8cb4f..38cf998e9e 100644
--- a/graphics/pstricks/contrib/pst-ode/pst-ode.pro
+++ b/graphics/pstricks/contrib/pst-ode/pst-ode.pro
@@ -1,7 +1,7 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PostScript prologue for pst-ode.tex.
-% Version 0.17, 2022/05/14
+% Version 0.18, 2022/11/24
%
% Alexander Grahn (C) 2012--today
%
@@ -49,16 +49,17 @@ ode@dict
/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
+ /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
@@ -69,50 +70,50 @@ end
ODESET
ode@dict
ddt mulvect /k1 exch def
- dup k1 b21 mulvect addvect
- tcur ddt 4 div add
+ dup k1 a21 mulvect addvect
+ tcur ddt c2 mul 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
+ 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 b41 mulvect addvect k2 b42 mulvect addvect k3 b43 mulvect addvect
- tcur ddt 12 mul 13 div add
+ 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 b51 mulvect addvect k2 b52 mulvect addvect k3 b53 mulvect addvect
- k4 b54 mulvect addvect
+ 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 b61 mulvect addvect k2 b62 mulvect addvect k3 b63 mulvect addvect
- k4 b64 mulvect addvect k5 b65 mulvect addvect
- tcur ddt 2 div add
+ 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 a41 mulvect k3 a43 mulvect addvect k4 a44 mulvect addvect
- k5 a45 mulvect addvect dup
+ 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 a51 mulvect k3 a53 mulvect addvect k4 a54 mulvect addvect
- k5 a55 mulvect addvect k6 a56 mulvect addvect subvect
+ 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)
@@ -122,52 +123,93 @@ end
% => [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
+%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
- tout tcur sub dup abs ddt abs lt {/ddt exch def}{pop} ifelse
+ ddt mulvect /k2 exch def
+ dup k2 0.5 mulvect addvect
end
- RKF45
- ode@tol div dup 1 gt {
- %failed step -> reduce step size
+ 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
- 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
+ /tcur tout def
+ /outStepCnt outStepCnt 1 add def
+ /tout tStart dt outStepCnt mul add def
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
+ %decrease overshooting step size
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
+ tout tcur sub dup abs ddt abs lt {/ddt exch def}{pop} ifelse
end
- 0 eq {
- (o) odeprint
+ RKF45
+ ode@tol div dup 1 gt {
+ %failed step -> reduce step size
ode@dict
- /tcur tout def
- /outStepCnt outStepCnt 1 add def
- /tout tStart dt outStepCnt mul add def
+ 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
}{
- (+) odeprint ODEINT %continue integration
+ %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
+ } ifelse
} bind def
/writeresult { %write output vector to file
/loopproc load forall