-1
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.

martineau
  • 119,623
  • 25
  • 170
  • 301
ElonMusk_codes
  • 1
  • Better than what way? – martineau Feb 26 '20 at 01:14
  • This code will not do it, since the position derivative is the velocity, which you do not pass as argument. You should indeed get an error in the second line of above code as `vx` is not defined. – Lutz Lehmann Feb 26 '20 at 09:13
  • You could be inspired by [Runge-Kutta code not converging with builtin method](https://stackoverflow.com/q/35071393/3088138), or if you want to shorten the solution of the second order ODE another way, you can eliminate the xposk and yposk variables in this situation. – Lutz Lehmann Feb 26 '20 at 09:26

1 Answers1

0

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
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51