%% %% This is file `pst-ode.tex', %% %% 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. %% %% `pst-ode' defines \pstODEsolve for integrating systems of first order %% ODEs using the Runge-Kutta-Fehlberg (RKF45) method with automatic %% step size adjustment %% \def\fileversion{0.2} \def\filedate{2012/09/14} \csname PSTODELoaded\endcsname \let\PSTODELoaded\endinput % % Requires some packages \ifx\PSTricksLoaded\endinput\else \input pstricks \fi % \message{`pst-ode' v\fileversion, \filedate\space (ag)} % \edef\PstAtCode{\the\catcode`\@} \catcode`\@=11\relax \pst@addfams{pst-ode} %% prologue for postcript \pstheader{pst-ode.pro}% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % pstODEsolve % % LaTeX command for integrating systems of first order ODEs using the Runge- % Kutta-Fehlberg method with automatic step size adjustment; % values of the integration parameter `t' as well as the solution (= state) % vectors `x(t)' at output points are stored as a long list in a Postscript % object; its content can be plotted using the listplot* functions % of pst-plot and pst-3dplot packages. % % Optionally, the result can be written to a file (ps2pdf -dNOSAFER ...) % % Usage: % % \pstODEsolve[Options] % {result}{output vector format}{ta}{tb}{number of output points} % {initial cond.}{function} % % options: % * append result appended to which must already exist (e. g. % from previous use of \pstODEsolve); usually the initial % condition vector argument is left empty in order to continue % integration from the last state % * saveData result is written to file .dat % * algebraic (system of) ODE given in infix notation % * algebraicIC initial condition vector given in infix notation % * silent suppress output of stepping information % * varsteptol relative tolerance for step size adjustment % % arguments: % % #1: Postscript identifier taking the result as a long list of state vectors % calculated at the output points % #2: output vector format, e. g. `(t) 0 1'; specifies which data to be written % to #1; (t) (parentheses required) puts value of integration parameter `t' % into the output list; 0, 1, 2, etc. specify the elements of the state % vector to be put into the output list % #3: start value of integration parameter (ta) % #4: end value of integration parameter (tb) % #5: number of output points, including ta and tb; must be >= 2 % #6: initial condition vector; if empty, continue integration from the last % state vector of the previous \pstODEsolve usage or from state vector set % by \pstODEsetRestorestate macro % #7: right hand side of ODE system; if given in Postscript notation (algebraic= % false), the function provided should pop (in reverse order!) the elements % of the current state vector from the operand stack and push the first % derivatives (right hand side of ODE system) back to it; inside the % function, the integration parameter can be accessed using `t' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \define@boolkey[psset]{pst-ode}[Pst@]{append}[true]{}% \define@boolkey[psset]{pst-ode}[Pst@]{saveData}[true]{}% \define@boolkey[psset]{pst-ode}[Pst@]{algebraicIC}[true]{}% \define@boolkey[psset]{pst-ode}[Pst@]{algebraic}[true]{}% \define@boolkey[psset]{pst-ode}[Pst@]{silent}[true]{}% \define@boolkey[psset]{pst-ode}[Pst@]{algebraicOutputFormat}[true]{}% \define@key[psset]{pst-ode}{varsteptol}{\def\ode@varsteptol{#1}}% \def\pstODEsolve{\def\pst@par{}\pst@object{pstODEsolve}}% \def\pstODEsolve@i#1#2#3#4#5#6#7{% \pst@killglue% \begingroup% \use@par% \def\filemode{w}% \ifPst@append\def\filemode{a}\fi% \edef\ode@init{#6{}}% \edef\ode@init{\expandafter\ode@zapspace\ode@init\@nil}% \edef\ode@arg{#7{}}% \ifPst@algebraicOutputFormat\edef\ode@fmt{#2{}}\fi% \pstVerb{ tx@odeDict begin \ifPst@silent userdict /odeprint systemdict /pop get put \else userdict /odeprint {print flush} put \fi /ode@tol \ode@varsteptol\space def % rel. tolerance for step size adjustment %process arguments \ifPst@saveData /statefile (#1.dat) (\filemode) file def \fi /tStart #3 def /tEnd #4 def /dt tEnd tStart sub #5\space 1 sub div def % output step size \ifx\@empty#6\@empty\else \ifPst@algebraicIC true setglobal globaldict /ode@laststate [ tx@Dict begin (\ode@init) AlgParser cvx exec end ] put false setglobal \else true setglobal globaldict /ode@laststate [tx@Dict begin 1 dict begin #6 end end] put false setglobal \fi \fi% /xlength ode@laststate length def % number of equations /xlength1 xlength 1 add def % number of equations plus 1 \ifPst@algebraic /ode@rpn (\expandafter\ode@zapspace\ode@arg\@nil) AlgParser cvx bind def /ODESET {%system of ODEs tx@Dict begin /x exch def /y x def ode@rpn end xlength array astore } bind def \else /ODESET { aload pop tx@Dict begin 0 begin #7 end end xlength array astore } bind def %ensure local scope of user defined variables in #7 /ODESET load 4 1 dict put \fi \ifPst@algebraicOutputFormat /ode@fmtrpn (\expandafter\ode@zapspace\ode@fmt\@nil) AlgParser cvx bind def /formatoutput {% tx@Dict begin ode@laststate /x exch def /y x def /t tcur def ode@fmtrpn end } bind def \else /formatoutput {[#2] assembleresult} def \fi% %perform ode integration (\string\n pstODEsolve RKF45 method; '-' failed step, '+' successful step, % 'o' output step : \string\n) odeprint /tcur tStart def % current parameter t value /tout tStart dt add def % next output t /ddt dt def % initial integration step size \ifPst@append\else [ [formatoutput] \ifPst@saveData dup writeresult \fi aload pop true setglobal ] globaldict exch /#1 exch cvx put false setglobal (o) odeprint \fi #5\space 1 sub { ode@laststate ODEINT [ exch aload pop true setglobal ] globaldict exch /ode@laststate exch put false setglobal [ #1 [formatoutput] \ifPst@saveData dup writeresult \fi aload pop true setglobal ] globaldict exch /#1 exch cvx put false setglobal } repeat \ifPst@saveData statefile closefile \fi end % tx@odeDict }% \endgroup% \ignorespaces% } \def\ode@zapspace#1#2\@nil{% helper for stripping spaces from argument \ifx#1 \relax\else#1\fi% \ifx\@empty#2\@empty\else\ode@zapspace#2\@nil\fi% } %macro for saving last state vector into PS object \def\pstODEsaveState#1{% #1: identifier \pstVerb{ true setglobal globaldict /#1 [ode@laststate cvx exec] cvx put false setglobal }% } %macro for restoring state vector from PS object or setting it from %literal vector \def\pstODEsetOrRestoreState#1{% #1: identifier or \pstVerb{ % space separated list of numbers true setglobal globaldict /ode@laststate [#1] put false setglobal }% } \psset[pst-ode]{append=false,saveData=false,algebraicIC=false, algebraic=false,silent=false,varsteptol=1e-6,algebraicOutputFormat=false, } \catcode`\@=\PstAtCode\relax %% END: pst-ode.tex \endinput