summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/dvips/pst-ode
diff options
context:
space:
mode:
authorKarl Berry <karl@freefriends.org>2012-08-16 22:48:03 +0000
committerKarl Berry <karl@freefriends.org>2012-08-16 22:48:03 +0000
commit9eafcfc80979b4c1fed5bdcf48308fccdbf8b56e (patch)
treee6a7979c4f1fe7af3b2c1596b9ca63a20d095efe /Master/texmf-dist/dvips/pst-ode
parenteb49231093d1336c3166b1b979ad69342430b282 (diff)
new pstricks package pst-ode (16aug12)
git-svn-id: svn://tug.org/texlive/trunk@27431 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf-dist/dvips/pst-ode')
-rw-r--r--Master/texmf-dist/dvips/pst-ode/pst-ode.pro136
1 files changed, 136 insertions, 0 deletions
diff --git a/Master/texmf-dist/dvips/pst-ode/pst-ode.pro b/Master/texmf-dist/dvips/pst-ode/pst-ode.pro
new file mode 100644
index 00000000000..22b347cd546
--- /dev/null
+++ b/Master/texmf-dist/dvips/pst-ode/pst-ode.pro
@@ -0,0 +1,136 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% PostScript prologue for pst-ode.tex.
+% Version 0.1, 2012/07/20
+%
+% 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 { %assembles state vector for output into file + writes it
+ outvect {
+ dup (t) eq {
+ pop tcur 20 string cvs statefile exch writestring
+ }{
+ ode@laststate exch get 20 string cvs statefile exch writestring
+ } ifelse
+ statefile ( ) writestring
+ } forall
+ statefile (\n) writestring
+} bind def
+/assembleresult { %assembles state vector for building table of results
+ outvect {
+ dup (t) eq {
+ pop tcur
+ }{
+ ode@laststate exch get
+ } ifelse
+ } forall
+} bind def
+end