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

[Alexander Grahn]

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