SIR传染病模型MATLAB实战：从微分方程建模到参数分析

1. 从一道例题开始：为什么微分方程是数学建模的“灵魂”

如果你参加过数学建模竞赛，或者处理过任何涉及动态变化、趋势预测的工程问题，那你一定绕不开“微分方程”这四个字。它听起来有点吓人，像是数学系高年级的专属课程，但实际上，它可能是连接你的问题描述和计算机求解之间最直接、最有力的桥梁。很多人学了一堆解法，背了一堆公式，但一到实际建模，还是无从下手——问题往往不在于数学本身，而在于如何把现实世界那个模糊的“变化关系”，翻译成微分方程这个精确的数学语言。

今天，我们不空谈理论，就从一道非常经典的建模例题出发，手把手带你走完“问题分析 -> 建立方程 -> MATLAB求解 -> 结果分析”的全过程。这道题就是**传染病模型（SIR模型）**的一个简化版。它几乎出现在每一本建模教材里，因为它完美地诠释了微分方程如何刻画一个系统中个体状态的动态转移。我们的目标很明确：不仅给你能直接运行的MATLAB源码，更要彻底讲清楚每一行代码背后的建模思想和数学原理，让你下次遇到类似问题，能自己“造轮子”，而不是只会“调包”。

你会发现，掌握微分方程建模，核心不是解方程（计算机可以代劳），而是定义方程。这就像编程，语法不难，难的是算法设计。我们通过SIR模型这个“麻雀”，把它彻底解剖。关键词很明确：数学建模、微分方程、MATLAB、源码。这篇文章就是为需要将理论付诸实践的理工科学生、工程师，以及所有希望用数学工具量化分析动态过程的朋友准备的。

2. SIR传染病模型：如何将一场疫情“翻译”成数学方程

在开始写代码之前，我们必须把模型搞清楚。SIR模型是Kermack和McKendrick在1927年提出的，至今仍是描述传染病动力学的基石。它的核心思想是将人群划分为三个“舱室”，并通过微分方程描述其间的流动。

2.1 模型假设与状态定义

任何建模的第一步都是做出合理且必要的简化假设。对于经典的SIR模型，我们假设：

1. 总人口数N恒定。不考虑出生、死亡、迁移等人口变动。这是一个很强的假设，但在疫情短期爆发期内（比如几个月），可以近似成立。
2. 均匀混合。人群中的个体充分接触，任何一个易感者（S）遇到任何一个感染者（I）的概率是相同的。这忽略了社交网络、空间结构的影响。
3. 康复者获得永久免疫。一旦从感染中恢复，将终身免疫，且不再具有传染性。
4. 疾病传播和康复过程是随机的，但在大群体中可以用平均速率来描述。

基于以上假设，我们定义三个核心状态变量，它们都是关于时间t的函数：

• S(t): 易感者数量。指未患病但有可能被感染的人群。
• I(t): 感染者数量。指已患病且具有传染性的人群。
• R(t): 康复者（或移出者）数量。指从感染中恢复并获得永久免疫的人群。

显然，在任何时刻，有S(t) + I(t) + R(t) = N（常数）。

2.2 构建微分方程：关键在于定义“流率”

微分方程描述的是变化率。所以，我们要问：每个群体数量的变化速度由什么决定？

• 易感者S的减少：S只会因为被感染而减少。一个新感染的发生，需要一个易感者（S）和一个感染者（I）接触。在均匀混合假设下，单位时间内一个感染者所能接触到的总人数是固定的，其中易感者的比例是S/N。因此，一个感染者单位时间内能新感染的人数正比于(S/N)。假设感染者的有效接触率为β（这是一个关键参数，表示传染力度），那么一个感染者单位时间新感染的人数为β * (S/N)。现在有I(t)个感染者，所以单位时间内总的新感染人数，即S的减少速率为β * (S/N) * I。由于S在减少，所以其变化率是负的：dS/dt = -β * (S/N) * I

• 感染者I的变化：I的人数增加来源于S被感染，减少来源于患者康复（或死亡，在本模型中并入R）。单位时间内，新增感染者就是上面计算的β * (S/N) * I。同时，假设感染者以固定的平均速率γ康复，即单位时间内康复的人数占感染者总数的比例是γ，所以康复人数为γ * I。因此，感染者的净变化率是：dI/dt = β * (S/N) * I - γ * I这个方程非常关键，它决定了疫情是爆发（dI/dt > 0）还是消退（dI/dt < 0）。

• 康复者R的增加：R只增加，不减少。其增加的速度就是感染者康复的速度：dR/dt = γ * I

至此，我们得到了经典的SIR模型方程组：

dS/dt = -β * (S/N) * I dI/dt = β * (S/N) * I - γ * I dR/dt = γ * I

其中，β称为感染率，γ称为康复率。它们的倒数1/γ具有明确的物理意义：平均感染期，即一个病人平均具有传染性的时间。

注意：这里β和γ的量纲都是1/时间（例如1/天）。β综合反映了病毒的传染能力和人群的接触频率。在后续的MATLAB实现中，我们将直接为这两个参数赋值。

2.3 模型的核心洞察：基本再生数 R0

在求解方程之前，我们可以直接从方程中获得一个极其重要的阈值参数——基本再生数 R0。它表示在一个完全易感的人群中（S ≈ N），一个感染者在其整个传染期内平均能感染多少人。

推导很简单：在疫情初期，S ≈ N， I很小。此时，一个感染者单位时间内能感染β * (N/N) = β个人。他的平均传染期是1/γ。所以，在整个传染期内，他能感染的人数R0 = β * (1/γ) = β / γ。

• 如果 R0 < 1：意味着每个感染者平均感染不到一个人，疫情将无法持续，会逐渐消失。从方程看，此时dI/dt ≈ (β - γ)I = γ(R0 - 1)I < 0，感染者数量会指数衰减。
• 如果 R0 > 1：意味着每个感染者平均能感染超过一个人，疫情可能爆发。初期dI/dt > 0，感染者数量会指数增长。

R0是流行病学中衡量传染病传播能力的最核心指标。我们的MATLAB模拟将直观地展示不同R0下的疫情发展轨迹。

3. MATLAB实战：从方程到代码的完整实现

理论模型建立后，接下来就是用MATLAB将它“复活”。我们将分步实现：首先定义微分方程系统，然后选择合适的数值方法求解，最后可视化结果并进行分析。

3.1 准备工作：定义方程函数（odefun）

MATLAB求解常微分方程组主要使用ode45等求解器，它们要求我们将方程组写成一个函数文件。这个函数返回的是每个状态变量的导数（变化率）。

我们创建一个名为sir_ode.m的函数文件。

function dydt = sir_ode(t, y, beta, gamma, N) % SIR模型微分方程函数 % 输入： % t: 时间（求解器自动传入，此处未显式使用，因为方程是自治的） % y: 状态向量 [S; I; R] % beta: 感染率 % gamma: 康复率 % N: 总人口 % 输出： % dydt: 导数向量 [dS/dt; dI/dt; dR/dt] S = y(1); I = y(2); % R = y(3); % 在计算导数时R不需要单独提取 % 计算变化率 dS_dt = -beta * (S / N) * I; dI_dt = beta * (S / N) * I - gamma * I; dR_dt = gamma * I; % 组装输出向量 dydt = [dS_dt; dI_dt; dR_dt]; end

关键点解析：

1. t是时间变量，虽然我们的方程不明显依赖于t（称为自治系统），但MATLAB的ODE求解器函数接口要求必须有前两个参数(t, y)。
2. y是一个列向量，包含了所有状态变量在某个时刻t的值。我们约定其顺序为[S; I; R]。
3. 我们将参数beta,gamma,N作为额外的输入参数传入，而不是在函数内部写死，这样使得函数更灵活，便于后续进行参数研究。
4. 方程的实现就是直接翻译我们上一节推导的数学公式。

3.2 主脚本：配置参数、调用求解器并绘图

接下来，我们编写主脚本run_sir_model.m来驱动整个模拟过程。

%% 清空环境 clear; close all; clc; %% 1. 设置模型参数 N = 1000; % 总人口 I0 = 1; % 初始感染者人数 R0_init = 0; % 初始康复者人数 S0 = N - I0 - R0_init; % 初始易感者人数 beta = 0.3; % 感染率 (每人每天) gamma = 0.1; % 康复率 (1/天) % 计算基本再生数 R0 R0 = beta / gamma; fprintf('基本再生数 R0 = %.2f\n', R0); %% 2. 设置时间跨度 tspan = [0, 150]; % 模拟从第0天到第150天 %% 3. 定义初始状态向量 y0 = [S0; I0; R0_init]; %% 4. 求解微分方程组 % 使用ode45求解器，通过匿名函数将额外参数(beta, gamma, N)传递给方程函数 [t, y] = ode45(@(t, y) sir_ode(t, y, beta, gamma, N), tspan, y0); % 提取结果 S = y(:, 1); I = y(:, 2); R = y(:, 3); %% 5. 可视化结果 figure('Position', [100, 100, 1200, 500]); % 设置图形窗口大小 % 子图1：三类人群数量随时间变化 subplot(1, 2, 1); plot(t, S, 'b-', 'LineWidth', 2); hold on; plot(t, I, 'r-', 'LineWidth', 2); plot(t, R, 'g-', 'LineWidth', 2); hold off; grid on; xlabel('时间 (天)'); ylabel('人数'); title(sprintf('SIR模型动态 (N=%d, \\beta=%.2f, \\gamma=%.2f, R_0=%.2f)', N, beta, gamma, R0)); legend('易感者 S(t)', '感染者 I(t)', '康复者 R(t)', 'Location', 'best'); % 子图2：感染者比例变化，更直观看到疫情峰值 subplot(1, 2, 2); plot(t, I/N, 'r-', 'LineWidth', 2); grid on; xlabel('时间 (天)'); ylabel('感染者比例'); title('感染者比例 I(t)/N 变化曲线'); % 标记峰值点 [maxI, idx] = max(I); peak_time = t(idx); peak_proportion = maxI / N; hold on; plot(peak_time, peak_proportion, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); text(peak_time+5, peak_proportion, sprintf('峰值: %.1f天, %.1f%%', peak_time, peak_proportion*100), ... 'VerticalAlignment', 'bottom'); hold off; %% 6. 输出关键指标 fprintf('模拟结果摘要：\n'); fprintf(' 最终易感者比例: %.2f%%\n', S(end)/N*100); fprintf(' 最终康复者比例: %.2f%%\n', R(end)/N*100); fprintf(' 疫情峰值发生在第 %.1f 天\n', peak_time); fprintf(' 峰值感染人数为 %.0f (占总人口 %.1f%%)\n', maxI, peak_proportion*100);

代码逐段解读：

1. 参数部分：我们设定了一个1000人的封闭社区，从1个感染者开始。beta=0.3,gamma=0.1，计算得R0=3，这是一个典型的高传染性疾病（如麻疹、COVID-19原始毒株的估计值），预示着疫情会爆发。
2. 求解器调用：ode45是MATLAB中最常用的非刚性常微分方程求解器，它采用Runge-Kutta方法，对于大多数问题精度和效率都很好。@(t,y) ...是匿名函数，它将我们预设的参数beta, gamma, N“绑定”到sir_ode函数上，以满足求解器对函数签名的要求（只接受t, y两个输入）。
3. 结果提取：ode45返回两个数组：t（时间点）和y（每个时间点对应的状态向量）。y是一个三列的矩阵，我们按列提取出S, I, R。
4. 可视化：我们创建了两个子图。左图展示S、I、R三类人群的绝对数量变化，右图专门展示感染者比例I/N，并标记出疫情峰值（感染人数最多的时刻），这对于评估医疗系统压力至关重要。
5. 关键指标输出：在命令行窗口打印出R0、最终人群分布和峰值信息，便于定量分析。

运行这个脚本，你将看到两条曲线生动地展示了疫情从发生、发展到消退的全过程。易感者S从初始值单调下降，感染者I先上升达到峰值后下降，康复者R单调上升直至疫情结束。最终，并非所有人都会被感染，会有一部分人始终是易感者（S(∞) > 0），这是SIR模型的一个重要结论。

4. 模型探究与参数敏感性分析：改变R0会发生什么？

一个模型的价值不仅在于能模拟一种情况，更在于能让我们进行“思想实验”，探究不同条件下的结果。我们最关心的参数就是R0。让我们修改主脚本，来模拟不同R0下的疫情发展。

我们在主脚本后面添加以下代码，或者新建一个分析脚本：

%% 参数敏感性分析：不同R0的影响 % 固定平均感染期 1/gamma = 10天 gamma_fixed = 0.1; beta_values = [0.05, 0.15, 0.25, 0.35]; % 对应 R0 = 0.5, 1.5, 2.5, 3.5 R0_values = beta_values / gamma_fixed; % 准备绘图 figure('Position', [100, 100, 1000, 800]); colors = lines(length(beta_values)); % 获取一组区分度好的颜色 % 循环模拟不同beta (R0) for i = 1:length(beta_values) beta_current = beta_values(i); R0_current = R0_values(i); % 求解方程 [t, y] = ode45(@(t, y) sir_ode(t, y, beta_current, gamma_fixed, N), tspan, y0); I_ratio = y(:, 2) / N; % 感染者比例 % 绘图 subplot(2, 2, 1); % 感染者比例对比 plot(t, I_ratio, '-', 'Color', colors(i, :), 'LineWidth', 1.5, 'DisplayName', sprintf('R_0=%.1f', R0_current)); hold on; subplot(2, 2, 2); % 最终感染规模对比（总康复者比例） final_R = y(end, 3) / N; bar(i, final_R*100, 'FaceColor', colors(i, :)); hold on; subplot(2, 2, 3); % 峰值感染比例对比 peak_I = max(y(:, 2)) / N; bar(i, peak_I*100, 'FaceColor', colors(i, :)); hold on; subplot(2, 2, 4); % 峰值到达时间对比 [~, idx] = max(y(:, 2)); peak_time = t(idx); bar(i, peak_time, 'FaceColor', colors(i, :)); hold on; fprintf('R0=%.1f: 最终感染比例=%.1f%%，峰值感染比例=%.1f%%，峰值时间=%.1f天\n', ... R0_current, final_R*100, peak_I*100, peak_time); end % 完善第一个子图（感染者比例曲线） subplot(2, 2, 1); hold off; grid on; xlabel('时间 (天)'); ylabel('感染者比例 I/N'); title('不同R0下的感染者比例变化'); legend('show', 'Location', 'northeast'); % 完善第二个子图（最终感染规模） subplot(2, 2, 2); hold off; grid on; set(gca, 'XTick', 1:length(R0_values), 'XTickLabel', arrayfun(@(x) sprintf('%.1f', x), R0_values, 'UniformOutput', false)); xlabel('基本再生数 R_0'); ylabel('最终感染比例 (%)'); title('总感染规模 vs R_0'); ylim([0 100]); % 完善第三个子图（峰值感染比例） subplot(2, 2, 3); hold off; grid on; set(gca, 'XTick', 1:length(R0_values), 'XTickLabel', arrayfun(@(x) sprintf('%.1f', x), R0_values, 'UniformOutput', false)); xlabel('基本再生数 R_0'); ylabel('峰值感染比例 (%)'); title('疫情峰值强度 vs R_0'); ylim([0 50]); % 根据模拟结果调整 % 完善第四个子图（峰值到达时间） subplot(2, 2, 4); hold off; grid on; set(gca, 'XTick', 1:length(R0_values), 'XTickLabel', arrayfun(@(x) sprintf('%.1f', x), R0_values, 'UniformOutput', false)); xlabel('基本再生数 R_0'); ylabel('峰值时间 (天)'); title('疫情达峰时间 vs R_0');

分析解读： 这段代码固定康复率gamma=0.1（即平均感染期10天），通过改变感染率beta来获得不同的R0（0.5, 1.5, 2.5, 3.5）。我们同时观察四个关键指标：

1. 感染者比例曲线：直观展示不同传播速度下的疫情发展轨迹。R0=0.5<1时，疫情很快熄灭；R0>1时，曲线出现明显的波峰，且R0越大，波峰越陡峭、越高。
2. 最终感染规模：即疫情结束后，总康复者比例R(∞)/N。可以看到，即使R0只是从2.5增加到3.5，最终感染比例也会有显著提升，这解释了为什么将R0控制在1以下是如此重要。
3. 峰值感染比例：医疗系统承受的最大压力。R0的微小增加会导致峰值比例的急剧上升，这对病床、医护人员等资源规划是致命性的。
4. 峰值到达时间：R0越大，疫情发展越快，峰值到来得越早，留给社会准备和响应的时间就越短。

通过这个简单的参数分析，微分方程模型的威力就显现出来了。它不再是黑箱，而是一个可以让我们定量评估干预措施（如戴口罩、减少接触可以降低β，加快检测隔离可以缩短平均感染期即增大γ）效果的“数字沙盘”。

5. 从SIR到更复杂的模型：建模思维的延伸

经典的SIR模型是一个强大的起点，但现实往往更复杂。建模的魅力就在于，你可以根据实际情况，像搭积木一样对模型进行扩展。这里介绍几种常见的变体，并说明如何在MATLAB框架下修改方程和代码。

5.1 SEIR模型：加入潜伏期

很多传染病（如流感、COVID-19）在感染后不会立即具有传染性，会有一段潜伏期。为此，我们在S和I之间增加一个暴露者（Exposed, E）舱室。E个体已感染病毒但尚未表现出症状且没有传染性，以速率α转化为感染者I。

方程变为：

dS/dt = -β * (S/N) * I dE/dt = β * (S/N) * I - α * E dI/dt = α * E - γ * I dR/dt = γ * I

MATLAB修改：你需要将状态向量y扩展为[S; E; I; R]，并相应修改odefun。参数α是潜伏期的倒数（1/α为平均潜伏期）。

5.2 SIRS模型：免疫力会衰减

对于一些疾病（如普通感冒），康复后获得的免疫力不是终身的，会随时间衰减，使个体重新变为易感者。我们在SIR基础上增加一个从R流回S的通道，速率为δ（免疫衰减率）。

方程变为：

dS/dt = -β * (S/N) * I + δ * R dI/dt = β * (S/N) * I - γ * I dR/dt = γ * I - δ * R

这个模型可以描述疾病的周期性流行。

5.3 带出生死亡的SIR模型（常微分方程版）

如果考虑时间跨度很长，需要加入人口动力学。假设有常数出生率Λ（直接进入易感者）和自然死亡率μ（从每个舱室中以相同速率移除）。

方程变为：

dS/dt = Λ - β * (S/N) * I - μ * S dI/dt = β * (S/N) * I - γ * I - μ * I dR/dt = γ * I - μ * R dN/dt = Λ - μ * N # 总人口数不再恒定

此时总人口N也是一个变量，由微分方程描述。

实现心得： 当你扩展模型时，最好的方法是先画出示意图，标出所有舱室和它们之间的流向及速率。然后对每个舱室，遵循“流入 - 流出 = 净变化”的原则列写方程。在MATLAB中，只需按照新的方程更新odefun函数，并相应调整初始向量y0和参数列表即可。求解器ode45对此是完全通用的。

6. 模型校准与实战中的注意事项

在真正的数学建模竞赛或科研中，我们通常不是凭空设定β和γ，而是需要利用真实数据来校准（Calibrate）这些参数。这通常转化为一个优化问题：寻找一组参数，使得模型输出的曲线（如每日新增感染数γ*I(t)）与真实数据最吻合。

6.1 最小二乘法参数估计

假设我们有一组真实数据，记录了从第0天到第T天的累计感染数（或每日新增）。我们可以用MATLAB的优化工具箱（如lsqcurvefit,fminsearch）来求解。

基本思路如下：

1. 定义一个损失函数，比如预测值与真实值之差的平方和。
2. 将SIR模型的求解过程封装成一个函数，输入是参数[beta, gamma]，输出是模型预测的时间序列。
3. 调用优化算法，自动搜索使损失函数最小的参数值。

这里给出一个简化的代码框架：

% 假设 real_cases 是一个向量，包含从第1天到第T天的每日新增病例数 % time_data 是对应的天数向量 [1, 2, ..., T] % 初始猜测参数 params_guess = [0.2, 0.1]; % [beta_guess, gamma_guess] % 定义目标函数：计算模型预测与真实数据的误差 function error = objective_function(params, time_data, real_cases, N, I0) beta = params(1); gamma = params(2); % 运行SIR模型，时间跨度要覆盖数据期 [t_model, y_model] = ode45(@(t,y) sir_ode(t,y,beta,gamma,N), [0, max(time_data)], [N-I0, I0, 0]); % 模型预测的每日新增病例：gamma * I(t) % 需要将模型输出的时间序列插值到与真实数据相同的时间点上 I_model = interp1(t_model, y_model(:,2), time_data); % 线性插值获取对应天数的I(t) predicted_new_cases = gamma * I_model; % 注意：这是连续近似，更精确的做法是计算积分 % 计算误差（这里用简单的平方和） error = sum((predicted_new_cases - real_cases).^2); end % 使用fminsearch进行无约束优化（需要Optimization Toolbox） options = optimset('Display', 'iter', 'MaxFunEvals', 1000); params_opt = fminsearch(@(p) objective_function(p, time_data, real_cases, N, I0), params_guess, options); fprintf('优化得到的参数：beta = %.4f, gamma = %.4f, R0 = %.2f\n', params_opt(1), params_opt(2), params_opt(1)/params_opt(2));

注意：这是一个高度简化的示例。实际数据拟合中，你需要考虑数据的噪声、模型初始条件的敏感性、是否需要对新增病例进行平滑处理、以及选择更合适的误差函数（如泊松似然）等问题。

6.2 建模与编码中的常见“坑”与技巧

1. 量纲与参数范围：确保所有参数（β,γ,α等）的量纲一致（通常是1/时间）。β和γ的值通常在0到1之间，但具体取决于时间单位（天、周）。不合理（如极大）的参数值可能导致求解器失败或结果无意义。
2. 初始条件敏感性与“零除”错误：在SIR方程中，如果初始I(0)=0，那么dS/dt和dI/dt始终为0，疫情不会发生。另外，在计算S/N时，虽然数学上没问题，但在某些极端参数下，数值计算可能导致问题。确保初始I(0)为一个小的正数（如1）。
3. 求解器选择：ode45适用于大多数非刚性问题。如果你的模型方程包含差异巨大的时间尺度（即“刚性”问题，例如某些化学反应模型），ode45会变得非常慢，此时应换用刚性求解器如ode15s或ode23s。
4. 结果验证：
   • 守恒量检查：对于经典SIR，总人口S+I+R应恒定。在代码中加入total = S+I+R; plot(t, total)来验证，应该是一条水平线。这是检查方程实现和求解是否正确的最快方法。
   • 长期行为：运行足够长的时间（tspan设大一些），观察感染者I(t)是否最终趋于0，以及易感者S(t)是否趋于一个大于0的稳定值。
   • 与解析结论对照：例如，最终易感者比例s∞ = S(∞)/N应满足方程s∞ = exp(-R0*(1 - s∞))。你可以用数值解验证这个隐式方程是否近似成立。
5. 性能与向量化：如果需要进行大规模参数扫描或长时间模拟，考虑将代码向量化。例如，如果要同时模拟100组不同参数，可以尝试在odefun内使用矩阵运算，或者使用parfor循环进行并行计算。

微分方程建模是一个“思考->翻译->实现->验证->分析”的循环。MATLAB是一个绝佳的平台，它将数学表达和数值计算无缝连接。从SIR这个经典模型入手，理解其每一个环节，你就掌握了用微分方程刻画动态世界的基本方法论。下次当你遇到种群竞争、谣言传播、化学反应、热量传导等问题时，不妨想想：能否划分几个“状态舱室”？它们之间如何“流动”？答案，往往就是一个微分方程组。