From 357637f662b6d1e080a92685dc01102ba935157b Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Thu, 24 Nov 2022 20:39:29 +0000 Subject: pst-ode (24nov22) git-svn-id: svn://tug.org/texlive/trunk@65096 c570f23f-e606-0410-a88d-b1316a301751 --- Master/texmf-dist/dvips/pst-ode/pst-ode.pro | 164 +++++++++++++++++----------- 1 file changed, 103 insertions(+), 61 deletions(-) (limited to 'Master/texmf-dist/dvips') diff --git a/Master/texmf-dist/dvips/pst-ode/pst-ode.pro b/Master/texmf-dist/dvips/pst-ode/pst-ode.pro index 7358e8cb4fb..38cf998e9ed 100644 --- a/Master/texmf-dist/dvips/pst-ode/pst-ode.pro +++ b/Master/texmf-dist/dvips/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 -- cgit v1.2.3