1. 项目概述从SEIRS模型看疫情传播的数学本质最近几年大家对于传染病传播模型应该都不陌生了。无论是新闻里看到的疫情预测曲线还是各种政策推演的背后都离不开这些数学模型的支持。今天我想和大家深入聊聊的就是这个在经典SIR模型基础上演化而来的SEIRS模型并且会手把手地带你用Matlab把它“跑”起来。这个模型之所以重要是因为它更贴近现实中许多传染病的特性比如流感、新冠等它们感染后获得的免疫力并非永久人会再次变为易感者同时疾病还存在潜伏期并非一接触就立刻发病。SEIRS模型正是通过引入“潜伏者E”和“免疫丧失从R变回S”这两个关键状态来刻画这些复杂动态。简单来说SEIRS模型把人群分成了四类易感者S、潜伏者E、感染者I和康复者R。易感者接触到病毒后不会马上发病而是先进入潜伏期E经过一段时间才具有传染性I康复后R获得暂时免疫力但这份免疫力可能随时间减弱使人重新变回易感者S。整个系统就像一个动态流转的池子而我们的目标就是用数学方程描述这个流转过程并通过计算机模拟来观察不同参数如接触率、潜伏期时长、免疫持续时间下疫情的发展轨迹。这对于理解疫情走势、评估防控措施如提高隔离率相当于降低接触率的效果具有非常直观的参考价值。接下来我将彻底拆解这个模型的每一个环节从微分方程的原理推导到Matlab代码的逐行实现再到如何调整参数模拟不同场景。无论你是数学建模的初学者还是想寻找一个可靠、可复现的传染病模拟案例这篇文章都能给你提供一套完整的“工具箱”。2. SEIRS模型的核心原理与方程拆解要真正玩转一个模型不能只当它是“黑箱”知道输入输出就完事。我们必须深入其数学核心理解每一个方程、每一个参数背后的物理意义和生物学假设。只有这样当模拟结果出现异常时你才知道该调整哪里当需要针对特定疾病定制模型时你也知道如何修改方程。2.1 模型状态定义与流转关系SEIRS模型是一个典型的仓室模型它基于几个核心假设总人口数N恒定不考虑出生、死亡或假设出生率等于自然死亡率即 S E I R N 为常数。这个假设简化了模型让我们专注于疾病本身的传播动力学。均匀混合人群充分混合任何一个易感者S接触任何一个感染者I的概率是相同的。这显然是对现实的简化但对于理解宏观趋势非常有效。疾病进程确定个体从潜伏期E到发病期I再到康复R其转移速率是常数这意味着潜伏期和感染期近似服从指数分布。基于此我们定义四个核心状态变量它们都是关于时间t的函数S(t): 时刻 t 易感者的数量。E(t): 时刻 t 处于潜伏期已感染但未发病、无传染性者的数量。I(t): 时刻 t 感染者已发病、有传染性的数量。R(t): 时刻 t 康复者已康复、具有暂时免疫力的数量。它们之间的流转关系可以用下面的流程图来直观理解注意这里用文字描述代替图表 易感者S以一定的速率被感染但并非直接变成感染者I而是先进入潜伏期E。潜伏者E经过平均潜伏期后转化为感染者I。感染者I经过平均感染期后康复并进入康复者R群体。康复者R获得的免疫力并非永久经过平均免疫期后会重新变为易感者S。此外模型中通常还会考虑因病死亡率即感染者I可能以一定概率死亡这会使得总人口N发生变化但在最基本的SEIRS模型中我们常先忽略这一点或假设死亡率包含在康复率中。2.2 微分方程组的推导与参数解读上述的流转关系可以用一组常微分方程ODE来精确描述。这是整个模型的心脏。我们定义几个关键参数β (Beta): 有效接触率或传播率。它表示一个感染者I在单位时间内有效接触并传染易感者S的人数。在实际中β 接触率 × 每次接触的传染概率。它通常与总人口N有关更常见的写法是 β * I / N表示一个易感者遇到感染者的概率。σ (Sigma): 潜伏期转化率。它是平均潜伏期d_E的倒数即 σ 1 / d_E。例如平均潜伏期为5天则 σ 0.2 /天。表示单位时间内潜伏者E转化为感染者I的比例。γ (Gamma): 康复率。它是平均感染期d_I的倒数即 γ 1 / d_I。例如平均感染期为7天则 γ ≈ 0.143 /天。表示单位时间内感染者I康复进入R的比例。ξ (Xi): 免疫丧失率。它是平均免疫期d_R的倒数即 ξ 1 / d_R。例如免疫力平均维持180天则 ξ ≈ 0.0056 /天。表示单位时间内康复者R丧失免疫、重新变为易感者S的比例。μ (Mu): 自然死亡率或出生率在总人口恒定假设下两者相等。这是一个可选参数用于模拟更长期的动态平衡。现在我们可以写出经典的SEIRS微分方程组dS/dt μ * N ξ * R - β * (I/N) * S - μ * S dE/dt β * (I/N) * S - σ * E - μ * E dI/dt σ * E - γ * I - μ * I dR/dt γ * I - ξ * R - μ * R方程逐行解读易感者S的变化率 dS/dt μ * N假设有常数出生率或输入新增人口全部为易感者。 ξ * R康复者R丧失免疫力重新变为易感者S。- β * (I/N) * S这是核心的感染项。(I/N)是随机接触一人恰好是感染者的概率。β * (I/N)就是一个易感者单位时间内被感染的概率力。再乘以易感者总数S就得到了单位时间内从S流向E的人数。- μ * S易感者的自然死亡。潜伏者E的变化率 dE/dt β * (I/N) * S从易感者S被感染而来是E的“流入项”。- σ * E潜伏者E以速率σ转化为感染者I是E的“流出项”。- μ * E潜伏者的自然死亡。感染者I的变化率 dI/dt σ * E从潜伏者E转化而来是I的“流入项”。- γ * I感染者I以速率γ康复进入R群体是I的“流出项”。- μ * I感染者的自然死亡如果考虑疾病致死率这里可以改为- (γα)*I其中α是因病死亡率。康复者R的变化率 dR/dt γ * I从感染者I康复而来是R的“流入项”。- ξ * R康复者R丧失免疫力变回易感者S是R的“流出项”。- μ * R康复者的自然死亡。注意在许多短期疫情模拟中为了简化会忽略自然出生死亡即 μ0因为疫情周期几个月到一两年远小于人的寿命。此时方程简化为dS/dt ξR - βIS/N,dE/dt βIS/N - σE,dI/dt σE - γI,dR/dt γI - ξR。同时SEIR N恒定。我们后续的Matlab实现将基于这个简化版本因为它更聚焦于疾病传播本身。关键衍生参数——基本再生数 R0 R0是一个极其重要的流行病学指标代表在完全易感人群中一个感染者在其整个传染期内平均能传染的人数。对于SEIRS模型R0的计算公式为R0 β / γ这是因为一个感染者的平均传染期是1/γ在此期间他单位时间传染β个人在完全易感环境下S/N≈1。如果R0 1疾病会流行开来如果R0 1疾病会逐渐消失。在SEIRS模型中由于存在免疫丧失即使R01疫情也可能表现为周期性爆发而非一次流行后终结。3. 基于Matlab的模型实现与代码精讲理论说得再多不如一行代码。下面我们就用Matlab将上述方程组“复活”通过数值求解来观察疫情动态。我将采用最清晰、最易于理解和修改的方式编写代码并解释每一个关键步骤。3.1 模型参数设置与初始化首先我们需要定义模型参数和初始条件。参数的选择需要参考真实疾病的数据这里我们以类似流感的参数为例进行演示。% SEIRS模型参数设置 % 时间参数 t_start 0; % 模拟开始时间天 t_end 500; % 模拟结束时间天 dt 0.1; % 时间步长天用于欧拉法值越小越精确但计算越慢 % 流行病学参数 beta 0.5; % 传播率 (1/天) - 一个感染者每天有效接触人数 sigma 0.2; % 潜伏期转化率 (1/天) - 平均潜伏期 1/sigma 5天 gamma 0.143; % 康复率 (1/天) - 平均感染期 1/gamma ≈ 7天 xi 0.005; % 免疫丧失率 (1/天) - 平均免疫期 1/xi 200天 % 人口参数 N 10000; % 总人口 I0 1; % 初始感染者数量 E0 0; % 初始潜伏者数量 R0 0; % 初始康复者数量 S0 N - I0 - E0 - R0; % 初始易感者数量 % 初始化数组存储结果 num_steps ceil((t_end - t_start) / dt) 1; % 计算总步数 time linspace(t_start, t_end, num_steps); % 时间向量 S zeros(1, num_steps); E zeros(1, num_steps); I zeros(1, num_steps); R zeros(1, num_steps); % 设置初始值 S(1) S0; E(1) E0; I(1) I0; R(1) R0;参数设置心得时间步长dt这里我们使用了最简单的欧拉法进行数值求解。欧拉法是一种显式方法计算简单但稳定性有条件。对于这个模型只要dt设置得足够小比如小于模型中最小时间常数的倒数这里最小的时间常数可能是1/σ5天所以dt0.1是安全的结果就是可靠的。如果你想追求更高的精度和稳定性可以使用Matlab内置的ODE求解器如ode45这部分后面会讲。参数beta这是最敏感、最需要校准的参数。β R0 * γ。如果你从文献中知道了某种疾病的R0可以通过这个关系反推beta。例如流感的R0大约在1.2-1.8之间若γ1/7≈0.143则beta的范围约为0.17-0.26。这里设为0.5是为了让疫情发展更迅速便于观察。初始条件通常假设开始时只有一个或少数几个感染者I0其余均为易感者。潜伏者E0和康复者R0通常设为0。3.2 使用欧拉法进行数值求解欧拉法的核心思想是用当前时刻的导数来近似下一个时刻的值。对于微分方程dy/dt f(y, t)其迭代公式为y_{n1} y_n f(y_n, t_n) * dt。% 使用欧拉法迭代求解SEIRS模型 for i 1:num_steps-1 % 当前时刻的各仓室人数 S_curr S(i); E_curr E(i); I_curr I(i); R_curr R(i); % 计算各仓室的变化率导数基于简化模型μ0 dS_dt xi * R_curr - (beta * I_curr / N) * S_curr; dE_dt (beta * I_curr / N) * S_curr - sigma * E_curr; dI_dt sigma * E_curr - gamma * I_curr; dR_dt gamma * I_curr - xi * R_curr; % 欧拉法更新下一时刻的人数 S(i1) S_curr dS_dt * dt; E(i1) E_curr dE_dt * dt; I(i1) I_curr dI_dt * dt; R(i1) R_curr dR_dt * dt; % 可选施加约束防止数值误差导致人数为负虽然欧拉法在dt很小时很少出现 S(i1) max(0, S(i1)); E(i1) max(0, E(i1)); I(i1) max(0, I(i1)); R(i1) max(0, R(i1)); end欧拉法实现的注意事项计算顺序在计算导数dS_dt时我们使用的是当前时刻i的S_curr,I_curr等值。这意味着所有状态是同步更新的。如果错误地先更新了S再用新的S去计算E的导数就会引入误差破坏模型的一致性。负值处理由于是离散近似在极端参数或大步长下可能会出现某个仓室人数为负的荒谬情况。虽然概率很低但加上max(0, ...)约束是一个良好的编程习惯能保证模拟的鲁棒性。更根本的解决方法是使用更小的时间步长dt或改用更稳定的求解器。精度验证一个简单的验证方法是将dt减半如从0.1变为0.05重新运行模拟。如果两条曲线几乎重合说明当前的dt已经足够精确。如果差异明显则需要继续减小dt。3.3 使用Matlab内置ODE求解器推荐对于常微分方程组Matlab提供了强大且稳定的求解器如ode45基于Runge-Kutta方法。使用它们不仅精度高而且无需手动设置步长代码也更简洁。我强烈推荐在实际建模中使用这种方法。% 使用ode45求解SEIRS模型 % 定义时间跨度 tspan [0, 500]; % 定义初始状态向量 [S0; E0; I0; R0] y0 [S0; E0; I0; R0]; % 调用ode45求解器 % seirs_ode 是定义微分方程组的函数句柄 % tspan 是时间范围 % y0 是初始条件 [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma, gamma, xi, N), tspan, y0); % 提取结果 S_ode y(:, 1); E_ode y(:, 2); I_ode y(:, 3); R_ode y(:, 4);这里我们需要单独定义一个函数文件seirs_ode.m来描述微分方程组function dydt seirs_ode(t, y, beta, sigma, gamma, xi, N) % SEIRS模型微分方程组 % 输入 % t: 时间未显式使用但ode45要求此参数 % y: 状态向量 [S; E; I; R] % beta, sigma, gamma, xi, N: 模型参数 % 输出 % dydt: 导数向量 [dS/dt; dE/dt; dI/dt; dR/dt] S y(1); E y(2); I y(3); R y(4); % 计算各导数 dS_dt xi * R - (beta * I / N) * S; dE_dt (beta * I / N) * S - sigma * E; dI_dt sigma * E - gamma * I; dR_dt gamma * I - xi * R; dydt [dS_dt; dE_dt; dI_dt; dR_dt]; end使用ODE求解器的优势自动变步长ode45会根据方程的特性自动调整积分步长在变化平缓处用大步长提高效率在变化剧烈处用小步长保证精度。高精度采用高阶Runge-Kutta方法数值误差远小于简单的欧拉法。代码简洁省去了手动迭代和步长管理的麻烦。稳定性好对于刚性方程系统中存在快变和慢变过程可以使用专门的求解器如ode15s。实操心得在数学建模竞赛或科研中优先使用ode45等内置求解器。自己写欧拉法虽然有助于理解原理但容易在步长选择上踩坑导致结果不准确或计算不稳定。把精力花在模型本身和参数分析上而不是数值积分的基础细节上。3.4 结果可视化与分析模拟完成后将结果可视化是分析和展示的关键。我们可以绘制各仓室人数随时间变化的曲线。% 绘制SEIRS模型模拟结果 figure(Position, [100, 100, 1200, 500]); % 设置图形窗口大小 % 子图1各仓室人数随时间变化 subplot(1, 2, 1); plot(t, S_ode, b-, LineWidth, 2, DisplayName, 易感者 S); hold on; plot(t, E_ode, m--, LineWidth, 1.5, DisplayName, 潜伏者 E); plot(t, I_ode, r-, LineWidth, 2, DisplayName, 感染者 I); plot(t, R_ode, g-., LineWidth, 1.5, DisplayName, 康复者 R); hold off; xlabel(时间 (天)); ylabel(人数); title(SEIRS模型动态模拟); legend(Location, best); grid on; grid minor; % 计算并标注关键点感染者峰值 [I_max, idx_max] max(I_ode); t_max t(idx_max); text(t_max, I_max, sprintf(峰值: %.0f人\n第%.1f天, I_max, t_max), ... VerticalAlignment, bottom, HorizontalAlignment, center, ... BackgroundColor, w, EdgeColor, r); % 子图2感染者比例流行曲线 subplot(1, 2, 2); plot(t, I_ode / N * 100, r-, LineWidth, 2); xlabel(时间 (天)); ylabel(感染者比例 (%)); title(疫情流行曲线感染者占比); grid on; grid minor; ylim([0, ceil(max(I_ode/N*100))]); % 设置y轴上限 % 在主图上添加模型参数信息 sgtitle(sprintf(SEIRS模型模拟 (N%d, \\beta%.2f, \\sigma%.3f, \\gamma%.3f, \\xi%.4f, R0%.2f), ... N, beta, sigma, gamma, xi, beta/gamma), FontSize, 12, FontWeight, bold);可视化技巧与解读多曲线对比将S、E、I、R画在同一张图上可以清晰看到它们之间的此消彼长和相位关系。通常易感者S单调下降或波动下降感染者I先上升后下降形成波峰康复者R单调上升或波动上升潜伏者E的曲线通常与I相似但略微超前。标注关键信息自动计算并标注感染者峰值I_max及其发生时间t_max能让读图者快速抓住疫情的核心特征。峰值高度和出现时间是评估疫情严重程度和速度的关键。流行曲线单独绘制感染者比例I/N的曲线这就是经典的“流行病曲线”。它的形状是单峰还是多峰、峰值的宽窄直接反映了传播动力学特征。参数标注使用sgtitle或text在图上标注模型参数特别是R0值使得图形包含完整信息便于报告和复现。运行上述代码你将得到一张清晰的图表直观展示疫情从开始、爆发、到达峰值再到消退甚至可能因免疫丧失而再次爆发的全过程。4. 模型参数的影响分析与情景模拟建好模型只是第一步更重要的是用它来回答问题“如果……会怎样”。通过调整参数我们可以模拟不同的公共卫生干预情景这是数学建模真正的威力所在。4.1 关键参数敏感性分析我们主要关注三个核心参数传播率β、康复率γ或感染期、免疫丧失率ξ。通过控制变量法一次只改变一个参数观察流行曲线的变化。% 参数敏感性分析比较不同传播率beta的影响 beta_values [0.3, 0.5, 0.7]; % 低、中、高传播率 gamma 0.143; % 固定康复率 sigma 0.2; xi 0.005; N 10000; I0 1; figure; hold on; colors lines(length(beta_values)); % 获取区分度高的颜色 for i 1:length(beta_values) beta beta_values(i); R0 beta / gamma; % 使用ode45求解 [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma, gamma, xi, N), [0, 300], [N-I0; 0; I0; 0]); I_sim y(:, 3); % 绘制感染者曲线 plot(t, I_sim, Color, colors(i, :), LineWidth, 2, ... DisplayName, sprintf(\\beta%.1f, R0%.2f, beta, R0)); end hold off; xlabel(时间 (天)); ylabel(感染者数量 I(t)); title(不同传播率(\beta)下的疫情发展); legend(Location, best); grid on;分析结果解读β (传播率) / R0这是影响最大的参数。如图所示R0越大β越大疫情峰值越高到达峰值的时间越早疫情发展越迅猛。当R0刚好大于1时疫情会缓慢发展当R0远大于1时疫情会呈爆炸式增长。这直观地说明了降低接触率如保持社交距离、戴口罩对于“拉平曲线”的重要性。γ (康复率)γ越大平均感染期1/γ越短。在其他条件不变时缩短感染期即γ增大会降低R0因为R0β/γ从而降低峰值和延缓疫情。这对应着通过有效治疗缩短患者传染期的措施。σ (潜伏期转化率)σ越大平均潜伏期1/σ越短。潜伏期变短会使潜伏者更快转为感染者导致疫情爆发更早、更集中但通常对最终感染总规模影响不大在SEIR模型中最终规模主要取决于R0。ξ (免疫丧失率)这是SEIRS区别于SEIR的关键。ξ越大免疫力丧失越快。当ξ0时疫情不会一次终结。在第一次流行波过后由于不断有康复者变回易感者当易感者积累到一定数量时就可能引发第二次、第三次流行波形成周期性爆发。模拟中延长模拟时间如2000天就能观察到这种周期现象。4.2 模拟公共卫生干预措施我们可以通过动态改变参数来模拟非药物干预措施。情景一在疫情中期实施“封控”降低β假设在第50天开始采取严格措施将传播率β从0.5降低到0.1持续60天后恢复。% 模拟动态干预中期降低传播率 beta0 0.5; % 初始传播率 beta_low 0.1; % 干预期间的传播率 intervention_start 50; intervention_duration 60; intervention_end intervention_start intervention_duration; % 定义随时间变化的beta函数 beta_func (t) beta0 * (t intervention_start) ... beta_low * (t intervention_start t intervention_end) ... beta0 * (t intervention_end); % 修改微分方程函数使beta成为时间的函数 odefun_dynamic (t, y) [ xi * y(4) - (beta_func(t) * y(3) / N) * y(1); % dS/dt (beta_func(t) * y(3) / N) * y(1) - sigma * y(2); % dE/dt sigma * y(2) - gamma * y(3); % dI/dt gamma * y(3) - xi * y(4) % dR/dt ]; % 求解 [t_dyn, y_dyn] ode45(odefun_dynamic, [0, 300], [N-I0; 0; I0; 0]); I_dyn y_dyn(:, 3); % 绘图并标注干预区间 figure; plot(t_dyn, I_dyn, b-, LineWidth, 2); xlabel(时间 (天)); ylabel(感染者数量 I(t)); title(模拟中期干预降低传播率); grid on; xline(intervention_start, r--, 干预开始, LabelVerticalAlignment, middle); xline(intervention_end, r--, 干预结束, LabelVerticalAlignment, middle); patch([intervention_start, intervention_end, intervention_end, intervention_start], ... [0, 0, max(ylim), max(ylim)], r, FaceAlpha, 0.1, EdgeColor, none); legend(感染者数量, Location, best);情景二提高检测与隔离率同时影响β和σ提高检测能力可以缩短从感染到被隔离的时间。这等效于同时减少了感染者的有效传播时间降低了β和缩短了潜伏期可能增加了σ因为更快被识别。我们可以创建一个更复杂的函数来模拟这种综合效果。模拟结果解读在“封控”情景的图中你会看到在干预区间内感染者数量的上升势头被迅速遏制曲线变得平缓甚至下降。干预结束后由于易感者比例仍然较高疫情可能会出现反弹。这个模拟清晰地展示了干预措施“拉平曲线、推迟峰值”的效果但也提示了如果措施不能持续或没有其他配合如疫苗接种疫情可能反复。注意事项在模拟动态干预时参数变化函数beta_func(t)在时间点intervention_start和intervention_end处是不连续的。ode45等求解器虽然能处理但可能会在间断点附近需要更多计算步骤。如果遇到收敛问题可以考虑使用ode15s或者将间断点作为事件明确告诉求解器。5. 模型校准、验证与常见问题排查一个模型如果无法用现实数据校准其预测能力就值得怀疑。此外在代码实现和模拟过程中你可能会遇到各种问题。5.1 如何利用真实数据校准模型参数校准的目标是找到一组模型参数使得模型的输出通常是每日新增感染数或累计感染数与观测数据最吻合。常用方法是最小二乘法或最大似然估计并配合优化算法如fminsearch寻找最优参数。假设我们有一组真实的每日新增感染数据I_data_new长度为n天的向量我们可以这样操作% 假设有真实数据这里用模拟数据加噪声代替 true_beta 0.55; true_gamma 0.14; [t_true, y_true] ode45((t,y) seirs_ode(t, y, true_beta, 0.2, true_gamma, 0.005, N), [0:1:100], [N-10; 0; 10; 0]); I_true_new -diff(y_true(:,1) y_true(:,2)); % 近似计算每日新增-(ΔSΔE) I_data_new I_true_new randn(size(I_true_new)) * 5; % 加入高斯噪声模拟真实数据 % 定义误差函数残差平方和 error_func (params) sum((simulate_new_cases(params) - I_data_new).^2); % 初始参数猜测 [beta, gamma] initial_guess [0.4, 0.1]; % 使用fminsearch寻找最小化误差的参数 options optimset(Display, iter, MaxFunEvals, 1000); estimated_params fminsearch(error_func, initial_guess, options); fprintf(真实参数: beta%.3f, gamma%.3f\n, true_beta, true_gamma); fprintf(估计参数: beta%.3f, gamma%.3f\n, estimated_params(1), estimated_params(2)); % 辅助函数给定参数模拟生成每日新增序列 function new_cases_sim simulate_new_cases(params) beta_est params(1); gamma_est params(2); sigma 0.2; xi 0.005; N 10000; [~, y_sim] ode45((t,y) seirs_ode(t, y, beta_est, sigma, gamma_est, xi, N), [0:1:100], [N-10; 0; 10; 0]); % 计算模拟的每日新增感染从易感者和潜伏者减少的总和 S_sim y_sim(:,1); E_sim y_sim(:,2); new_cases_sim -(diff(S_sim) diff(E_sim)); end校准心得与陷阱参数可识别性SEIRS模型参数较多可能存在“不同的参数组合产生相似的输出曲线”的情况即参数不可唯一识别。解决方法是利用更多数据如死亡数、血清学调查的抗体阳性率或引入先验知识如潜伏期、感染期来自医学文献来约束参数范围。数据质量真实数据存在报告延迟、检测能力变化、统计口径调整等问题。直接使用报告的“新增确诊”来拟合模型可能会产生偏差。通常需要对数据进行平滑、修正或使用更复杂的模型来刻画报告过程。优化算法选择fminsearch适用于参数较少、问题较简单的情况。对于更复杂的模型可能需要使用全局优化算法如particleswarm,ga来避免陷入局部最优。5.2 常见问题与调试技巧实录在实现和运行SEIRS模型时你可能会遇到以下典型问题问题1模拟结果中人数出现负值或总数不守恒。可能原因1时间步长dt太大欧拉法。欧拉法是条件稳定的如果dt大于系统最小时间常数的倒数解就会发散导致出现负值。解决显著减小dt例如从1改为0.1或0.01再试。或者直接改用ode45。可能原因2微分方程代码有误。检查导数dS_dt dE_dt dI_dt dR_dt是否等于0在μ0的简化模型中。如果不等于0总人口就不守恒。解决仔细核对方程代码确保正负号正确流入流出项平衡。可能原因3参数极端。例如beta值极大导致单步更新中S的减少量超过了S本身。解决检查参数取值的生物学合理性。beta通常与人口规模N有关确保beta*I/N是一个合理的概率小于1。问题2使用ode45时计算速度很慢或报错如“积分容差无法满足”。可能原因1模型是刚性的。当系统中不同状态变量的变化速率差异巨大时例如潜伏期很短σ很大而免疫期很长ξ很小就会产生刚性导致ode45需要极小的步长。解决换用适用于刚性问题的求解器如ode15s或ode23s。只需将ode45替换为ode15s即可。可能原因2参数或初始值导致函数值出现NaN或Inf。例如在计算(I/N)时如果初始N0会导致除零错误。解决在微分方程函数开头加入检查语句如if any(~isfinite(y))并设置合理的初始值。问题3模拟曲线与预期不符例如没有出现疫情高峰。可能原因1R0 1。检查你的beta和gamma参数计算R0 beta/gamma。如果R0小于或等于1疾病无法形成流行感染者数量会指数衰减。解决确保beta足够大使得R0 1。参考目标疾病的文献值设置R0。可能原因2初始感染者I0太少。在巨大的人口N中如果I01疫情起步会非常缓慢在有限的模拟时间内可能看不到明显高峰。解决可以适当增加I0例如设为总人口的0.1%或者延长模拟时间t_end。可能原因3免疫丧失率ξ太大。如果免疫力丧失太快易感者得到快速补充疫情可能表现为持续的低水平流行而非明显的单峰。解决根据目标疾病的免疫特性调整ξ。对于流感免疫期以月计对于新冠原始毒株免疫期可能以年计对于某些普通感冒冠状病毒免疫期可能只有几个月。问题4如何模拟疫苗接种疫苗接种可以看作是将易感者S直接转移到康复者R或部分免疫状态的过程。可以在微分方程中增加一项。例如假设以速率v为易感者接种疫苗接种后直接获得完全免疫进入R仓室dS/dt ... - v * S (增加一项) dR/dt ... v * S (增加一项)如果疫苗有效率不是100%或者免疫效果会衰减模型会变得更复杂可能需要引入新的仓室如部分免疫者。通过上述的代码实现、情景模拟和问题排查你应该已经能够驾驭这个SEIRS模型了。记住模型是对现实的简化它的价值不在于预测得多么精确而在于帮助我们理解系统内在的动力学机制并定量地比较不同干预策略的潜在效果。当你需要分析一种具有潜伏期和暂时免疫力的传染病时这个SEIRS模型及其Matlab实现就是一个非常有力的起点。 SEIRS传染病模型原理与Matlab实现：从微分方程到疫情模拟
1. 项目概述从SEIRS模型看疫情传播的数学本质最近几年大家对于传染病传播模型应该都不陌生了。无论是新闻里看到的疫情预测曲线还是各种政策推演的背后都离不开这些数学模型的支持。今天我想和大家深入聊聊的就是这个在经典SIR模型基础上演化而来的SEIRS模型并且会手把手地带你用Matlab把它“跑”起来。这个模型之所以重要是因为它更贴近现实中许多传染病的特性比如流感、新冠等它们感染后获得的免疫力并非永久人会再次变为易感者同时疾病还存在潜伏期并非一接触就立刻发病。SEIRS模型正是通过引入“潜伏者E”和“免疫丧失从R变回S”这两个关键状态来刻画这些复杂动态。简单来说SEIRS模型把人群分成了四类易感者S、潜伏者E、感染者I和康复者R。易感者接触到病毒后不会马上发病而是先进入潜伏期E经过一段时间才具有传染性I康复后R获得暂时免疫力但这份免疫力可能随时间减弱使人重新变回易感者S。整个系统就像一个动态流转的池子而我们的目标就是用数学方程描述这个流转过程并通过计算机模拟来观察不同参数如接触率、潜伏期时长、免疫持续时间下疫情的发展轨迹。这对于理解疫情走势、评估防控措施如提高隔离率相当于降低接触率的效果具有非常直观的参考价值。接下来我将彻底拆解这个模型的每一个环节从微分方程的原理推导到Matlab代码的逐行实现再到如何调整参数模拟不同场景。无论你是数学建模的初学者还是想寻找一个可靠、可复现的传染病模拟案例这篇文章都能给你提供一套完整的“工具箱”。2. SEIRS模型的核心原理与方程拆解要真正玩转一个模型不能只当它是“黑箱”知道输入输出就完事。我们必须深入其数学核心理解每一个方程、每一个参数背后的物理意义和生物学假设。只有这样当模拟结果出现异常时你才知道该调整哪里当需要针对特定疾病定制模型时你也知道如何修改方程。2.1 模型状态定义与流转关系SEIRS模型是一个典型的仓室模型它基于几个核心假设总人口数N恒定不考虑出生、死亡或假设出生率等于自然死亡率即 S E I R N 为常数。这个假设简化了模型让我们专注于疾病本身的传播动力学。均匀混合人群充分混合任何一个易感者S接触任何一个感染者I的概率是相同的。这显然是对现实的简化但对于理解宏观趋势非常有效。疾病进程确定个体从潜伏期E到发病期I再到康复R其转移速率是常数这意味着潜伏期和感染期近似服从指数分布。基于此我们定义四个核心状态变量它们都是关于时间t的函数S(t): 时刻 t 易感者的数量。E(t): 时刻 t 处于潜伏期已感染但未发病、无传染性者的数量。I(t): 时刻 t 感染者已发病、有传染性的数量。R(t): 时刻 t 康复者已康复、具有暂时免疫力的数量。它们之间的流转关系可以用下面的流程图来直观理解注意这里用文字描述代替图表 易感者S以一定的速率被感染但并非直接变成感染者I而是先进入潜伏期E。潜伏者E经过平均潜伏期后转化为感染者I。感染者I经过平均感染期后康复并进入康复者R群体。康复者R获得的免疫力并非永久经过平均免疫期后会重新变为易感者S。此外模型中通常还会考虑因病死亡率即感染者I可能以一定概率死亡这会使得总人口N发生变化但在最基本的SEIRS模型中我们常先忽略这一点或假设死亡率包含在康复率中。2.2 微分方程组的推导与参数解读上述的流转关系可以用一组常微分方程ODE来精确描述。这是整个模型的心脏。我们定义几个关键参数β (Beta): 有效接触率或传播率。它表示一个感染者I在单位时间内有效接触并传染易感者S的人数。在实际中β 接触率 × 每次接触的传染概率。它通常与总人口N有关更常见的写法是 β * I / N表示一个易感者遇到感染者的概率。σ (Sigma): 潜伏期转化率。它是平均潜伏期d_E的倒数即 σ 1 / d_E。例如平均潜伏期为5天则 σ 0.2 /天。表示单位时间内潜伏者E转化为感染者I的比例。γ (Gamma): 康复率。它是平均感染期d_I的倒数即 γ 1 / d_I。例如平均感染期为7天则 γ ≈ 0.143 /天。表示单位时间内感染者I康复进入R的比例。ξ (Xi): 免疫丧失率。它是平均免疫期d_R的倒数即 ξ 1 / d_R。例如免疫力平均维持180天则 ξ ≈ 0.0056 /天。表示单位时间内康复者R丧失免疫、重新变为易感者S的比例。μ (Mu): 自然死亡率或出生率在总人口恒定假设下两者相等。这是一个可选参数用于模拟更长期的动态平衡。现在我们可以写出经典的SEIRS微分方程组dS/dt μ * N ξ * R - β * (I/N) * S - μ * S dE/dt β * (I/N) * S - σ * E - μ * E dI/dt σ * E - γ * I - μ * I dR/dt γ * I - ξ * R - μ * R方程逐行解读易感者S的变化率 dS/dt μ * N假设有常数出生率或输入新增人口全部为易感者。 ξ * R康复者R丧失免疫力重新变为易感者S。- β * (I/N) * S这是核心的感染项。(I/N)是随机接触一人恰好是感染者的概率。β * (I/N)就是一个易感者单位时间内被感染的概率力。再乘以易感者总数S就得到了单位时间内从S流向E的人数。- μ * S易感者的自然死亡。潜伏者E的变化率 dE/dt β * (I/N) * S从易感者S被感染而来是E的“流入项”。- σ * E潜伏者E以速率σ转化为感染者I是E的“流出项”。- μ * E潜伏者的自然死亡。感染者I的变化率 dI/dt σ * E从潜伏者E转化而来是I的“流入项”。- γ * I感染者I以速率γ康复进入R群体是I的“流出项”。- μ * I感染者的自然死亡如果考虑疾病致死率这里可以改为- (γα)*I其中α是因病死亡率。康复者R的变化率 dR/dt γ * I从感染者I康复而来是R的“流入项”。- ξ * R康复者R丧失免疫力变回易感者S是R的“流出项”。- μ * R康复者的自然死亡。注意在许多短期疫情模拟中为了简化会忽略自然出生死亡即 μ0因为疫情周期几个月到一两年远小于人的寿命。此时方程简化为dS/dt ξR - βIS/N,dE/dt βIS/N - σE,dI/dt σE - γI,dR/dt γI - ξR。同时SEIR N恒定。我们后续的Matlab实现将基于这个简化版本因为它更聚焦于疾病传播本身。关键衍生参数——基本再生数 R0 R0是一个极其重要的流行病学指标代表在完全易感人群中一个感染者在其整个传染期内平均能传染的人数。对于SEIRS模型R0的计算公式为R0 β / γ这是因为一个感染者的平均传染期是1/γ在此期间他单位时间传染β个人在完全易感环境下S/N≈1。如果R0 1疾病会流行开来如果R0 1疾病会逐渐消失。在SEIRS模型中由于存在免疫丧失即使R01疫情也可能表现为周期性爆发而非一次流行后终结。3. 基于Matlab的模型实现与代码精讲理论说得再多不如一行代码。下面我们就用Matlab将上述方程组“复活”通过数值求解来观察疫情动态。我将采用最清晰、最易于理解和修改的方式编写代码并解释每一个关键步骤。3.1 模型参数设置与初始化首先我们需要定义模型参数和初始条件。参数的选择需要参考真实疾病的数据这里我们以类似流感的参数为例进行演示。% SEIRS模型参数设置 % 时间参数 t_start 0; % 模拟开始时间天 t_end 500; % 模拟结束时间天 dt 0.1; % 时间步长天用于欧拉法值越小越精确但计算越慢 % 流行病学参数 beta 0.5; % 传播率 (1/天) - 一个感染者每天有效接触人数 sigma 0.2; % 潜伏期转化率 (1/天) - 平均潜伏期 1/sigma 5天 gamma 0.143; % 康复率 (1/天) - 平均感染期 1/gamma ≈ 7天 xi 0.005; % 免疫丧失率 (1/天) - 平均免疫期 1/xi 200天 % 人口参数 N 10000; % 总人口 I0 1; % 初始感染者数量 E0 0; % 初始潜伏者数量 R0 0; % 初始康复者数量 S0 N - I0 - E0 - R0; % 初始易感者数量 % 初始化数组存储结果 num_steps ceil((t_end - t_start) / dt) 1; % 计算总步数 time linspace(t_start, t_end, num_steps); % 时间向量 S zeros(1, num_steps); E zeros(1, num_steps); I zeros(1, num_steps); R zeros(1, num_steps); % 设置初始值 S(1) S0; E(1) E0; I(1) I0; R(1) R0;参数设置心得时间步长dt这里我们使用了最简单的欧拉法进行数值求解。欧拉法是一种显式方法计算简单但稳定性有条件。对于这个模型只要dt设置得足够小比如小于模型中最小时间常数的倒数这里最小的时间常数可能是1/σ5天所以dt0.1是安全的结果就是可靠的。如果你想追求更高的精度和稳定性可以使用Matlab内置的ODE求解器如ode45这部分后面会讲。参数beta这是最敏感、最需要校准的参数。β R0 * γ。如果你从文献中知道了某种疾病的R0可以通过这个关系反推beta。例如流感的R0大约在1.2-1.8之间若γ1/7≈0.143则beta的范围约为0.17-0.26。这里设为0.5是为了让疫情发展更迅速便于观察。初始条件通常假设开始时只有一个或少数几个感染者I0其余均为易感者。潜伏者E0和康复者R0通常设为0。3.2 使用欧拉法进行数值求解欧拉法的核心思想是用当前时刻的导数来近似下一个时刻的值。对于微分方程dy/dt f(y, t)其迭代公式为y_{n1} y_n f(y_n, t_n) * dt。% 使用欧拉法迭代求解SEIRS模型 for i 1:num_steps-1 % 当前时刻的各仓室人数 S_curr S(i); E_curr E(i); I_curr I(i); R_curr R(i); % 计算各仓室的变化率导数基于简化模型μ0 dS_dt xi * R_curr - (beta * I_curr / N) * S_curr; dE_dt (beta * I_curr / N) * S_curr - sigma * E_curr; dI_dt sigma * E_curr - gamma * I_curr; dR_dt gamma * I_curr - xi * R_curr; % 欧拉法更新下一时刻的人数 S(i1) S_curr dS_dt * dt; E(i1) E_curr dE_dt * dt; I(i1) I_curr dI_dt * dt; R(i1) R_curr dR_dt * dt; % 可选施加约束防止数值误差导致人数为负虽然欧拉法在dt很小时很少出现 S(i1) max(0, S(i1)); E(i1) max(0, E(i1)); I(i1) max(0, I(i1)); R(i1) max(0, R(i1)); end欧拉法实现的注意事项计算顺序在计算导数dS_dt时我们使用的是当前时刻i的S_curr,I_curr等值。这意味着所有状态是同步更新的。如果错误地先更新了S再用新的S去计算E的导数就会引入误差破坏模型的一致性。负值处理由于是离散近似在极端参数或大步长下可能会出现某个仓室人数为负的荒谬情况。虽然概率很低但加上max(0, ...)约束是一个良好的编程习惯能保证模拟的鲁棒性。更根本的解决方法是使用更小的时间步长dt或改用更稳定的求解器。精度验证一个简单的验证方法是将dt减半如从0.1变为0.05重新运行模拟。如果两条曲线几乎重合说明当前的dt已经足够精确。如果差异明显则需要继续减小dt。3.3 使用Matlab内置ODE求解器推荐对于常微分方程组Matlab提供了强大且稳定的求解器如ode45基于Runge-Kutta方法。使用它们不仅精度高而且无需手动设置步长代码也更简洁。我强烈推荐在实际建模中使用这种方法。% 使用ode45求解SEIRS模型 % 定义时间跨度 tspan [0, 500]; % 定义初始状态向量 [S0; E0; I0; R0] y0 [S0; E0; I0; R0]; % 调用ode45求解器 % seirs_ode 是定义微分方程组的函数句柄 % tspan 是时间范围 % y0 是初始条件 [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma, gamma, xi, N), tspan, y0); % 提取结果 S_ode y(:, 1); E_ode y(:, 2); I_ode y(:, 3); R_ode y(:, 4);这里我们需要单独定义一个函数文件seirs_ode.m来描述微分方程组function dydt seirs_ode(t, y, beta, sigma, gamma, xi, N) % SEIRS模型微分方程组 % 输入 % t: 时间未显式使用但ode45要求此参数 % y: 状态向量 [S; E; I; R] % beta, sigma, gamma, xi, N: 模型参数 % 输出 % dydt: 导数向量 [dS/dt; dE/dt; dI/dt; dR/dt] S y(1); E y(2); I y(3); R y(4); % 计算各导数 dS_dt xi * R - (beta * I / N) * S; dE_dt (beta * I / N) * S - sigma * E; dI_dt sigma * E - gamma * I; dR_dt gamma * I - xi * R; dydt [dS_dt; dE_dt; dI_dt; dR_dt]; end使用ODE求解器的优势自动变步长ode45会根据方程的特性自动调整积分步长在变化平缓处用大步长提高效率在变化剧烈处用小步长保证精度。高精度采用高阶Runge-Kutta方法数值误差远小于简单的欧拉法。代码简洁省去了手动迭代和步长管理的麻烦。稳定性好对于刚性方程系统中存在快变和慢变过程可以使用专门的求解器如ode15s。实操心得在数学建模竞赛或科研中优先使用ode45等内置求解器。自己写欧拉法虽然有助于理解原理但容易在步长选择上踩坑导致结果不准确或计算不稳定。把精力花在模型本身和参数分析上而不是数值积分的基础细节上。3.4 结果可视化与分析模拟完成后将结果可视化是分析和展示的关键。我们可以绘制各仓室人数随时间变化的曲线。% 绘制SEIRS模型模拟结果 figure(Position, [100, 100, 1200, 500]); % 设置图形窗口大小 % 子图1各仓室人数随时间变化 subplot(1, 2, 1); plot(t, S_ode, b-, LineWidth, 2, DisplayName, 易感者 S); hold on; plot(t, E_ode, m--, LineWidth, 1.5, DisplayName, 潜伏者 E); plot(t, I_ode, r-, LineWidth, 2, DisplayName, 感染者 I); plot(t, R_ode, g-., LineWidth, 1.5, DisplayName, 康复者 R); hold off; xlabel(时间 (天)); ylabel(人数); title(SEIRS模型动态模拟); legend(Location, best); grid on; grid minor; % 计算并标注关键点感染者峰值 [I_max, idx_max] max(I_ode); t_max t(idx_max); text(t_max, I_max, sprintf(峰值: %.0f人\n第%.1f天, I_max, t_max), ... VerticalAlignment, bottom, HorizontalAlignment, center, ... BackgroundColor, w, EdgeColor, r); % 子图2感染者比例流行曲线 subplot(1, 2, 2); plot(t, I_ode / N * 100, r-, LineWidth, 2); xlabel(时间 (天)); ylabel(感染者比例 (%)); title(疫情流行曲线感染者占比); grid on; grid minor; ylim([0, ceil(max(I_ode/N*100))]); % 设置y轴上限 % 在主图上添加模型参数信息 sgtitle(sprintf(SEIRS模型模拟 (N%d, \\beta%.2f, \\sigma%.3f, \\gamma%.3f, \\xi%.4f, R0%.2f), ... N, beta, sigma, gamma, xi, beta/gamma), FontSize, 12, FontWeight, bold);可视化技巧与解读多曲线对比将S、E、I、R画在同一张图上可以清晰看到它们之间的此消彼长和相位关系。通常易感者S单调下降或波动下降感染者I先上升后下降形成波峰康复者R单调上升或波动上升潜伏者E的曲线通常与I相似但略微超前。标注关键信息自动计算并标注感染者峰值I_max及其发生时间t_max能让读图者快速抓住疫情的核心特征。峰值高度和出现时间是评估疫情严重程度和速度的关键。流行曲线单独绘制感染者比例I/N的曲线这就是经典的“流行病曲线”。它的形状是单峰还是多峰、峰值的宽窄直接反映了传播动力学特征。参数标注使用sgtitle或text在图上标注模型参数特别是R0值使得图形包含完整信息便于报告和复现。运行上述代码你将得到一张清晰的图表直观展示疫情从开始、爆发、到达峰值再到消退甚至可能因免疫丧失而再次爆发的全过程。4. 模型参数的影响分析与情景模拟建好模型只是第一步更重要的是用它来回答问题“如果……会怎样”。通过调整参数我们可以模拟不同的公共卫生干预情景这是数学建模真正的威力所在。4.1 关键参数敏感性分析我们主要关注三个核心参数传播率β、康复率γ或感染期、免疫丧失率ξ。通过控制变量法一次只改变一个参数观察流行曲线的变化。% 参数敏感性分析比较不同传播率beta的影响 beta_values [0.3, 0.5, 0.7]; % 低、中、高传播率 gamma 0.143; % 固定康复率 sigma 0.2; xi 0.005; N 10000; I0 1; figure; hold on; colors lines(length(beta_values)); % 获取区分度高的颜色 for i 1:length(beta_values) beta beta_values(i); R0 beta / gamma; % 使用ode45求解 [t, y] ode45((t,y) seirs_ode(t, y, beta, sigma, gamma, xi, N), [0, 300], [N-I0; 0; I0; 0]); I_sim y(:, 3); % 绘制感染者曲线 plot(t, I_sim, Color, colors(i, :), LineWidth, 2, ... DisplayName, sprintf(\\beta%.1f, R0%.2f, beta, R0)); end hold off; xlabel(时间 (天)); ylabel(感染者数量 I(t)); title(不同传播率(\beta)下的疫情发展); legend(Location, best); grid on;分析结果解读β (传播率) / R0这是影响最大的参数。如图所示R0越大β越大疫情峰值越高到达峰值的时间越早疫情发展越迅猛。当R0刚好大于1时疫情会缓慢发展当R0远大于1时疫情会呈爆炸式增长。这直观地说明了降低接触率如保持社交距离、戴口罩对于“拉平曲线”的重要性。γ (康复率)γ越大平均感染期1/γ越短。在其他条件不变时缩短感染期即γ增大会降低R0因为R0β/γ从而降低峰值和延缓疫情。这对应着通过有效治疗缩短患者传染期的措施。σ (潜伏期转化率)σ越大平均潜伏期1/σ越短。潜伏期变短会使潜伏者更快转为感染者导致疫情爆发更早、更集中但通常对最终感染总规模影响不大在SEIR模型中最终规模主要取决于R0。ξ (免疫丧失率)这是SEIRS区别于SEIR的关键。ξ越大免疫力丧失越快。当ξ0时疫情不会一次终结。在第一次流行波过后由于不断有康复者变回易感者当易感者积累到一定数量时就可能引发第二次、第三次流行波形成周期性爆发。模拟中延长模拟时间如2000天就能观察到这种周期现象。4.2 模拟公共卫生干预措施我们可以通过动态改变参数来模拟非药物干预措施。情景一在疫情中期实施“封控”降低β假设在第50天开始采取严格措施将传播率β从0.5降低到0.1持续60天后恢复。% 模拟动态干预中期降低传播率 beta0 0.5; % 初始传播率 beta_low 0.1; % 干预期间的传播率 intervention_start 50; intervention_duration 60; intervention_end intervention_start intervention_duration; % 定义随时间变化的beta函数 beta_func (t) beta0 * (t intervention_start) ... beta_low * (t intervention_start t intervention_end) ... beta0 * (t intervention_end); % 修改微分方程函数使beta成为时间的函数 odefun_dynamic (t, y) [ xi * y(4) - (beta_func(t) * y(3) / N) * y(1); % dS/dt (beta_func(t) * y(3) / N) * y(1) - sigma * y(2); % dE/dt sigma * y(2) - gamma * y(3); % dI/dt gamma * y(3) - xi * y(4) % dR/dt ]; % 求解 [t_dyn, y_dyn] ode45(odefun_dynamic, [0, 300], [N-I0; 0; I0; 0]); I_dyn y_dyn(:, 3); % 绘图并标注干预区间 figure; plot(t_dyn, I_dyn, b-, LineWidth, 2); xlabel(时间 (天)); ylabel(感染者数量 I(t)); title(模拟中期干预降低传播率); grid on; xline(intervention_start, r--, 干预开始, LabelVerticalAlignment, middle); xline(intervention_end, r--, 干预结束, LabelVerticalAlignment, middle); patch([intervention_start, intervention_end, intervention_end, intervention_start], ... [0, 0, max(ylim), max(ylim)], r, FaceAlpha, 0.1, EdgeColor, none); legend(感染者数量, Location, best);情景二提高检测与隔离率同时影响β和σ提高检测能力可以缩短从感染到被隔离的时间。这等效于同时减少了感染者的有效传播时间降低了β和缩短了潜伏期可能增加了σ因为更快被识别。我们可以创建一个更复杂的函数来模拟这种综合效果。模拟结果解读在“封控”情景的图中你会看到在干预区间内感染者数量的上升势头被迅速遏制曲线变得平缓甚至下降。干预结束后由于易感者比例仍然较高疫情可能会出现反弹。这个模拟清晰地展示了干预措施“拉平曲线、推迟峰值”的效果但也提示了如果措施不能持续或没有其他配合如疫苗接种疫情可能反复。注意事项在模拟动态干预时参数变化函数beta_func(t)在时间点intervention_start和intervention_end处是不连续的。ode45等求解器虽然能处理但可能会在间断点附近需要更多计算步骤。如果遇到收敛问题可以考虑使用ode15s或者将间断点作为事件明确告诉求解器。5. 模型校准、验证与常见问题排查一个模型如果无法用现实数据校准其预测能力就值得怀疑。此外在代码实现和模拟过程中你可能会遇到各种问题。5.1 如何利用真实数据校准模型参数校准的目标是找到一组模型参数使得模型的输出通常是每日新增感染数或累计感染数与观测数据最吻合。常用方法是最小二乘法或最大似然估计并配合优化算法如fminsearch寻找最优参数。假设我们有一组真实的每日新增感染数据I_data_new长度为n天的向量我们可以这样操作% 假设有真实数据这里用模拟数据加噪声代替 true_beta 0.55; true_gamma 0.14; [t_true, y_true] ode45((t,y) seirs_ode(t, y, true_beta, 0.2, true_gamma, 0.005, N), [0:1:100], [N-10; 0; 10; 0]); I_true_new -diff(y_true(:,1) y_true(:,2)); % 近似计算每日新增-(ΔSΔE) I_data_new I_true_new randn(size(I_true_new)) * 5; % 加入高斯噪声模拟真实数据 % 定义误差函数残差平方和 error_func (params) sum((simulate_new_cases(params) - I_data_new).^2); % 初始参数猜测 [beta, gamma] initial_guess [0.4, 0.1]; % 使用fminsearch寻找最小化误差的参数 options optimset(Display, iter, MaxFunEvals, 1000); estimated_params fminsearch(error_func, initial_guess, options); fprintf(真实参数: beta%.3f, gamma%.3f\n, true_beta, true_gamma); fprintf(估计参数: beta%.3f, gamma%.3f\n, estimated_params(1), estimated_params(2)); % 辅助函数给定参数模拟生成每日新增序列 function new_cases_sim simulate_new_cases(params) beta_est params(1); gamma_est params(2); sigma 0.2; xi 0.005; N 10000; [~, y_sim] ode45((t,y) seirs_ode(t, y, beta_est, sigma, gamma_est, xi, N), [0:1:100], [N-10; 0; 10; 0]); % 计算模拟的每日新增感染从易感者和潜伏者减少的总和 S_sim y_sim(:,1); E_sim y_sim(:,2); new_cases_sim -(diff(S_sim) diff(E_sim)); end校准心得与陷阱参数可识别性SEIRS模型参数较多可能存在“不同的参数组合产生相似的输出曲线”的情况即参数不可唯一识别。解决方法是利用更多数据如死亡数、血清学调查的抗体阳性率或引入先验知识如潜伏期、感染期来自医学文献来约束参数范围。数据质量真实数据存在报告延迟、检测能力变化、统计口径调整等问题。直接使用报告的“新增确诊”来拟合模型可能会产生偏差。通常需要对数据进行平滑、修正或使用更复杂的模型来刻画报告过程。优化算法选择fminsearch适用于参数较少、问题较简单的情况。对于更复杂的模型可能需要使用全局优化算法如particleswarm,ga来避免陷入局部最优。5.2 常见问题与调试技巧实录在实现和运行SEIRS模型时你可能会遇到以下典型问题问题1模拟结果中人数出现负值或总数不守恒。可能原因1时间步长dt太大欧拉法。欧拉法是条件稳定的如果dt大于系统最小时间常数的倒数解就会发散导致出现负值。解决显著减小dt例如从1改为0.1或0.01再试。或者直接改用ode45。可能原因2微分方程代码有误。检查导数dS_dt dE_dt dI_dt dR_dt是否等于0在μ0的简化模型中。如果不等于0总人口就不守恒。解决仔细核对方程代码确保正负号正确流入流出项平衡。可能原因3参数极端。例如beta值极大导致单步更新中S的减少量超过了S本身。解决检查参数取值的生物学合理性。beta通常与人口规模N有关确保beta*I/N是一个合理的概率小于1。问题2使用ode45时计算速度很慢或报错如“积分容差无法满足”。可能原因1模型是刚性的。当系统中不同状态变量的变化速率差异巨大时例如潜伏期很短σ很大而免疫期很长ξ很小就会产生刚性导致ode45需要极小的步长。解决换用适用于刚性问题的求解器如ode15s或ode23s。只需将ode45替换为ode15s即可。可能原因2参数或初始值导致函数值出现NaN或Inf。例如在计算(I/N)时如果初始N0会导致除零错误。解决在微分方程函数开头加入检查语句如if any(~isfinite(y))并设置合理的初始值。问题3模拟曲线与预期不符例如没有出现疫情高峰。可能原因1R0 1。检查你的beta和gamma参数计算R0 beta/gamma。如果R0小于或等于1疾病无法形成流行感染者数量会指数衰减。解决确保beta足够大使得R0 1。参考目标疾病的文献值设置R0。可能原因2初始感染者I0太少。在巨大的人口N中如果I01疫情起步会非常缓慢在有限的模拟时间内可能看不到明显高峰。解决可以适当增加I0例如设为总人口的0.1%或者延长模拟时间t_end。可能原因3免疫丧失率ξ太大。如果免疫力丧失太快易感者得到快速补充疫情可能表现为持续的低水平流行而非明显的单峰。解决根据目标疾病的免疫特性调整ξ。对于流感免疫期以月计对于新冠原始毒株免疫期可能以年计对于某些普通感冒冠状病毒免疫期可能只有几个月。问题4如何模拟疫苗接种疫苗接种可以看作是将易感者S直接转移到康复者R或部分免疫状态的过程。可以在微分方程中增加一项。例如假设以速率v为易感者接种疫苗接种后直接获得完全免疫进入R仓室dS/dt ... - v * S (增加一项) dR/dt ... v * S (增加一项)如果疫苗有效率不是100%或者免疫效果会衰减模型会变得更复杂可能需要引入新的仓室如部分免疫者。通过上述的代码实现、情景模拟和问题排查你应该已经能够驾驭这个SEIRS模型了。记住模型是对现实的简化它的价值不在于预测得多么精确而在于帮助我们理解系统内在的动力学机制并定量地比较不同干预策略的潜在效果。当你需要分析一种具有潜伏期和暂时免疫力的传染病时这个SEIRS模型及其Matlab实现就是一个非常有力的起点。