python - Is there a better way to tidy this chuck of code? It is the Runge-Kutta 4th order with 4 ODEs

Question

def xpos(x,y,t):    #dx/dt = v_x 
    return vx 

def ypos(x,y,t):  #dy/dt = v_y 
    return vy 

def xvel(x,y,t):   #dv_x/dt = -GMx/r^3 
    return -G*M*x/((x)**2 + (y)**2)**1.5

def yvel(x,y,t):    # dv_y/dt = -GMy/r^3
    return -G*M*y/((x)**2 + (y)**2)**1.5 

xposk1 = h*xpos(x,y,t)
yposk1 = h*ypos(x,y,t)
xvelk1 = h*xvel(x,y,t)
yvelk1 = h*yvel(x,y,t)

xposk2 = h*xpos(x+(0.5*xposk1),y+(0.5*yposk1),t+(0.5*h))
yposk2 = h*ypos(x+(0.5*xposk1),y+(0.5*yposk1),t+(0.5*h))
xvelk2 = h*xvel(x+(0.5*xvelk1),y+(0.5*yvelk1),t+(0.5*h))
yvelk2 = h*xvel(x+(0.5*xvelk1),y+(0.5*yvelk1),t+(0.5*h))

xposk3 = h*xpos(x+(0.5*xposk2),y+(0.5*yposk2),t+(0.5*h))
yposk3 = h*ypos(x+(0.5*xposk2),y+(0.5*yposk2),t+(0.5*h))
xvelk3 = h*xvel(x+(0.5*xvelk2),y+(0.5*yvelk2),t+(0.5*h))
yvelk3 = h*yvel(x+(0.5*xvelk2),y+(0.5*yvelk2),t+(0.5*h))

xposk4 = h*xpos(x+xposk3,y+yposk3,t+h)
yposk4 = h*ypos(x+xposk3,y+yposk3,t+h)
xvelk4 = h*xvel(x+xvelk3,y+yvelk3,t+h)
yvelk4 = h*yvel(x+xvelk3,y+yvelk3,t+h)

Here I have 4 ODEs (Ordinary Differential Equations) that I want to solve using the Runge-Kutta 4th order method. The code I've included above should be able to do it but I was curious if there is a much simpler and shorter way of doing it. I have only included the relevant (RK4) part of the code.

Answer 1

As there are not enough planetary orbit simulations using RK4 (if you stay on this topic, quickly switch to an adaptive time stepping method that follows the speed variations along an elliptical orbit).

The acceleration vector of the position vector x can be computed compactly as

def acc(x):
   r3 = sum(x**2)**1.5;
   return -G*M*x/r3

where it is supposed that x is always a numpy array.

Following the reduction of RK4 in the second order ODE case as computed here (and in many other places) the RK4 step can be implemented as

def RK4step(x,v,h):
    dv1 = h*acc(x)
    dv2 = h*acc(x+0.5*h*v)
    dv3 = h*acc(x+0.5*h*(v+0.5*dv1))
    dv4 = h*acc(x+h*(v+0.5*dv2))
    return x+h*(v+(dv1+dv2+dv3)/6), v+(dv1+2*(dv2+dv3)+dv4)/6