用MATLAB复现经典SEIR模型：从零开始搭建你的第一个疫情传播仿真（附完整代码）

用MATLAB构建SEIR模型零基础实现疫情传播动态仿真当第一次看到传染病传播曲线的陡峭上升时我被数学模型的预测能力震撼了。作为流行病学研究的基础工具SEIR模型用简洁的微分方程揭示了病毒扩散的内在规律。本文将带你从零开始用MATLAB构建这个经典模型无需深厚的数学背景只要跟着步骤操作90分钟内就能看到自己电脑上生成的疫情传播动态曲线。1. 准备工作与环境配置在开始建模之前我们需要确保MATLAB环境就绪。如果你尚未安装MATLAB官网提供30天试用版学生还可以申请教育版许可。安装时建议选择以下工具箱MATLAB基础组件必选Curve Fitting Toolbox用于后期数据分析Statistics and Machine Learning Toolbox可选安装完成后新建一个脚本文件快捷键CtrlN保存为SEIR_Model.m。我们将在这个文件中编写所有代码。对于初次接触MATLAB的用户了解几个基本操作会很有帮助% 这是单行注释 disp(Hello World) % 在命令窗口输出文本 % 常用工作区命令 clear all % 清空工作区变量 clc % 清空命令窗口 close all % 关闭所有图形窗口2. SEIR模型核心原理解析SEIR模型将人群划分为四个互斥状态易感者(Susceptible)可能被感染的健康人群潜伏者(Exposed)已感染但无症状且不具传染性感染者(Infectious)出现症状并具有传染性康复者(Recovered)痊愈并获得免疫力的人群模型的核心是一组耦合微分方程描述各状态间的转化关系dS/dt -βSI/N dE/dt βSI/N - σE dI/dt σE - γI dR/dt γI其中关键参数包括β感染率感染者每天接触并传染易感者的概率σ潜伏期倒数1/平均潜伏期天数γ康复率1/平均感染期天数这些参数决定了疫情发展的形态。例如COVID-19的典型参数值为平均潜伏期5.2天 → σ ≈ 0.19平均感染期10天 → γ ≈ 0.1基本传染数R₀≈2.5 → β ≈ R₀×γ 0.253. 基础SEIR模型的MATLAB实现现在我们将上述方程转化为可执行的MATLAB代码。首先定义模型参数和初始条件% 参数设置 N 1e6; % 总人口数 beta 0.25; % 感染率 sigma 0.19; % 潜伏期倒数 gamma 0.1; % 康复率 % 初始条件 I0 1; % 初始感染者 E0 0; % 初始潜伏者 R0 0; % 初始康复者 S0 N - I0; % 初始易感者接下来我们需要用欧拉法求解微分方程组。虽然MATLAB有内置的ODE求解器但为了教学目的我们手动实现% 时间设置 t_start 0; t_end 180; % 模拟180天 dt 1; % 时间步长(天) t t_start:dt:t_end; % 初始化数组 S zeros(size(t)); S(1) S0; E zeros(size(t)); E(1) E0; I zeros(size(t)); I(1) I0; R zeros(size(t)); R(1) R0; % 欧拉法迭代 for i 1:length(t)-1 S(i1) S(i) - beta*S(i)*I(i)/N * dt; E(i1) E(i) (beta*S(i)*I(i)/N - sigma*E(i)) * dt; I(i1) I(i) (sigma*E(i) - gamma*I(i)) * dt; R(i1) R(i) gamma*I(i) * dt; end最后我们绘制各人群随时间变化的曲线% 绘制结果 figure(Color,white,Position,[100,100,800,500]) plot(t,S,b, t,E,y, t,I,r, t,R,g,LineWidth,2) grid on xlabel(时间(天)) ylabel(人数) legend(易感者,潜伏者,感染者,康复者,Location,best) title(基础SEIR模型仿真结果) set(gca,FontSize,12)运行这段代码你将看到典型的疫情传播曲线感染者数量先上升达到峰值然后逐渐下降最终人群主要分布在康复者和少数未感染的易感者中。4. 模型优化与参数敏感性分析基础模型虽然能反映传播趋势但现实情况更为复杂。我们可以从几个方面进行改进4.1 加入干预措施假设在第60天实施社交距离政策使感染率β降低50%for i 1:length(t)-1 % 第60天实施干预 current_beta beta; if t(i) 60 current_beta beta * 0.5; end S(i1) S(i) - current_beta*S(i)*I(i)/N * dt; % 其余方程保持不变... end4.2 参数敏感性分析了解不同参数对结果的影响至关重要。我们可以创建参数扫描函数function peak_infections simulate_SEIR(beta, sigma, gamma) % 参数设置同上... % 运行模拟... peak_infections max(I); end % 测试不同R0值的影响 R0_values 1.5:0.5:3.5; peaks zeros(size(R0_values)); for i 1:length(R0_values) peaks(i) simulate_SEIR(R0_values(i)*gamma, sigma, gamma); end figure plot(R0_values, peaks/N*100, -o) xlabel(基本传染数 R_0) ylabel(感染峰值占总人口比例(%)) title(R0对疫情峰值的影响)4.3 可视化优化使用subplot创建多面板图表更全面地展示结果figure(Color,white,Position,[100,100,1000,700]) subplot(2,2,1) plot(t,S,b,LineWidth,2) title(易感者变化) grid on subplot(2,2,2) plot(t,E,y,LineWidth,2) title(潜伏者变化) grid on subplot(2,2,3) plot(t,I,r,LineWidth,2) title(感染者变化) grid on subplot(2,2,4) plot(t,R,g,LineWidth,2) title(康复者变化) grid on5. 进阶扩展思路掌握了基础SEIR模型后你可以尝试以下扩展方向5.1 模型变体SEIRS模型考虑免疫力随时间衰减康复者可能再次变为易感者SEIQRD模型加入隔离(Q)和死亡(D)人群年龄结构化模型考虑不同年龄组的接触模式差异5.2 实际数据拟合使用curve fitting工具箱将模型与实际疫情数据拟合% 假设有实际数据actual_cases ft fittype(A*exp(b*x),coefficients,{A,b}); fit_result fit(t(1:30), actual_cases(1:30), ft);5.3 空间扩展将模型扩展到空间维度模拟病毒在不同地区间的传播% 假设有5个相互连接的城市 num_cities 5; connectivity rand(num_cities); % 城市间连接强度 % 每个城市维护自己的SEIR状态 S zeros(length(t), num_cities); I zeros(length(t), num_cities); % 初始化... for i 1:length(t)-1 for city 1:num_cities % 计算来自其他城市的输入感染 imported_infections sum(connectivity(:,city).*I(i,:)); S(i1,city) S(i,city) - (beta*S(i,city)*I(i,city)/N 0.1*imported_infections) * dt; % 其余方程... end end6. 常见问题与调试技巧初学者在实现SEIR模型时常遇到以下问题数值不稳定时间步长dt过大导致解发散。尝试减小dt或改用ode45求解器[t,y] ode45(SEIR_equations, [0 180], [S0 E0 I0 R0]);参数设置不合理导致曲线形态异常。确保所有参数为非负值各人群初始值之和等于总人口Nβ/(γμ) ≈ R₀μ为死亡率单位不一致时间参数σ、γ的单位需统一通常为1/天可视化问题曲线重叠难以区分。可以使用不同线型和颜色添加图例和网格调整坐标轴范围当模型行为不符合预期时建议打印中间变量值检查计算过程简化模型排除干扰因素与已知正确结果进行对比% 调试示例检查人群总数是否守恒 total_population S E I R; if any(abs(diff(total_population)) 1e-6) warning(总人口不守恒) end通过本教程你不仅学会了SEIR模型的实现还掌握了传染病建模的基本思路。记得保存你的工作成果这些代码可以作为更复杂模型的基础。尝试改变参数观察曲线变化这是理解流行病动力学的绝佳方式。