用Python模拟疫情传播手把手教你用微分方程实现SIS模型附完整代码在数据科学和流行病学交叉领域数学建模正成为理解复杂传播现象的核心工具。SISSusceptible-Infectious-Susceptible模型作为经典传染病框架特别适合模拟新冠等治愈后仍可能重复感染的疾病。本文将抛开繁琐的理论推导带您用Python从零构建完整的疫情传播模拟器——通过不到100行代码您将掌握微分方程求解、参数调优和动态可视化的全流程实战技能。1. 环境配置与模型原理速览1.1 快速搭建科学计算环境推荐使用Anaconda创建专属建模环境conda create -n epidemic_model python3.8 conda activate epidemic_model conda install numpy scipy matplotlib jupyter1.2 SIS模型核心机制模型将人群划分为两类动态变化的群体易感者(S)可能被感染的健康人群感染者(I)具有传播能力的患病群体关键参数对传播趋势的影响参数物理意义典型取值范围对传播的影响λ日感染率0.1-0.5值越大传播越快μ日治愈率0.05-0.2值越大疫情消退越快提示实际建模时需通过历史数据拟合确定λ和μ的精确值2. 微分方程的数值解法实现2.1 从理论公式到Python函数SIS模型的微分方程表达式def sis_model(y, t, lambda_, mu): i y[0] didt lambda_ * i * (1 - i) - mu * i return [didt]2.2 使用SciPy求解微分方程from scipy.integrate import odeint import numpy as np # 参数设置 lambda_ 0.3 # 感染率 mu 0.1 # 治愈率 i0 0.01 # 初始感染比例 t np.linspace(0, 100, 1000) # 100天内的1000个时间点 # 求解方程 solution odeint(sis_model, [i0], t, args(lambda_, mu)) infected_curve solution[:, 0]3. 动态可视化与参数分析3.1 基础疫情曲线绘制import matplotlib.pyplot as plt plt.figure(figsize(10,6)) plt.plot(t, infected_curve, r, label感染比例) plt.xlabel(时间(天)) plt.ylabel(人群比例) plt.title(SIS模型传播趋势 (λ0.3, μ0.1)) plt.legend() plt.grid() plt.show()3.2 多参数对比实验通过修改λ和μ值观察不同传播模式params_combinations [ {lambda_: 0.2, mu: 0.1, color: b, label: 温和传播}, {lambda_: 0.4, mu: 0.1, color: r, label: 快速传播}, {lambda_: 0.3, mu: 0.2, color: g, label: 强控制} ] plt.figure(figsize(12,7)) for params in params_combinations: sol odeint(sis_model, [i0], t, args(params[lambda_], params[mu])) plt.plot(t, sol[:,0], params[color], labelparams[label]) plt.xlabel(时间(天)) plt.ylabel(感染比例) plt.title(不同参数下的传播趋势对比) plt.legend() plt.grid() plt.show()典型传播模式分析指数增长期感染人数初期快速增长拐点出现易感人群减少导致传播速度下降平衡状态新增感染与治愈人数达到动态平衡4. 模型进阶与实战技巧4.1 关键指标计算计算基本再生数R0和平衡点R0 lambda_ / mu equilibrium max(0, 1 - 1/R0) print(f基本再生数R0: {R0:.2f}) print(f理论平衡点: {equilibrium:.2%})4.2 实时交互式模拟使用IPython widgets创建参数调节面板from ipywidgets import interact interact( lambda_(0.1, 0.5, 0.05), mu(0.01, 0.3, 0.01), i0(0.001, 0.2, 0.01) ) def interactive_sis(lambda_0.3, mu0.1, i00.01): sol odeint(sis_model, [i0], t, args(lambda_, mu)) plt.figure(figsize(10,6)) plt.plot(t, sol[:,0], b) plt.ylim(0, 1) plt.title(fSIS模型动态模拟 (λ{lambda_}, μ{mu})) plt.show()4.3 实际应用中的注意事项数据校准通过真实疫情数据反推参数空间因素考虑地理分布的异质性干预策略模拟隔离、疫苗接种等措施效果# 完整代码示例 def full_sis_simulation(lambda_0.3, mu0.1, i00.01, days100): t np.linspace(0, days, days*10) solution odeint(sis_model, [i0], t, args(lambda_, mu)) plt.figure(figsize(12,6)) plt.plot(t, solution[:,0], b-, linewidth2) plt.xlabel(Time (days)) plt.ylabel(Fraction Infected) plt.title(fSIS Model: λ{lambda_}, μ{mu}, Initial{i0}) plt.grid() plt.show() print(fFinal infected rate: {solution[-1,0]:.2%}) print(fPredicted equilibrium: {max(0, 1 - mu/lambda_):.2%})在多次项目实践中发现当R0接近1时模型对参数变化极为敏感。曾遇到λ仅变化0.02就导致预测结果差异超过30%的情况这提示我们在实际应用中需要高频更新参数估计。
