在 ODEINT 中求解微分方程

Question

我快疯了，我不知道我的代码有什么问题。我正在尝试使用 ODEINT 求解微分方程，但由于某些原因它无法正常工作，问题是，它是 5 的矩阵： equation

def model(y,t):
    
    #constant
    u = 7
    
    #l(t)
    l = 8.924 \
    - 1.584 * cos(math.radians((pi*t) / 1.51))\
    + 7.897 * sin(math.radians((pi*t) / 3.02))\
    - 10.434 * cos(math.radians((pi*t) / 4.53))\
    + 4.293 * cos(math.radians((pi*t) / 6.04))
    
    p0 = y[0]
    p1 = y[1]
    p2 = y[2]
    p3 = y[3]
    p4 = y[4]
    
    #Differential equations
    dp0dt = -l*p0 + u*p1
    dp1dt = l*p0 - (l+u)*p1 + u*p2
    dp2dt = l*p1 - (l+u)*p2 + u*p3
    dp3dt = l*p2 - (l+u)*p3 + u*p4
    dp4dt = l*p3 - u*p4
    
    return dp0dt, dp1dt, dp2dt, dp3dt, dp4dt

这是 ODEINT 和绘图代码：

#initial condition
y0 = [1,0,0,0,0]

#time
time = np.linspace(0,8)

#solve ode
y = odeint(model,y0,time)

p0 = y[:,0]
p1 = y[:,1]
p2 = y[:,2]
p3 = y[:,3]
p4 = y[:,4]

#plot
plt.plot(time,p4)
plt.xlabel('time')
plt.ylabel('p4')
plt.show()

它应该这样绘制： p4

Answer 1

主要问题是您在定义函数 l 中使用了 math.radians。在数学中，cos 等参数几乎总是以弧度表示。如果它们涉及 pi 那么它们肯定是弧度的。

所以我将那部分重写为

    #l(t)
    l = 8.924 \
    - 1.584 * math.cos(math.pi*t / 1.51)\
    + 7.897 * math.sin(math.pi*t / 3.02)\
    - 10.434* math.cos(math.pi*t / 4.53)\
    + 4.293 * math.cos(math.pi*t / 6.04)

你在 model 函数中的身份也是错误的——我冒昧地在问题中修正了它们，因为我认为这是一个复制和粘贴问题，而不是你真正的代码问题

通过这个修复，我从你的代码中得到了这张图

在我看来，这正是您想要的

在 ODEINT 中求解微分方程

Solving a Differential equation in ODEINT

python

differential-equations

odeint