правильный способ кодирования ступени Кутта 4-й шаг - PullRequest
Новая
       65

правильный способ кодирования ступени Кутта 4-й шаг

0 голосов
/ 19 апреля 2020

xdot возвращает [dp_dt, dphi_dt], но второй шаг метода RK4 не позволяет вызывать k1, поскольку k1 является списком ... 

def EOM(s)

    "Equations of Motion"
    p = s[0]
    phi = s[1]
    p_hat = p*b/(2*V)
    C = -0.06*phi+0.033*p_hat+0.073*p_hat**3-0.36*p_hat*phi**2+1.47*p_hat**2*phi
    dp_dt = 1/(2*Ixx)*rho*V**2*S*b*C
    dphi_dt = p 
    sdot = []
    sdot.append(dp_dt)
    sdot.append(dphi_dt)
    return sdot



def rk4(xold):


    xnew = []
    xdot = EOM(xold)
    for i, state in enumerate(xold):
        k1=xdot[i]
        #this is a list of [p,phi] ... how 
        k2=EOM(xold(1)+dt/2,xold(2)+dt/2*k1)
        k3=EOM(xold(1)+dt/2,xold(2)+dt/2*k2)
        k4=EOM(xold(1)+dt,xold(2)+dt*k3)
        xnew.append(state+dt/6*(k1+2*k2+2*k3+k4))
    return xnew

1 Ответ

0 голосов
/ 19 апреля 2020

Потому что это список. python не выполняет арифметические c векторные операции над списками, для этого используйте numpy для получения массивов, которые действуют как векторы.

Для реализации RK4 без numpy см. Решение Модель Лоренца с использованием Runge Kutta 4-го порядка в Python без пакета . Это не совсем то же самое, что использование списка скалярных функций ODE, но принципы должны быть видны.

Ваш код может быть исправлен для работы без numpy, но это решение только для 2 или 3 переменные, с большим он становится подверженным ошибкам 

def rk4(xold):
    k1 = EOM(xold)
    #this is a list of [p_dot,phi_dot]
    k2=EOM(xold[0]+dt/2*k1[0],xold[1]+dt/2*k1[1],) # trailing comma generates list
    k3=EOM([xold[0]+dt/2*k2[0],xold[1]+dt/2*k2[1]]) # or do it directly as list
    k4=EOM([xold[i]+dt*k3[i] for i in [0,1]]) # or use list generation
    return [ xold[i]+dt/6*(k1[i]+2*k2[i]+2*k3[i]+k4[i] for i in [0,1] ]))
...