用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-, labelPopulation) 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-, labelSusceptible) plt.plot(t, I, r-, labelInfected) plt.plot(t, R, g-, labelRecovered) 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β/γ是关键参数当R01时疫情会扩散R01时疫情将消退4. 常见问题与调试技巧在实现微分方程模型时经常会遇到以下问题4.1 求解失败原因排查初值设置不当确保初值在合理范围内如人口不为负对于刚性方程可能需要使用ode15s(MATLAB)或LSODA(Python)等专用求解器参数范围不合理感染率β应在(0,1)之间时间步长不宜过大特别是变化剧烈阶段方程刚性(stiff)问题# Python中使用LSODA方法 solution odeint(model, y0, t, args(params,), mxstep5000)% 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, labelSusceptible) ax2.semilogy(t, I, labelInfected) # 对数坐标动态参数调整界面% 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, argsparams) 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, colork, 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]