从零到一用Python微分方程模拟传染病传播以SIR模型为例在公共卫生领域传染病传播模型一直是预测疫情发展趋势的重要工具。SIR模型作为经典的传染病动力学模型通过微分方程组描述了易感者(S)、感染者(I)和康复者(R)三类人群的动态变化。本文将带你从零开始用Python实现SIR模型的数值求解与可视化让你不仅能理解模型背后的数学原理还能通过代码直观地观察不同参数下疫情发展的趋势。1. 环境准备与模型基础在开始编码之前我们需要先理解SIR模型的基本假设和数学表达。SIR模型将人群分为三类易感者(Susceptible)可能被感染的健康人群感染者(Infectious)已经感染并具有传染性的人群康复者(Recovered)已经康复并获得免疫力的人群模型的基本微分方程组如下ds/dt -β * s * i di/dt β * s * i - γ * i dr/dt γ * i其中β(感染率) λ(接触率) × p(传染概率)γ(康复率) 1/传染期s, i, r分别表示三类人群占总人口的比例(s i r 1)要运行本文的代码示例需要安装以下Python库pip install numpy scipy matplotlib2. 构建SIR模型的微分方程系统我们将使用SciPy库中的odeint函数来求解这个微分方程组。首先定义模型方程import numpy as np from scipy.integrate import odeint def sir_model(y, t, beta, gamma): s, i, r y dsdt -beta * s * i didt beta * s * i - gamma * i drdt gamma * i return [dsdt, didt, drdt]这里我们创建了一个sir_model函数它接受当前状态y(包含s,i,r值)、时间t和参数beta、gamma返回微分方程组的导数。提示在定义微分方程时参数的顺序必须与odeint要求的顺序一致即先状态变量再时间最后是其他参数。3. 设置初始条件与参数接下来我们需要设置模拟的初始条件和参数值# 初始条件假设总人口为1比例表示 s0 0.99 # 初始易感者比例 i0 0.01 # 初始感染者比例 r0 0.0 # 初始康复者比例 y0 [s0, i0, r0] # 模型参数 beta 0.3 # 感染率 gamma 0.1 # 康复率 # 时间点单位天 t np.linspace(0, 200, 200)参数的选择对模型行为有重要影响参数典型范围含义影响因素β0.1-0.5感染率接触频率、传染概率γ0.05-0.2康复率疾病持续时间4. 求解微分方程并可视化现在我们可以求解微分方程并绘制结果了# 求解微分方程 solution odeint(sir_model, y0, t, args(beta, gamma)) s, i, r solution.T # 绘制结果 import matplotlib.pyplot as plt plt.figure(figsize(10,6)) plt.plot(t, s, labelSusceptible) plt.plot(t, i, labelInfected) plt.plot(t, r, labelRecovered) plt.xlabel(Days) plt.ylabel(Proportion of Population) plt.title(SIR Model Simulation) plt.legend() plt.grid() plt.show()运行这段代码你将看到经典的SIR模型曲线图展示了三类人群随时间的变化趋势。5. 参数敏感性分析理解参数如何影响模型行为至关重要。我们可以创建一组不同参数的模拟进行比较# 定义不同参数组合 param_sets [ {beta: 0.2, gamma: 0.1, color: blue, label: Low spread}, {beta: 0.3, gamma: 0.1, color: green, label: Moderate spread}, {beta: 0.4, gamma: 0.1, color: red, label: High spread} ] plt.figure(figsize(12,8)) for params in param_sets: sol odeint(sir_model, y0, t, args(params[beta], params[gamma])) plt.plot(t, sol[:,1], colorparams[color], labelf{params[label]} (β{params[beta]})) plt.xlabel(Days) plt.ylabel(Infected Proportion) plt.title(Infection Curves with Different β Values) plt.legend() plt.grid() plt.show()这个分析展示了感染率β如何影响疫情发展β值越大感染峰值越高且来得越早β值越小疫情发展越平缓基本再生数R0β/γ决定疫情是否会爆发6. 模型扩展与实际应用基础SIR模型可以进一步扩展以适应更复杂的场景。以下是几个常见的扩展方向加入出生和死亡率def sir_birth_death(y, t, beta, gamma, mu): s, i, r y dsdt mu - beta*s*i - mu*s didt beta*s*i - gamma*i - mu*i drdt gamma*i - mu*r return [dsdt, didt, drdt]加入疫苗接种def sir_vaccination(y, t, beta, gamma, p): s, i, r y dsdt -beta*s*i - p*s didt beta*s*i - gamma*i drdt gamma*i p*s return [dsdt, didt, drdt]拟合实际数据 可以使用scipy.optimize.curve_fit来估计模型参数from scipy.optimize import curve_fit def fit_sir(params, t, beta, gamma): return odeint(sir_model, params[0:3], t, args(beta, gamma))[:,1] # 假设real_data是实际感染数据 popt, pcov curve_fit(fit_sir, t, real_data, p0[0.99,0.01,0.0,0.3,0.1])注意在实际应用中模型参数需要通过流行病学数据估计不同疾病的参数差异可能很大。7. 交互式模拟与参数探索为了更直观地理解参数影响我们可以创建一个交互式模拟from ipywidgets import interact def plot_sir_interactive(beta0.3, gamma0.1, s00.99, i00.01): y0 [s0, i0, 1-s0-i0] sol odeint(sir_model, y0, t, args(beta, gamma)) plt.figure(figsize(10,6)) plt.plot(t, sol[:,0], labelSusceptible) plt.plot(t, sol[:,1], labelInfected) plt.plot(t, sol[:,2], labelRecovered) plt.ylim(0,1) plt.legend() plt.show() interact(plot_sir_interactive, beta(0.1,0.5,0.01), gamma(0.01,0.3,0.01), s0(0.8,0.999,0.001), i0(0.001,0.2,0.001))这个交互式工具允许你实时调整参数并观察曲线变化非常适合教学和探索性分析。8. 模型验证与局限性虽然SIR模型提供了传染病传播的有用洞见但它也有若干局限性简化假设均匀混合人群实际中存在社交网络结构常数参数实际中干预措施会改变β值忽略年龄结构、空间因素等验证方法比较模型预测与实际疫情数据进行参数敏感性分析使用AIC/BIC等准则比较不同模型在实际项目中我通常会先运行基础SIR模型建立直觉然后根据具体需求逐步引入更复杂的因素如潜伏期、年龄分层或空间扩散等。
