Runge-kutta fourth order for 3 coupled second order equations.

nonlinear dynamicsnumerical methodsrunge-kutta-methodsstability-in-odes

Someone, please help me by checking whether the steps of the application of RK4 in my calculation is correct or not. I could not find any paper/books/write with similar problems or examples. Calculations for $x_1$, $\dot x_1$, $y_1$, $\dot y_1$,$\ z_1$, $\dot z_1$ is given below with initial conditions : At, $t_0 =0 sec, r_0 = 3.5 * 10^8 km, \dot r_0 =x_0=0, $$\theta_0 =(π/2)°$,$\dot \theta_0 =y_0 = 0, \phi_0=0^∘, $$\dot \phi_0 = z_0 =0 $,

if, \begin{align} \dot r= x,\, \ddot r = \dot x,\,
\dot \theta&=y,\ddot \theta =\dot y,\,
\dot \phi =z,\,
\ddot \phi=\dot z
\end{align}
then, the six first-order differential equations:
\begin{alignat}1
\dot r& = x &= f_1(t,r, \theta,\phi,x,y,z)
\\
\dot x &= r[y^2 +(z +Ω)^2(\sin\theta)^2 – \beta z\sin\theta B_\theta)] &= f_2(t,r,\theta,\phi,x,y,z)
\\
\dot \theta &= y &= f_3(t,r,\theta, \phi, x,y,z)
\\
\dot y &= \frac{1}{r}[- 2xy – r(z+ Ω)^2 \sin\theta \cos\theta + \beta r z \sin\theta B_r] &=f_4
\\
\dot \theta &=z &= f_5(t,r,\theta,\phi,x,y,z)
\\
\dot z &= \frac{1}{r \sin\theta}[-2x(z+Ω)\sin\theta – 2ry(z+ Ω)\cos\theta + \beta(xB_\theta – ryB_r)]&= f_6
\end{alignat}
The solutions are:
\begin{align}
\ x_1 &=\ x_0 + \frac{1}{6} (k_0 + 2k_1 + 2k_2+ k_3)
\\
\dot x_1 &= \dot x_0 + \frac{1}{6} (l_0 + 2l_1 + 2l_2+ l_3)
\\
y_1 &= y_0 + \frac{1}{6} (m_0 + 2m_1 + 2m_2+ m_3)
\\
\dot y_1 &= \dot y_0 + \frac{1}{6} (n_0 + 2n_1 + 2n_2+ n_3)
\\
z_1 &= z_0 + \frac{1}{6} (p_0 + 2p_1 + 2p_2+ p_3)
\\
\dot z_1 &= \dot z_0 + \frac{1}{6} (q_0 + 2q_1 + 2q_2+ q_3)
\end{align}
To find x_1, we can find $${k_0 = h* f_1(t_0,x_0,y_0,z_0,r_0,\theta_0, \phi_0)}$$ etc.
But to find x_2, what would be the values of $${r_1,\theta_1,\phi_1\quad in\quad f_1(t_1,x_1,y_1,z_!,r_1,\theta_1, \phi_1)} $$ ?

https://drive.google.com/file/d/1VVzkxPaIX7m6zQ5bP8hZVtImPZwz6XA-/view?usp=sharing

Soving 3 coupled second order equations using RK4

calculation
enter image description here

Best Answer

You are right to be suspect of this code, they made a common beginner error, and some others in the course of calculation and transcription.

  • The order of the variables is different in the argument tuples and the derivative functions. The arguments are $(t,x,y,z,r,θ,ϕ)$, the functions $(f_1,...,f_6)$ are for $(\dot r,\dot x=\ddot r,\dot θ,\dot y, \dot ϕ,\dot z)$. With $k_j=hf_1, l_j=hf_2,m_j,n_j,p_j,q_j=hf_6$ the first updated point should be $$(t_0+h/2,x_0+l_0/2,y_0+n_j/2,z_0+q_0/2,r_0+k_0/2,θ+m_0/2,ϕ_0+p_0/2).$$ Compare how the arguments that are actually used mix the two argument orders.

  • In the formulation used the first equation that is actually solved is $\dot x=f_1(t,..)=x$, giving an exponential function $x(t)=x_0e^t$ as exact solution.

  • In $l_0$ the factor $r_0$ is for some reason distributed to the terms in the sum factor, but not to the only non-zero middle term.

  • The formula for $n_0$ is sign-incompatible with the formula for $\dot y$, the value is zero anyway.

  • In the last 3 or 4 equations of each stage, the updated values are not used?

  • There may be some mismatch of degrees and radians?

  • On the first page in the last line one finds $p_2=...=0.5·0.06·10^{-12}=0.3·10^{-12}$ which is wrong by a factor of $10$.

  • The method used is not RK4, it has 5 stages that are all except the first evaluated at $t+h/2$. RK4 has 4 stages that are evaluated at $t,t+h/2, t+h/2, t+h$.

  • The collection formulas for the next point are completely corrupt.

In python one can arrange the step computation as

from math import sin, cos, pi
Ω=9.49e-7
β=3.12e-18
def acc(u,v):
    r,θ,ϕ = u
    x,y,z = v
    B_θ=-8.6e-6*cos(θ)
    B_r=25893.2e-9*sin(θ)
    ax = r*(y**2 + (z+Ω)**2 * sin(θ)**2 - β*z*sin(θ)*B_θ)
    ay = (-2*x*y-r*(z+Ω)**2*sin(θ)*cos(θ)+β*r*z*sin(θ)*B_r)/r
    az = (-2*x*(z+Ω)*sin(θ)-2*r*y*(z+Ω)*cos(θ)+β*(x*B_θ-r*y*B_r))/(r*sin(θ))
    return np.array([ax,ay,az])

h=0.5

1st stage

r0,θ0,ϕ0 = 0.7e+8, 0.5*pi, 0
u0, v0 = np.array([r0,θ0,ϕ0]), np.zeros(3)

a0 = acc(u0,v0)
print("acc0 = ",a0)
#>>> acc0 =  [ 6.30420700e-05 -5.51459066e-29 -0.00000000e+00]

2nd stage

u1 = u0+0.5*h*v0
v1 = v0+0.5*h*a0
print("[r1,θ1,ϕ1] = ",u1)
print("[x1,y1,z1] = ",v1)
#>>> [r1,θ1,ϕ1] =  [7.00000000e+07 1.57079633e+00 0.00000000e+00]
#>>> [x1,y1,z1] =  [ 1.57605175e-05 -1.37864766e-29  0.00000000e+00]

a1 = acc(u1,v1)
print("acc1 = ",a1)
#>>> acc1 =  [ 6.30420700e-05 -5.51459066e-29 -4.27335175e-19]

3rd stage

u2 = u0+0.5*h*v1
v2 = v0+0.5*h*a1
print("[r2,θ2,ϕ2] = ",u2)
print("[x2,y2,z2] = ",v2)
#>>> [r2,θ2,ϕ2] =  [7.00000000e+07 1.57079633e+00 0.00000000e+00]
#>>> [x2,y2,z2] =  [ 1.57605175e-05 -1.37864766e-29 -1.06833794e-19]

a2 = acc(u2,v2)
print("acc2 = ",a2)
#>>> acc2 =  [ 6.30420700e-05 -5.51459066e-29 -4.27335174e-19]

4th stage

u3 = u0+h*v2
v3 = v0+h*a2
print("[r3,θ3,ϕ3] = ",u3)
print("[x3,y3,z3] = ",v3)
#>>> [r3,θ3,ϕ3] =  [ 7.00000000e+07  1.57079633e+00 -5.34168968e-20]
#>>> [x3,y3,z3] =  [ 3.15210350e-05 -2.75729533e-29 -2.13667587e-19]

a3 = acc(u3,v3)
print("acc3 = ",a3)
#>>> acc3 =  [ 6.30420700e-05 -5.51459066e-29 -4.27335174e-19]

next point

unext = u0+h/6*(v0+2*v1+2*v2+v3)
vnext = v0+h/6*(a0+2*a1+2*a2+a3)

print("[rn,θn,ϕn] = ",unext)
print("[xn,yn,zn] = ",vnext)
#>>> [rn,θn,ϕn] =  [ 7.00000000e+07  1.57079633e+00 -3.56112645e-20]
#>>> [xn,yn,zn] =  [ 3.15210350e-05 -2.75729533e-29 -2.13667587e-19]
Related Question