简介SIRS传染病动力学模型的Matlab源码与数据包面向生物数学、流行病学建模领域的学习者与科研人员适用于数学建模竞赛、课程实验和传染病趋势仿真。模型将人群分为易感、感染、康复、免疫四类通过四个微分方程描述状态转化借助该代码可直观观察传染率、康复率、免疫率对感染峰值和持续时间的影响也能对比不同公共卫生干预策略的效果。压缩包共包含四个文件含两个m文件和两张示例图片整体仅80KB轻量紧凑便于快速下载运行和二次修改其中脚本承担方程定义和主流程模拟图片可直接用于结果对照或报告插图。目前已有199人学习下载非常适合作为SIRS入门到进阶的起点代码稍加调整即可扩展为SEIR等拓展模型是完成建模作业与科研复现的实用工具。1. 为什么SIRS模型比SIR更贴近真实疫情SIRS模型Susceptible-Infectious-Recovered-Susceptible的Matlab实现源码是流行病建模里性价比最高的一类脚本。相比SIR模型SIRS多了一条从R回到S的箭头把康复者免疫力消退后再次易感的真实场景纳入微分方程。源码包里的run_SIRS.m承担参数配置、求解与作图Fsirs.m则是状态方程右端函数另附一张SIRS_model.jpg用来对照状态流转关系。做传染病数据复现、控制措施评估甚至把人群传播模型平移到网络扩散分析的人拿到这份代码后改几个参数就能出图。它不要求你把每个方程从头推导但前提是看懂参数弹性和初值约束这正是下面几章要展开的内容。2. SIRS状态方程与Fsirs.m中的右端函数实现2.1 三状态加上“免疫衰减”的耦合逻辑SIRS虽然常被画成S、I、R三个框但摘要里提到的“免疫”并不是独立微分方程而是R状态里“已经康复但尚未失去免疫力”的群体。真正参与微分方程的是S、I、R三个变量NSIR在无出生和死亡假设下是常数。康复者以μ的速率重新变回易感者这一项写成μR它让模型具备“第二波疫情”的表达能力当μ足够大时I曲线会出现明显的多峰振荡而不是SIR那种单调走完一波就结束的单峰形态。用微分方程表达右端函数需要计算三项变化速率dS/dt-β·S·I/Nμ·RdI/dtβ·S·I/N-γ·IdR/dtγ·I-μ·R。这里的βSI/N不是简单的两个人数相乘而是质量作用律下的有效接触项。除以N是为了把接触率换算为比例避免人口规模直接放大传播项。γ是康复率等于平均感染期的倒数例如平均10天转阴γ≈0.1/天。μ是免疫丧失率如果某种病原体感染后抗体平均维持半年μ≈1/180≈0.0056/天。这几个速率都习惯用“1/天”作量纲看参数表时最容易出错的是β。R0是SIRS的关键判据R0β/γ。当R0小于1时感染人数单调下降R0大于1时初期会出现指数增长段I曲线抬升到峰值后回落。与SIR不同R0大于1的SIRS不会恢复成全体易感而是收敛到地方病平衡点S的终值落在γN/β附近I的终值由μ、γ共同决定。理解这层关系后面做参数扫描时就能预判曲线的走向而不是等图画出来才发现参数方向设反了。2.2 Fsirs.m右端函数怎么组织最稳妥源码包里的Fsirs.m是求解器每次迭代都要调用的函数常见实现是把时间t放在第一个参数位即使方程里对t无引用也要保留占位。下面这份写法在多数Matlab版本上都能直接运行适合当作基本模板function dydt Fsirs(t, y, beta, gamma, mu, N) % SIRS 模型的右端函数 % y(1) 易感者 Sy(2) 感染者 Iy(3) 康复者 R S y(1); I y(2); R y(3); dSdt -beta * S * I / N mu * R; dIdt beta * S * I / N - gamma * I; dRdt gamma * I - mu * R; dydt [dSdt; dIdt; dRdt]; end这段代码里dydt必须按列向量拼接因为ode45要求右端函数返回的导数与状态向量y维度一致且方向相同。很多人在把SIR模板改成SIRS时会在dRdt这一行漏掉-mu*R结果就是模型退化成SIRR曲线只升不降第二波疫情完全消失。调用时不要直接在Fsirs里引用全局变量Matlab全局变量的调试成本偏高。常见做法是让run_SIRS.m用匿名函数把参数绑进去(t,y) Fsirs(t,y,beta,gamma,mu,N)。这样Fsirs保持纯函数换参数不用改方程文件多个工况循环时也不会因为共享变量产生脏数据。2.3 参数表与最容易踩的量纲陷阱参数含义常用初值量纲说明beta有效传染率0.4~1.01/天等于有效接触次数乘以单次传染概率gamma康复率0.05~0.31/天等于1除以平均感染期天数mu免疫丧失率0~0.051/天等于1除以平均免疫持续天数N总人口10000人SIR恒定初值必须与之一致最容易踩的量纲陷阱是把beta写成人均接触人数而不是每日速率。假如beta2它不代表一个感染者能传染两个人而是每人每天的有效接触次数乘概率后的综合值。另一个常见错误是初值不守恒比如N10000初值却写S099990此时S/N约等于10传播项被人为放大一个数量级I峰值和峰位都会失真。SIRS模型不要求S0I0R0严格等于整数N但误差应在浮点精度范围内。建议在run_SIRS.m开头加一行assert(abs(S0I0R0-N)1e-6)把错误挡在求解之前。检查量纲时还要注意mu与gamma不要写成百分比gamma0.1代表每天10%的感染者康复而不是整个感染期康复10%理解错会造成曲线拖尾长到几乎不下降。3. run_SIRS.m主脚本参数扫描与情景仿真3.1 脚本初始化与时间网格设计run_SIRS.m的典型结构分为参数区、求解区、作图区三块。参数区不要把所有数字散写在各行我一般习惯把beta、gamma、mu、N、S0、I0、R0集中放在脚本顶部并用注释标注来源这样后续用matlab优化工具箱做拟合时只需要替换这一整块。求解区的核心是选择时间跨度常见两种写法tspan[0 365]和tspan0:1:365。前者把自适应步长完全交给ode45后者强制求解器在整点输出结果便于和多天粒度的真实疫情数据对齐。这里有一个容易被忽略的点打开Matlab后双击run_SIRS.m不一定会运行很多人以为是matlab安装环境问题其实只是当前文件夹没有切到源码目录。确认左上方当前文件夹路径位于SIRS.zip解压后的目录Fsirs.m文件必须和run_SIRS.m在同一级路径下否则会报Undefined function Fsirs。基础环境里不需要额外工具箱R2020a以后的版本跑这份代码都没有问题。如果你是在新装的Matlab上第一次跑脚本先把当前文件夹设置对再处理参数扫描。3.2 多组参数循环模拟与数据收集单次求解只需要一次ode45调用但真正有用的是把参数扫描结果整理成可对比的数据结构。下面这段脚本把beta从0.3扫到0.9记录每组参数下的感染峰值和峰值时间并预留了存放完整曲线的cell数组beta_list [0.3, 0.6, 0.9]; gamma 0.1; mu 0.02; N 10000; S0 N - 1; I0 1; R0 0; tspan 0:1:365; y0 [S0, I0, R0]; peak_info zeros(length(beta_list), 3); % 行存 beta, 峰值, 峰位 curves cell(length(beta_list), 1); for k 1:length(beta_list) beta beta_list(k); [t, y] ode45((t,y) Fsirs(t,y,beta,gamma,mu,N), tspan, y0); [I_peak, idx] max(y(:,2)); peak_info(k,:) [beta, I_peak, t(idx)]; curves{k} [t, y]; end循环里用k作为索引而不是直接遍历beta_list中的值再用浮点数比对避免0.3在二进制表示下不等价于列表元素的问题。峰值和峰位用max的第二个输出得到如果要评估“多大beta会导致医疗资源击穿”这组数据可以直接画成曲线peak_info(:,1)为横轴peak_info(:,2)为纵轴。curves里保存了每次完整求解结果后续做动画或拖尾分析时不用重算代价是365乘3的double矩阵乘以工况数量内存开销极小不需要预分配大数组。3.3 用SIRS_model.jpg核对输出图像源码包里的SIRS_model.jpg画的是状态转移示意而不是仿真结果图所以用它核对位置要谨慎。真正的核对点是仿真曲线的特征形状当beta大于gamma时I曲线应在模拟期前四分之一段出现明显峰值R曲线同步上升后缓慢回落S曲线降到最低点后略有回升这是因为mu让一部分R重新回到S。如果S曲线单调下降没有回升先检查Fsirs.m里dSdt是否带上了muR项如果R曲线完全不下落检查dRdt是否少了-muR。用Matlab的绘图工具直接看这三条曲线的相对顺序比盯着一堆数字更直观。下面这张表总结了最需要看的三个特征检查项曲线正常特征常见参数错误I曲线先升后降存在明显峰值单调下降beta小于gammaR曲线上升至高位后缓慢回落只升不降mu漏项或mu0S曲线先降后出现小幅回升持续下降dSdt漏掉加mu*R运行参数扫描后画多子图也很顺手常见做法是subplot(3,1,k)配合hold on绘制三条状态曲线legend放在每个子图顶部。如果你习惯的是把外部观测数据从csv导入到matlab中进行fft仿真那一套流程要特别记住这里y矩阵的行是时间点、列是状态变量读y(:,2)得到感染人数和fft流程里行对应通道的约定完全不同不要套用。另外有人把tspan写成linspace(0,365,1000)并期望它带来更高精度其实tspan只是输出采样点不影响ode45求解步长真正的精度控制靠odeset里的RelTol和AbsTol。下一章讲的就是这些数值控制参数。4. 数值稳定性、初值敏感性与常见报错排查4.1 刚性问题为什么有时必须换ode15sSIRS的三条方程在参数悬殊时会出现明显的刚性。比如把mu设成0.001gamma设为0.5I和R的动力学在几天内完成而S的变化要持续数百天直接调用ode45不是报错而是求解步长会被压到极小仿真一年需要几十万步跑起来明显卡顿有些版本还会提示Solver unable to meet integration tolerances without reducing the step size below the smallest value allowed。这不是代码逻辑错误是数值方法选型不对。常见做法是换成ode15s它对刚性系统使用后向差分格式稳定域更大步长可以放宽opts odeset(RelTol, 1e-6, AbsTol, 1e-8, NonNegative, 1:3); [t, y] ode15s((t,y) Fsirs(t,y,beta,gamma,mu,N), tspan, y0, opts);RelTol控制每个分量相对误差AbsTol控制接近零时的绝对误差。SIRS状态变量都是人数在疫情末期I可能降到10的负二次方量级AbsTol设1e-8足够保证不出现负人数。NonNegative的作用是强制S、I、R非负防止数值振荡把I算成负值后又通过传播项把整个系统带偏。相比ode45ode15s在参数正常时不明显更慢但能安抚刚性参数组合建议在参数扫描代码里直接默认使用它省一次返工。4.2 初值设置违反人口守恒的后果SIRS的初值有三个硬约束总和等于N、I0大于0、各分量非负。违反第一项时模型不会报错只会给出荒谬的传播规模和峰值时间因为传播项的分母始终是固定的N而不是当前总人口。我之前排查过一个案例S09999、I01、R00但N误写成1000导致S/N约等于10初始再生数被放大了十倍I峰值对应到实际数据完全对不上。这种错误用assert检查最有效assert(abs(S0 I0 R0 - N) 1e-6, 初值总和不等于N); assert(I0 0, 初始感染者必须大于0);如果初值全为零N不会被改变求解器会直接返回全零曲线看起来像没有疫情发生。另一个敏感点是I0相对于N的比例当I01、N1000万时早期增长段被压缩到初始几步曲线会显得平坦但这其实是正常的实际疫情数据也经常需要把潜伏期过程包含在初始阶段。做参数扫描时建议给每种工况打印一行日志把S0、I0、R0和N四个值都显示出来肉眼扫一遍比事后核对曲线快得多。4.3 报错信息与修复手段报错或异常表现常见原因处理方法Undefined function Fsirs文件名大小写或路径不在当前目录运行which Fsirs确认路径切换目录Dimensions of matrices being concatenated are not consistenty0写成行向量或Fsirs返回值方向不一致把y0和dydt统一为列向量Solver unable to meet integration tolerances...参数刚性组合改用ode15s并设置NonNegativeI曲线初始段为0I0被误设为0检查y0中第二个数大于0beta很大时曲线依然无峰值beta_list里误写0打印beta_list检查参数录入第一行是最常见的“假报错”很多新装Matlab的机器双击脚本不运行其实是当前文件夹不对。注意如果从旧版本Matlab迁移代码Fsirs.m头部的function签名不要改动函数名大小写在这种调用中是敏感的。最后一行提到的退化情况也值得留意beta设为0时方程线性化只会衰减不会传播仿真很快但没有任何探究意义参数扫描前先检查beta_list里没有0。如果碰到了表格之外的报错优先看错误栈指向的是Fsirs.m里的哪一行因为绝大多数维度问题都出在状态向量的拼接逻辑上。5. 把SIRS扩充成带干预措施的时变参数模型5.1 把beta改成时变函数静态beta只能描述无干预传播真正做政策评估时需要让beta随时间变化。常见做法是不改动Fsirs.m而是在运行时构造一个beta(t)函数并传入beta_t (t) 0.8 * (t 30) 0.3 * (t 30 t 90) 0.6 * (t 90); [t, y] ode45((t,y) Fsirs(t, y, beta_t(t), gamma, mu, N), tspan, y0);这里第30天到第90天把传播率压到0.3模拟封锁措施第90天放松到0.6。由于beta_t(t)在每一个求解步内都被重新计算变化时刻附近可能出现步长加密这是正常的。用这种方式不需要新增单独的M文件也方便比较多种措施组合。5.2 量化隔离措施的相对效果把上一节的代码放进参数扫描框架里可以算出每组干预策略下的累计感染负担。累计感染不是直接取y(:,2)的最大值而是用trapz(t, y(:,2))计算感染人日数以30天为界的前后对照能算出隔离措施减少了多少感染人日。对于需要在报告里写量化结论的场景这个数字比“曲线变得更平”更有说服力。5.3 用matlab优化工具箱反推参数如果手里有真实疫情数据可以用matlab优化工具箱里的lsqcurvefit或fminsearch反推beta和gamma。核心是把ode45的求解过程包成一个可被优化器调用的函数输入参数向量p输出感染人数的时间序列function I_pred sim_model(p, t_obs, N, mu) beta p(1); gamma p(2); [~, y] ode45((t,y) Fsirs(t,y,beta,gamma,mu,N), t_obs, [N-1, 1, 0]); I_pred y(:,2); end目标函数计算模拟感染人数与真实感染人数之差的平方和然后交给fminsearch做无梯度搜索。注意sim_model里的t_obs必须是观测数据的时间点而真实数据中潜伏期造成的滞后会让拟合出来的beta系统性偏小。拟合出的参数序列还有一个用处如果后续要做深度学习matlab时间序列预测把beta、gamma随时间的变化当作特征输入比直接用原始I曲线更贴近传播机制。将拟合结果代回run_SIRS.m顶部替换beta、gamma后如果残差仍有周期性波动下一步就该检查是否需要在beta里加入季节项而不是继续堆参数。本文还有配套的精品资源点击获取 SIRS模型Matlab实现：从微分方程到参数扫描的完整源码解析
简介SIRS传染病动力学模型的Matlab源码与数据包面向生物数学、流行病学建模领域的学习者与科研人员适用于数学建模竞赛、课程实验和传染病趋势仿真。模型将人群分为易感、感染、康复、免疫四类通过四个微分方程描述状态转化借助该代码可直观观察传染率、康复率、免疫率对感染峰值和持续时间的影响也能对比不同公共卫生干预策略的效果。压缩包共包含四个文件含两个m文件和两张示例图片整体仅80KB轻量紧凑便于快速下载运行和二次修改其中脚本承担方程定义和主流程模拟图片可直接用于结果对照或报告插图。目前已有199人学习下载非常适合作为SIRS入门到进阶的起点代码稍加调整即可扩展为SEIR等拓展模型是完成建模作业与科研复现的实用工具。1. 为什么SIRS模型比SIR更贴近真实疫情SIRS模型Susceptible-Infectious-Recovered-Susceptible的Matlab实现源码是流行病建模里性价比最高的一类脚本。相比SIR模型SIRS多了一条从R回到S的箭头把康复者免疫力消退后再次易感的真实场景纳入微分方程。源码包里的run_SIRS.m承担参数配置、求解与作图Fsirs.m则是状态方程右端函数另附一张SIRS_model.jpg用来对照状态流转关系。做传染病数据复现、控制措施评估甚至把人群传播模型平移到网络扩散分析的人拿到这份代码后改几个参数就能出图。它不要求你把每个方程从头推导但前提是看懂参数弹性和初值约束这正是下面几章要展开的内容。2. SIRS状态方程与Fsirs.m中的右端函数实现2.1 三状态加上“免疫衰减”的耦合逻辑SIRS虽然常被画成S、I、R三个框但摘要里提到的“免疫”并不是独立微分方程而是R状态里“已经康复但尚未失去免疫力”的群体。真正参与微分方程的是S、I、R三个变量NSIR在无出生和死亡假设下是常数。康复者以μ的速率重新变回易感者这一项写成μR它让模型具备“第二波疫情”的表达能力当μ足够大时I曲线会出现明显的多峰振荡而不是SIR那种单调走完一波就结束的单峰形态。用微分方程表达右端函数需要计算三项变化速率dS/dt-β·S·I/Nμ·RdI/dtβ·S·I/N-γ·IdR/dtγ·I-μ·R。这里的βSI/N不是简单的两个人数相乘而是质量作用律下的有效接触项。除以N是为了把接触率换算为比例避免人口规模直接放大传播项。γ是康复率等于平均感染期的倒数例如平均10天转阴γ≈0.1/天。μ是免疫丧失率如果某种病原体感染后抗体平均维持半年μ≈1/180≈0.0056/天。这几个速率都习惯用“1/天”作量纲看参数表时最容易出错的是β。R0是SIRS的关键判据R0β/γ。当R0小于1时感染人数单调下降R0大于1时初期会出现指数增长段I曲线抬升到峰值后回落。与SIR不同R0大于1的SIRS不会恢复成全体易感而是收敛到地方病平衡点S的终值落在γN/β附近I的终值由μ、γ共同决定。理解这层关系后面做参数扫描时就能预判曲线的走向而不是等图画出来才发现参数方向设反了。2.2 Fsirs.m右端函数怎么组织最稳妥源码包里的Fsirs.m是求解器每次迭代都要调用的函数常见实现是把时间t放在第一个参数位即使方程里对t无引用也要保留占位。下面这份写法在多数Matlab版本上都能直接运行适合当作基本模板function dydt Fsirs(t, y, beta, gamma, mu, N) % SIRS 模型的右端函数 % y(1) 易感者 Sy(2) 感染者 Iy(3) 康复者 R S y(1); I y(2); R y(3); dSdt -beta * S * I / N mu * R; dIdt beta * S * I / N - gamma * I; dRdt gamma * I - mu * R; dydt [dSdt; dIdt; dRdt]; end这段代码里dydt必须按列向量拼接因为ode45要求右端函数返回的导数与状态向量y维度一致且方向相同。很多人在把SIR模板改成SIRS时会在dRdt这一行漏掉-mu*R结果就是模型退化成SIRR曲线只升不降第二波疫情完全消失。调用时不要直接在Fsirs里引用全局变量Matlab全局变量的调试成本偏高。常见做法是让run_SIRS.m用匿名函数把参数绑进去(t,y) Fsirs(t,y,beta,gamma,mu,N)。这样Fsirs保持纯函数换参数不用改方程文件多个工况循环时也不会因为共享变量产生脏数据。2.3 参数表与最容易踩的量纲陷阱参数含义常用初值量纲说明beta有效传染率0.4~1.01/天等于有效接触次数乘以单次传染概率gamma康复率0.05~0.31/天等于1除以平均感染期天数mu免疫丧失率0~0.051/天等于1除以平均免疫持续天数N总人口10000人SIR恒定初值必须与之一致最容易踩的量纲陷阱是把beta写成人均接触人数而不是每日速率。假如beta2它不代表一个感染者能传染两个人而是每人每天的有效接触次数乘概率后的综合值。另一个常见错误是初值不守恒比如N10000初值却写S099990此时S/N约等于10传播项被人为放大一个数量级I峰值和峰位都会失真。SIRS模型不要求S0I0R0严格等于整数N但误差应在浮点精度范围内。建议在run_SIRS.m开头加一行assert(abs(S0I0R0-N)1e-6)把错误挡在求解之前。检查量纲时还要注意mu与gamma不要写成百分比gamma0.1代表每天10%的感染者康复而不是整个感染期康复10%理解错会造成曲线拖尾长到几乎不下降。3. run_SIRS.m主脚本参数扫描与情景仿真3.1 脚本初始化与时间网格设计run_SIRS.m的典型结构分为参数区、求解区、作图区三块。参数区不要把所有数字散写在各行我一般习惯把beta、gamma、mu、N、S0、I0、R0集中放在脚本顶部并用注释标注来源这样后续用matlab优化工具箱做拟合时只需要替换这一整块。求解区的核心是选择时间跨度常见两种写法tspan[0 365]和tspan0:1:365。前者把自适应步长完全交给ode45后者强制求解器在整点输出结果便于和多天粒度的真实疫情数据对齐。这里有一个容易被忽略的点打开Matlab后双击run_SIRS.m不一定会运行很多人以为是matlab安装环境问题其实只是当前文件夹没有切到源码目录。确认左上方当前文件夹路径位于SIRS.zip解压后的目录Fsirs.m文件必须和run_SIRS.m在同一级路径下否则会报Undefined function Fsirs。基础环境里不需要额外工具箱R2020a以后的版本跑这份代码都没有问题。如果你是在新装的Matlab上第一次跑脚本先把当前文件夹设置对再处理参数扫描。3.2 多组参数循环模拟与数据收集单次求解只需要一次ode45调用但真正有用的是把参数扫描结果整理成可对比的数据结构。下面这段脚本把beta从0.3扫到0.9记录每组参数下的感染峰值和峰值时间并预留了存放完整曲线的cell数组beta_list [0.3, 0.6, 0.9]; gamma 0.1; mu 0.02; N 10000; S0 N - 1; I0 1; R0 0; tspan 0:1:365; y0 [S0, I0, R0]; peak_info zeros(length(beta_list), 3); % 行存 beta, 峰值, 峰位 curves cell(length(beta_list), 1); for k 1:length(beta_list) beta beta_list(k); [t, y] ode45((t,y) Fsirs(t,y,beta,gamma,mu,N), tspan, y0); [I_peak, idx] max(y(:,2)); peak_info(k,:) [beta, I_peak, t(idx)]; curves{k} [t, y]; end循环里用k作为索引而不是直接遍历beta_list中的值再用浮点数比对避免0.3在二进制表示下不等价于列表元素的问题。峰值和峰位用max的第二个输出得到如果要评估“多大beta会导致医疗资源击穿”这组数据可以直接画成曲线peak_info(:,1)为横轴peak_info(:,2)为纵轴。curves里保存了每次完整求解结果后续做动画或拖尾分析时不用重算代价是365乘3的double矩阵乘以工况数量内存开销极小不需要预分配大数组。3.3 用SIRS_model.jpg核对输出图像源码包里的SIRS_model.jpg画的是状态转移示意而不是仿真结果图所以用它核对位置要谨慎。真正的核对点是仿真曲线的特征形状当beta大于gamma时I曲线应在模拟期前四分之一段出现明显峰值R曲线同步上升后缓慢回落S曲线降到最低点后略有回升这是因为mu让一部分R重新回到S。如果S曲线单调下降没有回升先检查Fsirs.m里dSdt是否带上了muR项如果R曲线完全不下落检查dRdt是否少了-muR。用Matlab的绘图工具直接看这三条曲线的相对顺序比盯着一堆数字更直观。下面这张表总结了最需要看的三个特征检查项曲线正常特征常见参数错误I曲线先升后降存在明显峰值单调下降beta小于gammaR曲线上升至高位后缓慢回落只升不降mu漏项或mu0S曲线先降后出现小幅回升持续下降dSdt漏掉加mu*R运行参数扫描后画多子图也很顺手常见做法是subplot(3,1,k)配合hold on绘制三条状态曲线legend放在每个子图顶部。如果你习惯的是把外部观测数据从csv导入到matlab中进行fft仿真那一套流程要特别记住这里y矩阵的行是时间点、列是状态变量读y(:,2)得到感染人数和fft流程里行对应通道的约定完全不同不要套用。另外有人把tspan写成linspace(0,365,1000)并期望它带来更高精度其实tspan只是输出采样点不影响ode45求解步长真正的精度控制靠odeset里的RelTol和AbsTol。下一章讲的就是这些数值控制参数。4. 数值稳定性、初值敏感性与常见报错排查4.1 刚性问题为什么有时必须换ode15sSIRS的三条方程在参数悬殊时会出现明显的刚性。比如把mu设成0.001gamma设为0.5I和R的动力学在几天内完成而S的变化要持续数百天直接调用ode45不是报错而是求解步长会被压到极小仿真一年需要几十万步跑起来明显卡顿有些版本还会提示Solver unable to meet integration tolerances without reducing the step size below the smallest value allowed。这不是代码逻辑错误是数值方法选型不对。常见做法是换成ode15s它对刚性系统使用后向差分格式稳定域更大步长可以放宽opts odeset(RelTol, 1e-6, AbsTol, 1e-8, NonNegative, 1:3); [t, y] ode15s((t,y) Fsirs(t,y,beta,gamma,mu,N), tspan, y0, opts);RelTol控制每个分量相对误差AbsTol控制接近零时的绝对误差。SIRS状态变量都是人数在疫情末期I可能降到10的负二次方量级AbsTol设1e-8足够保证不出现负人数。NonNegative的作用是强制S、I、R非负防止数值振荡把I算成负值后又通过传播项把整个系统带偏。相比ode45ode15s在参数正常时不明显更慢但能安抚刚性参数组合建议在参数扫描代码里直接默认使用它省一次返工。4.2 初值设置违反人口守恒的后果SIRS的初值有三个硬约束总和等于N、I0大于0、各分量非负。违反第一项时模型不会报错只会给出荒谬的传播规模和峰值时间因为传播项的分母始终是固定的N而不是当前总人口。我之前排查过一个案例S09999、I01、R00但N误写成1000导致S/N约等于10初始再生数被放大了十倍I峰值对应到实际数据完全对不上。这种错误用assert检查最有效assert(abs(S0 I0 R0 - N) 1e-6, 初值总和不等于N); assert(I0 0, 初始感染者必须大于0);如果初值全为零N不会被改变求解器会直接返回全零曲线看起来像没有疫情发生。另一个敏感点是I0相对于N的比例当I01、N1000万时早期增长段被压缩到初始几步曲线会显得平坦但这其实是正常的实际疫情数据也经常需要把潜伏期过程包含在初始阶段。做参数扫描时建议给每种工况打印一行日志把S0、I0、R0和N四个值都显示出来肉眼扫一遍比事后核对曲线快得多。4.3 报错信息与修复手段报错或异常表现常见原因处理方法Undefined function Fsirs文件名大小写或路径不在当前目录运行which Fsirs确认路径切换目录Dimensions of matrices being concatenated are not consistenty0写成行向量或Fsirs返回值方向不一致把y0和dydt统一为列向量Solver unable to meet integration tolerances...参数刚性组合改用ode15s并设置NonNegativeI曲线初始段为0I0被误设为0检查y0中第二个数大于0beta很大时曲线依然无峰值beta_list里误写0打印beta_list检查参数录入第一行是最常见的“假报错”很多新装Matlab的机器双击脚本不运行其实是当前文件夹不对。注意如果从旧版本Matlab迁移代码Fsirs.m头部的function签名不要改动函数名大小写在这种调用中是敏感的。最后一行提到的退化情况也值得留意beta设为0时方程线性化只会衰减不会传播仿真很快但没有任何探究意义参数扫描前先检查beta_list里没有0。如果碰到了表格之外的报错优先看错误栈指向的是Fsirs.m里的哪一行因为绝大多数维度问题都出在状态向量的拼接逻辑上。5. 把SIRS扩充成带干预措施的时变参数模型5.1 把beta改成时变函数静态beta只能描述无干预传播真正做政策评估时需要让beta随时间变化。常见做法是不改动Fsirs.m而是在运行时构造一个beta(t)函数并传入beta_t (t) 0.8 * (t 30) 0.3 * (t 30 t 90) 0.6 * (t 90); [t, y] ode45((t,y) Fsirs(t, y, beta_t(t), gamma, mu, N), tspan, y0);这里第30天到第90天把传播率压到0.3模拟封锁措施第90天放松到0.6。由于beta_t(t)在每一个求解步内都被重新计算变化时刻附近可能出现步长加密这是正常的。用这种方式不需要新增单独的M文件也方便比较多种措施组合。5.2 量化隔离措施的相对效果把上一节的代码放进参数扫描框架里可以算出每组干预策略下的累计感染负担。累计感染不是直接取y(:,2)的最大值而是用trapz(t, y(:,2))计算感染人日数以30天为界的前后对照能算出隔离措施减少了多少感染人日。对于需要在报告里写量化结论的场景这个数字比“曲线变得更平”更有说服力。5.3 用matlab优化工具箱反推参数如果手里有真实疫情数据可以用matlab优化工具箱里的lsqcurvefit或fminsearch反推beta和gamma。核心是把ode45的求解过程包成一个可被优化器调用的函数输入参数向量p输出感染人数的时间序列function I_pred sim_model(p, t_obs, N, mu) beta p(1); gamma p(2); [~, y] ode45((t,y) Fsirs(t,y,beta,gamma,mu,N), t_obs, [N-1, 1, 0]); I_pred y(:,2); end目标函数计算模拟感染人数与真实感染人数之差的平方和然后交给fminsearch做无梯度搜索。注意sim_model里的t_obs必须是观测数据的时间点而真实数据中潜伏期造成的滞后会让拟合出来的beta系统性偏小。拟合出的参数序列还有一个用处如果后续要做深度学习matlab时间序列预测把beta、gamma随时间的变化当作特征输入比直接用原始I曲线更贴近传播机制。将拟合结果代回run_SIRS.m顶部替换beta、gamma后如果残差仍有周期性波动下一步就该检查是否需要在beta里加入季节项而不是继续堆参数。本文还有配套的精品资源点击获取