企业官网建设流程全解析

1. 项目概述：当SIS模型遇上新冠疫情

如果你是刚开始接触数学建模的Python新手，面对“新冠疫情SIS模型”这个题目，可能会觉得既熟悉又陌生。熟悉的是“新冠疫情”，陌生的是“SIS模型”和如何用Python把它实现出来。这个项目本质上，是让我们用一个经典的传染病动力学模型，去模拟和分析新冠疫情中“感染-康复-再感染”这一可能存在的动态过程。SIS模型是“易感者-感染者-易感者”模型的简称，它描述的是个体在感染康复后，并不能获得永久免疫力，而是会再次变为易感者，重新面临感染风险。这在一些细菌性感染或病毒变异性强的场景中（例如某些冠状病毒的反复感染）有很强的现实意义。

通过这个项目，你将不仅仅学会敲几行Python代码画出一条曲线，更重要的是理解如何将一个现实世界的复杂问题（疫情传播），抽象成一个包含核心参数的数学模型，并利用编程工具进行仿真、分析和可视化。这恰恰是数学建模的核心思想：用数学的语言描述世界，用计算的力量预测趋势。无论你是为了准备数学建模竞赛，还是单纯想提升自己用Python解决实际问题的能力，这个从理论到代码的完整实践过程，都将是一次极佳的练手机会。接下来，我会带你从零开始，拆解SIS模型的每一个细节，并用清晰的Python代码将其实现，同时分享我在建模过程中踩过的坑和总结出的实用技巧。

2. SIS模型的核心原理与参数解读

在动手写代码之前，我们必须先把模型本身的“数学骨架”吃透。SIS模型是一个典型的仓室模型，它把总人口（N）划分为两个互斥的“仓室”：易感者（Susceptible，记为S）和感染者（Infected，记为I）。显然，在任何时刻，都有 S(t) + I(t) = N。模型的核心是描述这两个群体数量如何随时间t变化。

2.1 模型微分方程及其含义

SIS模型的基本微分方程组通常写作：

1. dS/dt = -β * S * I / N + γ * I
2. dI/dt = β * S * I / N - γ * I

这两个方程是模型的心脏，每一部分都有明确的流行病学含义：

• β（感染率）：这是一个关键参数。β * S * I / N 被称为“感染项”或“新发病例数”。它表示单位时间内，新感染发生的数量。为什么是S * I / N呢？这基于一个“均匀混合”的假设：一个感染者单位时间内接触的人数是固定的，其中与易感者接触的比例是 S/N，因此有效接触率是 β * (S/N)。再乘以感染者总数I，就得到了总的新感染数。β综合反映了病毒的传染力（一个感染者能传染多少人）和人群的接触频率。
• γ（康复率）：γ * I 被称为“移出项”（这里是从感染者仓室移出）。它的倒数 1/γ 具有重要的现实意义：它代表平均感染期（Duration of Infection）。例如，如果 γ=0.2/天，那么平均感染期就是 1/0.2 = 5天。这意味着一个感染者平均经过5天会康复（并失去免疫力，变回易感者）。

注意：这里有一个非常重要的点，也是新手容易混淆的地方。在SIR模型中，康复者会进入“移除者（R）”仓室并获得永久免疫。而在SIS模型中，康复者直接回到了“易感者（S）”仓室。因此，方程(1)中有一个“+γI”项，表示康复者又重新加入了易感者大军。这是SIS与SIR最本质的区别。

2.2 基本再生数R0与模型动力学

在传染病模型中，有一个灵魂参数叫基本再生数（Basic Reproduction Number, R0）。对于SIS模型，R0 = β / γ。它的流行病学意义是：在一个全部是易感者的人群中，一个感染者在整个传染期内平均能感染的人数。

• R0 > 1：这意味着一个感染者在其传染期内，能感染超过一个人。感染项（βSI/N）的力量大于康复项（γI），疫情会扩散，最终系统会达到一个地方病平衡点（Endemic Equilibrium），即感染者和易感者以一定的比例长期共存，疫情不会消失。这也是SIS模型常用于描述慢性、反复流行疾病的原因。
• R0 < 1：这意味着疾病无法在人群中持续传播，感染项弱于康复项，疫情会逐渐衰减至消失。

在SIS模型中，地方病平衡点时感染者的比例 I*/N 可以直接由R0计算出来：I*/N = 1 - 1/R0。这个公式非常直观：R0越大，最终稳态的感染者比例就越高。例如，若R0=2，则最终会有 1 - 1/2 = 50% 的人口处于感染状态（在模型假设下）。这个公式为我们后续验证代码的正确性提供了一个黄金标准。

2.3 针对新冠疫情的参数考量

将SIS模型应用于新冠疫情分析时，我们需要对参数进行审慎的思考和估计，这本身也是建模的一部分。

• 感染率 β：新冠的β值变化极大，它强烈依赖于防控措施（如口罩、社交距离、封锁）。在疫情初期，没有干预的情况下，估算的R0大约在2.5-3.5之间。我们可以通过假设一个平均感染期来反推β。例如，假设平均感染期1/γ=10天（即γ=0.1），若取R0=3.0，则 β = R0 * γ = 3.0 * 0.1 = 0.3/天。
• 康复率 γ：如前所述，γ=1/(平均感染期)。对于新冠，从感染到具有传染性再到康复/不再排毒的时间窗较为复杂，通常简化估计为7-14天。我们常取γ在0.07~0.14/天之间。
• 模型局限性：必须清醒认识到，经典的SIS模型对新冠的模拟是高度简化的。它忽略了潜伏期、无症状感染、年龄结构、空间异质性、病毒变异导致的免疫逃逸（这恰恰是SIS适用的点）、以及更复杂的防控政策等因素。因此，我们的模拟更多是原理性展示和教学目的，揭示“感染-康复-再感染”这一循环的动态特性，而非精确预测。

3. Python实现：从方程到动态模拟

理解了模型原理，我们就可以用Python这个强大的工具，让静态的方程“动起来”。我们将使用SciPy库来数值求解微分方程，用Matplotlib进行可视化。这是数学建模中非常标准且高效的流程。

3.1 环境准备与库安装

首先确保你的Python环境已经安装了必要的科学计算库。打开你的终端或命令提示符，使用pip进行安装：

pip install numpy scipy matplotlib

• numpy：提供高效的数组运算，是科学计算的基础。
• scipy：特别是其中的integrate模块，提供了求解微分方程的函数，如odeint或solve_ivp，我们将使用后者，因为它更现代、功能更丰富。
• matplotlib：绘图库，用于将模拟结果可视化。

实操心得：建议使用Jupyter Notebook或VS Code等支持交互和分步执行的开发环境来做这种探索性建模。你可以方便地修改参数，立即看到图形结果，非常适合调试和优化模型。

3.2 定义模型微分方程函数

这是最关键的一步，我们将微分方程组翻译成Python函数。我们使用solve_ivp函数，它要求我们定义一个函数，输入是时间t和状态变量y，输出是导数dy/dt。

import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt def sis_model(t, y, beta, gamma, N): """ 定义SIS模型的微分方程。 参数: t: 时间（solve_ivp自动传入，此处未显式使用，但函数签名需要） y: 状态变量数组，y[0]=S, y[1]=I beta: 感染率 gamma: 康复率 N: 总人口 返回: dydt: 导数数组，[dS/dt, dI/dt] """ S, I = y # 计算导数 dS_dt = -beta * S * I / N + gamma * I dI_dt = beta * S * I / N - gamma * I return [dS_dt, dI_dt]

代码解读：

1. 函数sis_model严格对应了我们之前的微分方程。
2. 参数y是一个包含两个元素的列表（或数组），分别代表当前时刻的S和I。
3. 计算dS_dt和dI_dt时，我们直接套用了公式。注意beta * S * I / N这一项在两个方程中符号相反。
4. 函数返回的也是一个列表[dS_dt, dI_dt]，告诉求解器每个状态变量的变化速度。

3.3 设置参数与初始条件并进行模拟

接下来，我们设定具体的参数、初始条件，并调用solve_ivp进行求解。

# 模型参数 N = 1000 # 总人口 beta = 0.3 # 感染率，/天 gamma = 0.1 # 康复率，/天 R0 = beta / gamma # 基本再生数 print(f"基本再生数 R0 = {R0:.2f}") # 初始条件：假设初始有10个感染者，其余为易感者 I0 = 10 S0 = N - I0 y0 = [S0, I0] # 求解器需要的初始状态向量 # 模拟时间范围 (天) t_start = 0 t_end = 150 # 模拟150天 t_eval = np.linspace(t_start, t_end, 301) # 在0-150天内均匀取301个时间点输出结果 # 使用solve_ivp求解微分方程 solution = solve_ivp( sis_model, [t_start, t_end], y0, args=(beta, gamma, N), # 传递给sis_model的额外参数 t_eval=t_eval, # 指定输出时间点 method='RK45', # 龙格-库塔法，默认且通常足够精确 rtol=1e-6, atol=1e-9 # 相对和绝对误差容限，控制精度 ) # 提取结果 t = solution.t S = solution.y[0] I = solution.y[1]

参数设置详解：

• N=1000：我们模拟一个1000人的封闭群体。
• beta=0.3, gamma=0.1：如前所述，这给出了R0=3.0。我们预期疫情会爆发并达到地方病平衡。
• I0=10：模拟疫情初期引入少量感染者的情况。
• t_eval：我们并不需要每一天的解，np.linspace生成了一个从0到150共301个点的数组（包含首尾），意味着大约每0.5天输出一个结果。这让我们绘制的曲线更平滑。
• solve_ivp参数：args用于传递模型函数所需的额外参数（beta, gamma, N）。RK45是经典的四阶/五阶龙格-库塔法，适用于大多数常微分方程。rtol和atol是精度控制参数，在大多数情况下默认值即可，对于学术研究或需要高精度对比时，可以适当收紧。

3.4 结果可视化与分析

得到数据后，可视化是理解模型行为最直观的方式。

# 创建图形和坐标轴 plt.figure(figsize=(12, 8)) # 绘制易感者(S)和感染者(I)数量随时间变化曲线 plt.plot(t, S, 'b-', linewidth=2, label=f'Susceptible (S)') plt.plot(t, I, 'r-', linewidth=2, label=f'Infected (I)') # 计算并绘制理论上的地方病平衡点（稳态） I_endemic = N * (1 - 1/R0) if R0 > 1 else 0 S_endemic = N - I_endemic plt.axhline(y=S_endemic, color='b', linestyle='--', alpha=0.7, label=f'S Endemic ({S_endemic:.0f})') plt.axhline(y=I_endemic, color='r', linestyle='--', alpha=0.7, label=f'I Endemic ({I_endemic:.0f})') # 图表装饰 plt.xlabel('Time (days)', fontsize=12) plt.ylabel('Number of People', fontsize=12) plt.title(f'SIS Model Simulation for COVID-19-like Disease\n$\\beta$={beta}, $\\gamma$={gamma}, $R_0$={R0:.2f}, N={N}, I0={I0}', fontsize=14) plt.legend(loc='best', fontsize=11) plt.grid(True, which='both', linestyle='--', alpha=0.5) plt.xlim([t_start, t_end]) # 可以添加第二个子图，展示感染者比例 plt.figure(figsize=(12, 4)) plt.plot(t, I/N, 'g-', linewidth=2, label='Infected Fraction (I/N)') plt.axhline(y=1 - 1/R0, color='g', linestyle='--', alpha=0.7, label=f'Theoretical Equilibrium ({1-1/R0:.3f})') plt.xlabel('Time (days)', fontsize=12) plt.ylabel('Fraction of Population', fontsize=12) plt.title('Infected Fraction Over Time', fontsize=14) plt.legend(loc='best', fontsize=11) plt.grid(True, which='both', linestyle='--', alpha=0.5) plt.xlim([t_start, t_end]) plt.tight_layout() plt.show()

图形解读： 运行代码后，你会看到两张图。第一张图展示了S和I的绝对数量变化。你可以清晰地看到：

1. 疫情爆发期：初期，感染者I迅速上升，易感者S快速下降。
2. 达到峰值：感染人数达到一个峰值。
3. 趋于平衡：随后，两条曲线不再大幅波动，而是逐渐趋于水平，并最终稳定在两条虚线标注的理论平衡点附近。S和I都不再为零，这就是地方病平衡。在我们的参数下（R0=3），最终大约有667人易感，333人感染（因为 I*/N = 1 - 1/3 ≈ 0.333）。

第二张图专门展示感染者比例I/N，它更直观地反映了疫情的严重程度，并且其平衡值直接等于1 - 1/R0。图中绿色虚线标出了这个理论值，模拟曲线最终与之重合，这很好地验证了我们代码实现的正确性。

4. 深入分析与模型探索

一个基本的模拟跑通了，但建模的魅力在于探索。我们可以通过改变参数，来回答一系列“如果...那么...”的问题。

4.1 探究不同R0对疫情的影响

R0是决定疫情走向的“开关”。我们可以通过固定康复率γ，改变感染率β来获得不同的R0，进行批量模拟。

# 定义一组R0值 R0_values = [0.8, 1.5, 2.5, 4.0] gamma_fixed = 0.1 beta_values = [R0 * gamma_fixed for R0 in R0_values] # 根据R0计算对应的beta # 设置统一的初始条件和时间范围 N = 1000 I0 = 10 t_end = 150 plt.figure(figsize=(14, 8)) for i, (R0, beta) in enumerate(zip(R0_values, beta_values)): # 求解模型 sol = solve_ivp(sis_model, [0, t_end], [N-I0, I0], args=(beta, gamma_fixed, N), t_eval=np.linspace(0, t_end, 301)) # 绘制感染者比例曲线 plt.plot(sol.t, sol.y[1]/N, linewidth=2, label=f'$R_0$={R0} ($\\beta$={beta:.2f})') # 标注理论平衡点 if R0 > 1: eq_point = 1 - 1/R0 plt.axhline(y=eq_point, color=plt.gca().lines[-1].get_color(), linestyle=':', alpha=0.5) plt.xlabel('Time (days)', fontsize=12) plt.ylabel('Fraction Infected (I/N)', fontsize=12) plt.title('Impact of Basic Reproduction Number $R_0$ on SIS Model Dynamics', fontsize=14) plt.legend(loc='upper right', fontsize=11) plt.grid(True, linestyle='--', alpha=0.5) plt.ylim([0, 0.6]) # 限制y轴范围以便观察 plt.show()

分析结果：

• R0=0.8 (<1)：感染者比例单调下降至0，疾病最终消失。这是我们可以通过强有力干预（降低β）达到的理想状态。
• R0=1.5, 2.5, 4.0 (>1)：疫情都会爆发并达到一个稳定的地方病平衡。R0越大，疫情峰值越高，达到峰值越快，并且最终的稳态感染比例也越高。从图中可以清晰看到，R0=4.0的曲线不仅峰值最高，其平衡点也接近0.75（即75%的人口感染），远高于R0=1.5时的平衡点（约33%）。

4.2 模拟干预措施：动态变化的β

现实中的防控措施（如封控、戴口罩、接种疫苗）会改变感染率β。我们可以让β随时间变化，来模拟干预的引入和解除。

def sis_model_dynamic_beta(t, y, gamma, N): """SIS模型，但感染率beta是时间的函数""" S, I = y # 定义随时间变化的beta：前50天无干预，第50天起强力干预，第100天干预放松 if t < 50: beta = 0.3 # 原始高传播 elif t < 100: beta = 0.1 # 强力干预，降低传播 else: beta = 0.2 # 干预放松，传播部分回升 dS_dt = -beta * S * I / N + gamma * I dI_dt = beta * S * I / N - gamma * I return [dS_dt, dI_dt] # 参数 N = 1000 I0 = 10 gamma = 0.1 t_end = 200 # 求解 sol_dynamic = solve_ivp(sis_model_dynamic_beta, [0, t_end], [N-I0, I0], args=(gamma, N), t_eval=np.linspace(0, t_end, 401), max_step=0.5) # 可视化 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10), sharex=True) # 子图1：感染者比例 ax1.plot(sol_dynamic.t, sol_dynamic.y[1]/N, 'r-', linewidth=2, label='Infected Fraction (I/N)') ax1.set_ylabel('Fraction Infected', fontsize=12) ax1.set_title('SIS Model with Time-Varying Intervention (Dynamic $\\beta$)', fontsize=14) ax1.legend(loc='upper right') ax1.grid(True, linestyle='--', alpha=0.5) # 在图上标注干预阶段 ax1.axvspan(0, 50, color='yellow', alpha=0.1, label='No Intervention') ax1.axvspan(50, 100, color='green', alpha=0.1, label='Strong Intervention') ax1.axvspan(100, 200, color='orange', alpha=0.1, label='Relaxed Intervention') ax1.legend(loc='upper left') # 子图2：动态的beta值 beta_t = [] for t in sol_dynamic.t: if t < 50: beta_t.append(0.3) elif t < 100: beta_t.append(0.1) else: beta_t.append(0.2) ax2.plot(sol_dynamic.t, beta_t, 'b-', linewidth=2, drawstyle='steps-post') # steps-post能更好显示阶跃变化 ax2.set_xlabel('Time (days)', fontsize=12) ax2.set_ylabel('Infection Rate ($\\beta$)', fontsize=12) ax2.set_ylim([0, 0.35]) ax2.grid(True, linestyle='--', alpha=0.5) plt.tight_layout() plt.show()

模拟解读： 这张图生动地展示了干预措施的效果。

1. 0-50天（无干预）：β=0.3 (R0=3)，疫情迅速爆发，感染者比例冲向平衡点（约33%）。
2. 50-100天（强力干预）：β骤降至0.1 (R0=1)，此时感染项与康复项力量发生逆转。你可以看到感染者比例开始快速下降，因为新的感染速度已经赶不上康复速度。
3. 100-200天（干预放松）：β回升至0.2 (R0=2)。由于人群中仍有相当数量的易感者，疫情再次抬头，但上升速度比第一阶段慢，并最终在一个新的、更低的平衡点（I*/N = 1 - 1/2 = 0.5）附近稳定下来。

这个模拟清晰地告诉我们，非药物干预措施（NPIs）通过降低β（即R0）来控制疫情是立竿见影的。但一旦放松，疫情可能反弹。这也部分解释了现实中疫情为何会呈现“波浪形”发展。

5. 模型扩展与思考

基础的SIS模型是一个强大的起点，但现实世界更复杂。作为建模者，我们可以尝试扩展它，使其更贴近实际。

5.1 考虑出生与死亡：SI模型

在经典SIS模型中，我们通常假设总人口N不变。但若考虑一个更长期的视角，可以引入自然出生率和死亡率。假设出生率Λ，自然死亡率μ（所有仓室相同），则模型变为： dS/dt = ΛN - βSI/N - μS + γI dI/dt = βSI/N - γI - μI 这个模型被称为有出生死亡的SIS模型，或简称SI模型（在一些文献中）。它的平衡态分析会更复杂一些，但基本结论不变：当R0 > 1时，疾病会形成地方性流行。在Python实现上，只需修改sis_model函数中的导数计算公式即可。

5.2 离散时间模型与随机性

我们目前使用的是连续时间的确定性微分方程模型。另一种思路是建立离散时间随机模型。例如，我们可以用“天”作为时间步长，根据概率来模拟每个个体每天的状态转移（易感->感染，感染->康复易感）。这需要用到随机模拟（如蒙特卡洛方法）。虽然计算量更大，但它能更好地捕捉小群体中的随机波动现象，例如疾病在传播早期由于随机性而偶然消失的可能性。这可以通过Python的numpy.random模块来实现状态转移的随机抽样。

5.3 结合真实数据进行参数估计

一个更高级、也更有挑战性的方向是参数估计。如果我们有某地区新冠疫情的时间序列数据（如每日新增感染数），我们可以利用这些数据来反推模型的参数β和γ。这通常通过优化算法来实现：定义一个目标函数（如模型预测的新增感染数与实际新增感染数之间的误差平方和），然后使用SciPy的optimize模块（例如curve_fit或minimize函数）来寻找使误差最小的β和γ值。这个过程能让你真正将模型与数据连接起来，但需要注意模型假设与数据真实性的匹配度。

6. 常见问题、调试技巧与心得

在实际动手编码和模拟的过程中，你几乎一定会遇到下面这些问题。这里我把自己踩过的坑和解决方法总结一下。

6.1 数值求解器不收敛或结果异常

• 问题表现：运行solve_ivp后，得到的解出现NaN（非数字），或者曲线剧烈震荡、发散。
• 排查步骤：
  1. 检查微分方程函数：这是最常见的问题源。仔细核对sis_model函数中的每一个符号（正负号）和运算顺序。确保dS/dt + dI/dt = 0（在无出生死亡的情况下），这是一个很好的快速验证方法。打印几个时间点的导数看看是否合理。
  2. 检查参数范围：β和γ通常应为正数。总人口N应为正。初始感染数I0应小于N。
  3. 调整求解器参数：尝试使用不同的求解方法（method），对于刚性问题（某些参数下方程变化速率差异极大），RK45可能不稳定，可以尝试Radau或BDF方法。适当减小最大步长（max_step），或收紧误差容限（rtol,atol）。
  4. 简化问题：先尝试用一组非常简单的参数（如β很小，γ很大，使R0<1）测试，看模型是否能正确模拟疾病消失。然后再逐步调整到目标参数。

6.2 结果与理论平衡点不符

• 问题表现：模拟曲线稳定后的值，与计算出的理论平衡点I* = N*(1-1/R0)有较大差距。
• 可能原因与解决：
  1. 模拟时间不够长：系统可能尚未达到稳态。尝试延长t_end，比如模拟500天或1000天。
  2. 理论公式应用错误：确保你计算理论平衡点时使用的是正确的公式，并且R0 = β/γ。再次确认你的β和γ值。
  3. 模型实现有误：这是根本原因。回到6.1的第一步，彻底检查微分方程代码。一个有效的debug方法是：在模拟结束后，计算最后时刻的β * S * I / N和γ * I。在平衡点附近，这两个值应该近似相等（因为dI/dt ≈ 0）。如果不相等，说明你的方程或参数有问题。

6.3 如何提高代码的可复用性和可读性

对于初学者，把一切写在一个脚本里没问题。但当项目复杂后，良好的代码习惯至关重要。

• 将模型定义、参数设置、模拟运行、可视化分离：可以写成不同的函数或放在不同的代码块中。例如：

  def run_sis_simulation(beta, gamma, N, I0, days): # 封装求解过程 ... return t, S, I def plot_results(t, S, I, beta, gamma, N): # 封装绘图过程 ...

• 使用字典管理参数：将所有参数放在一个字典里，便于管理和传递。

  params = { 'beta': 0.3, 'gamma': 0.1, 'N': 1000, 'I0': 10, 't_end': 150 } # 调用函数时使用 **params 解包 t, S, I = run_sis_simulation(**params)

• 为函数和变量添加清晰的注释和文档字符串：就像我在示例代码中做的那样。几个月后回头看，你会感谢自己的。

6.4 模型局限性与应用思考

在完成这个项目后，务必清醒地认识到这个简单SIS模型的局限性，并思考如何改进：

• 忽略潜伏期：新冠有显著的潜伏期（Exposed阶段），SEIS或SEIR模型更合适。
• 忽略无症状感染：无症状感染者具有传染性但行为模式不同，需要更复杂的仓室划分。
• 同质混合假设：模型假设人群完全均匀混合，这与现实中的社交网络结构不符。网络模型能更好地描述这一点。
• 忽略空间因素：疫情传播有地理空间上的扩散过程，需要元胞自动机或反应扩散方程。
• 忽略免疫逃逸：虽然SIS假设康复后无免疫力，但新冠的免疫逃逸是复杂的，并非简单的“康复即易感”，免疫力会随时间衰减或对变异株无效。

因此，这个SIS项目更像是一个教学沙盒和思维起点。它教会你传染病建模的基本范式：定义仓室、建立方程、设置参数、数值求解、分析结果。当你掌握了这个流程，再去学习更复杂的SEIR、网络模型、甚至基于智能体的模拟（ABM）时，就会有一个坚实的认知基础。真正的数学建模，是从这个简单的“玩具模型”出发，根据面对的具体问题，不断引入新的因素，权衡复杂性与准确性，一步步构建起更有解释力和预测力的模型的过程。