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

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 ОДУ (обыкновенные дифференциальные уравнения), которые я хочу решить, используя метод 4-го порядка Рунге-Кутты. Код, который я включил выше, должен уметь это делать, но мне было любопытно, есть ли более простой и короткий способ сделать это. Я включил только соответствующую (RK4) часть кода.


python runge-kutta
person ElonMusk_codes    schedule 26.02.2020    source источник
comment
Лучше чем каким образом?   -  person martineau    schedule 26.02.2020
comment
Этот код не будет этого делать, поскольку производная положения - это скорость, которую вы не передаете в качестве аргумента. Вы действительно должны получить ошибку во второй строке приведенного выше кода, поскольку vx не определен.   -  person Lutz Lehmann    schedule 26.02.2020
comment
Вас может вдохновить код Рунге-Кутты, не сходящийся со встроенным методом, или если вы хотите сократить решение второго Закажите ODE другим способом, вы можете исключить переменные xposk и yposk в этой ситуации.   -  person Lutz Lehmann    schedule 26.02.2020

Ответы (1)


arrow_upward
0
arrow_downward

Поскольку недостаточно моделирования планетарной орбиты с использованием 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
person Lutz Lehmann    schedule 26.02.2020