[pstricks] Searching for a numeric solution for coupled differential equation systems of order 2

Alexander Grahn A.Grahn at hzdr.de
Fri May 25 14:37:37 CEST 2012

On Fri, May 25, 2012 at 01:02:23PM +0200, Alexander Grahn wrote:
>On Fri, May 25, 2012 at 01:52:18PM +0300, Manfred Braun wrote:
>>On 25.05.2012, at 12:15, Alexander Grahn wrote:
>>
>>> On Thu, May 24, 2012 at 10:21:06PM +0200, Juergen Gilg wrote:
>>>> Dear PSTricks list,
>>>>
>>>> is there anybody out there, knowing about how to generate a PS code for
>>>> "Runge-Kutta 4" to solve a "coupled differential equation system of
>>>> order 2"?
>>>>
>>>> There is the need for such a code to animate a double-pendulum in real
>>>> time, following the coupled differential equations with the variables
>>>> (\varphi_1, \varphi_2) and their derivatives like:
>>>>
>>>> (1)
>>>> l_1\ddot{\varphi}_1+\frac{m_2}{m_1+m_2}l_2\ddot{\varphi}_2\cos(\varphi_1-\varphi_2)-\frac{m_2}{m_1+m_2}l_2\dot{\varphi}_2^2\sin(\varphi_1-\varphi_2)+g\sin\varphi_1=0
>>>>
>>>> (2)
>>>> l_2\ddot{\varphi}_2+l_1\ddot{\varphi}_1\cos(\varphi_1-\varphi_2)-l_1\dot{\varphi}_1^2\sin(\varphi_1-\varphi_2)+g\sin\varphi_2=0
>>>
>>> I'd first transform this set of DEs of second order into a set of four
>>> DEs of first order,
>>
>>Yes, this is the right way.
>>>
>>> However, it may still not be amenable to Runge-Kutta, because the right
>>> hand sides of the resulting four equation set of DEs, put into normal
>>> form, still contain derivatives:
>>
>>No! Take the following steps:
>>
>>1. Solve the two equations for \ddot{\varphi_1} and \ddot{\varphi_2}:
>>    \ddot{\varphi}_1 = f_1(\varphi_1, \dot{\varphi}_1, \varphi_2, \dot{\varphi}_2)
>>    \ddot{\varphi}_2 = f_2(\varphi_1, \dot{\varphi}_1, \varphi_2, \dot{\varphi}_2)
>>
>>2. Define  x_1 = \varphi_1,  x_2 = \dot{\varphi}_1,  x_3 = \varphi_2,  x_4 = \dot{\varphi}_2
>>
>>3. The new variables satisfy the system of four 1st-order equations
>>    \dot{x}_1 = x_2
>>    \dot{x}_2 = f_1(x_1, x_2, x_3, x_4)
>>    \dot{x}_3 = x_4
>>    \dot{x}_4 = f_2(x_1, x_2, x_3, x_4)
>>
>>This can be solved by Runge-Kutta or any other procedure.
>>
>>Manfred
>
>Thanks, Manfred! In the meantime I tried to write the set of equations
>down. Here is the result, without errors I hope:
>
>The equations can be solved further to eliminate derivatives on the right
>hand side. The resulting set of first order ODEs can be integrated by
>\odesolve from the animate package documentation, provided the ODEs are
>not stiff. Initial conditions for $\varphi_1$, $\dot\varphi_1$,
>$\varphi_2$, $\dot\varphi_2$ are required:
>
>\documentclass{article}
>\usepackage{amsmath}
>\parindent0pt
>\begin{document}
>
>\section{Original set of DEs of second order}
>
>\begin{gather*}
>l_1\ddot{\varphi}_1+\frac{m_2}{m_1+m_2}l_2\ddot{\varphi}_2\cos(\varphi_1-\varphi_2)-\frac{m_2}{m_1+m_2}l_2\dot{\varphi}_2^2\sin(\varphi_1-\varphi_2)+g\sin\varphi_1=0\\
>l_2\ddot{\varphi}_2+l_1\ddot{\varphi}_1\cos(\varphi_1-\varphi_2)-l_1\dot{\varphi}_1^2\sin(\varphi_1-\varphi_2)+g\sin\varphi_2=0
>\end{gather*}
>
>\section{Set of first order ODEs in normal form, amenable to Runge-Kutta}
>
>\begin{align*}
>\dot\varphi_{11}& = \varphi_{12}\\
>\dot\varphi_{12}& = \frac{\frac{a_1-b_1}{l_1}-\frac{c_1}{l_1}\frac{a_2-b_2}{l_2}}{1-\frac{c_1 c_2}{l_1 l_2}}\\
>\dot\varphi_{21}& = \varphi_{22}\\
>\dot\varphi_{22}& = \frac{\frac{a_2-b_2}{l_2}-\frac{c_2}{l_2}\frac{a_1-b_1}{l_1}}{1-\frac{c_1 c_2}{l_1 l_2}}\\
>\end{align*}
>with
>\begin{align*}
>a_1& = \frac{m_2}{m_1+m_2} l_2 \varphi_{22}^2 \sin(\varphi_{11}-\varphi_{21})\\
>b_1& = g \sin\varphi_{11}\\
>c_1& = \frac{m_2}{m_1+m_2} l_2 \cos(\varphi_{11}-\varphi_{21})\\
>a_2& = l_1 \varphi_{12} \sin(\varphi_{11}-\varphi_{21})\\
>b_2& = g \sin\varphi_{21}\\
>c_2& = l_1 \cos(\varphi_{11}-\varphi_{21})
>\end{align*}
>\end{document}

What follows is the (untested) PS code for the set of ODEs to be
integrated with \odesolve from the animate package documentation. Anything
marked as ?.*?  must be replaced by meaningful values. Note that angles
are measured in degrees. ps2pdf must be invoked as ps2pdf -dNOSAFER,
because the result is written into a file.

%define m_1, m_2, l_1, l_2, g
\pstVerb{
/m1 ?? def
/m2 ?? def
/l1 ?? def
/l2 ?? def
/g  ?? def
}

%implement ODE set in RP notation as a PS function
\pstVerb{
/setofodes {
%pop elements of current solution vector from operand stack, note the
%reverse order of popping
/phi22 exch def  %\dot\varphi_2
/phi21 exch def  %\varphi2
/phi12 exch def  %\dot\varphi_1
/phi11 exch def  %\varphi1

%auxiliary variables a_1, b_1, c_1, a_2, b_2, c_2
/a1 m2 m1 m2 add div l2 mul phi22 phi22 mul mul phi11 phi21 sub sin mul def
/b1 g phi11 sin mul def
/c1 m2 m1 m2 add div l2 mul phi11 phi21 sub cos mul def

/a2 l1 phi12 mul phi11 phi21 sub sin mul def
/b2 g phi21 sin mul def
/c2 l1 phi11 phi21 sub cos mul def

%compute right hand side of odes (4 values), leave them on the
%operand stack
phi12 %\dot\varphi_1

%\ddot\varphi_1
a1 b1 sub l1 div c1 l1 div a2 b2 sub l2 div mul sub 1 c1 c2 mul l1 l2 mul
div sub

phi22 %\dot\varphi_2

%\ddot\varphi_2
a2 b2 sub l2 div c2 l2 div a1 b1 sub l1 div mul sub 1 c1 c2 mul l1 l2 mul
div sub
} def
}

\odesolve{result.dat}{% output file name
%output format: integr. param. t' & complete solution vector (4 elements)
%on one table row
(t) 0 1 2 3
}{?ta?}{?tb?}{% start and end of integration parameter t'
?number of output points?
}{
% initial values \varphi1, \dot\varphi_1, \varphi2, \dot\varphi_2 at t=0
?? ?? ?? ??
}{setofodes} % name of the PS function for right hand side of ODE set

Alexander