# PYTHON Runge Kutta 4 algorithm for 2nd order ODE

I have to solve the ODE system with RK4, the problem is, my Python gives small errors each step, and on long interval the result becomes incorrect.

Here is my Phyton code:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

myu=398600.4418E+9
J2=1.08262668E-3
req=6378137

t0=86400*2.3567000000000000E+04
tN= 86400*2.3567250000000000E+04

def f(t, y):
r = (y[0]**2 + y[1]**2 + y[2]**2)**0.5
r2=r**2
r3=r**3

w=1+1.5*J2*(req*req/r2)*(1-5*y[2]*y[2]/r2)

return np.array([
y[3],
y[4],
y[5],
-myu*y[0]*w/r3,
-myu*y[1]*w/r3,
-myu*y[2]*w/r3,
])

def rk4(f, h, y0, t0):
k1 = f(t0, y0)
k2 = f(t0 + h / 2, y0 + h / 2 * k1)
k3 = f(t0 + h / 2, y0 + h / 2 * k2)
k4 = f(t0 + h, y0 + h * k3)
return y0 + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6

y0 = np.array([-9.0656779992979474E+05, -4.8397431127596954E+06, -5.0408120071376814E+06, -7.6805804020022015E+02, 5.4710987127502622E+03, -5.1022193482389539E+03])

N = 1000
h = (tN - t0) / N

fig = plt.figure()

for i in range(N):
t = t0 + h
y = rk4(f, h, y0, t0)
#    ax.scatter(y[0], y[1], y[2], c='k', s=2)
t0 = t
y0 = y

#plt.show()

print(y[0])

• Please check if I rewrote your ODEs correctly. Apr 12, 2017 at 12:46
• Welcome to math.SE. Please see math.meta.stackexchange.com/questions/588/… for information on how to attract quality answers. Proper formatting is expected; for information about writing math at this site see e.g. here, here, here and here.
– mlc
Apr 12, 2017 at 12:46
• What is the question? Any numerical algorithm will produce small errors which become large on longer intervals. Numerics is never exact. Do you think there is a mistake in your code? Or are you looking for a different algorithm? These would be two very different questions. Apr 12, 2017 at 12:48
• I have to solve the ODE and get exact solution. Can you suggest any solution? Apr 12, 2017 at 12:49
• @Simon that's not entirely true. Unless the method is unconditionally unstable, you can always perform a stability analysis and find integration steps that keep approximation errors bounded. Apr 12, 2017 at 12:53

As noted in the last section of the accepted answer at https://scicomp.stackexchange.com/questions/26620/small-errors-accumulate-while-solving-ode-of-motion, the actual error is not numerical, but an omission in transcribing the ODE system. The acceleration of the $$z$$ component has a slightly different correction factor than the $$x$$ and $$y$$ components, see the (transcribed with glitches) documentation for GLONASS orbit calculations.

To demonstrate that this was really the only error, independent of the integration method, correct the system to

def f(t, y):
r = (y[0]**2 + y[1]**2 + y[2]**2)**0.5
r2=r**2
r3=r**3

w  = 1+1.5*J2*(req*req/r2)*(1-5*y[2]*y[2]/r2)
w2 = 1+1.5*J2*(req*req/r2)*(3-5*y[2]*y[2]/r2) # note the 3 here

return np.array([
y[3],
y[4],
y[5],
-myu*y[0]*w /r3,
-myu*y[1]*w /r3,
-myu*y[2]*w2/r3,
])


and leave the remaining code for the RK4 integration unchanged. This immediately gives the vector at the end of the integration interval for N=1000 steps as

['1092060.0210', '-1971932.7224', '6685288.7526', '-408.3387', '-7208.3118', '-2057.9943']


Increasing the step number by a factor of 4 to N=4000 gives the result

['1092060.1195', '-1971931.2089', '6685289.2632', '-408.3384', '-7208.3122', '-2057.9927']


and with N=10000 steps

['1092060.1198', '-1971931.2044', '6685289.2646', '-408.3384', '-7208.3122', '-2057.9927']


This compares to what you in the comments indicated the result should be:

[1092060.1198177, -1971931.204341, 6685289.264684,  -408.338412, -7208.312193, -2057.992705]