Python 3.x Runge Kutta simple orbit

Question

I am in the early stages of creating a program to plot orbits using the Runge-Kutta method, and would like to plot the orbit in 2D, however, no matter what the initial conditions are, i get a straight line. I have seen a similar question but it didn't solve my problem. Why is this happening?

import numpy as np
import matplotlib.pyplot as mpl
def derX(vx):
    return vx

def derY(vy):
    return vy

def derVx(x,y):
    return -(G*M*x)/((x**2 + y**2)**(3/2))

def timestep(x,k1,k2,k3,k4):
    return x + (step/6)*(k1 + 2*k2 +2*k3 + k4)

G=6.67408E-11 #m^3/kg s^2
M=5.972E24 #kg, mass of Earth
step=100 #seconds
x=4596194 #initial conditions in m and m/s
y=4596194
vx=-6646
vy=6646
t=0
T=3600
bodyx = 444 #stationary body position metres
bodyy = 444
tarray=[]
xarray=[]
yarray=[]
vxarray=[]
vyarray=[]
while t<T:
    k1 = np.zeros(4)
    k2 = np.zeros(4)
    k3 = np.zeros(4)
    k4 = np.zeros(4)
    tarray.append(t)
    xarray.append(x)
    yarray.append(y)
    vxarray.append(vx)
    vyarray.append(vy)
    x = bodyx - x
    y = bodyy - y
    k1[0]=derX(vx)
    k1[1]=derY(vy)
    k1[2]=derVx(x,y)
    k1[3]=derVx(y,x)


    k2[0]=derX(vx+(step/2)*k1[2])
    k2[1]=derY(vy+(step/2)*k1[3])
    k2[2]=derVx(x+(step/2)*k1[0],y+(step/2)*k1[1])
    k2[3]=derVx(y+(step/2)*k1[1],x+(step/2)*k1[0])


    k3[0]=derX(vx+(step/2)*k2[2])
    k3[1]=derY(vy+(step/2)*k2[3])
    k3[2]=derVx(x+(step/2)*k2[0],y+(step/2)*k2[1])
    k3[3]=derVx(y+(step/2)*k2[1],x+(step/2)*k2[0])


    k4[0]=derX(vx+step*k3[2])
    k4[1]=derY(vy+step*k3[3])
    k4[2]=derVx(x+step*k3[0],y+step*k3[1])
    k4[3]=derVx(y+step*k3[1],vx+step*k3[0])

    t=t+step
    x=timestep(x,k1[0],k2[0],k3[0],k4[0])
    y=timestep(x,k1[1],k2[1],k3[1],k4[1])
    vx=timestep(x,k1[2],k2[2],k3[2],k4[2])
    vy=timestep(x,k1[3],k2[3],k3[3],k4[3])

mpl.plot(xarray, yarray)

Answer 1

There is a spurious v in the computation of k4[3].

The call of timestep has x as argument where it should be y, vx, vy.

And another error seems to be that in the difference computation

x = bodyx - x
y = bodyy - y

you also change the absolute position. Also the force direction becomes reversed.

Change that to something like

diffx = x - bodyx
diffy = y - bodyy

and use these relative positions in the force computation.

---

To compare, the built-in procedures produce scipy.integrate.odeint with data

G=6.67408E-11 #m^3/kg s^2
M=5.972E24 #kg, mass of Earth

bodyx = 444 #stationary body position metres
bodyy = 444

def system(u,t):
    x,y,vx,vy = u
    x -= bodyx
    y -= bodyy
    f = -(G*M)/((x**2 + y**2)**(1.5))
    return [ vx, vy, f*x, f*y ]

x0=4596194 #initial conditions in m and m/s
y0=4596194
vx0=-6646
vy0=6646
u0 = [ x0, y0, vx0, vy0 ]

T= np.linspace(0,3600,36+1)

sol = odeint(system, u0, T)

mpl.plot(sol[:,0], sol[:,1]); mpl.show()

gives a nicely curved bow, about 1/4 of a full orbit.