PYTHON Runge Kutta 4 algorithm for 2nd order ODE

Question

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()
ax = fig.add_subplot(111, projection='3d')

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])

Answer 1

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]