编写runge kutta 4阶步骤的正确方法答案

【问题标题】：proper way to code runge kutta 4th order steps编写runge kutta 4阶步骤的正确方法
【发布时间】：2020-04-19 03:43:07
【问题描述】：

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

【问题讨论】：

标签： python runge-kutta

【解决方案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] ]))

【讨论】：

因此，如果我使用 np.array 将 xdot 更改为数组，我将如何将数组的每个项调用到 rk4 中。我尝试在将 k1[0] 更改为数组后调用它，但它说“列表”对象不可调用
你没有。另外，再次检查什么参数是标量时间，什么是状态向量，EOM 只有状态向量作为参数。所以你会调用k2 = EOM(xold+0.5*dt*k1)，前提是所有元组都转换为np.array。