无法复现教材季节性驱动SEIR模型图（疑似数值不稳定）

问题概述 #

• 目标：复现Keeling 2008《Modeling Infectious Diseases in Humans and Animals》中图5.6，验证季节性驱动SEIR模型实现正确性
• 异常现象：仅第0行第0列、第1行第0列子图生成正常，其余子图因感染比例数值解变为0，计算ln(0)产生负无穷而无输出
• 触发警告：
  • ODEintWarning: Excess work done on this call
  • RuntimeWarning: divide by zero encountered in log
  • RuntimeWarning: invalid value encountered in log
• 已知背景：教材官方代码提及大Beta1值可能引发数值错误，但图中所用Beta1值理论上不应触发该问题；短期可修改季节性参数生成有效数据，但需严格复现原图

错误原因分析 #

1. ODE求解器数值精度不足：当感染比例趋近于极小时，求解器可能因步长设置或稳定性限制，直接将其截断为0，而非保留极小非零值
2. 对数计算的数学奇点：一旦感染比例变为0，ln(0)直接产生无效值，导致后续绘图逻辑失效
3. 模型实现细节偏差：季节性驱动项的相位、幅度计算错误，或SEIR仓室转换速率与教材公式不符，导致部分参数组合下感染无法维持传播

可行解决方案 #

1. 优化ODE求解器参数 #

• 提高odeint的数值精度阈值，例如设置：

  from scipy.integrate import odeint
  odeint(model, y0, t, args=(params,), rtol=1e-10, atol=1e-12)
  

• 尝试替换求解器为scipy.integrate.solve_ivp，该求解器在处理低幅值、稀疏解时稳定性更优

2. 规避对数计算的奇点 #

• 在计算自然对数前添加极小偏移量，避免直接计算ln(0)：

  import numpy as np
  log_I = np.log(np.maximum(I, 1e-15))  # 1e-15为可调整的极小值
  

3. 核对模型实现细节 #

• 确认季节性Beta的表达式与教材一致，典型形式为：

  def beta(t, beta0, beta1, phi):
      return beta0 * (1 + beta1 * np.cos(2 * np.pi * t + phi))
  

  检查相位phi、周期设置是否与图5.6参数匹配
• 逐一验证SEIR各仓室的微分方程，确保易感者、暴露者、感染者、康复者的转换速率完全符合教材公式

4. 微调初始条件 #

• 针对出错的子图参数组合，将初始感染比例从0调整为极小值（如1e-6），确保感染有足够的初始基数维持传播

验证建议 #

• 先复现无季节性的基础SEIR模型，确认数值求解无错误后再引入季节性驱动
• 对出错的参数组合单独模拟，输出各仓室的时间序列，观察感染比例是真实衰减至0，还是求解器数值截断导致的假零值

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