diff options
Diffstat (limited to 'Master/texmf-dist/tex/generic/pgf/math/pgfmathode.code.tex')
-rw-r--r-- | Master/texmf-dist/tex/generic/pgf/math/pgfmathode.code.tex | 178 |
1 files changed, 178 insertions, 0 deletions
diff --git a/Master/texmf-dist/tex/generic/pgf/math/pgfmathode.code.tex b/Master/texmf-dist/tex/generic/pgf/math/pgfmathode.code.tex new file mode 100644 index 00000000000..d8837897558 --- /dev/null +++ b/Master/texmf-dist/tex/generic/pgf/math/pgfmathode.code.tex @@ -0,0 +1,178 @@ +% Copyright 2011 by Christophe Jorssen +% +% This file may be distributed and/or modified +% +% 1. under the LaTeX Project Public License and/or +% 2. under the GNU Public License. +% +% See the file doc/generic/pgf/licenses/LICENSE for more details. +\def\pgfmathode@stopadd{\pgfmathode@stopadd} + +\def\pgfmathode@addindextofunc[#1]#2{% + \begingroup + \gdef\pgfmathode@func{}% + \pgfmathode@addindextofunc@i[#1]#2\pgfmathode@stopadd} + +\def\pgfmathode@addindextofunc@i[#1]#2#3{% + \ifx#1#2 + \expandafter\gdef\expandafter\pgfmathode@func\expandafter{% + \pgfmathode@func #2[\pgfmathode@timeindex]}% + \else + \expandafter\gdef\expandafter\pgfmathode@func\expandafter{% + \pgfmathode@func #2}% + \fi + \ifx#3\pgfmathode@stopadd + \expandafter\gdef\expandafter\pgfmathode@func\expandafter{% + \expandafter{\pgfmathode@func}}% + \let\next\endgroup + \else + \def\next{\pgfmathode@addindextofunc@i[#1]#3}% + \fi + \next} + +%********************************************************************** +% Know limitation: fk cannot be a time function. Should be possible to +% use \Sol[0] as time in the set of equations. +%********************************************************************** + +% Numerically solve a differential system with 1st order Runge-Kutta +% (Euler method) +% dq1/dt = f1(t,q1,q2,...,qN) +% dq2/dt = f2(t,q1,q2,...,qN) +% ... +% dqN/dt = fN(t,q1,q2,...,qN) +% +% #1: A pgf array where the solution will be stored. +% At the end: +% #1 -> {{t(0),q1(t(0)),q2(t(0)),...,qN(t(0))}, +% {t(1),q1(t(1)),q2(t(1)),...,qN(t(1))}, +% ... +% {t(Nstep),q1(t(Nstep)),q2(t(Nstep)),...,qN(t(Nstep))}} +% #2: A pgf array with f1,f2,...,fN +% {f1} +% #3: A pgf array with initial conditions. +% {t(0),q1(t(0)),q2(t(0)),...,qN(t(0))} +% #4: t(Nstep) +% #5: Nstep (number of steps) +% +% * First order example: dq/dt = -q + 5 with IC=(t=0,q(0)=0) +% \pgfmathodeRKI{\Sol}{{-\Sol[1]+5}}{{0,0}}{3.5}{100} +% * Second order example: {dq1/dt=q2, dq2/dt=-q2/2-sin(q1)} (damped +% pendulum) +% \pgfmathodeRKI{\Sol}{{\Sol[2],-\Sol[2]/2-sin(deg(\Sol[1]))}}% +% {{0,-1,5}}{15}{50} + +\def\pgfmathode@generic@init@RK{% + \begingroup + % Sol array (#1) first element is given by the set of initial + % conditions (#3). We need to enclose it in brace so that is a now + % a "matrix" array. + \gdef#1{{#3}}% + % The set of rhs {f1(t,q1,q2,...,qN),...,fN(t,q1,q2,...,qN)} is + % given in terms of #1[0],#1[1],...,#1[N]. We need to transform to + % #1[\pgfmathode@timeindex][0],...,#1[\pgfmathode@timeindex][N], that is + % add [\pgfmathode@timeindex] after every occurence of #1 in #2. + % No need to brace #2 since it is already surounded by braces. The + % result is stored in \pgfmathode@func + \pgfmathode@addindextofunc[#1]#2% + % Set time step. Actually the number of time steps is #5-1, so we + % divide by (#5-1)+1 = #5. + \pgfmathsetmacro\pgfmathode@Sol@t{#1[0][0]}% + \pgfmathsetmacro\pgfmathode@Deltat{(#4-\pgfmathode@Sol@t)/#5}% + % Compute the order of the system: it is equal to the number of IC + % minus one (time). + \pgfmathsetmacro\pgfmathode@dimfunc{dim(#3)-1}% + \pgfmathsetmacro\pgfmathode@numstep{#5-1}} + +\def\pgfmathodeRKI#1#2#3#4#5{% + \pgfmathode@generic@init@RK + % Start the time loop for Runge-Kutta 1st order + \foreach \pgfmathode@timeindex in {0,...,\pgfmathode@numstep} {% + \message{Step \pgfmathode@timeindex...}% + % Start the function loop + \foreach \pgfmathode@funcindex in {1,...,\pgfmathode@dimfunc} {% + % The RKI algorithm + % qk[t(i)] = Deltat*fk(t(i-1),q1(t(i-1)),...qN(t(i-1))) + + % qk[t(i-1)] + % Note: this part can be really slow due to array managment in + % pgfmath. + % TODO: find a better way (at least inside this loop) + \pgfmathparse{% + \pgfmathode@Deltat*(\pgfmathode@func[\pgfmathode@funcindex-1])+ + #1[\pgfmathode@timeindex][\pgfmathode@funcindex]}% + % Store the result in a helper macro + \expandafter\xdef\csname pgfmathode@Sol@\pgfmathode@funcindex + \endcsname{\pgfmathresult}% + } + % Store step i + \pgfmathparse{\pgfmathode@Sol@t+\pgfmathode@Deltat}% + \xdef\pgfmathode@Sol@t{\pgfmathresult}% + \xdef\pgfmathode@temp{}% + \foreach \pgfmathode@funcindex in {1,...,\pgfmathode@dimfunc} {% + \ifx\pgfmathode@temp\pgfutil@empty + \xdef\pgfmathode@temp{% + \csname pgfmathode@Sol@\pgfmathode@funcindex\endcsname}% + \else + \xdef\pgfmathode@temp{% + \pgfmathode@temp, + \csname pgfmathode@Sol@\pgfmathode@funcindex\endcsname}% + \fi} + \xdef#1{% + {\expandafter\pgfutil@firstofone#1,% + {\pgfmathode@Sol@t,\pgfmathode@temp}}}% + } + \endgroup} + +% RKIV: far from working yet! +\def\pgfmathodeRKIV#1#2#3#4#5{% + \pgfmathode@generic@init@RK + % Start the time loop for Runge-Kutta 4th order + \foreach \pgfmathode@timeindex in {0,...,\pgfmathode@numstep} {% + \message{Step \pgfmathode@timeindex...}% + % Start the function loop + \foreach \pgfmathode@funcindex in {1,...,\pgfmathode@dimfunc} {% + % The RKIV algorithm + % qk[t(i)] = qk[t(i-1)] + + % (1/6) * (j1k + 2j2k + 2j3k + j4k) + % where j1k = fk(t(i-1),q1(t(i-1)),...qN(t(i-1))) * Deltat + % j2k = fk(t(i-1)+.5*Deltat, + % q1(t(i-1))+.5*j11,...qN(t(i-1))+.5*j1N) * Deltat + % j3k = fk(t(i-1)+.5*Deltat, + % q1(t(i-1))+.5*j21,...qN(t(i-1))+.5*j2N) * Deltat + % j4k = fk(t(i-1)+Deltat, + % q1(t(i-1))+j31,...qN(t(i-1))+j3N) * Deltat + % + % Compute j1k + \pgfmathparse{% + \pgfmathode@Deltat*(\pgfmathode@func[\pgfmathode@funcindex-1])+ + #1[\pgfmathode@timeindex][\pgfmathode@funcindex]}% + \expandafter\xdef\csname pgfmathode@RKIV@1@\pgfmathode@funcindex + \endcsname{\pgfmathresult}% + % Compute j2k + \pgfmathparse{% + \pgfmathode@Deltat*(\pgfmathode@func[\pgfmathode@funcindex-1])+ + #1[\pgfmathode@timeindex][\pgfmathode@funcindex]}% + \expandafter\xdef\csname pgfmathode@RKIV@1@\pgfmathode@funcindex + \endcsname{\pgfmathresult}% + } + % Store step i + \pgfmathparse{\pgfmathode@Sol@t+\pgfmathode@Deltat}% + \xdef\pgfmathode@Sol@t{\pgfmathresult}% + \xdef\pgfmathode@temp{}% + \foreach \pgfmathode@funcindex in {1,...,\pgfmathode@dimfunc} {% + \ifx\pgfmathode@temp\pgfutil@empty + \xdef\pgfmathode@temp{% + \csname pgfmathode@Sol@\pgfmathode@funcindex\endcsname}% + \else + \xdef\pgfmathode@temp{% + \pgfmathode@temp, + \csname pgfmathode@Sol@\pgfmathode@funcindex\endcsname}% + \fi} + \xdef#1{% + {\expandafter\pgfutil@firstofone#1,% + {\pgfmathode@Sol@t,\pgfmathode@temp}}}% + } + \endgroup} + + +\endinput
\ No newline at end of file |