python - 编写runge kutta 4阶步骤的正确方法

问题描述

xdot返回[dp_dt, dphi_dt]但不允许k1调用 RK4 方法的第二步，因为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

标签： pythonrunge-kutta

因为它是一个列表。python不对列表进行算术向量运算，因为它用于numpy获取充当向量的数组。

对于 RK4 的实现，numpy请参阅在没有包的情况下在 Python 中使用 Runge Kutta 4th Order 求解 Lorentz 模型。它与使用标量 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] ]))

python - 编写runge kutta 4阶步骤的正确方法

问题描述

解决方案

推荐阅读