带时滞的SIR模型ddeint求解绘图异常问题排查

核心错误点 #

• 初始条件定义错误：ddeint求解延迟微分方程（DDE）需要传入历史函数（定义t≤0时系统的状态），而非ODE那样的单个初始值。你写的g = lambda t : array([1000,0.3,0.1])完全不符合要求，应该返回初始时刻的S、I、R值，且求解时必须传入这个历史函数，而非(S0,I0,R0)。
• 时滞项提取错误：Y(t-τ)返回的是完整的状态数组[S(t-τ), I(t-τ), R(t-τ)]，你直接把Is = Y(t-10)是把整个数组当成了感染人数，应该提取第二个元素（索引1），即Is = Y(t-10)[1]，同理Im = Y(t-2)[1]。
• 代码语法错误：pip install ddeint不能写在import语句中间，这是终端命令，需单独执行（Jupyter环境可用%pip install ddeint）。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from ddeint import ddeint

# 模型参数
N = 1000
beta = 0.3
gamma = 1/10
I0 = 1
R0 = 0
S0 = N - I0 - R0
t = np.linspace(0, 160, 160)

# 定义历史函数：t≤0时，系统处于初始状态
def history(t):
    return np.array([S0, I0, R0])

# 带时滞的SIR模型
def model(Y, t):
    S, I, R = Y(t)  # 当前时刻的状态
    # 提取时滞时刻的感染人数I(t-10)和I(t-2)
    I_t10 = Y(t - 10)[1]
    I_t2 = Y(t - 2)[1]
    
    # 若保留你原代码的自定义公式，用下面这行dSdt
    dSdt = -beta * S * I_t10 / N + gamma * I_t2
    # 若为标准时滞感染SIR模型，用这行dSdt：dSdt = -beta * S * I_t10 / N
    dIdt = beta * S * I_t10 / N - gamma * I
    dRdt = gamma * I
    
    return np.array([dSdt, dIdt, dRdt])

# 求解DDE
ret = ddeint(model, history, t)
S, I, R = ret.T

# 绘图（补充你确认可行的绘图逻辑）
plt.figure(figsize=(10,6))
plt.plot(t, S, label='易感者S')
plt.plot(t, I, label='感染者I')
plt.plot(t, R, label='康复者R')
plt.xlabel('时间')
plt.ylabel('人数')
plt.legend()
plt.title('带时滞的SIR模型仿真')
plt.show()

额外说明 #

你原代码中dSdt的+ gamma*Im项不属于标准SIR模型，如果你是要自定义“康复者回归易感者”的时滞机制，修正时滞项提取后即可正常运行；如果是误加的，换成注释里的标准公式即可。

内容的提问来源于stack exchange，提问作者Galbotrix