Есть ли способ лучше привести в порядок этот кусок кода? Это Рунге-Кутта 4-го порядка с 4 ODE.

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)

Здесь у меня есть 4 ODE (обыкновенные дифференциальные уравнения), которые я хочу решить, используя метод 4-го порядка Рунге-Кутты. Код, который я включил выше, должен уметь это делать, но мне было любопытно, есть ли более простой и короткий способ сделать это. Я включил только соответствующую (RK4) часть кода.

Answer 1

Поскольку моделирования планетарной орбиты с использованием RK4 недостаточно (если вы останетесь в этой теме, быстро переключитесь на метод адаптивного шага по времени, который следует за изменениями скорости вдоль эллиптической орбиты).

Вектор ускорения вектора положения x можно компактно вычислить как

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

где предполагается, что x всегда является массивом numpy.

После уменьшения RK4 в случае ODE второго порядка, вычисленного здесь (и во многих других местах), шаг RK4 может быть реализован как

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