简介这份资源面向生物数学、传染病动力学方向的学习者与科研人员提供一套基于Matlab的传染病模型代码与配套分析材料帮助理解SIR、SEIR等模型的构建思路并完成无病平衡点与地方病平衡点的稳定性分析。压缩包共13个文件包含8个m脚本、4个fig图形文件和1份pdf文档整体约2.48MBm文件主要用于定义模型参数、求解微分方程及稳定性计算fig文件保存了不同平衡点情形下的仿真图形pdf则提供理论参考。资源中涉及后向分支、E0/E1/E2平衡点局部与全局稳定性等典型分析场景并借助dde23求解器处理延迟微分方程可帮助读者对照代码复现疾病传播动态、验证基本再生数R0对长期行为的影响。目前已有3887人学习下载适合作为生物数学建模入门与课程实践的参考素材。1. 传染病模型从纸面到代码为什么我劝你先跑通 SEIR 再谈参数拟合很多人第一次接触传染病动力学是在某次公共卫生事件的新闻里看到“基本传染数”这个词然后兴冲冲打开 Matlab想用一组真实数据反推参数。结果往往是曲线怎么调都对不上参数一改结果就飞最后怀疑模型本身有问题。其实问题不在模型而在于跳过了“先跑通标准模型、再谈拟合”这一步。Matlab 传染病模型这条路径核心不是某个神奇函数而是把 SIR、SEIR 这类常微分方程组用 ode45 稳定地解出来再逐步加入潜伏期、干预项和年龄分层。它适合做课程设计的学生、做应急推演的工程师以及需要快速验证防控策略的研究者。这一章先把“模型到底在算什么”讲清楚后面几章再落到代码、参数和踩坑上。传染病模型本质上是把人群按状态分箱再写每个箱子的人数随时间怎么变。SIR 只有易感、感染、移除三个箱子SEIR 多了一个潜伏期箱子。箱子之间的流动速率就是参数接触率、潜伏期倒数、恢复率。Matlab 的优势在于 ode45 这类自适应步长求解器对刚性不强的传染病方程足够稳而且画图、调参、批量扫描都在一个环境里完成。但很多人一上来就把真实确诊数塞进去做最小二乘忽略了报告延迟和检测能力变化导致拟合出的参数没有物理意义。我的建议是先用一组人为设定的参数生成“干净”曲线确认求解器和初值没问题再引入真实数据。这一步花不了半小时却能省掉后面几天的玄学调参。2. 用 ode45 跑通 SEIR从方程到可复现曲线2.1 状态方程怎么写才不容易出错SEIR 的标准形式是四个状态S、E、I、R。写成微分方程时最容易翻车的地方是符号和单位。我一般把总人数 N 固定方程里所有项都除以 N这样接触率 β 的量纲就是“每人每天有效接触次数”而不是“总接触次数”。下面这段代码把方程写成一个独立函数输入是时间 t、状态向量 y 和参数结构体 p。function dydt seir_eq(t, y, p) % y [S; E; I; R]均为人数 S y(1); E y(2); I y(3); R y(4); N S E I R; % 总人数动态计算防止初值不守恒 beta p.beta; % 接触率 sigma p.sigma; % 潜伏期倒数 gamma p.gamma; % 恢复率 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逻辑说明dS 是易感者被感染的速率用 S*I/N 表示随机混合下的有效接触dE 是新增潜伏减去转为感染dI 是潜伏转出减去恢复dR 只增不减。参数说明beta 越大传播越快sigma 决定潜伏期长短gamma 决定感染期长短。注意 N 用动态求和而不是外部传入这样即使初值有微小舍入误差也不会出现总人数漂移。2.2 主脚本初值、时间跨度与求解器选项有了方程函数主脚本负责设初值、调 ode45、画图。初值我习惯用“1 个感染者 若干易感者”的设定总人数取 10000便于观察比例。时间跨度先给 160 天因为大多数急性呼吸道传染病的流行周期在这个量级。% 主脚本run_seir.m clear; clc; p.beta 0.35; % 接触率 p.sigma 1/5.2; % 潜伏期约 5.2 天 p.gamma 1/7.0; % 感染期约 7 天 N0 10000; y0 [N0-1; 0; 1; 0]; % 初始1 个感染者其余易感 tspan [0 160]; opts odeset(RelTol,1e-6,AbsTol,1e-8); [t, y] ode45((t,y) seir_eq(t,y,p), tspan, y0, opts); % 画图 figure; plot(t, y(:,1), LineWidth, 1.5); hold on; plot(t, y(:,2), LineWidth, 1.5); plot(t, y(:,3), LineWidth, 1.5); plot(t, y(:,4), LineWidth, 1.5); xlabel(时间天); ylabel(人数); legend(S,E,I,R,Location,best); grid on;逻辑说明odeset 里把相对误差和绝对误差收紧到 1e-6 和 1e-8是因为传染病曲线在峰值附近变化快默认精度有时会让峰值位置偏移一两天。参数说明beta0.35 对应基本传染数 R0 beta/gamma ≈ 2.45属于中等偏高的呼吸道传播水平sigma 和 gamma 分别控制潜伏和恢复节奏。跑完这张图你应该能看到 I 曲线先升后降R 单调上升S 单调下降四条线加起来始终等于 10000。2.3 基本传染数 R0 的两种算法与验证R0 是传染病模型里最常被问到的量。对于 SEIRR0 beta / gamma注意它和 sigma 无关因为潜伏期不改变最终感染比例只改变时间尺度。但很多人会误以为 R0 要乘上 sigma这是常见误解。验证方法很简单把 beta 和 gamma 代入算出 R0然后看 I 曲线是否在 R01 时出现明显峰值R01 时单调衰减。下面这段代码做批量扫描把不同 beta 下的峰值画出来。beta_list 0.1:0.05:0.6; peak_I zeros(size(beta_list)); for k 1:length(beta_list) p.beta beta_list(k); [~, y] ode45((t,y) seir_eq(t,y,p), [0 300], y0, opts); peak_I(k) max(y(:,3)); end plot(beta_list, peak_I, o-); xlabel(beta); ylabel(感染峰值); grid on;逻辑说明固定 sigma 和 gamma只改 beta观察峰值变化。参数说明beta 从 0.1 到 0.6对应 R0 从 0.7 到 4.2。你会看到 R0 略大于 1 时峰值很小R0 超过 2 后峰值迅速抬升。这个扫描能帮你建立“参数改多少、结果变多少”的直觉后面拟合真实数据时就不会盲目。3. 参数拟合与干预项把模型拉到真实数据附近3.1 最小二乘拟合的 Matlab 实现与初值陷阱有了正向模型下一步是用真实或半真实的累计感染数据反推 beta 和 gamma。Matlab 里常用 lsqcurvefit 或 fminsearch。我一般先用 fminsearch 做粗调再用 lsqcurvefit 精修。关键陷阱是初值如果 beta 初值给得离真值太远fminsearch 容易陷到局部极小拟合出的曲线“看起来还行”但参数没有物理意义。% 假设 obs_t 和 obs_I 是观测时间与感染人数 obj (p) sum((simulate_I(p, obs_t) - obs_I).^2); p0 [0.3, 1/5, 1/7]; % [beta, sigma, gamma] 初值 p_est fminsearch(obj, p0); function I_sim simulate_I(p, t_query) y0 [10000-1; 0; 1; 0]; [t, y] ode45((t,y) seir_eq(t,y,struct(beta,p(1),sigma,p(2),gamma,p(3))), ... [0 max(t_query)], y0); I_sim interp1(t, y(:,3), t_query, pchip); end逻辑说明obj 是残差平方和simulate_I 把 ode45 结果插值到观测时间点。参数说明p0 的三个初值分别对应中等传播、5 天潜伏、7 天恢复。注意 interp1 用 pchip 而不是 linear因为感染曲线在峰值附近曲率大线性插值会引入偏差。如果拟合结果里 sigma 跑到 0.5 以上潜伏期不到 2 天基本可以判定过拟合或数据质量问题。3.2 加入干预项时变 beta 的两种写法真实疫情里防控措施会改变接触率。最常见的做法是把 beta 写成时间的阶跃函数或指数衰减函数。阶跃函数适合“某天开始封控”的场景指数衰减适合“逐渐减少接触”的场景。下面给出阶跃版本的代码。function dydt seir_interv(t, y, p) S y(1); E y(2); I y(3); R y(4); N S E I R; if t p.t_lock beta p.beta0; else beta p.beta0 * p.lock_ratio; % 例如降到 30% end dS -beta * S * I / N; dE beta * S * I / N - p.sigma * E; dI p.sigma * E - p.gamma * I; dR p.gamma * I; dydt [dS; dE; dI; dR]; end逻辑说明t_lock 是干预开始时间lock_ratio 是接触率保留比例。参数说明lock_ratio0.3 表示接触减少 70%。注意阶跃函数在 t_lock 处不连续ode45 会自动缩小步长但如果你把 RelTol 设得太松可能会跳过阶跃点导致曲线出现不自然的折角。我的习惯是在 t_lock 前后各加一个输出点或者直接用 odeset 的 MaxStep 限制最大步长。3.3 拟合结果怎么判断“能信”拟合完不要只看残差平方和。我一般看三个东西一是参数是否落在合理区间beta 对应 R0 在 1 到 6 之间sigma 对应潜伏期 2 到 10 天gamma 对应感染期 3 到 14 天二是峰值时间和峰值大小是否和观测一致三是把拟合参数代回模型看 S 和 R 曲线是否出现违反常识的形态比如 R 超过总人数。如果三项都过才认为这组参数可以拿去做情景推演。4. 避坑与排查传染病模型里最容易翻车的五件事4.1 曲线爆炸或出现负人数现象ode45 跑出来的 I 曲线冲到 10 的几十次方或者 S 变成负数。原因beta 给得过大或者时间跨度太长数值误差累积。解决先把 beta 降到 0.1 量级试跑确认求解器正常如果必须用大 beta把 RelTol 收紧到 1e-8并检查方程里有没有漏掉除以 N。4.2 总人数不守恒现象SEIR 随时间缓慢漂移。原因方程里用了外部传入的固定 N而初值之和与 N 有舍入差异。解决像 2.1 那样在方程内部动态求和或者把初值归一化到总和为 1用比例而不是人数做状态变量。4.3 拟合参数没有物理意义现象拟合出的 gamma 大于 1对应感染期不到 1 天。原因观测数据里混入了报告延迟或者初值给得太离谱。解决先对观测数据做滑动平均再拟合同时给参数加边界比如 gamma 限制在 1/14 到 1/3 之间。4.4 干预项导致求解器卡死现象加入阶跃 beta 后ode45 运行时间暴涨或报错。原因不连续点让自适应步长反复试探。解决把阶跃改成平滑过渡比如用 1/(1exp(-k*(t-t_lock))) 做软切换k 取 0.5 到 2 之间。4.5 批量扫描时内存溢出现象参数扫描循环里每次保存完整时间序列几百组参数就把内存吃满。原因ode45 输出点数随参数变化累积起来很大。解决每次只保存峰值、峰值时间、最终感染比例这几个标量不要存完整 y 矩阵。5. 从单组参数到情景推演一个可复用的批量框架前面几章把单组参数的模型跑通了但实际做防控推演时你需要同时比较多种干预强度、多种开始时间。我一般写一个批量框架输入是参数网格输出是每个组合的峰值和总感染比例。下面这段代码把 beta0、lock_ratio、t_lock 三个维度组合起来跑完直接出表格。beta0_list [0.25 0.35 0.45]; ratio_list [0.2 0.4 0.6]; tlock_list [20 40 60]; results []; for b beta0_list for r ratio_list for tl tlock_list p.beta0 b; p.lock_ratio r; p.t_lock tl; p.sigma 1/5.2; p.gamma 1/7.0; y0 [9999; 0; 1; 0]; [t, y] ode45((t,y) seir_interv(t,y,p), [0 300], y0, ... odeset(RelTol,1e-6,MaxStep,1)); peak max(y(:,3)); total_inf y(end,4) y(end,3) y(end,2); results [results; b, r, tl, peak, total_inf]; end end end T array2table(results, VariableNames, ... {beta0,lock_ratio,t_lock,peak_I,total_infected}); disp(T);逻辑说明三层循环遍历所有组合每次跑完只记录峰值和最终感染规模。参数说明MaxStep 设为 1 天是为了在阶跃点附近不跳过干预开始时间。跑完你会得到一张 27 行的表按 total_infected 排序就能看出哪种组合最有效。我的血泪经验是不要一上来就做三维网格先把 beta0 固定只扫 lock_ratio 和 t_lock确认趋势符合直觉后再加维度。否则一旦结果反直觉你很难判断是模型问题还是代码问题。最后说一个我自己的习惯每次改完参数或方程先跑一组“基准场景”并保存曲线后面任何改动都和基准对比。这样即使中间某次调参把结果带偏也能快速定位是哪一步引入的。传染病模型不是越复杂越好能把 SEIR 加一个时变 beta 跑稳、拟合合理、推演可解释已经能覆盖大多数课程设计和应急推演的需求。希望帮到你。本文还有配套的精品资源点击获取