对runge-kutta四阶步骤进行编码的正确方法

Question

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

Answer 1

因为它是一个列表python不会对列表执行算术向量运算，因此使用numpy来获取用作向量的数组

有关不带numpy的RK4的实现，请参见Solving the Lorentz model using Runge Kutta 4th Order in Python without a package。它与使用标量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] ]))