用MATLAB复现经典SEIR模型：从零开始搭建你的第一个疫情传播仿真（附完整代码）
用MATLAB构建SEIR模型：零基础实现疫情传播动态仿真
当第一次看到传染病传播曲线的陡峭上升时，我被数学模型的预测能力震撼了。作为流行病学研究的基础工具，SEIR模型用简洁的微分方程揭示了病毒扩散的内在规律。本文将带你从零开始，用MATLAB构建这个经典模型，无需深厚的数学背景，只要跟着步骤操作，90分钟内就能看到自己电脑上生成的疫情传播动态曲线。
1. 准备工作与环境配置
在开始建模之前，我们需要确保MATLAB环境就绪。如果你尚未安装，MATLAB官网提供30天试用版，学生还可以申请教育版许可。安装时建议选择以下工具箱：
- MATLAB基础组件（必选）
- Curve Fitting Toolbox（用于后期数据分析）
- Statistics and Machine Learning Toolbox（可选）
安装完成后，新建一个脚本文件（快捷键Ctrl+N），保存为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.25
3. 基础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(i+1) = S(i) - beta*S(i)*I(i)/N * dt; E(i+1) = E(i) + (beta*S(i)*I(i)/N - sigma*E(i)) * dt; I(i+1) = I(i) + (sigma*E(i) - gamma*I(i)) * dt; R(i+1) = 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(i+1) = 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(i+1,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模型的实现，还掌握了传染病建模的基本思路。记得保存你的工作成果，这些代码可以作为更复杂模型的基础。尝试改变参数观察曲线变化，这是理解流行病动力学的绝佳方式。