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

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)}$$ ?

Soving 3 coupled second order equations using RK4

• This will be much easier to read when it is properly formatted. Start here to find out how: math.stackexchange.com/help/notation Commented Sep 16, 2022 at 20:04
• RK = Runge-Kutta. Please say it. Commented Sep 16, 2022 at 20:26
• Use $\dot x$ or $\ddot x$ for $\dot x$ or $\ddot x$ resp. Commented Sep 16, 2022 at 20:29
• I just edited your last group of formulas in order you have a model. Commented Sep 16, 2022 at 20:41
• Almost no-one will implement RK4 in 6 dimensions component-wise. Use the system formulation, hints in the links in math.stackexchange.com/a/4177218/115115 Commented Sep 17, 2022 at 10:31

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]

• Sir, the code gives: File ~\.spyder-py3\temp.py:22 in <module> u0, v0 = np.array([r0,θ0,ϕ0]), np.zeros(3) NameError: name 'np' is not defined. Cannot import numpy in math Commented Sep 26, 2022 at 17:56
• This import is so common I leave it often out. import numpy as np and also import matplotlib.pyplot as plt. If you can not import numpy at all, this is a more serious problem. Commented Sep 26, 2022 at 17:58
• Done as your suggestion: import numpy as np import matplotlib.pyplot as plt from math import sin, cos, pi. It now works fine. Commented Sep 26, 2022 at 18:02
• sir, Since none of the terms in the equations contains ''t', it appears that the value of say y1 evaluated at t = 0 seconds or y1 at t= 60 seconds, with RK4. For example, t0+h/2 or 60 +h/2 does not seem to affect the values of r, x1, y1, z3 etc. If I want to see the values of r,x, theta, y etc at each time ( 1 sec, 2s, 3s, .... 60 s) how would I approach it, or what changes should be made in the codes you have provided? Commented Sep 30, 2022 at 13:22
• You build a function table $(t,y)$. That the equations are autonomous just says that the absolute time or the start time has no importance, but time differences or offsets still retain their importance. Commented Sep 30, 2022 at 13:25