RK4 for comet trajectory in Python

Question

I'm trying to write a python program which simulates the trajectory of a comet using the Runge-Kutta 4th degree method. I'm given the acceleration formula

$$\ddot x = -\frac{GM}{r^3} x$$

I believe I've created the correct system of ODEs and implemented hem correctly; but when I run the program the numbers just explode dramatically.

from math import hypot, sqrt
from pylab import plot, show

#define constants
G = 6.674 * (10**-17)   #the gravitational 
constant N*m2/kg2
M = 1.989 * (10**30)    #mass of the Sun kg

c = G*M                 #useful constant

def a(x,r):
    return -(c*x)/(r**3)

h = 3600*24

rx = [4000000000]
vx = [0]
ry = [0]
vy = [.5]

for i in range(100000):
    r = sqrt(rx[-1]**2+ry[-1]**2)
    # print(r)
    #input()


    k1vx = a(rx[-1],r)
    k1rx = vx[-1]

    k2vx = a(rx[-1]+k1rx*h*.5,r)
    k2rx = vx[-1]*k1vx*h*.5

    k3vx = a(rx[-1]+k2rx*h*.5,r)
    k3rx = vx[-1]*k2vx*h*.5

    k4vx = a(rx[-1]+k3rx*h,r)
    k4rx = vx[-1]*k3vx*h

    vx.append(vx[-1]+(h/6)*(k1vx+2*k2vx+2*k3vx+k4vx))
    rx.append(rx[-1]+(h/6)*(k1rx+2*k2rx+2*k3rx+k4rx))

    k1vy = a(ry[-1],r)
    k1ry = vy[-1]

    k2vy = a(ry[-1]+k1ry*h*.5,r)
    k2ry = vy[-1]*k1vy*h*.5

    k3vy = a(ry[-1]+k2ry*h*.5,r)
    k3ry = vy[-1]*k2vy*h*.5

    k4vy = a(ry[-1]+k3ry*h,r)
    k4ry = vy[-1]*k3vy*h

    vy.append(vy[-1]+(h/6)*(k1vy+2*k2vy+2*k3vy+k4vy))
    ry.append(ry[-1]+(h/6)*(k1ry+2*k2ry+2*k3ry+k4ry))

    #print(rx[-1])
    #print(ry[-1])
    #input()

print(ry)
plot(rx,ry)
show()

Answer 1

You need to solve a coupled system as a coupled system, you can not separate the x updates from the y updates.

Change the acceleration computation to compute the radius directly from the actual position

def a(x): r = sum(x**2)**0.5; return -c*x/r**3

and correct the k2xr etc. computation to have an additive update

    k2rx = vx[-1]+k1vx*h*.5

where rx and vx` are now lists of vectors, numpy arrays.

See other examples here and on stackoverflow for solar system simulations.

---

Fixed-point RK4 is not a good method for planet simulations and even less for comet simulations, as the speed and curvature of the orbit changes drastically during each cycle. Close to the sun one would need extremely small step sizes, far away one can use longer time steps.