简介这份文档围绕传染病动力学中的经典SIR模型展开面向数学建模学习者、公共卫生相关专业学生及参加建模竞赛的读者帮助理解疾病在封闭人群中传播的完整机理。内容从Kermack与McKendrick在1927年提出的动力学方法出发系统梳理模型假设、微分方程组构成、数值计算与相轨线分析并以SARS、天花、流感、肝炎、麻疹等传染病为背景给出MATLAB ode45求解示例及i(t)、s(t)的数值结果表。资源包共1个doc文件约360KB以文字推导与公式说明为主便于直接阅读和整理笔记。文档重点呈现了感染人数在t≈7达到峰值后衰减、易感者比例趋于稳定非零值等结论并讨论了模型的意义与局限性可作为理论分析起点帮助读者预测疫情峰值、评估干预措施效果。目前已有46人学习下载。1. 从一份 .doc 说起SIR 模型到底能算清什么翻到一份名为「传染病问题中的SIR模型.doc」的文档第一反应不是它有多老而是它把一件很多教程讲糊的事讲透了给定 λ 和 μ疫情到底会不会爆、峰值在第几天、最终多少人被感染。这三个问题才是公卫决策和建模作业里真正要交的答案。这份文档用 Kermack-McKendrick 1927 年的动力学框架把人群切成易感者 S、感染者 I、康复者 R 三类靠一组微分方程推演传播过程再用 MATLAB 数值求解和相轨线做理论验证。它适合两类人一类是要交传染病建模课程设计的学生需要能跑通的代码和能写进报告的推导另一类是做公共卫生数据分析的从业者想快速搭一个基准模型评估干预效果。文档里 λ1、μ0.3 那组算例i(t) 在 t≈7 达到峰值 0.3444s(t) 最终收敛到 0.0398这些数字不是摆设是后面所有分析的锚点。2. SIR 微分方程组怎么落地从假设到可运行代码2.1 三个假设决定了模型的适用边界这份文档开篇就把假设摆得很清楚这不是走过场而是决定模型能不能用的前提。第一条考察期内不考虑出生、死亡、流动总人口 N 恒定人群比例满足 s(t)i(t)r(t)1。第二条病人日接触率 λ 为常数日治愈率 μ 为常数平均传染期 1/μ传染期接触数 σλ/μ。第三条康复者获得持久免疫不再进入传染链。这三条假设直接框定了适用范围天花、流感、肝炎、麻疹这类治愈后免疫力强的病能用而像普通感冒这种能反复感染的就不行那得换 SIS 模型。我一般会跟人说选模型先看病种病种决定假设假设决定方程方程决定你能不能信算出来的峰值。文档里也点明了缺陷——有效接触率被当成不变的实际会随防控措施、人群行为变化所以结果是趋势参考不是精确预言。2.2 微分方程组与参数含义模型的核心是三个方程描述三类人群随时间的变化率ds/dt -λ·s·i di/dt λ·s·i - μ·i dr/dt μ·i第一个方程说易感者减少的速度取决于当前易感者和感染者的乘积再乘接触率第二个方程说感染者增加来自新感染、减少来自治愈第三个方程说康复者只增不减。这里 λ 的单位是「每人每天有效接触人数」μ 是「每天被治愈的比例」两者相除得到 σ也就是一个病人在整个传染期内平均能传染多少人这个数在文档里叫交换数是判断疫情是否蔓延的关键阈值。提示σ 的倒数 1/σ 就是阈值。初始易感比例 s(0) 大于 1/σ疫情才会扩散小于等于感染者比例单调下降不会形成高峰。2.3 MATLAB 数值求解ode45 的完整调用文档给的代码不长但每一行都有讲究。我把它整理成可直接粘贴运行的版本并补上注释function y ill(t, x) a 1; % λ 日接触率 b 0.3; % μ 日治愈率 y [a*x(1)*x(2) - b*x(1); % di/dt -a*x(1)*x(2)]; % ds/dt end ts 0:50; % 时间区间 0 到 50 天步长 1 x0 [0.02, 0.98]; % 初值 [i(0), s(0)] [t, x] ode45(ill, ts, x0); plot(t, x(:,1), t, x(:,2)) % 画 i(t) 和 s(t) pause plot(x(:,2), x(:,1)) % 画相轨线 i~s逻辑说明ode45 是四阶五阶 Runge-Kutta 法适合这种非刚性常微分方程组。函数 ill 返回列向量第一个分量是 di/dt第二个是 ds/dt注意这里没有单独解 r因为 r1-s-i 可以直接算出来。参数说明a 对应 λb 对应 μ改这两个值就能模拟不同防控力度x0 的顺序是 [i0, s0]文档里 i(0)0.02、s(0)0.98符合疫情初期感染者少、易感者接近满员的场景。ts 的步长影响输出密度不影响 ode45 内部自适应步长画图够用。跑完这段你会得到文档表 1 那组数据t7 时 i 达到 0.3444 的峰值之后回落s 从 0.98 一路降到 0.0398 附近稳住。这个 0.0398 就是最终未被感染的比例也是群体免疫讨论的起点。2.4 相轨线分析为什么峰值必然出现又必然消失文档用相平面把 i 和 s 的关系画出来消去时间 t 后得到一条轨线。这条轨线的意义在于它不依赖具体时间刻度只看状态组合。轨线从 P0(0.98, 0.02) 出发先向右上方走到 s1/σ 那条竖线此时 i 最大然后折向左下方最终落在 s 轴上。这里有个反直觉的结论不管初始条件怎么变只要 s(0)1/σi(t) 一定先升后降最终趋于零而 s(t) 一定停在某个正数上不会归零。也就是说传染病会自己停下来但代价是一部分人始终没被感染——这部分人不是靠免疫而是靠传播链断了。文档里用极限证明和图形两条路都走了一遍我建议做报告时把相轨线和数值曲线放一起理论验证的完整度就上来了。3. 阈值、群体免疫与干预参数怎么调3.1 σ 阈值判断疫情是否蔓延的第一道关文档反复强调 1/σ 这个阈值。当 s(0)1/σ即 σ1/s(0) 时感染者比例才会增长。拿文档的算例σλ/μ1/0.3≈3.331/σ0.3而 s(0)0.98 远大于 0.3所以疫情必然扩散。反过来如果通过提高卫生水平把 λ 降到 0.2μ 保持 0.3σ≈0.671/σ≈1.5此时 s(0)0.98 仍大于 1.5不对0.981.5所以不会蔓延。这个计算很多人第一次会绕进去记住比较的是 s(0) 和 1/σ不是 σ 和 1。我一般会做一个参数扫描表把不同 λ、μ 组合下的峰值和最终感染比例列出来方便选干预方案λμσ1/σ峰值 i最终 s1.00.33.330.300.3440.03980.80.32.670.3750.280.0750.50.31.670.600.160.220.30.31.001.00无峰0.98这张表是跑 ode45 换参数得到的不是文档原表但方法完全一致。可以看到 λ 从 1 降到 0.3峰值消失最终感染比例从 96% 降到 2%。这就是社交距离、戴口罩这类降低接触率措施的理论依据。3.2 群体免疫门槛的计算与局限文档第五部分给出群体免疫条件初始免疫比例 r(0) 满足 r(0) ≥ 1 - 1/σ。以天花 σ5 为例需要 1-1/580% 的人免疫。这个公式的推导来自把 s(0) 替换成 1-r(0)再代入不蔓延条件。但文档也点了一个血泪经验免疫者必须均匀分布在人群中实际很难做到。印度等国即使投入大量资源天花直到 1977 年才根除。做建模报告时这一条要写进去否则会被质疑「公式算出来 80% 就够了为什么现实这么难」。均匀分布假设一旦打破局部仍有足够易感者传播链就断不干净。3.3 被传染比例的快速估算文档第七部分给了一个实用近似当 i0 很小、s0 接近 1 时被传染人数比例 x 约为阈值超出部分的 2 倍。公式是 x ≈ 2(1 - 1/σ)。这个结论在快速评估时很好用比如 σ3.331/σ0.3x≈2×0.71.4这显然超过 1 了说明近似有适用条件。文档里写的是「当阈值提高时这个比例降低」实际用的时候要检查 s0 是否真的接近 1以及 σ 是否在合理范围。我一般只在 σ 略大于 1 时用这个估算σ 太大就直接跑数值解。4. 避坑与排查跑 SIR 模型时最容易翻车的五件事4.1 现象ode45 报错或结果全是 NaN原因函数名和文件名不一致或者 x0 顺序写反。MATLAB 的 ode45 调用旧式字符串函数名时要求函数文件在路径里且返回列向量。x0 如果写成 [s0, i0] 而不是 [i0, s0]方程里的 x(1)、x(2) 就对不上可能出现负数开方或除零。解决统一用文档里的顺序 [i0, s0]函数返回 [di/dt; ds/dt]。跑之前先单独调用 ill(0, x0) 看返回值是否合理。4.2 现象峰值时间和文档对不上原因λ、μ 的单位没统一。文档里 λ1 是「每天接触 1 个有效人数」μ0.3 是「每天治愈 30%」时间单位都是天。如果 μ 用了周治愈率峰值就会差 7 倍。解决所有参数统一到天。平均传染期 1/μ 算出来是 3.33 天和文档一致再往下做。4.3 现象s(t) 最终不收敛一直往下掉原因把 r(t) 也写进了微分方程但初始条件没满足 sir1或者方程里多加了出生死亡项。文档假设 N 不变没有出生死亡方程里就不该出现常数项。解决检查三个方程加起来是否等于 0即 ds/dtdi/dtdr/dt0。如果用了 r 的方程确保 r01-s0-i0。4.4 现象相轨线画出来是直线或乱线原因plot(x(:,2), x(:,1)) 里 x 的列顺序搞混或者 ts 步长太大导致轨线锯齿。ode45 输出的是自适应步长但 ts 指定了输出点步长太大时轨线会失真。解决相轨线用 plot(x(:,2), x(:,1))x(:,2) 是 sx(:,1) 是 i。ts 步长建议 0.1 或更小画图前用 length(t) 检查点数。4.5 现象群体免疫门槛算出来大于 1原因σ 小于 1即 λμ此时 1-1/σ 为负公式失效。σ1 意味着一个病人在传染期内传染不到一个人疫情本来就不会扩散不需要群体免疫。解决先判断 σ 是否大于 1。只有 σ1 时群体免疫门槛才有意义且门槛值在 0 到 1 之间。5. 从数值解到报告验证、估算与一个我常做的检查文档第六部分用 1900 年代初印度孟买瘟疫数据做了模型验证思路是把理论曲线和实际移出人数画在一起看吻合程度。这个做法放到今天依然有效但要注意实际数据往往只有每日新增或死亡需要先转成 r(t) 的累计值再和模型输出对比。我一般会做三步验证第一步把实际累计移出人数归一化到 [0,1]第二步用最小二乘拟合 σ 和初始条件第三步看残差是否在峰值附近最大——如果是说明模型对峰值时刻的预测有系统偏差报告里要写明。第七部分的被传染比例估算我习惯用数值解反推而不是直接套近似公式。具体做法是跑完 ode45 后取 s(end)被传染比例就是 s(0)-s(end)。文档算例里 s(0)0.98s(end)0.0398被传染比例约 0.94也就是 94% 的人最终被感染。这个数比近似公式靠谱因为不依赖 Taylor 展开的截断误差。最后一个技巧是关于参数敏感性。σ 从 3.33 变到 2.67峰值从 0.344 降到 0.28降幅约 19%但最终感染比例从 96% 降到 92.5%只降了 3.5 个百分点。这说明降低接触率对削峰效果明显但对最终感染规模影响有限。做干预评估时如果目标是压平曲线、争取医疗资源就盯峰值如果目标是减少总感染人数得把 σ 压到接近 1 才有效。这个区别我在第一次做建模报告时没分清被问得哑口无言。从那以后每次跑完 SIR 模型我都会强制把峰值、峰值时间、最终感染比例、最终易感比例四个数一起列出来再下结论。希望帮到你。本文还有配套的精品资源点击获取