尧图网站建设 
简介本资源是一套面向生物信息学、数学建模与公共卫生研究初学者的MATLAB传染病动力学仿真源码集聚焦SI、SIS、SIR三类经典微分方程模型解决传染病传播机制理解、参数敏感性分析及防控策略模拟等核心问题。压缩包共含多个.m主程序文件如si_sim.m、sis_sim.m、sir_sim.m及配套参数配置脚本全部为可直接运行的MATLAB函数与ODE求解脚本用于构建微分方程组、调用ode45数值求解器、生成S/I/R时间演化曲线并实现动态可视化包体大小仅75KB轻量易部署。已有4784人学习下载适合高校数学建模课程实践、流行病学仿真实验及科研入门训练。读者可直接复现三种模型的完整建模流程——从方程定义、参数设置、数值求解到结果绘图并通过对比不同β感染率与γ恢复率下的传播曲线深入掌握模型差异与现实映射逻辑。1. 项目概述从代码集锦到建模实战手头有一套别人分享的“MATLAB源码集锦”里面打包了传染病SI、SIS、SIR三种经典模型的代码这大概是很多数学建模新手入门时都会遇到的“宝藏资源”。但说实话我见过太多人拿到这种源码包后只是机械地运行一下看到几条曲线跳出来就觉得“哦我懂了”然后就把代码扔在一边。这其实浪费了这些代码最大的价值——它们不仅是拿来“跑”的更是用来“拆”和“改”的是理解传染病动力学建模思想最直接的脚手架。SI、SIS、SIR这三个模型可以说是流行病学数学建模的“三原色”。它们用极其简洁的微分方程刻画了人群在疾病传播过程中的状态流转。SI模型描述的是像某些无法治愈的慢性感染或早期理论中的理想情况SIS模型则加入了康复后可能再次被感染的特性比如普通感冒而SIR模型引入了“移出者”的概念刻画了那些感染后获得持久免疫力或死亡而不再参与传播的过程例如麻疹、天花。掌握它们你不仅是在学几个公式更是在掌握一种用数学语言量化社会现象、预测事态发展的思维方式。这套源码的价值就在于它把抽象的方程变成了可视化的动态过程让你能亲手调整参数直观感受“传播率”、“恢复率”这些关键因子如何左右一场“虚拟疫情”的走向。接下来我不会仅仅带大家浏览一遍代码。我们将以这套源码为蓝本深入这三个模型的数学核心并手把手教你如何从“会用”进阶到“会改”。你会学到如何将代码结构拆解成清晰的模块如何根据实际问题调整模型假设比如考虑隔离措施、疫苗接种以及如何对输出结果进行专业的分析和可视化呈现。无论你是正在备战数学建模竞赛的学生还是对复杂系统建模感兴趣的爱好者这篇内容都将帮你把这一套“源码集锦”真正转化为你工具箱里一件趁手的兵器。2. 模型核心思想与数学原理拆解在打开MATLAB代码之前我们必须先吃透模型背后的数学思想。代码只是方程的解算器和画笔方程本身才是灵魂。这三个模型都基于一个核心框架仓室模型。我们把总人口N划分为几个互斥的“仓室”个体在不同仓室间流动流动的速率由微分方程描述。2.1 SI模型最简单的传播逻辑SI模型是所有传染病模型的起点。它假设人群只分为两类易感者和感染者。一旦感染终身患病且具有传染性不会康复也不会死亡或者说模型不关心感染后的其他结局。这听起来很极端但它适用于描述某些无法治愈的慢性病在封闭群体中早期的传播动态或者在理论研究中作为基准。它的微分方程极其简洁易感者 dS/dt -β * S * I / N感染者 dI/dt β * S * I / N这里S是易感者数量I是感染者数量N是总人口恒定。β是整个模型的关键称为有效接触率或传播系数。方程 dS/dt -β * S * I / N 是理解所有后续模型的基础单位时间内新感染的人数正比于当前易感者数量S、当前感染者数量I以及他们之间发生有效接触的概率由β/N体现。这个“有效接触”是指足以导致疾病传播的接触。由于总人口N不变且S I N所以实际上只需要解一个关于I的方程即可。注意很多初学者会混淆β和另一个常见参数“基本再生数R0”。在SI模型中由于没有康复所有感染者最终都会感染所有人R0的概念不那么突出但β直接决定了疫情发展的速度。β越大曲线上升得越陡峭。2.2 SIS模型考虑康复与再感染SIS模型在SI的基础上增加了一点现实性感染者可以康复。但康复后他们并没有获得免疫力而是立刻重新变为易感者。这就形成了一个S - I - S的循环。很多细菌性感染如淋病、普通感冒针对特定毒株可以用SIS模型近似描述。它的微分方程变为易感者 dS/dt -β * S * I / N γ * I感染者 dI/dt β * S * I / N - γ * I这里引入了一个新参数γ称为恢复率或移出率。它的倒数 1/γ 可以理解为平均感染期。例如γ0.2/天意味着平均感染期是5天。方程中-γ * I 项表示感染者以速率γ康复而γ * I 项则表示这些康复者重新加入了易感者行列。SIS模型的一个重要特征是可能存在一个地方病平衡点。当传播和恢复达到平衡时感染人数会稳定在一个非零的水平而不是像SI模型那样感染所有人。这个平衡点取决于 β 和 γ 的比值。这里就引出了流行病学中至关重要的概念——基本再生数 R0。在SIS模型中R0 β / γ。当 R0 1 时疾病会流行并最终达到地方病平衡当 R0 1 时疾病会逐渐消失。2.3 SIR模型经典的“免疫者”模型SIR模型是三个模型中最著名、应用最广的。它在SIS的基础上增加了一个新的仓室移出者。个体从感染者状态移出后不再回到易感者状态而是进入R仓室。R可以代表获得持久免疫力的人也可以代表因病死亡的人——总之他们不再参与疾病的传播过程。麻疹、天花、水痘等疾病符合这个特征。它的微分方程系统是易感者 dS/dt -β * S * I / N感染者 dI/dt β * S * I / N - γ * I移出者 dR/dt γ * I注意dS/dt的方程和SI模型一样因为易感者只会减少。dI/dt的方程和SIS模型一样但移出项γ*I不再流回S而是流向了R。总人口 N S I R 依然守恒。SIR模型的行为非常经典疫情从少数感染者开始易感者数量减少感染者数量先上升后下降最终所有感染者都移出疫情结束留下部分从未感染的易感者和大量具有免疫力的移出者。最终未感染的人口比例取决于R0。基本再生数 R0 β / γ在这里同样至关重要。它表示在一个全部为易感者的人群中一个感染者在其整个传染期内平均能感染的人数。R0 1是疾病能够流行的阈值。实操心得理解参数的单位是避免错误的关键。β的单位通常是“每人每天”表示一个感染者每天能有效感染的人数在完全易感人群中。γ的单位是“每天”表示每天康复的比例。在编程时确保时间步长如ode45求解器中的时间向量与参数的单位通常是“天”匹配否则模拟出的疫情速度会失真。3. MATLAB源码结构解析与关键函数实现拿到一套打包好的源码最忌讳的就是直接点“运行”。我们先要像解剖一样看清它的代码结构、数据流和关键函数。一套良好的SI/SIS/SIR源码集锦通常会包含以下几个部分主脚本、模型方程函数、参数设置区、求解调用和绘图部分。我们来逐一拆解。3.1 主脚本框架与参数初始化一个典型的main.m或run_models.m脚本会像下面这样组织。它的核心作用是设置环境、定义参数、调用不同的模型函数并展示结果。% 清空环境关闭所有图形窗口确保干净的开始 clear; close all; clc; % 参数设置区域 % 这是你需要反复修改和实验的地方 N 1000; % 总人口 I0 1; % 初始感染者人数 S0 N - I0; % 初始易感者人数 R0 0; % 初始移出者人数 (SIR模型用) % 流行病学参数 beta 0.3; % 传播系数 (每人每天) gamma 0.1; % 恢复率 (每天) % 计算基本再生数 R0 R0_basic beta / gamma; fprintf(基本再生数 R0 %.2f\n, R0_basic); % 时间设置 t_start 0; t_end 200; % 模拟总时长 (天) tspan [t_start, t_end]; % 时间区间供ODE求解器使用 % 调用不同模型进行模拟 % 通常会有三个独立的代码块或子函数调用 % 1. 模拟SI模型 % 2. 模拟SIS模型 % 3. 模拟SIR模型 % 绘制结果对比图 % 将三个模型的结果绘制在一张或几张图中进行对比注意事项参数初始化部分至关重要。beta和gamma的值需要你根据实际疾病的文献或题目要求进行设定。例如流感的R0大约在1.2-1.8之间麻疹则高达12-18。通过调整这些参数你可以模拟不同传染力的疾病。初始感染者I0通常设为1但也可以研究输入性病例I0稍大的影响。3.2 核心微分方程组的ODE函数定义模型的灵魂在于微分方程组。在MATLAB中我们使用一个独立的函数文件例如odefun_SIR.m来定义这些方程供ode45这类求解器调用。这是将数学公式转化为代码的关键一步。以SIR模型为例我们创建一个函数文件odefun_SIR.mfunction dydt odefun_SIR(t, y, beta, gamma, N) % SIR模型的微分方程函数 % 输入: % t: 时间 (未直接使用但ODE求解器要求此参数) % y: 状态向量 [S; I; R] % beta, gamma, N: 参数 % 输出: % dydt: 导数向量 [dS/dt; dI/dt; dR/dt] S y(1); I y(2); R y(3); % 虽然方程中用不到dR/dt里的R但这里解出以保持向量完整 dS_dt -beta * S * I / N; dI_dt beta * S * I / N - gamma * I; dR_dt gamma * I; dydt [dS_dt; dI_dt; dR_dt]; end对于SIS模型odefun_SIS.mfunction dydt odefun_SIS(t, y, beta, gamma, N) % SIS模型的微分方程函数 % 状态向量 y [S; I] S y(1); I y(2); dS_dt -beta * S * I / N gamma * I; dI_dt beta * S * I / N - gamma * I; dydt [dS_dt; dI_dt]; end对于SI模型就更简单了由于SIN我们通常只对I建立方程或者依然用两个变量但令gamma0。关键细节函数头function dydt odefun_SIR(t, y, beta, gamma, N)中的参数传递方式。我们通过匿名函数或在调用ode45时传递额外参数。推荐使用匿名函数方式因为它更清晰ode45((t,y) odefun_SIR(t,y,beta,gamma,N), tspan, y0)这样betagammaN这些参数就被“绑定”到了微分方程函数中。3.3 使用ODE求解器进行数值积分定义了微分方程我们需要用数值方法求解它。MATLAB的ode45龙格-库塔法是解决非刚性常微分方程的首选对于这些传染病模型完全够用。在主脚本中调用求解器的代码块如下% --- SIR模型求解 --- y0_SIR [S0; I0; R0]; % 初始条件向量 % 使用匿名函数传递参数 [t_SIR, y_SIR] ode45((t,y) odefun_SIR(t,y,beta,gamma,N), tspan, y0_SIR); % 提取结果 S_SIR y_SIR(:, 1); I_SIR y_SIR(:, 2); R_SIR y_SIR(:, 3); % --- SIS模型求解 --- y0_SIS [S0; I0]; [t_SIS, y_SIS] ode45((t,y) odefun_SIS(t,y,beta,gamma,N), tspan, y0_SIS); S_SIS y_SIS(:, 1); I_SIS y_SIS(:, 2); % --- SI模型求解 --- % SI模型可以看作gamma0的SIR或SIS特例这里用解析解或简单ODE求解 % 方法1使用ODE求解与SIR类似但gamma0 beta_SI 0.5; % SI模型的传播率可能设得不同 gamma_SI 0; y0_SI [S0; I0]; [t_SI, y_SI] ode45((t,y) odefun_SIR(t,y,beta_SI,gamma_SI,N), tspan, y0_SI); S_SI y_SI(:, 1); I_SI y_SI(:, 2);ode45返回两个数组t是时间点向量y是对应时间点的状态向量每一列是一个状态变量。我们需要把这些数据提取出来以便绘图和分析。实操心得ode45是自适应步长的它返回的时间点t不是均匀的。如果你需要均匀时间序列上的结果比如每天的数据可以在tspan中指定更多输出点例如tspan linspace(0, 200, 201)这样就会输出0,1,2,...,200共201个时间点的解。但注意这并不改变求解精度只是增加了输出密度。3.4 结果可视化与对比分析绘图是将数学结果直观化的关键。一个好的对比图能清晰展示三个模型的本质差异。figure(Position, [100, 100, 1200, 400]); % 设置图形窗口位置和大小 % 子图1三类人群数量随时间变化以SIR为例 subplot(1, 3, 1); plot(t_SIR, S_SIR, b-, LineWidth, 2); hold on; plot(t_SIR, I_SIR, r-, LineWidth, 2); plot(t_SIR, R_SIR, g-, LineWidth, 2); hold off; xlabel(时间 (天)); ylabel(人口数量); title(SIR模型动态); legend(易感者 S, 感染者 I, 移出者 R, Location, best); grid on; % 子图2感染者比例对比三个模型 subplot(1, 3, 2); plot(t_SI, I_SI/N, k--, LineWidth, 1.5); hold on; plot(t_SIS, I_SIS/N, m-., LineWidth, 1.5); plot(t_SIR, I_SIR/N, r-, LineWidth, 2); hold off; xlabel(时间 (天)); ylabel(感染者比例 (I/N)); title(感染者比例对比); legend(SI模型, SIS模型, SIR模型, Location, best); grid on; % 添加R0信息标注 text(0.05*t_end, 0.9, sprintf(R_0 %.2f, R0_basic), Units, normalized, FontSize, 10); % 子图3相平面图SIR模型展示S-I关系 subplot(1, 3, 3); plot(S_SIR/N, I_SIR/N, b-, LineWidth, 1.5); xlabel(易感者比例 (S/N)); ylabel(感染者比例 (I/N)); title(SIR模型相平面 (S-I)); grid on; % 添加阈值线 S 1/R0 if R0_basic 0 threshold 1/R0_basic; line([threshold, threshold], ylim, Color, r, LineStyle, --, LineWidth, 1); text(threshold0.02, 0.05, S1/R_0, Color, r); end这段代码生成一个包含三个子图的综合视图。第一个子图展示SIR模型中三类人群的完整动态第二个子图直接对比三个模型中感染者比例的变化这是最核心的观察第三个子图是SIR模型的相平面图它展示了感染者比例随易感者比例变化的轨迹那条红色的阈值线S 1/R0是理论上的“疫情峰值线”当易感者比例降到这条线以下时感染者比例开始下降。避坑技巧绘图时将数量转化为比例除以总人口N常常更有意义因为它消除了人口规模的影响便于不同规模人群间的比较。另外使用hold on和hold off管理图形叠加使用legend和xlabel等函数完善图表信息这些都是生成专业图表的基本功。在竞赛或报告中清晰的图表能极大提升作品质量。4. 从理解到应用模型拓展与实战技巧掌握了基础模型的实现我们就可以开始“魔改”了。数学建模的魅力在于根据实际问题调整模型。下面介绍几个常见的拓展方向和实现技巧。4.1 引入出生与死亡SEIR模型雏形基础SIR模型假设人口恒定但长期模拟中出生和自然死亡是必须考虑的。我们可以引入一个常数出生率Λ和自然死亡率μ。此时微分方程变为 dS/dt Λ - β S I / N - μ S dI/dt β S I / N - γ I - μ I dR/dt γ I - μ R 总人口N不再恒定但通常假设出生率等于死亡率Λ μ N来保持人口稳定。在MATLAB中实现只需修改对应的ODE函数即可。这种模型更适用于研究疾病的地方性流行。4.2 加入潜伏期SEIR模型很多传染病如流感、COVID-19感染后不会立即具有传染性而是有一段潜伏期。这时需要在S和I之间增加一个暴露者仓室。这就是SEIR模型 dS/dt -β S I / N dE/dt β S I / N - σ E dI/dt σ E - γ I dR/dt γ I 其中σ是潜伏期转发病率平均潜伏期是1/σ。实现时状态向量变为[S; E; I; R]初始条件通常设E00。这个模型能模拟出疫情峰值的延迟。4.3 模拟干预措施动态参数β模型最直接的应用之一是评估公共卫生干预措施的效果。例如在疫情爆发后第50天实施社交隔离使传播率β降低一半。这可以通过在ODE函数中让β成为时间t的函数来实现。function dydt odefun_SIR_intervention(t, y, beta0, gamma, N, intervention_day, reduction_factor) % 带干预措施的SIR模型 if t intervention_day beta beta0; else beta beta0 * reduction_factor; % 例如 reduction_factor 0.5 end S y(1); I y(2); R y(3); dS_dt -beta * S * I / N; dI_dt beta * S * I / N - gamma * I; dR_dt gamma * I; dydt [dS_dt; dI_dt; dR_dt]; end在调用时需要传递干预时间和降低因子。通过对比干预前后模拟曲线的差异可以直观展示隔离措施“压平曲线”的效果。4.4 参数估计与模型拟合拿到一套源码除了做模拟实验更高阶的用法是利用它来拟合真实数据从而估计参数β和γ。假设你有一段时间内每天的新增感染病例数据你可以构建一个最小化模型输出与实际数据差异的优化问题。基本思路是定义一个损失函数如误差平方和衡量模型预测的新增感染人数dI/dt γI注意新增感染是βSI/N与实际数据的差距。使用MATLAB的优化工具箱函数如fminsearch或lsqcurvefit寻找使损失函数最小的β和γ。将估计出的参数代回模型观察拟合效果。这是一个进阶话题涉及最优化理论但却是将理论模型应用于实际问题的关键一步。在源码基础上添加一个参数估计模块能极大提升代码的实用价值。常见问题拟合时可能会遇到参数不唯一或过拟合的问题。解决方法是一是利用先验知识约束参数范围如平均感染期1/γ应在合理的医学范围内二是使用更复杂但更真实的模型如SEIR三是确保有足够多的数据点。5. 常见问题排查与调试心得实录在实际运行和修改这些模型代码时你肯定会遇到各种问题。下面是我在多次使用和教学中总结的一些典型“坑”和解决方法。5.1 模型结果与理论预期不符问题表现比如SIR模型中感染者曲线没有出现先上升后下降的经典峰形而是一直上升或一直下降。可能原因1参数设置导致R01。检查你的β和γ。如果R0 β/γ 1疾病无法流行感染者会指数衰减。确保你设置的R0 1例如β0.3 γ0.1 R03。可能原因2初始易感者比例太低。如果你设置的初始感染者I0太大导致S0/N很小可能一开始就跨过了流行阈值。尝试用I01 S0N-1。可能原因3时间范围太短。疫情发展需要时间如果t_end设得太小比如只有10天可能还没看到峰值就结束了。尝试延长模拟时间到100或200天。调试方法首先打印出计算出的R0值。然后单独运行模型并绘制S和I的绝对数值和比例。观察初期前几个时间步I是否在增长。5.2 ODE求解器报错或结果异常问题表现ode45返回NaN非数字或报错“积分容差无法满足”。可能原因1微分方程中存在除零操作。在SIR模型的方程dS/dt -β * S * I / N中如果N设置为0会导致问题。确保N0。更隐蔽的情况是在模拟过程中S和I由于数值误差变为负数导致后续计算异常。ODE求解器可能会产生微小的负值。可能原因2参数值极端。如果β或γ的值非常大如1000方程会变得非常“刚性”ode45可能失效。可以考虑使用适合刚性方程的求解器ode15s。可能原因3时间跨度或初始条件设置不当。调试方法在ODE函数内部开头添加一段保护性代码强制状态变量为非负但这在物理上可能不精确仅用于调试function dydt odefun_SIR(t, y, beta, gamma, N) y(y0) 0; % 防止负值仅用于调试 S y(1); I y(2); R y(3); ... end如果问题解决说明确实出现了负值。更严谨的做法是检查模型假设或使用更稳定的求解器。5.3 图形显示问题或对比不清晰问题表现曲线挤在一起图例覆盖曲线或子图布局混乱。可能原因1没有使用hold on和hold off正确管理绘图。在同一坐标系画多条线时必须在第一条线之前或之后用hold on画完后用hold off。可能原因2线型、颜色区分度不够。使用-实线、--虚线、:点线、-.点划线结合不同颜色b,r,g,m,k来区分曲线。可能原因3图例位置不佳。使用legend(Name1,Name2,..., Location, best)让MATLAB自动选择最佳位置或手动指定如northwest左上角。调试方法养成好的绘图习惯。在绘图代码块前后使用figure创建新窗口使用subplot进行多图排列时注意三个数字的含义行数、列数、当前子图索引为每个子图添加清晰的xlabelylabel和title。5.4 性能优化与代码整洁当模型变得复杂如多群体、网络模型或需要进行大量参数扫描时代码速度可能成为瓶颈。技巧1向量化参数扫描。如果需要研究不同β值的影响避免在for循环内调用ode45。可以尝试将参数作为向量处理但更常见且清晰的做法还是使用循环优先保证代码可读性。技巧2预分配数组。如果在循环中存储结果务必预先分配好存储数组如results zeros(numSimulations, length(t))这能显著提升MATLAB循环速度。技巧3将代码模块化。将模型定义、求解、绘图、分析分别写成独立的函数或脚本段。主脚本只负责调用和高级控制。这样不仅调试方便也便于代码复用。最后分享一个最深的体会理解模型最好的方式不是只看而是动手去“破坏”它。尝试把SIR模型中的gamma设为0看看是不是变成了SI模型尝试在SIS模型中让beta变得非常小观察感染者是否最终消失通过这些实验你对模型参数敏感性和疾病动力学的理解会深刻得多。这套MATLAB源码集锦就是你开始这些实验的完美沙盒。本文还有配套的精品资源点击获取 从SI到SIR：MATLAB传染病模型源码拆解与建模实战指南
简介本资源是一套面向生物信息学、数学建模与公共卫生研究初学者的MATLAB传染病动力学仿真源码集聚焦SI、SIS、SIR三类经典微分方程模型解决传染病传播机制理解、参数敏感性分析及防控策略模拟等核心问题。压缩包共含多个.m主程序文件如si_sim.m、sis_sim.m、sir_sim.m及配套参数配置脚本全部为可直接运行的MATLAB函数与ODE求解脚本用于构建微分方程组、调用ode45数值求解器、生成S/I/R时间演化曲线并实现动态可视化包体大小仅75KB轻量易部署。已有4784人学习下载适合高校数学建模课程实践、流行病学仿真实验及科研入门训练。读者可直接复现三种模型的完整建模流程——从方程定义、参数设置、数值求解到结果绘图并通过对比不同β感染率与γ恢复率下的传播曲线深入掌握模型差异与现实映射逻辑。1. 项目概述从代码集锦到建模实战手头有一套别人分享的“MATLAB源码集锦”里面打包了传染病SI、SIS、SIR三种经典模型的代码这大概是很多数学建模新手入门时都会遇到的“宝藏资源”。但说实话我见过太多人拿到这种源码包后只是机械地运行一下看到几条曲线跳出来就觉得“哦我懂了”然后就把代码扔在一边。这其实浪费了这些代码最大的价值——它们不仅是拿来“跑”的更是用来“拆”和“改”的是理解传染病动力学建模思想最直接的脚手架。SI、SIS、SIR这三个模型可以说是流行病学数学建模的“三原色”。它们用极其简洁的微分方程刻画了人群在疾病传播过程中的状态流转。SI模型描述的是像某些无法治愈的慢性感染或早期理论中的理想情况SIS模型则加入了康复后可能再次被感染的特性比如普通感冒而SIR模型引入了“移出者”的概念刻画了那些感染后获得持久免疫力或死亡而不再参与传播的过程例如麻疹、天花。掌握它们你不仅是在学几个公式更是在掌握一种用数学语言量化社会现象、预测事态发展的思维方式。这套源码的价值就在于它把抽象的方程变成了可视化的动态过程让你能亲手调整参数直观感受“传播率”、“恢复率”这些关键因子如何左右一场“虚拟疫情”的走向。接下来我不会仅仅带大家浏览一遍代码。我们将以这套源码为蓝本深入这三个模型的数学核心并手把手教你如何从“会用”进阶到“会改”。你会学到如何将代码结构拆解成清晰的模块如何根据实际问题调整模型假设比如考虑隔离措施、疫苗接种以及如何对输出结果进行专业的分析和可视化呈现。无论你是正在备战数学建模竞赛的学生还是对复杂系统建模感兴趣的爱好者这篇内容都将帮你把这一套“源码集锦”真正转化为你工具箱里一件趁手的兵器。2. 模型核心思想与数学原理拆解在打开MATLAB代码之前我们必须先吃透模型背后的数学思想。代码只是方程的解算器和画笔方程本身才是灵魂。这三个模型都基于一个核心框架仓室模型。我们把总人口N划分为几个互斥的“仓室”个体在不同仓室间流动流动的速率由微分方程描述。2.1 SI模型最简单的传播逻辑SI模型是所有传染病模型的起点。它假设人群只分为两类易感者和感染者。一旦感染终身患病且具有传染性不会康复也不会死亡或者说模型不关心感染后的其他结局。这听起来很极端但它适用于描述某些无法治愈的慢性病在封闭群体中早期的传播动态或者在理论研究中作为基准。它的微分方程极其简洁易感者 dS/dt -β * S * I / N感染者 dI/dt β * S * I / N这里S是易感者数量I是感染者数量N是总人口恒定。β是整个模型的关键称为有效接触率或传播系数。方程 dS/dt -β * S * I / N 是理解所有后续模型的基础单位时间内新感染的人数正比于当前易感者数量S、当前感染者数量I以及他们之间发生有效接触的概率由β/N体现。这个“有效接触”是指足以导致疾病传播的接触。由于总人口N不变且S I N所以实际上只需要解一个关于I的方程即可。注意很多初学者会混淆β和另一个常见参数“基本再生数R0”。在SI模型中由于没有康复所有感染者最终都会感染所有人R0的概念不那么突出但β直接决定了疫情发展的速度。β越大曲线上升得越陡峭。2.2 SIS模型考虑康复与再感染SIS模型在SI的基础上增加了一点现实性感染者可以康复。但康复后他们并没有获得免疫力而是立刻重新变为易感者。这就形成了一个S - I - S的循环。很多细菌性感染如淋病、普通感冒针对特定毒株可以用SIS模型近似描述。它的微分方程变为易感者 dS/dt -β * S * I / N γ * I感染者 dI/dt β * S * I / N - γ * I这里引入了一个新参数γ称为恢复率或移出率。它的倒数 1/γ 可以理解为平均感染期。例如γ0.2/天意味着平均感染期是5天。方程中-γ * I 项表示感染者以速率γ康复而γ * I 项则表示这些康复者重新加入了易感者行列。SIS模型的一个重要特征是可能存在一个地方病平衡点。当传播和恢复达到平衡时感染人数会稳定在一个非零的水平而不是像SI模型那样感染所有人。这个平衡点取决于 β 和 γ 的比值。这里就引出了流行病学中至关重要的概念——基本再生数 R0。在SIS模型中R0 β / γ。当 R0 1 时疾病会流行并最终达到地方病平衡当 R0 1 时疾病会逐渐消失。2.3 SIR模型经典的“免疫者”模型SIR模型是三个模型中最著名、应用最广的。它在SIS的基础上增加了一个新的仓室移出者。个体从感染者状态移出后不再回到易感者状态而是进入R仓室。R可以代表获得持久免疫力的人也可以代表因病死亡的人——总之他们不再参与疾病的传播过程。麻疹、天花、水痘等疾病符合这个特征。它的微分方程系统是易感者 dS/dt -β * S * I / N感染者 dI/dt β * S * I / N - γ * I移出者 dR/dt γ * I注意dS/dt的方程和SI模型一样因为易感者只会减少。dI/dt的方程和SIS模型一样但移出项γ*I不再流回S而是流向了R。总人口 N S I R 依然守恒。SIR模型的行为非常经典疫情从少数感染者开始易感者数量减少感染者数量先上升后下降最终所有感染者都移出疫情结束留下部分从未感染的易感者和大量具有免疫力的移出者。最终未感染的人口比例取决于R0。基本再生数 R0 β / γ在这里同样至关重要。它表示在一个全部为易感者的人群中一个感染者在其整个传染期内平均能感染的人数。R0 1是疾病能够流行的阈值。实操心得理解参数的单位是避免错误的关键。β的单位通常是“每人每天”表示一个感染者每天能有效感染的人数在完全易感人群中。γ的单位是“每天”表示每天康复的比例。在编程时确保时间步长如ode45求解器中的时间向量与参数的单位通常是“天”匹配否则模拟出的疫情速度会失真。3. MATLAB源码结构解析与关键函数实现拿到一套打包好的源码最忌讳的就是直接点“运行”。我们先要像解剖一样看清它的代码结构、数据流和关键函数。一套良好的SI/SIS/SIR源码集锦通常会包含以下几个部分主脚本、模型方程函数、参数设置区、求解调用和绘图部分。我们来逐一拆解。3.1 主脚本框架与参数初始化一个典型的main.m或run_models.m脚本会像下面这样组织。它的核心作用是设置环境、定义参数、调用不同的模型函数并展示结果。% 清空环境关闭所有图形窗口确保干净的开始 clear; close all; clc; % 参数设置区域 % 这是你需要反复修改和实验的地方 N 1000; % 总人口 I0 1; % 初始感染者人数 S0 N - I0; % 初始易感者人数 R0 0; % 初始移出者人数 (SIR模型用) % 流行病学参数 beta 0.3; % 传播系数 (每人每天) gamma 0.1; % 恢复率 (每天) % 计算基本再生数 R0 R0_basic beta / gamma; fprintf(基本再生数 R0 %.2f\n, R0_basic); % 时间设置 t_start 0; t_end 200; % 模拟总时长 (天) tspan [t_start, t_end]; % 时间区间供ODE求解器使用 % 调用不同模型进行模拟 % 通常会有三个独立的代码块或子函数调用 % 1. 模拟SI模型 % 2. 模拟SIS模型 % 3. 模拟SIR模型 % 绘制结果对比图 % 将三个模型的结果绘制在一张或几张图中进行对比注意事项参数初始化部分至关重要。beta和gamma的值需要你根据实际疾病的文献或题目要求进行设定。例如流感的R0大约在1.2-1.8之间麻疹则高达12-18。通过调整这些参数你可以模拟不同传染力的疾病。初始感染者I0通常设为1但也可以研究输入性病例I0稍大的影响。3.2 核心微分方程组的ODE函数定义模型的灵魂在于微分方程组。在MATLAB中我们使用一个独立的函数文件例如odefun_SIR.m来定义这些方程供ode45这类求解器调用。这是将数学公式转化为代码的关键一步。以SIR模型为例我们创建一个函数文件odefun_SIR.mfunction dydt odefun_SIR(t, y, beta, gamma, N) % SIR模型的微分方程函数 % 输入: % t: 时间 (未直接使用但ODE求解器要求此参数) % y: 状态向量 [S; I; R] % beta, gamma, N: 参数 % 输出: % dydt: 导数向量 [dS/dt; dI/dt; dR/dt] S y(1); I y(2); R y(3); % 虽然方程中用不到dR/dt里的R但这里解出以保持向量完整 dS_dt -beta * S * I / N; dI_dt beta * S * I / N - gamma * I; dR_dt gamma * I; dydt [dS_dt; dI_dt; dR_dt]; end对于SIS模型odefun_SIS.mfunction dydt odefun_SIS(t, y, beta, gamma, N) % SIS模型的微分方程函数 % 状态向量 y [S; I] S y(1); I y(2); dS_dt -beta * S * I / N gamma * I; dI_dt beta * S * I / N - gamma * I; dydt [dS_dt; dI_dt]; end对于SI模型就更简单了由于SIN我们通常只对I建立方程或者依然用两个变量但令gamma0。关键细节函数头function dydt odefun_SIR(t, y, beta, gamma, N)中的参数传递方式。我们通过匿名函数或在调用ode45时传递额外参数。推荐使用匿名函数方式因为它更清晰ode45((t,y) odefun_SIR(t,y,beta,gamma,N), tspan, y0)这样betagammaN这些参数就被“绑定”到了微分方程函数中。3.3 使用ODE求解器进行数值积分定义了微分方程我们需要用数值方法求解它。MATLAB的ode45龙格-库塔法是解决非刚性常微分方程的首选对于这些传染病模型完全够用。在主脚本中调用求解器的代码块如下% --- SIR模型求解 --- y0_SIR [S0; I0; R0]; % 初始条件向量 % 使用匿名函数传递参数 [t_SIR, y_SIR] ode45((t,y) odefun_SIR(t,y,beta,gamma,N), tspan, y0_SIR); % 提取结果 S_SIR y_SIR(:, 1); I_SIR y_SIR(:, 2); R_SIR y_SIR(:, 3); % --- SIS模型求解 --- y0_SIS [S0; I0]; [t_SIS, y_SIS] ode45((t,y) odefun_SIS(t,y,beta,gamma,N), tspan, y0_SIS); S_SIS y_SIS(:, 1); I_SIS y_SIS(:, 2); % --- SI模型求解 --- % SI模型可以看作gamma0的SIR或SIS特例这里用解析解或简单ODE求解 % 方法1使用ODE求解与SIR类似但gamma0 beta_SI 0.5; % SI模型的传播率可能设得不同 gamma_SI 0; y0_SI [S0; I0]; [t_SI, y_SI] ode45((t,y) odefun_SIR(t,y,beta_SI,gamma_SI,N), tspan, y0_SI); S_SI y_SI(:, 1); I_SI y_SI(:, 2);ode45返回两个数组t是时间点向量y是对应时间点的状态向量每一列是一个状态变量。我们需要把这些数据提取出来以便绘图和分析。实操心得ode45是自适应步长的它返回的时间点t不是均匀的。如果你需要均匀时间序列上的结果比如每天的数据可以在tspan中指定更多输出点例如tspan linspace(0, 200, 201)这样就会输出0,1,2,...,200共201个时间点的解。但注意这并不改变求解精度只是增加了输出密度。3.4 结果可视化与对比分析绘图是将数学结果直观化的关键。一个好的对比图能清晰展示三个模型的本质差异。figure(Position, [100, 100, 1200, 400]); % 设置图形窗口位置和大小 % 子图1三类人群数量随时间变化以SIR为例 subplot(1, 3, 1); plot(t_SIR, S_SIR, b-, LineWidth, 2); hold on; plot(t_SIR, I_SIR, r-, LineWidth, 2); plot(t_SIR, R_SIR, g-, LineWidth, 2); hold off; xlabel(时间 (天)); ylabel(人口数量); title(SIR模型动态); legend(易感者 S, 感染者 I, 移出者 R, Location, best); grid on; % 子图2感染者比例对比三个模型 subplot(1, 3, 2); plot(t_SI, I_SI/N, k--, LineWidth, 1.5); hold on; plot(t_SIS, I_SIS/N, m-., LineWidth, 1.5); plot(t_SIR, I_SIR/N, r-, LineWidth, 2); hold off; xlabel(时间 (天)); ylabel(感染者比例 (I/N)); title(感染者比例对比); legend(SI模型, SIS模型, SIR模型, Location, best); grid on; % 添加R0信息标注 text(0.05*t_end, 0.9, sprintf(R_0 %.2f, R0_basic), Units, normalized, FontSize, 10); % 子图3相平面图SIR模型展示S-I关系 subplot(1, 3, 3); plot(S_SIR/N, I_SIR/N, b-, LineWidth, 1.5); xlabel(易感者比例 (S/N)); ylabel(感染者比例 (I/N)); title(SIR模型相平面 (S-I)); grid on; % 添加阈值线 S 1/R0 if R0_basic 0 threshold 1/R0_basic; line([threshold, threshold], ylim, Color, r, LineStyle, --, LineWidth, 1); text(threshold0.02, 0.05, S1/R_0, Color, r); end这段代码生成一个包含三个子图的综合视图。第一个子图展示SIR模型中三类人群的完整动态第二个子图直接对比三个模型中感染者比例的变化这是最核心的观察第三个子图是SIR模型的相平面图它展示了感染者比例随易感者比例变化的轨迹那条红色的阈值线S 1/R0是理论上的“疫情峰值线”当易感者比例降到这条线以下时感染者比例开始下降。避坑技巧绘图时将数量转化为比例除以总人口N常常更有意义因为它消除了人口规模的影响便于不同规模人群间的比较。另外使用hold on和hold off管理图形叠加使用legend和xlabel等函数完善图表信息这些都是生成专业图表的基本功。在竞赛或报告中清晰的图表能极大提升作品质量。4. 从理解到应用模型拓展与实战技巧掌握了基础模型的实现我们就可以开始“魔改”了。数学建模的魅力在于根据实际问题调整模型。下面介绍几个常见的拓展方向和实现技巧。4.1 引入出生与死亡SEIR模型雏形基础SIR模型假设人口恒定但长期模拟中出生和自然死亡是必须考虑的。我们可以引入一个常数出生率Λ和自然死亡率μ。此时微分方程变为 dS/dt Λ - β S I / N - μ S dI/dt β S I / N - γ I - μ I dR/dt γ I - μ R 总人口N不再恒定但通常假设出生率等于死亡率Λ μ N来保持人口稳定。在MATLAB中实现只需修改对应的ODE函数即可。这种模型更适用于研究疾病的地方性流行。4.2 加入潜伏期SEIR模型很多传染病如流感、COVID-19感染后不会立即具有传染性而是有一段潜伏期。这时需要在S和I之间增加一个暴露者仓室。这就是SEIR模型 dS/dt -β S I / N dE/dt β S I / N - σ E dI/dt σ E - γ I dR/dt γ I 其中σ是潜伏期转发病率平均潜伏期是1/σ。实现时状态向量变为[S; E; I; R]初始条件通常设E00。这个模型能模拟出疫情峰值的延迟。4.3 模拟干预措施动态参数β模型最直接的应用之一是评估公共卫生干预措施的效果。例如在疫情爆发后第50天实施社交隔离使传播率β降低一半。这可以通过在ODE函数中让β成为时间t的函数来实现。function dydt odefun_SIR_intervention(t, y, beta0, gamma, N, intervention_day, reduction_factor) % 带干预措施的SIR模型 if t intervention_day beta beta0; else beta beta0 * reduction_factor; % 例如 reduction_factor 0.5 end S y(1); I y(2); R y(3); dS_dt -beta * S * I / N; dI_dt beta * S * I / N - gamma * I; dR_dt gamma * I; dydt [dS_dt; dI_dt; dR_dt]; end在调用时需要传递干预时间和降低因子。通过对比干预前后模拟曲线的差异可以直观展示隔离措施“压平曲线”的效果。4.4 参数估计与模型拟合拿到一套源码除了做模拟实验更高阶的用法是利用它来拟合真实数据从而估计参数β和γ。假设你有一段时间内每天的新增感染病例数据你可以构建一个最小化模型输出与实际数据差异的优化问题。基本思路是定义一个损失函数如误差平方和衡量模型预测的新增感染人数dI/dt γI注意新增感染是βSI/N与实际数据的差距。使用MATLAB的优化工具箱函数如fminsearch或lsqcurvefit寻找使损失函数最小的β和γ。将估计出的参数代回模型观察拟合效果。这是一个进阶话题涉及最优化理论但却是将理论模型应用于实际问题的关键一步。在源码基础上添加一个参数估计模块能极大提升代码的实用价值。常见问题拟合时可能会遇到参数不唯一或过拟合的问题。解决方法是一是利用先验知识约束参数范围如平均感染期1/γ应在合理的医学范围内二是使用更复杂但更真实的模型如SEIR三是确保有足够多的数据点。5. 常见问题排查与调试心得实录在实际运行和修改这些模型代码时你肯定会遇到各种问题。下面是我在多次使用和教学中总结的一些典型“坑”和解决方法。5.1 模型结果与理论预期不符问题表现比如SIR模型中感染者曲线没有出现先上升后下降的经典峰形而是一直上升或一直下降。可能原因1参数设置导致R01。检查你的β和γ。如果R0 β/γ 1疾病无法流行感染者会指数衰减。确保你设置的R0 1例如β0.3 γ0.1 R03。可能原因2初始易感者比例太低。如果你设置的初始感染者I0太大导致S0/N很小可能一开始就跨过了流行阈值。尝试用I01 S0N-1。可能原因3时间范围太短。疫情发展需要时间如果t_end设得太小比如只有10天可能还没看到峰值就结束了。尝试延长模拟时间到100或200天。调试方法首先打印出计算出的R0值。然后单独运行模型并绘制S和I的绝对数值和比例。观察初期前几个时间步I是否在增长。5.2 ODE求解器报错或结果异常问题表现ode45返回NaN非数字或报错“积分容差无法满足”。可能原因1微分方程中存在除零操作。在SIR模型的方程dS/dt -β * S * I / N中如果N设置为0会导致问题。确保N0。更隐蔽的情况是在模拟过程中S和I由于数值误差变为负数导致后续计算异常。ODE求解器可能会产生微小的负值。可能原因2参数值极端。如果β或γ的值非常大如1000方程会变得非常“刚性”ode45可能失效。可以考虑使用适合刚性方程的求解器ode15s。可能原因3时间跨度或初始条件设置不当。调试方法在ODE函数内部开头添加一段保护性代码强制状态变量为非负但这在物理上可能不精确仅用于调试function dydt odefun_SIR(t, y, beta, gamma, N) y(y0) 0; % 防止负值仅用于调试 S y(1); I y(2); R y(3); ... end如果问题解决说明确实出现了负值。更严谨的做法是检查模型假设或使用更稳定的求解器。5.3 图形显示问题或对比不清晰问题表现曲线挤在一起图例覆盖曲线或子图布局混乱。可能原因1没有使用hold on和hold off正确管理绘图。在同一坐标系画多条线时必须在第一条线之前或之后用hold on画完后用hold off。可能原因2线型、颜色区分度不够。使用-实线、--虚线、:点线、-.点划线结合不同颜色b,r,g,m,k来区分曲线。可能原因3图例位置不佳。使用legend(Name1,Name2,..., Location, best)让MATLAB自动选择最佳位置或手动指定如northwest左上角。调试方法养成好的绘图习惯。在绘图代码块前后使用figure创建新窗口使用subplot进行多图排列时注意三个数字的含义行数、列数、当前子图索引为每个子图添加清晰的xlabelylabel和title。5.4 性能优化与代码整洁当模型变得复杂如多群体、网络模型或需要进行大量参数扫描时代码速度可能成为瓶颈。技巧1向量化参数扫描。如果需要研究不同β值的影响避免在for循环内调用ode45。可以尝试将参数作为向量处理但更常见且清晰的做法还是使用循环优先保证代码可读性。技巧2预分配数组。如果在循环中存储结果务必预先分配好存储数组如results zeros(numSimulations, length(t))这能显著提升MATLAB循环速度。技巧3将代码模块化。将模型定义、求解、绘图、分析分别写成独立的函数或脚本段。主脚本只负责调用和高级控制。这样不仅调试方便也便于代码复用。最后分享一个最深的体会理解模型最好的方式不是只看而是动手去“破坏”它。尝试把SIR模型中的gamma设为0看看是不是变成了SI模型尝试在SIS模型中让beta变得非常小观察感染者是否最终消失通过这些实验你对模型参数敏感性和疾病动力学的理解会深刻得多。这套MATLAB源码集锦就是你开始这些实验的完美沙盒。本文还有配套的精品资源点击获取