Python使用RK4计算二阶ODE

Question

下面粘贴的是我的python代码。它是一个四阶龙格库塔，用于计算二阶常微分方程：y'+4y'+2y=0，初始条件y（0）=1，y'（0）=3

我需要帮助修复它。当我运行我的代码时，我的解析解与我的数值解不匹配，我的教授说它们应该是相同的。我试着编辑了一大堆，但似乎无法找出问题所在。谢谢大家!

    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()

Answer 1

我假设您运行的代码中的缩进是正确的，因为数值解在其他地方看起来是完全错误的。所以我猜你的精确解是错误的，或者精确解是一个与数值代码不同的方程

对于特征多项式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时间循环的缩进，并添加精确解，即可绘制出图