用Python和MATLAB搞定数学建模：从人口预测到传染病模型实战
用Python和MATLAB搞定数学建模：从人口预测到传染病模型实战
数学建模竞赛中，微分方程模型一直是解决实际问题的利器。但对于许多参赛者来说，如何将抽象的数学公式转化为可运行的代码，往往成为拦路虎。本文将带你用Python和MATLAB两大工具，从零实现经典的人口预测和传染病模型，通过代码直观理解模型行为，掌握参数调整技巧，并对比两种语言在科学计算中的差异。
1. 环境准备与工具对比
在开始建模前，我们需要配置好开发环境。Python方面推荐使用Anaconda发行版，它集成了我们需要的SciPy、NumPy和Matplotlib库。MATLAB则需要安装基础模块和Symbolic Math Toolbox。
两种工具在微分方程求解上各有优势：
| 特性 | Python (SciPy) | MATLAB (ode45) |
|---|---|---|
| 语法简洁度 | 中等（需导入多个库） | 高（内置函数直接调用） |
| 可视化便捷性 | 优秀（Matplotlib灵活） | 良好（内置绘图函数） |
| 求解速度 | 较快（依赖实现） | 极快（高度优化） |
| 学习曲线 | 平缓（通用语言） | 陡峭（专用语言） |
| 开源/商业 | 完全开源 | 商业软件 |
安装必要的Python库：
pip install numpy scipy matplotlibMATLAB用户只需确认已安装以下工具箱：
ver % 查看已安装工具箱2. Logistic人口增长模型实现
Logistic模型是描述资源有限环境下人口增长的经典模型，其微分方程为：
$$ \frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right) $$
其中P为人口数量，r为固有增长率，K为环境承载力。
2.1 Python实现
使用SciPy的odeint求解器：
import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def logistic_model(P, t, r, K): return r * P * (1 - P/K) # 参数设置 r = 0.1 # 增长率 K = 1000 # 环境容量 P0 = 10 # 初始人口 t = np.linspace(0, 100, 1000) # 时间范围 # 求解微分方程 solution = odeint(logistic_model, P0, t, args=(r, K)) # 可视化 plt.figure(figsize=(10,6)) plt.plot(t, solution, 'b-', label='Population') plt.xlabel('Time') plt.ylabel('Population') plt.title('Logistic Growth Model') plt.grid() plt.legend() plt.show()2.2 MATLAB实现
使用ode45求解器：
function logistic_growth % 定义参数 r = 0.1; % 增长率 K = 1000; % 环境容量 P0 = 10; % 初始值 tspan = [0 100]; % 时间范围 % 定义微分方程 function dPdt = ode_logistic(t, P) dPdt = r * P * (1 - P/K); end % 求解 [t, P] = ode45(@ode_logistic, tspan, P0); % 绘图 figure plot(t, P, 'LineWidth', 2) xlabel('Time') ylabel('Population') title('Logistic Growth Model') grid on end2.3 参数敏感性分析
调整r和K值会显著影响模型行为：
- 增长率r：值越大，人口达到稳定状态越快
- 承载力K：决定最终稳定人口规模
- 初始值P0：不影响最终结果，只改变达到稳定的时间
提示：在实际应用中，可通过最小二乘法拟合历史数据来估计r和K的真实值
3. SIR传染病模型实战
SIR模型将人群分为易感者(S)、感染者(I)和康复者(R)三类，其微分方程组为：
$$ \begin{cases} \frac{dS}{dt} = -\beta SI/N \ \frac{dI}{dt} = \beta SI/N - \gamma I \ \frac{dR}{dt} = \gamma I \end{cases} $$
3.1 Python实现
def sir_model(y, t, beta, gamma): S, I, R = y N = S + I + R dSdt = -beta * S * I / N dIdt = beta * S * I / N - gamma * I dRdt = gamma * I return [dSdt, dIdt, dRdt] # 参数设置 beta = 0.3 # 感染率 gamma = 0.1 # 康复率 N = 1000 # 总人口 I0 = 1 # 初始感染者 S0 = N - I0 # 初始易感者 R0 = 0 # 初始康复者 t = np.linspace(0, 200, 1000) # 求解 solution = odeint(sir_model, [S0, I0, R0], t, args=(beta, gamma)) S, I, R = solution.T # 可视化 plt.figure(figsize=(12,8)) plt.plot(t, S, 'b-', label='Susceptible') plt.plot(t, I, 'r-', label='Infected') plt.plot(t, R, 'g-', label='Recovered') plt.xlabel('Days') plt.ylabel('Population') plt.title('SIR Model') plt.legend() plt.grid() plt.show()3.2 MATLAB实现
function sir_model % 参数设置 beta = 0.3; % 感染率 gamma = 0.1; % 康复率 N = 1000; % 总人口 I0 = 1; % 初始感染者 S0 = N - I0; % 初始易感者 R0 = 0; % 初始康复者 tspan = [0 200]; % 定义微分方程组 function dydt = ode_sir(t, y) S = y(1); I = y(2); R = y(3); dSdt = -beta * S * I / N; dIdt = beta * S * I / N - gamma * I; dRdt = gamma * I; dydt = [dSdt; dIdt; dRdt]; end % 求解 [t, y] = ode45(@ode_sir, tspan, [S0; I0; R0]); % 绘图 figure plot(t, y(:,1), 'b-', 'LineWidth', 2); hold on plot(t, y(:,2), 'r-', 'LineWidth', 2) plot(t, y(:,3), 'g-', 'LineWidth', 2) xlabel('Days') ylabel('Population') title('SIR Model') legend('Susceptible','Infected','Recovered') grid on end3.3 模型扩展与参数优化
实际应用中，我们常需要调整基础SIR模型：
加入出生率和死亡率：
def sir_birth_death(y, t, beta, gamma, mu): S, I, R = y N = S + I + R dSdt = mu*N - beta*S*I/N - mu*S dIdt = beta*S*I/N - gamma*I - mu*I dRdt = gamma*I - mu*R return [dSdt, dIdt, dRdt]考虑疫苗接种：
function dydt = ode_sir_vaccine(t, y) S = y(1); I = y(2); R = y(3); vacc_rate = 0.01; % 每日疫苗接种率 dSdt = -beta*S*I/N - vacc_rate*S; dIdt = beta*S*I/N - gamma*I; dRdt = gamma*I + vacc_rate*S; dydt = [dSdt; dIdt; dRdt]; end
注意：基本再生数R0=β/γ是关键参数，当R0>1时疫情会扩散，R0<1时疫情将消退
4. 常见问题与调试技巧
在实现微分方程模型时，经常会遇到以下问题：
4.1 求解失败原因排查
初值设置不当：
- 确保初值在合理范围内（如人口不为负）
- 对于刚性方程，可能需要使用
ode15s(MATLAB)或LSODA(Python)等专用求解器
参数范围不合理：
- 感染率β应在(0,1)之间
- 时间步长不宜过大，特别是变化剧烈阶段
方程刚性(stiff)问题：
# Python中使用LSODA方法 solution = odeint(model, y0, t, args=(params,), mxstep=5000)% MATLAB中使用ode15s [t,y] = ode15s(@ode_func, tspan, y0);
4.2 可视化优化技巧
多子图对比：
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,5)) ax1.plot(t, S, label='Susceptible') ax2.semilogy(t, I, label='Infected') # 对数坐标动态参数调整界面：
% MATLAB中使用uicontrol创建交互界面 f = figure; uicontrol('Style','slider','Min',0.1,'Max',0.5,... 'Position',[100 20 120 20],'Callback',@update_plot);
4.3 性能优化建议
向量化运算：
# 避免循环，使用NumPy向量运算 dSdt = -beta * S * I / N # 优于循环计算每个点Jacobian矩阵提供（对刚性方程）：
options = odeset('Jacobian',@jacobian_func); [t,y] = ode15s(@ode_func, tspan, y0, options);并行计算（参数扫描时）：
from multiprocessing import Pool def simulate(params): return odeint(model, y0, t, args=params) with Pool(4) as p: results = p.map(simulate, param_list)
5. 进阶应用：模型拟合与预测
有了基础模型后，我们通常需要用真实数据来校准模型参数。
5.1 参数估计方法
最小二乘法拟合：
from scipy.optimize import curve_fit def model_wrapper(t, beta, gamma): sol = odeint(sir_model, [S0,I0,R0], t, args=(beta,gamma)) return sol[:,1] # 返回感染人数I popt, pcov = curve_fit(model_wrapper, t_data, I_data, p0=[0.3, 0.1], bounds=(0, [1,1]))MCMC方法（不确定性量化）：
import emcee def log_likelihood(theta, t, data): beta, gamma = theta model = odeint(sir_model, [S0,I0,R0], t, args=(beta,gamma)) return -0.5*np.sum((model[:,1]-data)**2) sampler = emcee.EnsembleSampler(nwalkers, ndim, log_likelihood, args=(t_data, I_data)) sampler.run_mcmc(p0, nsteps)
5.2 预测与干预分析
通过模型可以进行政策模拟：
% 模拟社交距离措施（降低beta） beta_normal = 0.3; beta_lockdown = 0.1; lockdown_start = 30; lockdown_end = 60; function dydt = ode_sir_intervention(t, y) if t >= lockdown_start && t <= lockdown_end beta = beta_lockdown; else beta = beta_normal; end % 其余部分与标准SIR相同 end5.3 模型验证技术
残差分析：
residuals = I_data - model_prediction plt.plot(t_data, residuals) plt.axhline(0, color='k', linestyle='--')交叉验证：
- 将数据分为训练集和测试集
- 用训练集估计参数
- 在测试集上验证预测效果
敏感性分析：
% 计算基本再生数R0对参数变化的敏感度 beta_range = linspace(0.1, 0.5, 20); gamma_range = linspace(0.05, 0.2, 20); [B,G] = meshgrid(beta_range, gamma_range); R0 = B./G; surf(B,G,R0)
6. 扩展模型与应用场景
基础模型可以扩展应用于更复杂的场景：
6.1 SEIR模型（考虑潜伏期）
def seir_model(y, t, beta, sigma, gamma): S,E,I,R = y N = S + E + I + R dSdt = -beta * S * I / N dEdt = beta * S * I / N - sigma * E dIdt = sigma * E - gamma * I dRdt = gamma * I return [dSdt, dEdt, dIdt, dRdt]6.2 空间扩散模型
% 使用偏微分方程工具箱求解空间扩散 model = createpde(); geometryFromEdges(model,@circleg); applyBoundaryCondition(model,'dirichlet','Edge',1:4,'u',0); specifyCoefficients(model,'m',0,'d',1,'c',1,'a',0,'f',0); setInitialConditions(model,@initialfun); generateMesh(model); tlist = linspace(0,1,20); results = solvepde(model,tlist);6.3 经济-流行病耦合模型
def economic_sir(y, t, beta, gamma, alpha): S,I,R,G = y # G代表经济活动水平 N = S + I + R dSdt = -beta*S*I/N - alpha*G*S dIdt = beta*S*I/N - gamma*I dRdt = gamma*I dGdt = -k1*I + k2*(1-G) return [dSdt,dIdt,dRdt,dGdt]