简介这份资源面向生物数学、传染病动力学方向的学习者与研究者提供一套基于Matlab的传染病模型代码与配套分析材料帮助理解SIR、SEIR等模型的构建思路并完成无病平衡点与地方病平衡点的稳定性分析。压缩包共13个文件约2.48MB以8个m脚本、4个fig图形文件和1个pdf文档为主m文件承担模型参数定义、微分方程求解与稳定性计算fig文件保存了不同平衡点情形下的动态曲线pdf则提供理论参考。内容涉及后向分支、基本再生数R0判断以及E0、E1、E2等平衡点的局部与全局稳定性讨论并借助dde23求解器处理延迟微分方程。已有3887人学习适合希望把理论推导落到可运行代码、并借助图形结果验证结论的读者参考。1. 从一张“SIR 曲线跑飞了”的截图说起群里有人丢出一张图Matlab 里 SIR 模型的感染人数曲线先冲顶再断崖评论区第一反应是“参数调得真好看”。但把代码拉下来一跑同样的 beta、gamma换一组初值曲线直接贴地或者爆表。问题不在模型在数值解法和参数边界没对齐。这份 Matlab 传染病模型资源核心就是几套能直接跑的 SIR、SEIR、SEIRS 脚本外加参数扫描和拟合的入口函数。它解决的不是“生物数学是什么”而是“拿到一组真实或模拟的感染数据怎么在 Matlab 里把模型跑通、把参数估出来、把图出成能放进报告的样子”。适合做课程设计、数模竞赛、公共卫生方向初探的从业者也适合已经会解 ODE 但没系统写过传染病模型代码的人。下面按“先跑通再调参再避坑”的顺序拆。2. SIR 与 SEIR 的 Matlab 实现从方程到 ode45 的映射2.1 状态变量与方程组的代码化SIR 三个状态S 易感、I 感染、R 移除。方程本身不复杂但写成 Matlab 函数时很多人会在向量索引上翻车。资源里常见的写法是把状态拼成一个列向量 y然后按位置取。下面这段是 SIR 右端函数的标准骨架我一般会把它单独存成sir_rhs.m方便 ode45 调用。function dydt sir_rhs(t, y, beta, gamma) % y(1)S, y(2)I, y(3)R % 总人数 N 在外部固定这里不重新计算 S y(1); I y(2); R y(3); dydt zeros(3,1); dydt(1) -beta * S * I; % 易感者被感染 dydt(2) beta * S * I - gamma * I; % 感染者转移除 dydt(3) gamma * I; % 移除者累积 end逻辑说明这里用的是质量作用项 betaSI没有除以 N。如果你的数据是比例形式SIR1那 beta 的含义会变后面拟合时初值也要跟着改。参数 beta 是有效接触率gamma 是移除率1/gamma 就是平均感染期。常见做法是把 N 归一化到 1这样数值更稳ode45 的容差也好设。2.2 ode45 调用与初值、时间向量的设置有了右端函数主脚本里就是设初值、时间跨度、调 ode45。资源里给的例子通常用[0, 160]天初值 I0 取 1e-3 到 1e-2 之间。下面这段可以直接抄。N 1; % 归一化总人数 I0 0.001; % 初始感染比例 S0 N - I0; R0 0; y0 [S0; I0; R0]; beta 0.4; % 有效接触率 gamma 0.15; % 移除率 tspan [0, 160]; [t, y] ode45((t,y) sir_rhs(t, y, beta, gamma), 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(比例);参数说明beta 和 gamma 的比值就是基本再生数 R0 beta/gamma。上面这组 beta0.4、gamma0.15R0 约 2.67曲线会明显冲顶。如果你把 gamma 改成 0.4R0 降到 1感染曲线就单调下降这是判断参数是否合理的快速检查。ode45 默认相对容差 1e-3对传染病模型通常够用但如果初值 I0 特别小比如 1e-6建议把odeset(RelTol,1e-6,AbsTol,1e-8)加上否则早期增长段会失真。2.3 SEIR 扩展多一个潜伏期状态怎么加SEIR 在 S 和 I 之间插一个 E暴露但未传染。方程多一项代码改动很小但索引容易乱。资源里 SEIR 的右端函数一般长这样function dydt seir_rhs(t, y, beta, sigma, gamma) % y(1)S, y(2)E, y(3)I, y(4)R S y(1); E y(2); I y(3); R y(4); dydt zeros(4,1); dydt(1) -beta * S * I; dydt(2) beta * S * I - sigma * E; dydt(3) sigma * E - gamma * I; dydt(4) gamma * I; endsigma 是潜伏期转为感染期的速率1/sigma 是平均潜伏期。调用方式和 SIR 一样只是 y0 变成四维初值 E0 通常设 0 或很小的值。注意SEIR 的 R0 表达式是 beta/gamma和 sigma 无关但潜伏期会拉长疫情到达峰值的时间。如果你用同一组 beta、gamma 跑 SIR 和 SEIR会发现 SEIR 的峰值更低、更靠后这是模型结构带来的差异不是代码写错了。3. 参数估计与拟合怎么从数据反推 beta 和 gamma3.1 目标函数与 lsqcurvefit 的接法拿到实际感染数据后最直接的做法是最小二乘拟合。Matlab 里lsqcurvefit比fminsearch更稳因为它能设参数上下界。资源里通常把 ode45 包一层输出 I 的模拟值再和观测数据做差。下面是一个可复用的拟合脚本骨架。function I_sim simulate_I(params, tdata, y0) beta params(1); gamma params(2); [~, y] ode45((t,y) sir_rhs(t, y, beta, gamma), tdata, y0); I_sim y(:,2); end % 主脚本 tdata (1:1:60); % 观测时间点 I_obs ...; % 实际感染比例列向量 y0 [0.999; 0.001; 0]; lb [0.01, 0.01]; % 参数下界 ub [2.0, 1.0]; % 参数上界 p0 [0.3, 0.1]; % 初值猜测 obj (p) simulate_I(p, tdata, y0) - I_obs; params_hat lsqcurvefit(obj, p0, [], [], lb, ub);逻辑说明simulate_I里每次迭代都重新解一次 ODE计算量不小但 60 个点、两个参数秒级能出结果。参数说明lb 和 ub 必须给否则 beta 可能被拟合成负数ODE 直接报错。p0 的选择有讲究一般先用 R0 的粗略估计反推 beta比如你估计 R0 在 2 到 3 之间gamma 取 0.1那 beta 初值就设 0.2 到 0.3。3.2 多起点与参数可辨识性单次 lsqcurvefit 容易掉进局部极小。资源里如果带了多起点循环别跳过。常见做法是在 lb 和 ub 之间随机撒 20 到 50 个初值每个跑一遍取残差最小的那组。下面这段是简化版。best_res inf; best_p []; for k 1:50 p0 lb (ub - lb) .* rand(1,2); [p_try, res] lsqcurvefit(obj, p0, [], [], lb, ub); if res best_res best_res res; best_p p_try; end end注意如果不同起点收敛到差别很大的 beta、gamma但残差差不多说明参数不可辨识。这在 SIR 里很常见因为 beta 和 gamma 以比值形式影响曲线。解决办法是固定其中一个或者用先验信息约束。我一般会先固定 gamma 为文献常见值比如 0.1 到 0.2只拟合 beta看残差是否可接受再决定要不要放开。3.3 拟合结果的置信区间与残差检查拟合完不能只看曲线重合。资源里如果有nlparci的调用那是用来算参数置信区间的。更实用的做法是画残差图plot(tdata, I_obs - simulate_I(best_p, tdata, y0))。如果残差有系统性的正负交替说明模型结构不对比如该用 SEIR 却用了 SIR。另外R0 的置信区间可以用 bootstrap对残差重采样生成多组伪数据每组重新拟合取 beta/gamma 的 2.5% 和 97.5% 分位数。这个流程资源里不一定有现成函数但值得自己补上。4. 避坑与排查那些让曲线跑飞的操作4.1 初值 I0 设成 0 导致 ode45 输出全零现象跑出来的 I 曲线从头到尾都是 0S 也不变。原因I0 设了 0右端函数里 betaSI 恒为 0系统没有动力学。解决I0 至少给 1e-6归一化模型里常用 1e-3 或 1e-4。如果数据里第一天感染数就是 0那说明观测窗口选晚了往前推几天再拟合。4.2 时间向量步长太大ode45 返回的点不够密现象plot 出来的曲线是折线峰值被削平。原因tspan 只给了起点和终点ode45 自适应步长虽然内部算得密但默认只返回部分点。解决用tspan linspace(0, 160, 500)强制输出 500 个点或者用[t, y] ode45(..., tspan)后自己插值。注意强制输出点不影响计算精度只影响返回密度。4.3 beta 或 gamma 为负ODE 报错或结果无意义现象lsqcurvefit 迭代过程中报 “函数值出现复数” 或 “积分失败”。原因参数越界变成负数S 或 I 出现负值下一步 sqrt 或 log 就崩。解决lb 必须设正数下界比如 1e-4。如果拟合需要允许参数接近 0用lb [1e-4, 1e-4]不要用 0。4.4 归一化与未归一化混用R0 算错现象同一组数据两次拟合得到的 R0 差一倍。原因一次用 SIR1一次用 SIRNbeta 的量纲变了。解决在脚本开头明确 N 是 1 还是实际人口数所有初值和参数都按同一套来。我一般统一用归一化最后报告 R0 时只报比值不报绝对 beta。4.5 数据里的报告延迟没考虑拟合残差偏大现象拟合曲线总是比实际数据晚几天到峰值。原因实际感染数据有报告延迟SIR 假设感染即报告。解决在目标函数里加一个时间平移参数 tau模拟值取simulate_I(p, tdata - tau, y0)把 tau 也放进拟合参数。这个改动小但能明显改善残差。5. 进阶技巧用参数扫描看 R0 的相图与峰值关系跑通单条曲线后最有价值的操作是参数扫描。固定 gamma让 beta 从 0.05 扫到 0.8步长 0.01每次解一次 SIR记录感染峰值和峰值时间。下面这段可以直接跑出图是峰值随 R0 变化的曲线。gamma 0.15; beta_vec 0.05:0.01:0.8; peak_I zeros(size(beta_vec)); peak_t zeros(size(beta_vec)); y0 [0.999; 0.001; 0]; for k 1:length(beta_vec) beta beta_vec(k); [t, y] ode45((t,y) sir_rhs(t, y, beta, gamma), [0, 300], y0); [peak_I(k), idx] max(y(:,2)); peak_t(k) t(idx); end R0 beta_vec / gamma; yyaxis left; plot(R0, peak_I, LineWidth, 1.5); ylabel(感染峰值); yyaxis right; plot(R0, peak_t, --, LineWidth, 1.5); ylabel(峰值时间); xlabel(R0); grid on;逻辑说明R0 小于 1 时峰值就是初值 I0曲线单调下降R0 刚过 1 时峰值时间很晚峰值也低R0 继续增大峰值升高、峰值时间提前。这个相图能帮你快速判断一组拟合参数是否落在合理区间。参数说明tspan 上限给 300 天是为了让 R0 接近 1 时曲线也能跑完否则峰值还没出现就截断了。另一个实用技巧是画 (S, I) 相平面。把 ode45 输出的 S 和 I 直接 plot能看到轨线从 (S0, I0) 出发穿过 I 的零增长线 S gamma/beta然后回到 I0。这条零增长线的位置就是 R0 的倒数比看时间序列更直观。我每次拟合完都会补这张图用来检查拟合参数下的轨线是否和观测数据的相轨迹一致。如果观测数据点落在相平面上的位置和模拟轨线偏离说明模型结构或参数有问题光看时间序列可能被掩盖。从那以后我每次跑传染病模型都强制先画相平面再画时间序列顺序反了就容易只盯着曲线形状调参忽略动力学约束。希望帮到你。本文还有配套的精品资源点击获取