简介：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三个变量，N=S+I+R在无出生和死亡假设下是常数。康复者以μ的速率重新变回易感者，这一项写成μR，它让模型具备“第二波疫情”的表达能力：当μ足够大时，I曲线会出现明显的多峰振荡，而不是SIR那种单调走完一波就结束的单峰形态。
用微分方程表达，右端函数需要计算三项变化速率：dS/dt=-β·S·I/N+μ·R，dI/dt=β·S·I/N-γ·I，dR/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) 易感者 S，y(2) 感染者 I，y(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，结果就是模型退化成SIR，R曲线只升不降，第二波疫情完全消失。调用时不要直接在Fsirs里引用全局变量，Matlab全局变量的调试成本偏高。常见做法是让run_SIRS.m用匿名函数把参数绑进去：@(t,y) Fsirs(t,y,beta,gamma,mu,N)。这样Fsirs保持纯函数，换参数不用改方程文件，多个工况循环时也不会因为共享变量产生脏数据。
2.3 参数表与最容易踩的量纲陷阱
| 参数 | 含义 | 常用初值 | 量纲 | 说明 |
|---|---|---|---|---|
| beta | 有效传染率 | 0.4~1.0 | 1/天 | 等于有效接触次数乘以单次传染概率 |
| gamma | 康复率 | 0.05~0.3 | 1/天 | 等于1除以平均感染期天数 |
| mu | 免疫丧失率 | 0~0.05 | 1/天 | 等于1除以平均免疫持续天数 |
| N | 总人口 | 10000 | 人 | S+I+R恒定，初值必须与之一致 |
最容易踩的量纲陷阱是把beta写成人均接触人数而不是每日速率。假如beta=2，它不代表一个感染者能传染两个人，而是每人每天的有效接触次数乘概率后的综合值。另一个常见错误是初值不守恒，比如N=10000，初值却写S0=99990，此时S/N约等于10，传播项被人为放大一个数量级，I峰值和峰位都会失真。SIRS模型不要求S0+I0+R0严格等于整数N，但误差应在浮点精度范围内。建议在run_SIRS.m开头加一行assert(abs(S0+I0+R0-N)<1e-6)，把错误挡在求解之前。检查量纲时还要注意mu与gamma不要写成百分比，gamma=0.1代表每天10%的感染者康复，而不是整个感染期康复10%，理解错会造成曲线拖尾长到几乎不下降。
3. run_SIRS.m主脚本：参数扫描与情景仿真
3.1 脚本初始化与时间网格设计
run_SIRS.m的典型结构分为参数区、求解区、作图区三块。参数区不要把所有数字散写在各行，我一般习惯把beta、gamma、mu、N、S0、I0、R0集中放在脚本顶部，并用注释标注来源，这样后续用matlab优化工具箱做拟合时，只需要替换这一整块。求解区的核心是选择时间跨度，常见两种写法：tspan=[0 365]和tspan=0: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小于gamma |
| R曲线 | 上升至高位后缓慢回落 | 只升不降：mu漏项或mu=0 |
| S曲线 | 先降后出现小幅回升 | 持续下降：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 刚性问题：为什么有时必须换ode15s
SIRS的三条方程在参数悬殊时会出现明显的刚性。比如把mu设成0.001，gamma设为0.5，I和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算成负值后又通过传播项把整个系统带偏。相比ode45，ode15s在参数正常时不明显更慢，但能安抚刚性参数组合，建议在参数扫描代码里直接默认使用它，省一次返工。
4.2 初值设置违反人口守恒的后果
SIRS的初值有三个硬约束：总和等于N、I0大于0、各分量非负。违反第一项时模型不会报错，只会给出荒谬的传播规模和峰值时间，因为传播项的分母始终是固定的N而不是当前总人口。我之前排查过一个案例：S0=9999、I0=1、R0=0，但N误写成1000，导致S/N约等于10，初始再生数被放大了十倍，I峰值对应到实际数据完全对不上。这种错误用assert检查最有效：
assert(abs(S0 + I0 + R0 - N) < 1e-6, '初值总和不等于N'); assert(I0 > 0, '初始感染者必须大于0');如果初值全为零，N不会被改变，求解器会直接返回全零曲线，看起来像没有疫情发生。另一个敏感点是I0相对于N的比例，当I0=1、N=1000万时，早期增长段被压缩到初始几步，曲线会显得平坦，但这其实是正常的，实际疫情数据也经常需要把潜伏期过程包含在初始阶段。做参数扫描时，建议给每种工况打印一行日志，把S0、I0、R0和N四个值都显示出来，肉眼扫一遍比事后核对曲线快得多。
4.3 报错信息与修复手段
| 报错或异常表现 | 常见原因 | 处理方法 |
|---|---|---|
| Undefined function 'Fsirs' | 文件名大小写或路径不在当前目录 | 运行which Fsirs确认路径，切换目录 |
| Dimensions of matrices being concatenated are not consistent | y0写成行向量或Fsirs返回值方向不一致 | 把y0和dydt统一为列向量 |
| Solver unable to meet integration tolerances... | 参数刚性组合 | 改用ode15s并设置NonNegative |
| I曲线初始段为0 | I0被误设为0 | 检查y0中第二个数大于0 |
| beta很大时曲线依然无峰值 | 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里加入季节项，而不是继续堆参数。
本文还有配套的精品资源，点击获取
)
