Python使用RK4评估二阶ODE

如何解决Python使用RK4评估二阶ODE

下面粘贴的是我的python代码。它是四阶Runge-Kutta，它以初始条件y（0）= 1，y'（0）= 3评估二阶ODE：y''+ 4y'+ 2y = 0。

我需要修复它的帮助。我的教授说，当我运行代码时，我的解析解与数值解不符，我的教授说它们应该相同。我试过编辑很多，似乎无法找出问题所在。谢谢！

    import numpy as np
    import matplotlib.pyplot as plt

    def ode(y):
        return np.array([y[1],(-2*y[0]-4*y[1])])

    tStart=0

    tEnd=5

    h=.1

    t=np.arange(tStart,tEnd+h,h)

    y=np.zeros((len(t),2))

    y[0,:]=[1,3]

    for i in range(1,len(t)):
        k1=ode(y[i-1,:])
        k2=ode(y[i-1,:]+h*k1/2)
        k3=ode(y[i-1,:]+h*k2/2)
        k4=ode(y[i-1,:]+h*k3)
    
    y[i,:]=y[i-1,:]+(h/6)*(k1+2*k3+2*k3+k4)

    plt.plot(t,y[:,0])
    plt.plot(t,1-t)
    plt.grid()
    plt.gca().legend(('y(t)',"y'(t)"))
    plt.show()

解决方法

我认为缩进在您运行的代码中是正确的，因为数值解在其他情况下看起来完全错误。因此，我的猜测是您的精确解是错误的，或者确切的解是与数字代码不同的方程式。

使用特征多项式q^2+4q+2=(q+2)^2-2=0的精确解是

y(t)=exp(-2t)*(A*cosh(sqrt(2)*t)+B*sinh(sqrt(2)*t))

根据初始条件A=1和-2*A+sqrt(2)*B=3

更正RK4时间循环的缩进并添加精确的解决方案即可得出图