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

calculation enter image description here

  • 2
    $\begingroup$ This will be much easier to read when it is properly formatted. Start here to find out how: math.stackexchange.com/help/notation $\endgroup$
    – David K
    Commented Sep 16, 2022 at 20:04
  • 1
    $\begingroup$ RK = Runge-Kutta. Please say it. $\endgroup$
    – Jean Marie
    Commented Sep 16, 2022 at 20:26
  • 2
    $\begingroup$ Use $\dot x$ or $\ddot x$ for $\dot x$ or $\ddot x$ resp. $\endgroup$
    – Jean Marie
    Commented Sep 16, 2022 at 20:29
  • $\begingroup$ I just edited your last group of formulas in order you have a model. $\endgroup$
    – Jean Marie
    Commented Sep 16, 2022 at 20:41
  • 1
    $\begingroup$ 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 $\endgroup$ Commented Sep 17, 2022 at 10:31

1 Answer 1


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
def acc(u,v):
    r,θ,ϕ = u
    x,y,z = v
    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])


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]
  • $\begingroup$ 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 $\endgroup$ Commented Sep 26, 2022 at 17:56
  • $\begingroup$ 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. $\endgroup$ Commented Sep 26, 2022 at 17:58
  • $\begingroup$ Done as your suggestion: import numpy as np import matplotlib.pyplot as plt from math import sin, cos, pi. It now works fine. $\endgroup$ Commented Sep 26, 2022 at 18:02
  • $\begingroup$ 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? $\endgroup$ Commented Sep 30, 2022 at 13:22
  • $\begingroup$ 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. $\endgroup$ Commented Sep 30, 2022 at 13:25

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .