summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/tex/generic/pgf/math/pgfmathode.code.tex
diff options
context:
space:
mode:
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.tex178
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