简介：这份资源面向生物数学、传染病动力学方向的学习者与科研人员，围绕Matlab在传染病建模中的应用展开，帮助读者理解SIR、SEIR等模型的构建思路，并借助代码完成无病平衡点与地方病平衡点的稳定性分析。压缩包共13个文件，包含8个m脚本、4个fig图形文件和1份pdf文档，整体约2.48MB；m文件用于定义模型参数、方程与求解流程，fig文件保存了不同参数下的动态仿真结果，pdf则提供理论参考。资源中涉及后向分支、E0至E2各平衡点的局部与全局稳定性判定，以及基于dde23求解器的延迟微分方程模拟，可帮助读者对照代码复现疾病传播规律、验证基本再生数R0对长期行为的影响。目前已有3887人学习下载，适合希望将数值计算与传染病防控分析结合的中高级读者参考。

1. 用 Matlab 搭传染病模型：从 SIR 曲线到参数反推，一次讲透

很多人第一次在 Matlab 里跑传染病模型，都是被一条“S 型曲线”骗进来的——以为把三个微分方程丢进ode45就完事了，结果参数一改，曲线要么直接躺平，要么炸成一条直线，完全对不上现实里的感染波峰。传染病模型在 Matlab 里真正难的不是写方程，而是把“感染率、恢复率、初始条件”这几个旋钮拧到能复现真实传播节奏的位置，并且知道每个旋钮拧过头会发生什么。这篇笔记面向两类人：一类是刚接触 SIR/SEIR、想用 Matlab 把模型跑起来的学生和初级工程师；另一类是做数据分析、需要快速验证一套传播假设是否站得住脚的从业者。下面从方程怎么落地、参数怎么标定、坑在哪，一路写到怎么用模型做情景推演，全部给可复现的代码和参数说明。

2. SIR 与 SEIR 的方程拆解：先搞清楚每个变量代表谁

2.1 三个仓室和四个仓室的区别到底在哪

SIR 把人群分成三类：易感者 S、感染者 I、移除者 R。它的核心假设是“感染后经过一个平均周期就永久移除”，适合病程短、康复后免疫明显的场景。SEIR 多了一个暴露者 E，也就是已经感染但还没传染性的阶段，适合有潜伏期的传播过程。选哪个不是看哪个高级，而是看你的数据里有没有“潜伏期”这个可观测的滞后。如果感染曲线相对移除曲线有明显的时间延迟，SIR 会把这个延迟硬塞进恢复率里，导致参数失真，这时候就该上 SEIR。

在 Matlab 里，SIR 是一个三维状态向量，SEIR 是四维。维度一变，后面所有参数标定、初值设置、画图逻辑都要跟着改。我一般建议先把 SIR 跑通，确认自己对ode45的返回结构、时间向量、参数传递方式都熟了，再扩到 SEIR。直接上 SEIR 很容易在“E 的初值给多少”这个问题上卡住，因为 E 通常不可直接观测。

2.2 把微分方程写成 Matlab 能吃的函数句柄

SIR 的标准形式是：dS/dt = -βSI/N，dI/dt = βSI/N - γI，dR/dt = γI。其中 β 是有效接触率，γ 是恢复率，N 是总人口。下面是一个可以直接保存为sir_model.m的函数：

function dydt = sir_model(t, y, beta, gamma, N) % y(1)=S, y(2)=I, y(3)=R S = y(1); I = y(2); R = y(3); dS = -beta * S * I / N; dI = beta * S * I / N - gamma * I; dR = gamma * I; dydt = [dS; dI; dR]; end

逻辑说明：这里把beta、gamma、N作为额外参数传入，而不是写死在函数里，是为了后面做参数扫描时不用反复改文件。N放在分母上是把“接触人数”归一化成比例，避免人口规模不同导致 β 不可比。注意y是列向量，返回的dydt也必须是列向量，否则ode45会报维度不一致。

参数说明：beta的单位是“每人每天有效接触次数”，典型范围在 0.1 到 1.0 之间；gamma是恢复率的倒数对应平均感染期，比如平均 7 天恢复，γ≈0.143。N如果设成 100 万，初值 I 给 1，那早期增长会非常慢，因为 S/N 接近 1 但 I 太小，这是正常的指数增长起步阶段。

2.3 用 ode45 求解并画出三条曲线

N = 1e6; beta = 0.35; gamma = 0.14; I0 = 100; S0 = N - I0; R0 = 0; y0 = [S0; I0; R0]; tspan = [0 200]; [t, y] = ode45(@(t,y) sir_model(t, y, beta, gamma, N), tspan, y0); S = y(:,1); I = y(:,2); R = y(:,3); plot(t, S, 'b', t, I, 'r', t, R, 'g', 'LineWidth', 1.5); legend('S','I','R'); xlabel('天数'); ylabel('人数'); title('SIR 模型基本再生数 R0 = \beta/\gamma'); grid on;

逻辑说明：ode45是变步长求解器，适合这种非刚性方程。匿名函数@(t,y) sir_model(...)把额外参数包进去，是 Matlab 里传参的标准做法。画图时直接对y的列做索引，S 是第一列，I 是第二列，R 是第三列，顺序不能乱。

参数说明：tspan给 200 天是为了看到完整的感染波峰和回落。如果beta/gamma小于 1，感染曲线不会出现明显波峰，I 会单调下降，这时候要么调大 β，要么检查初值 I0 是不是太小导致还没涨起来就结束了。LineWidth设 1.5 是为了在论文或报告里截图清晰，不是必须的。

3. 参数标定与情景推演：让模型对得上数据

3.1 基本再生数 R0 怎么算、怎么用

R0 = β/γ，这是判断传播是否扩散的核心指标。R0 > 1 时感染会增长，R0 < 1 时感染会自然消退。在 Matlab 里可以直接算：

R0 = beta / gamma; fprintf('基本再生数 R0 = %.2f\n', R0);

逻辑说明：这个值不是模型跑出来的，而是参数组合出来的。它的意义在于给你一个快速判断：如果 R0 只有 1.2，那感染波峰不会太高，但拖得久；如果 R0 到 3 以上，波峰又高又陡，对医疗资源冲击大。

参数说明：β 和 γ 的估计通常来自数据拟合。如果手头有每日新增病例，可以用lsqcurvefit或fminsearch去最小化模型输出和实际数据的残差。我一般先用fminsearch粗调，再用lsqcurvefit精修，因为前者不需要工具箱，后者对初值敏感但收敛快。

3.2 用 fminsearch 做参数反推的最小示例

假设你有一组观测到的 I 数据，存在I_obs里，时间点存在t_obs里：

function err = fit_err(params, t_obs, I_obs, N) beta = params(1); gamma = params(2); I0 = I_obs(1); S0 = N - I0; y0 = [S0; I0; 0]; [~, y] = ode45(@(t,y) sir_model(t, y, beta, gamma, N), t_obs, y0); I_model = y(:,2); err = sum((I_model - I_obs).^2); end % 调用 params0 = [0.3, 0.1]; params_hat = fminsearch(@(p) fit_err(p, t_obs, I_obs, N), params0);

逻辑说明：fit_err返回的是残差平方和，fminsearch会不断调整params让这个值变小。注意每次迭代都要重新跑ode45，所以如果数据点很多、时间跨度很长，这个过程会比较慢。可以把t_obs限制在感染上升期和波峰附近，减少计算量。

参数说明：params0是初值猜测，给得太离谱会导致收敛到局部极小。我一般会先用网格搜索扫一遍 β 和 γ 的粗略范围，比如 β 从 0.1 到 0.6，γ 从 0.05 到 0.3，找到残差最小的格子，再把这个格子中心作为fminsearch的初值。N在整个拟合过程中保持不变，不要把它也当成待拟合参数，否则会出现参数不可辨识的问题。

3.3 情景推演：改一个参数看波峰怎么变

beta_list = [0.2, 0.35, 0.5]; gamma = 0.14; figure; hold on; for i = 1:length(beta_list) beta = beta_list(i); [t, y] = ode45(@(t,y) sir_model(t, y, beta, gamma, N), [0 200], y0); plot(t, y(:,2), 'LineWidth', 1.5, 'DisplayName', sprintf('\\beta=%.2f', beta)); end legend('show'); xlabel('天数'); ylabel('感染者人数'); title('不同 \beta 下的感染曲线对比'); grid on;

逻辑说明：循环里每次改 β 重新求解，把 I 曲线叠在一起。DisplayName配合legend('show')可以自动生成图例，不用手动写字符串数组。

参数说明：β 从 0.2 到 0.5，R0 从 1.43 到 3.57，波峰高度和出现时间都会明显变化。如果 β 低于 γ，也就是 R0 < 1，曲线会直接下降，这时候图上几乎看不到波峰，不要以为是代码写错了。

4. 避坑与排查：这五个地方最容易翻车

4.1 初值 I0 给太小导致曲线“没反应”

现象：跑出来的 I 曲线几乎贴着零轴，看起来像模型没生效。原因：I0 相对于 N 太小，比如 N=1e6 时 I0=1，早期增长是指数级的，但在线性坐标下前几十天都看不出变化。解决：要么把 I0 调大到 100 或 1000，要么把 y 轴改成对数坐标semilogy，这样早期增长就能看清。

4.2 把 beta 和 gamma 的单位搞混

现象：拟合出来的 γ 大于 1，或者 β 小得离谱。原因：β 的单位是“每人每天”，γ 的单位是“每天”，如果数据的时间单位是周而不是天，γ 会差 7 倍。解决：统一时间单位，所有参数都按天算。如果原始数据是周报，先把时间轴乘以 7 转成天，再拟合。

4.3 ode45 报“维度不一致”或“输入参数不足”

现象：运行时报错，提示Not enough input arguments或Dimensions of arrays being concatenated are not consistent。原因：sir_model的签名和匿名函数里的参数数量对不上，或者y0给成了行向量而dydt返回列向量。解决：检查sir_model定义里的参数个数，确保匿名函数@(t,y) sir_model(t, y, beta, gamma, N)传了 5 个参数；y0用分号分隔成列向量。

4.4 拟合时参数跑到负值

现象：fminsearch返回的 β 或 γ 是负数，物理上无意义。原因：fminsearch是无约束优化，参数可以取任意实数。解决：在fit_err里加惩罚项，比如if beta <= 0 || gamma <= 0, err = 1e10; return; end，把负参数直接判为无效。或者改用lsqcurvefit并设置上下界lb和ub。

4.5 把 R0 当成固定常数到处用

现象：不同情景下都用同一个 R0 去解释，结论互相矛盾。原因：R0 是 β 和 γ 的组合，β 会随接触行为变化，γ 会随医疗条件变化，R0 本身不是常数。解决：每次情景推演都重新算 R0，并且在报告里写清楚“在当前 β 和 γ 假设下 R0 等于多少”，而不是笼统说“R0 是 2.5”。

5. 进阶技巧：用 SEIR 处理潜伏期，用事件函数处理干预

5.1 SEIR 的方程和初值设置

SEIR 在 SIR 基础上加了 E 仓室：dS/dt = -βSI/N，dE/dt = βSI/N - σE，dI/dt = σE - γI，dR/dt = γI。σ 是潜伏期转为感染期的速率，1/σ 是平均潜伏期。代码结构：

function dydt = seir_model(t, y, beta, sigma, gamma, N) S = y(1); E = y(2); I = y(3); R = y(4); dS = -beta * S * I / N; dE = beta * S * I / N - sigma * E; dI = sigma * E - gamma * I; dR = gamma * I; dydt = [dS; dE; dI; dR]; end

逻辑说明：E 的初值通常给 0 或者给一个很小的数，因为暴露者往往不可观测。如果 E0=0 且 I0>0，早期 E 会从零开始上升，感染曲线会比 SIR 多一个滞后。

参数说明：σ 的典型值对应潜伏期 3 到 7 天，即 σ 在 0.14 到 0.33 之间。β 和 γ 的含义与 SIR 一致。注意 SEIR 的 R0 仍然是 β/γ，σ 不影响 R0，但影响波峰出现的时间。

5.2 用事件函数模拟“封控”或“干预”的突然生效

Matlab 的ode45支持事件函数，可以在某个条件触发时改变参数或停止求解。比如模拟第 30 天开始干预，β 从 0.35 降到 0.15：

function [value, isterminal, direction] = intervention_event(t, y) value = t - 30; % 第30天触发 isterminal = 1; % 停止求解 direction = 1; % 只在上穿时触发 end % 第一段：0到30天 [t1, y1] = ode45(@(t,y) sir_model(t, y, 0.35, 0.14, N), [0 30], y0, ... odeset('Events', @intervention_event)); % 第二段：30天之后，用第一段的末状态作为初值 y0_new = y1(end, :); [t2, y2] = ode45(@(t,y) sir_model(t, y, 0.15, 0.14, N), [30 200], y0_new); % 合并 t_all = [t1; t2(2:end)]; y_all = [y1; y2(2:end,:)];

逻辑说明：事件函数返回的value在 t=30 时过零，isterminal=1让求解器停下。然后手动把第一段的末状态取出来，作为第二段的初值，换一个 β 继续跑。合并时去掉t2的第一个点，避免时间重复。

参数说明：direction=1表示只在value从负变正时触发，如果写成 0 则双向触发。odeset('Events', ...)是固定写法。第二段的 β 从 0.35 降到 0.15，R0 从 2.5 降到 1.07，感染曲线会在干预后继续上升一小段然后回落，这个“惯性”是 SEIR 和 SIR 都会有的，不是 bug。

5.3 用热图看 β 和 γ 同时变化时的波峰高度

beta_grid = 0.1:0.05:0.6; gamma_grid = 0.05:0.02:0.3; peak_I = zeros(length(beta_grid), length(gamma_grid)); for i = 1:length(beta_grid) for j = 1:length(gamma_grid) [~, y] = ode45(@(t,y) sir_model(t, y, beta_grid(i), gamma_grid(j), N), [0 300], y0); peak_I(i,j) = max(y(:,2)); end end imagesc(gamma_grid, beta_grid, peak_I); colorbar; xlabel('\gamma'); ylabel('\beta'); title('不同 \beta 和 \gamma 下的感染波峰高度');

逻辑说明：双重循环遍历参数网格，每次求解后取 I 的最大值存进矩阵。imagesc用颜色表示波峰高度，横轴 γ，纵轴 β。颜色越亮表示波峰越高。

参数说明：网格步长决定分辨率，0.05 和 0.02 是折中，再细会明显变慢。tspan给 300 天是为了确保波峰已经出现，如果 R0 接近 1，波峰可能很晚才来，需要再拉长时间。这个热图能帮你快速找到“波峰高度可接受”的参数区间，而不是一个个试。

5.4 我自己的习惯：先固定 N 和 I0，再调 β 和 γ

跑了这么多轮，我最大的教训是：不要同时改所有参数。N 和 I0 一旦确定就不要再动，否则每次改 β 或 γ 都要重新理解曲线为什么变。我一般先把 N 设成实际人口规模，I0 设成 100 左右，然后固定 γ 在 0.1 到 0.2 之间，只调 β 看波峰怎么变。等 β 的范围摸清了，再回头微调 γ。这样每一步只有一个变量在动，出了问题能立刻定位。希望帮到你。

本文还有配套的精品资源，点击获取