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

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

0 голосов
/ 26 февраля 2020 
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) часть кода.

1 Ответ

0 голосов
/ 26 февраля 2020

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

вектор ускорения вектора положения 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
...