当前位置：首页 > news >正文

用Python模拟SIS模型：从微分方程到代码实现，可视化疫情传播全过程

news 2026/5/27 11:19:01

用Python模拟SIS模型：从微分方程到代码实现，可视化疫情传播全过程

在流行病学研究中，SIS（Susceptible-Infectious-Susceptible）模型是描述传染病传播动态的基础工具之一。与SIR模型不同，SIS模型假设康复后的个体不会获得永久免疫力，而是重新变为易感状态——这种特性使其特别适合模拟流感、普通感冒等反复感染的疾病。本文将带您从零开始，用Python构建完整的SIS模型仿真系统，通过代码直观展现λ（感染率）与μ（康复率）如何共同塑造疫情发展曲线。

1. 环境配置与基础准备

开始前需要确保已安装科学计算三件套：NumPy用于数值运算，SciPy提供微分方程求解器，Matplotlib实现可视化。推荐使用Anaconda环境管理工具一键安装：

conda install numpy scipy matplotlib

对于动态可视化，我们额外引入matplotlib.animation模块。以下是基础配置代码：

import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation

定义模型核心参数时，建议采用字典结构便于管理：

params = { 'lambda': 0.3, # 感染率 'mu': 0.1, # 康复率 'i0': 0.01, # 初始感染比例 't_max': 100, # 模拟时长 'dt': 0.1 # 时间步长 }

提示：λ/μ的比值决定疫情是否爆发，当λ/μ > 1时称为基本再生数R0>1，此时疾病会持续传播

2. 微分方程的数值解法

SIS模型的核心微分方程为：

di/dt = λ*i*(1-i) - μ*i

2.1 使用SciPy求解

SciPy的odeint函数能高效求解常微分方程。首先定义方程函数：

def sis_model(i, t, lambda_, mu): return lambda_ * i * (1 - i) - mu * i

接着设置时间点并求解：

t = np.linspace(0, params['t_max'], 1000) solution = odeint(sis_model, params['i0'], t, args=(params['lambda'], params['mu']))

2.2 手动实现欧拉法

为深入理解数值计算过程，我们可以手动实现欧拉方法：

def euler_sis(params): i_values = [params['i0']] t_values = np.arange(0, params['t_max'], params['dt']) for _ in range(1, len(t_values)): di = (params['lambda'] * i_values[-1] * (1 - i_values[-1]) - params['mu'] * i_values[-1]) i_values.append(i_values[-1] + di * params['dt']) return t_values, i_values

两种方法的对比结果可通过下表呈现：

方法	计算速度	精度	适用场景
SciPy odeint	快	高	生产环境、快速原型
欧拉法	慢	中等	教学理解、自定义算法

3. 动态可视化实现

3.1 静态曲线绘制

基础可视化展示感染比例随时间变化：

plt.plot(t, solution, label='Infected proportion') plt.xlabel('Time (days)') plt.ylabel('Proportion of population') plt.title('SIS Model Dynamics') plt.grid(True) plt.legend()

3.2 创建传播动画

动态演示更能直观展示传播过程：

fig, ax = plt.subplots() line, = ax.plot([], [], lw=2) ax.set_xlim(0, params['t_max']) ax.set_ylim(0, 0.5) def init(): line.set_data([], []) return line, def animate(frame): x = t[:frame] y = solution[:frame] line.set_data(x, y) return line, ani = FuncAnimation(fig, animate, frames=len(t), init_func=init, blit=True) plt.close()

将动画保存为GIF或HTML：

ani.save('sis_simulation.gif', writer='pillow', fps=30)

4. 参数敏感性分析

4.1 关键参数影响

通过调整λ和μ观察不同传播场景：

scenarios = [ {'lambda': 0.4, 'mu': 0.2}, # R0=2, 地方性流行 {'lambda': 0.2, 'mu': 0.3}, # R0<1, 疾病消失 {'lambda': 0.5, 'mu': 0.1} # R0=5, 大范围流行 ] plt.figure(figsize=(10,6)) for scenario in scenarios: sol = odeint(sis_model, params['i0'], t, args=(scenario['lambda'], scenario['mu'])) plt.plot(t, sol, label=f"λ={scenario['lambda']}, μ={scenario['mu']}")

4.2 相空间分析

绘制di/dt随i变化的相位图，揭示系统平衡点：

i_range = np.linspace(0, 1, 100) di_dt = params['lambda'] * i_range * (1 - i_range) - params['mu'] * i_range plt.plot(i_range, di_dt) plt.axhline(0, color='gray', linestyle='--') plt.xlabel('Infected proportion (i)') plt.ylabel('Rate of change (di/dt)')

平衡点出现在di/dt=0处，解析解为：

i* = 1 - μ/λ (当λ > μ时)

5. 模型扩展与实践应用

5.1 时变参数处理

现实中的防控措施会使λ随时间变化。修改模型为：

def time_varying_sis(i, t, mu): lambda_t = 0.5 * np.exp(-0.05 * t) # 随时间递减的感染率 return lambda_t * i * (1 - i) - mu * i

5.2 网络结构集成

在复杂网络中的SIS模拟需要networkx库：

import networkx as nx def network_sis(G, lambda_, mu, steps): infected = set(np.random.choice(G.nodes(), int(0.1*G.number_of_nodes()))) for _ in range(steps): new_infected = set() for node in infected: for neighbor in G.neighbors(node): if neighbor not in infected and np.random.rand() < lambda_: new_infected.add(neighbor) recovered = set() for node in infected: if np.random.rand() < mu: recovered.add(node) infected = (infected - recovered) | new_infected yield len(infected) / G.number_of_nodes()