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

1. 从一道例题开始为什么微分方程是数学建模的“灵魂”如果你参加过数学建模竞赛或者处理过任何涉及动态变化、趋势预测的工程问题那你一定绕不开“微分方程”这四个字。它听起来有点吓人像是数学系高年级的专属课程但实际上它可能是连接你的问题描述和计算机求解之间最直接、最有力的桥梁。很多人学了一堆解法背了一堆公式但一到实际建模还是无从下手——问题往往不在于数学本身而在于如何把现实世界那个模糊的“变化关系”翻译成微分方程这个精确的数学语言。今天我们不空谈理论就从一道非常经典的建模例题出发手把手带你走完“问题分析 - 建立方程 - MATLAB求解 - 结果分析”的全过程。这道题就是**传染病模型SIR模型**的一个简化版。它几乎出现在每一本建模教材里因为它完美地诠释了微分方程如何刻画一个系统中个体状态的动态转移。我们的目标很明确不仅给你能直接运行的MATLAB源码更要彻底讲清楚每一行代码背后的建模思想和数学原理让你下次遇到类似问题能自己“造轮子”而不是只会“调包”。你会发现掌握微分方程建模核心不是解方程计算机可以代劳而是定义方程。这就像编程语法不难难的是算法设计。我们通过SIR模型这个“麻雀”把它彻底解剖。关键词很明确数学建模、微分方程、MATLAB、源码。这篇文章就是为需要将理论付诸实践的理工科学生、工程师以及所有希望用数学工具量化分析动态过程的朋友准备的。2. SIR传染病模型如何将一场疫情“翻译”成数学方程在开始写代码之前我们必须把模型搞清楚。SIR模型是Kermack和McKendrick在1927年提出的至今仍是描述传染病动力学的基石。它的核心思想是将人群划分为三个“舱室”并通过微分方程描述其间的流动。2.1 模型假设与状态定义任何建模的第一步都是做出合理且必要的简化假设。对于经典的SIR模型我们假设总人口数N恒定。不考虑出生、死亡、迁移等人口变动。这是一个很强的假设但在疫情短期爆发期内比如几个月可以近似成立。均匀混合。人群中的个体充分接触任何一个易感者S遇到任何一个感染者I的概率是相同的。这忽略了社交网络、空间结构的影响。康复者获得永久免疫。一旦从感染中恢复将终身免疫且不再具有传染性。疾病传播和康复过程是随机的但在大群体中可以用平均速率来描述。基于以上假设我们定义三个核心状态变量它们都是关于时间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 准备工作定义方程函数odefunMATLAB求解常微分方程组主要使用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关键点解析t是时间变量虽然我们的方程不明显依赖于t称为自治系统但MATLAB的ODE求解器函数接口要求必须有前两个参数(t, y)。y是一个列向量包含了所有状态变量在某个时刻t的值。我们约定其顺序为[S; I; R]。我们将参数beta,gamma,N作为额外的输入参数传入而不是在函数内部写死这样使得函数更灵活便于后续进行参数研究。方程的实现就是直接翻译我们上一节推导的数学公式。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_time5, 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);代码逐段解读参数部分我们设定了一个1000人的封闭社区从1个感染者开始。beta0.3,gamma0.1计算得R03这是一个典型的高传染性疾病如麻疹、COVID-19原始毒株的估计值预示着疫情会爆发。求解器调用ode45是MATLAB中最常用的非刚性常微分方程求解器它采用Runge-Kutta方法对于大多数问题精度和效率都很好。(t,y) ...是匿名函数它将我们预设的参数beta, gamma, N“绑定”到sir_ode函数上以满足求解器对函数签名的要求只接受t, y两个输入。结果提取ode45返回两个数组t时间点和y每个时间点对应的状态向量。y是一个三列的矩阵我们按列提取出S, I, R。可视化我们创建了两个子图。左图展示S、I、R三类人群的绝对数量变化右图专门展示感染者比例I/N并标记出疫情峰值感染人数最多的时刻这对于评估医疗系统压力至关重要。关键指标输出在命令行窗口打印出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);分析解读 这段代码固定康复率gamma0.1即平均感染期10天通过改变感染率beta来获得不同的R00.5, 1.5, 2.5, 3.5。我们同时观察四个关键指标感染者比例曲线直观展示不同传播速度下的疫情发展轨迹。R00.51时疫情很快熄灭R01时曲线出现明显的波峰且R0越大波峰越陡峭、越高。最终感染规模即疫情结束后总康复者比例R(∞)/N。可以看到即使R0只是从2.5增加到3.5最终感染比例也会有显著提升这解释了为什么将R0控制在1以下是如此重要。峰值感染比例医疗系统承受的最大压力。R0的微小增加会导致峰值比例的急剧上升这对病床、医护人员等资源规划是致命性的。峰值到达时间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 γ * IMATLAB修改你需要将状态向量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来求解。基本思路如下定义一个损失函数比如预测值与真实值之差的平方和。将SIR模型的求解过程封装成一个函数输入是参数[beta, gamma]输出是模型预测的时间序列。调用优化算法自动搜索使损失函数最小的参数值。这里给出一个简化的代码框架% 假设 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/时间。β和γ的值通常在0到1之间但具体取决于时间单位天、周。不合理如极大的参数值可能导致求解器失败或结果无意义。初始条件敏感性与“零除”错误在SIR方程中如果初始I(0)0那么dS/dt和dI/dt始终为0疫情不会发生。另外在计算S/N时虽然数学上没问题但在某些极端参数下数值计算可能导致问题。确保初始I(0)为一个小的正数如1。求解器选择ode45适用于大多数非刚性问题。如果你的模型方程包含差异巨大的时间尺度即“刚性”问题例如某些化学反应模型ode45会变得非常慢此时应换用刚性求解器如ode15s或ode23s。结果验证守恒量检查对于经典SIR总人口SIR应恒定。在代码中加入total SIR; plot(t, total)来验证应该是一条水平线。这是检查方程实现和求解是否正确的最快方法。长期行为运行足够长的时间tspan设大一些观察感染者I(t)是否最终趋于0以及易感者S(t)是否趋于一个大于0的稳定值。与解析结论对照例如最终易感者比例s∞ S(∞)/N应满足方程s∞ exp(-R0*(1 - s∞))。你可以用数值解验证这个隐式方程是否近似成立。性能与向量化如果需要进行大规模参数扫描或长时间模拟考虑将代码向量化。例如如果要同时模拟100组不同参数可以尝试在odefun内使用矩阵运算或者使用parfor循环进行并行计算。微分方程建模是一个“思考-翻译-实现-验证-分析”的循环。MATLAB是一个绝佳的平台它将数学表达和数值计算无缝连接。从SIR这个经典模型入手理解其每一个环节你就掌握了用微分方程刻画动态世界的基本方法论。下次当你遇到种群竞争、谣言传播、化学反应、热量传导等问题时不妨想想能否划分几个“状态舱室”它们之间如何“流动”答案往往就是一个微分方程组。