尧图网站建设 
1. 项目概述从一次直播课看数学建模的实战入门最近看到不少同学在讨论“数学建模清风第一次直播”的内容核心聚焦在传染病模型和微分方程拟合上。这让我想起了自己早年接触数学建模的经历当时也是从一个经典的传染病问题入门的。这次直播的主题选得非常精准它几乎是所有数模新手通往实战的必经之路。传染病模型特别是SIR及其变体不仅仅是公共卫生领域的工具更是一个绝佳的数学建模教学案例。它完美串联了问题背景分析、模型假设、微分方程建立、参数求解以及结果分析这一整套建模流程。而微分方程拟合则是将抽象的数学方程与现实数据连接起来的关键桥梁是检验模型、预测趋势的核心技能。无论你是数学、计算机还是其他理工科的学生掌握这套从微分方程到参数拟合的完整链条都能为你解决各类动态系统问题如种群增长、谣言传播、化学反应动力学等打下坚实的基础。接下来我就结合这次直播可能涉及的内容以及我个人的一些实战经验为大家拆解一下这个经典课题里的门道。2. 核心思路拆解为什么是传染病模型和微分方程拟合2.1 传染病模型作为“第一课”的不可替代性选择传染病模型作为数学建模的入门案例绝非偶然。首先它的背景极其直观每个人都对“传染”有基本认知这降低了理解门槛。其次其核心动态——健康者、感染者、康复者之间的转化——可以用清晰的“仓室”来划分对应到建模中的状态变量。例如经典的SIR模型就是将总人口分为易感者(S)、感染者(I)、移除者(R)三类。这种划分方式逻辑清晰是建立微分方程组的完美起点。更重要的是传染病模型具备良好的可扩展性。你可以从最简单的SI模型开始逐步增加考虑免疫的SIR再到考虑潜伏期的SEIR甚至加入出生死亡、疫苗接种、空间扩散等因素。这种由简入繁的递进能让学习者逐步体会到如何根据实际问题复杂度来调整模型结构这是建模思维的核心训练。2.2 微分方程拟合连接理论与现实的“标尺”建立了微分方程模型里面往往包含一些未知参数比如传染率、恢复率。这些参数无法直接从理论推导必须通过实际观测数据来“反推”这个过程就是参数估计或拟合。微分方程拟合是数学建模中技术含量很高的一环。它不同于简单的曲线拟合如多项式拟合你的拟合对象不是一个显式函数而是一个微分方程组的解。这意味着你需要利用数值方法如龙格-库塔法先求解微分方程组得到一个“数值解”再将这个数值解与真实数据进行比较通过优化算法如最小二乘法调整参数使得数值解尽可能贴近真实数据。这个过程深刻地体现了“用数学模型描述世界并用数据校准模型”的完整科学范式。掌握它你就掌握了量化分析大多数动态过程的能力。3. 从SIR模型出发基础框架与方程建立3.1 SIR模型的基本假设与仓室划分我们以最经典的SIR模型为例来具体看看它是如何构建的。首先我们必须做出明确的假设这是所有建模的起点总人口恒定不考虑出生、死亡和迁移总人口数N S(t) I(t) R(t) 为常数。均匀混合人群充分混合任何一个易感者与任何一个感染者接触的机会均等。传染率与恢复率恒定单位时间内一个感染者能使β个易感者被感染感染者以固定速率γ恢复并具有永久免疫力。疾病传播瞬时完成没有潜伏期。基于这些假设三个仓室间的转移关系就非常清晰了易感者(S)通过接触感染者(I)而变成感染者感染者(I)以一定速率恢复变成移除者(R)。移除者不再参与疾病传播过程。3.2 微分方程组的推导与解释根据上述转移关系我们可以用微分方程来描述每个仓室人数随时间的变化率易感者(S)的变化率 dS/dt易感者只会减少减少的速度取决于当前感染者数量I和易感者数量S以及传染接触的有效性β。并且由于总人口均匀混合一个易感者被感染的概率与感染者占总人口的比例(I/N)成正比。因此单位时间内新增的感染人数为 β * S * (I/N)。由于S在减少所以方程是负的dS/dt -β * S * I / N。感染者(I)的变化率 dI/dt感染者一方面从易感者中补充进来即β * S * I / N另一方面以速率γ恢复移出。所以dI/dt β * S * I / N - γ * I。移除者(R)的变化率 dR/dt移除者只从感染者中恢复而来所以dR/dt γ * I。注意这里β是有效接触率它包含了接触频率和传染概率的综合影响。有时也会看到形式dS/dt -β S I这里的β含义不同是β / N的关系使用时务必明确参数定义这是初期容易混淆的地方。这三个方程就构成了封闭的SIR模型微分方程组。给定初始时刻的S(0), I(0), R(0)和参数β, γ我们就可以通过数值计算预测未来任意时刻各仓室的人数变化。4. 模型求解数值方法入门与编程实现4.1 为什么需要数值求解我们得到了SIR模型的微分方程组但除了极特殊情况这类非线性方程组是找不到解析解即用初等函数公式表达的解的。因此我们必须依赖数值方法来获得方程在离散时间点上的近似解。这就好比我们无法写出一个复杂函数曲线具体的公式但可以用计算机算出这条曲线上成千上万个点的坐标然后用线连起来就能无限逼近真实曲线。4.2 龙格-库塔法最常用的“发动机”在科学计算中四阶龙格-库塔法RK4是求解常微分方程初值问题最经典、最常用的方法之一。它精度和稳定性兼顾对于像SIR这样的非刚性方程非常有效。其核心思想是利用当前点的导数信息通过一个巧妙的加权平均来预测下一个点的函数值。虽然其推导涉及泰勒展开但对于使用者而言我们可以把它理解为一个高度精确的“递推公式”。以SIR模型为例我们的状态变量是一个向量Y [S, I, R]微分方程右边定义了导数dY/dt f(t, Y) [-β*S*I/N, β*S*I/N - γ*I, γ*I]。RK4的每一步计算如下k1 f(t, Y) k2 f(t dt/2, Y dt*k1/2) k3 f(t dt/2, Y dt*k2/2) k4 f(t dt, Y dt*k3) Y_new Y (dt/6) * (k1 2*k2 2*k3 k4)其中dt是我们选择的时间步长。通过循环迭代我们就可以从初始时刻t0算到结束时刻t_end得到一系列[S, I, R]的数值。4.3 编程实现示例Python在实际操作中我们通常不需要自己编写RK4的循环。利用Python的SciPy库可以非常便捷地实现。下面是一个完整的示例import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 1. 定义SIR模型的微分方程 def sir_model(t, y, beta, gamma, N): S, I, R y dSdt -beta * S * I / N dIdt beta * S * I / N - gamma * I dRdt gamma * I return [dSdt, dIdt, dRdt] # 2. 设置参数和初始条件 N 1000 # 总人口 I0, R0 10, 0 # 初始感染者和移除者 S0 N - I0 - R0 # 初始易感者 beta 0.3 # 传染率 gamma 0.1 # 恢复率平均感染期 1/gamma 10天 y0 [S0, I0, R0] # 初始条件向量 # 3. 定义时间范围0到200天 t_span [0, 200] t_eval np.linspace(0, 200, 200) # 希望输出的时间点 # 4. 数值求解 solution solve_ivp(sir_model, t_span, y0, args(beta, gamma, N), t_evalt_eval, methodRK45) # 5. 提取结果并绘图 S solution.y[0] I solution.y[1] R solution.y[2] t solution.t plt.figure(figsize(10,6)) plt.plot(t, S, labelSusceptible, linewidth2) plt.plot(t, I, labelInfected, linewidth2) plt.plot(t, R, labelRecovered, linewidth2) plt.xlabel(Time (days)) plt.ylabel(Number of people) plt.title(SIR Model Simulation (β0.3, γ0.1)) plt.legend() plt.grid(True) plt.show()这段代码清晰地展示了从定义模型到求解可视化的全流程。solve_ivp函数内部默认使用的就是类似RK45自适应步长的龙格-库塔法的算法我们无需操心具体实现细节。实操心得在初次编程时最容易出错的地方是参数单位的一致性。比如如果我们的数据是“每日新增感染数”那么beta和gamma的单位必须是“每天”。如果时间步长dt或求解器输出的时间序列单位是“周”那么参数也需要相应调整为“每周”。务必在代码注释中明确所有参数的单位这是保证结果有意义的前提。5. 核心环节基于真实数据的参数拟合5.1 问题定义与损失函数构建模拟运行只是第一步真正的挑战在于当一场真实的传染病发生时我们如何确定模型中的参数β和γ这就是拟合要解决的问题。假设我们手头有一组时间序列数据例如每天报告的累计感染人数通常对应I(t)R(t)或每日新增感染人数。 我们的目标是找到一组参数(beta, gamma)使得SIR模型在这些参数下运行得到的数值解I_model(t) R_model(t)与真实的累计感染数据C_data(t)之间的差距最小。这个“差距”需要用数学语言量化即构建一个损失函数。最常用的是最小二乘法损失函数定义为所有时间点上模型值与观测值之差的平方和Loss(β, γ) Σ [ (I_model(t_i) R_model(t_i) - C_data(t_i))^2 ]我们的任务就是寻找使Loss最小的(β, γ)组合。这是一个典型的无约束非线性优化问题。5.2 拟合流程与优化算法选择整个拟合流程可以概括为以下步骤数据预处理清洗数据确保时间序列连续处理可能的异常值或报告延迟。将数据整理成(时间, 累计感染数)的数组。定义带参数的模型函数编写一个函数输入参数(beta, gamma)和初始条件输出模型模拟的I(t)R(t)时间序列。定义损失函数如上所述计算模型输出与真实数据之间的误差平方和。选择优化器并求解调用优化算法自动搜索最小化损失函数的参数。对于SIR模型这类规模较小、参数较少2-3个的问题SciPy中的curve_fit或minimize函数就足够强大。curve_fit专为曲线拟合设计接口简单minimize功能更通用可以尝试不同的优化算法如Nelder-Mead, BFGS, L-BFGS-B等。5.3 实战代码使用curve_fit进行拟合假设我们有一份模拟的“真实数据”它是由参数(β0.25, γ0.05)生成的并加上了一些随机噪声。我们现在假装不知道这些参数尝试从数据中反推出来。from scipy.optimize import curve_fit import pandas as pd # 生成带噪声的模拟“真实数据”假设我们不知道真实参数 true_beta, true_gamma 0.25, 0.05 solution_true solve_ivp(sir_model, [0, 150], y0, args(true_beta, true_gamma, N), t_evalnp.linspace(0, 150, 151), dense_outputTrue) C_true solution_true.y[1] solution_true.y[2] # IR即累计感染 np.random.seed(42) C_data C_true * (1 0.05 * np.random.randn(len(C_true))) # 添加5%的高斯噪声 C_data np.maximum(C_data, 0) # 确保非负 # 定义一个用于拟合的包装函数 def fit_function(t, beta, gamma): 对于给定的参数beta, gamma返回模型预测的累计感染数 IR sol solve_ivp(sir_model, [t[0], t[-1]], y0, args(beta, gamma, N), t_evalt, methodRK45) # 返回累计感染数 return sol.y[1] sol.y[2] # 执行拟合 time_points np.linspace(0, 150, 151) # 提供参数的初始猜测值这对收敛很重要 initial_guess [0.2, 0.1] # 设置参数边界例如β和γ都应大于0 bounds ([0, 0], [np.inf, np.inf]) try: popt, pcov curve_fit(fit_function, time_points, C_data, p0initial_guess, boundsbounds, maxfev5000) fitted_beta, fitted_gamma popt print(f真实参数: β{true_beta:.4f}, γ{true_gamma:.4f}) print(f拟合参数: β{fitted_beta:.4f}, γ{fitted_gamma:.4f}) print(f基本再生数 R0 (拟合) {fitted_beta/fitted_gamma:.4f}) except Exception as e: print(f拟合失败: {e}) # 可视化拟合效果 C_fitted fit_function(time_points, fitted_beta, fitted_gamma) plt.figure(figsize(10,6)) plt.scatter(time_points, C_data, alpha0.5, labelNoisy Data (Simulated), s10) plt.plot(time_points, C_true, k--, labelTrue Model (Noise-Free), linewidth2) plt.plot(time_points, C_fitted, r-, labelFitted Model, linewidth2) plt.xlabel(Time (days)) plt.ylabel(Cumulative Infections) plt.title(SIR Model Parameter Fitting) plt.legend() plt.grid(True) plt.show()这段代码演示了完整的拟合过程。curve_fit会调用fit_function多次不断调整beta和gamma直到模型输出的曲线与散点数据C_data的差距最小。最终输出的popt就是最优参数估计值。6. 关键细节与进阶考量6.1 初始条件的敏感性拟合结果不仅依赖于参数(β, γ)也强烈依赖于初始条件[S0, I0, R0]。在实际疫情早期I0和R0往往难以准确估计。一个常见的处理方法是将I0也作为一个待拟合的参数与β, γ一同优化。但这会增加优化问题的复杂度从2维搜索变为3维搜索并且可能使得解的不确定性增大。更稳健的做法是利用疫情早期指数增长阶段的数据通过线性回归粗略估计初始增长率和I0将其作为优化起点。6.2 模型选择与复杂度权衡SIR模型是基础但现实往往更复杂。当拟合效果不佳时可能需要考虑更复杂的模型SEIR模型在S和I之间增加一个潜伏期仓室(E)。这需要多拟合一个参数潜伏期的倒数σ。考虑时变参数在疫情中干预措施如封控、戴口罩会导致传染率β随时间下降。此时可以假设β是一个随时间变化的函数如分段常数或指数衰减这大大增加了拟合难度通常需要更专业的贝叶斯推断方法。考虑人口动力学对于长周期传染病需加入出生和自然死亡项。注意事项切勿盲目追求复杂模型。奥卡姆剃刀原理告诉我们在能同等解释数据的情况下应选择最简单的模型。增加参数虽然可能让曲线拟合得更好损失函数更小但极易导致“过拟合”——即模型过分迎合数据中的噪声反而失去了预测未来或解释现象的能力。通常先用简单模型尝试如果系统性地偏离数据如峰值时间、峰值高度、曲线对称性再考虑增加复杂度。6.3 基本再生数R0的计算与意义在拟合得到β和γ后一个至关重要的衍生指标是基本再生数R0。在SIR模型中R0 β / γ。它的流行病学含义是在一个全部是易感者的人群中一个感染者在其整个传染期内平均能传染的人数。R0 1意味着疾病会蔓延R0 1则疾病会逐渐消失。通过拟合估计R0可以为公共卫生决策如需要将传染率降低多少才能控制疫情提供关键定量依据。7. 常见问题与排查技巧实录在实际操作中从模型建立到成功拟合你会遇到各种“坑”。下面是我总结的一些典型问题及解决方法。7.1 拟合失败或结果不合理问题表现优化算法不收敛或收敛到明显错误的参数值如负值或数量级离谱。排查思路检查参数初始猜测值优化算法像爬山初始点如果离山顶太远可能困在局部洼地或根本找不到路。尝试多个不同的初始猜测值观察结果是否稳定。可以参考文献中类似疾病的参数范围进行设置。检查参数边界像传染率β、恢复率γ这些参数从物理意义上必须大于0。在curve_fit或minimize中设置合理的下界如bounds(0, [np.inf, np.inf])可以防止算法搜索到无意义的区域。缩放问题如果状态变量S,I,R的数量级是百万而参数β、γ在0.1左右数值计算可能会产生精度问题。考虑将人口数归一化例如令N1S,I,R代表比例这样参数和变量的数量级更接近有助于优化器稳定工作。数据本身的问题检查数据是否存在大量零值、突变或平台期。早期数据太少或噪声太大也可能导致拟合失败。有时需要对数据进行平滑处理或选取疫情发展较为明显的阶段进行拟合。7.2 模型模拟曲线与数据形状不匹配问题表现即使拟合出参数模拟曲线的上升速度、峰值位置或下降趋势也与数据有肉眼可见的系统性偏差。可能原因与对策现象可能原因对策建议模型上升比数据慢低估了传染率β或高估了初始易感者S0尝试增大β的初始猜测值检查并修正初始条件I0早期数据可能漏报严重。模型峰值过早到达传染率β估计过高或恢复率γ估计过低尝试引入时变β考虑防控措施的影响。或检查模型是否忽略了潜伏期考虑SEIR。模型下降比数据慢恢复率γ估计过低增大γ的初始猜测值。或考虑是否有部分感染者未被报告即实际I比数据大导致模型中的“感染者池”清空得慢。曲线无法重现数据的“长尾”或“平台”SIR模型假设永久免疫疫情后会结束。长尾可能源于1. 不断有新的输入病例。2. 免疫力非永久如SIRS模型。3. 数据报告延迟或积累效应。根据实际背景判断。如果是输入病例可在模型中加入一个很小的外部输入项。如果怀疑免疫力减弱可改用SIRS模型。7.3 数值求解不稳定问题表现在求解微分方程时解出现剧烈振荡或溢出变成NaN。解决方法减小时间步长对于solve_ivp可以尝试更严格的最大步长max_step或使用更适合刚性问题的求解器如methodRadau或methodBDF。检查方程实现仔细核对微分方程代码确保正负号、分母不为零等。一个常见的错误是在分母中直接使用了I当I0时会导致除零错误。在实际编程中可以给分母加一个极小值epsilon或确保初始I0不为0。归一化变量如前所述将变量归一化到[0,1]区间能有效提升数值稳定性。7.4 关于“过拟合”的再提醒当你使用复杂模型如SEIR或带时变参数的SIR并拟合了大量参数后即使得到了与历史数据完美贴合的曲线也切勿过于乐观。一定要进行交叉验证或后验预测检查。简单来说就是用一部分数据如前80%的时间点来拟合参数然后用拟合的模型去预测剩余20%的数据看预测效果如何。如果预测效果很差说明模型可能只是记住了历史数据的噪声泛化能力不足。在数学建模竞赛中清晰地展示模型验证过程是获得高分的关键。最后我想分享的一点个人体会是传染病模型拟合是一个需要耐心和反复调试的过程。它不像解数学题有标准答案更像是一个与数据对话、不断修正对疾病传播认知的科学探索。第一次尝试可能结果不理想但这正是学习的价值所在——通过分析失败的原因你对模型假设、参数意义和数据特性的理解会深刻得多。不妨从这次直播介绍的经典SIR模型开始亲手运行一遍代码尝试拟合一份公开的疫情数据如某次流感的早期数据你会对整个数学建模的流程有一个非常扎实的感性认识。 从SIR传染病模型到微分方程拟合：数学建模实战入门指南
1. 项目概述从一次直播课看数学建模的实战入门最近看到不少同学在讨论“数学建模清风第一次直播”的内容核心聚焦在传染病模型和微分方程拟合上。这让我想起了自己早年接触数学建模的经历当时也是从一个经典的传染病问题入门的。这次直播的主题选得非常精准它几乎是所有数模新手通往实战的必经之路。传染病模型特别是SIR及其变体不仅仅是公共卫生领域的工具更是一个绝佳的数学建模教学案例。它完美串联了问题背景分析、模型假设、微分方程建立、参数求解以及结果分析这一整套建模流程。而微分方程拟合则是将抽象的数学方程与现实数据连接起来的关键桥梁是检验模型、预测趋势的核心技能。无论你是数学、计算机还是其他理工科的学生掌握这套从微分方程到参数拟合的完整链条都能为你解决各类动态系统问题如种群增长、谣言传播、化学反应动力学等打下坚实的基础。接下来我就结合这次直播可能涉及的内容以及我个人的一些实战经验为大家拆解一下这个经典课题里的门道。2. 核心思路拆解为什么是传染病模型和微分方程拟合2.1 传染病模型作为“第一课”的不可替代性选择传染病模型作为数学建模的入门案例绝非偶然。首先它的背景极其直观每个人都对“传染”有基本认知这降低了理解门槛。其次其核心动态——健康者、感染者、康复者之间的转化——可以用清晰的“仓室”来划分对应到建模中的状态变量。例如经典的SIR模型就是将总人口分为易感者(S)、感染者(I)、移除者(R)三类。这种划分方式逻辑清晰是建立微分方程组的完美起点。更重要的是传染病模型具备良好的可扩展性。你可以从最简单的SI模型开始逐步增加考虑免疫的SIR再到考虑潜伏期的SEIR甚至加入出生死亡、疫苗接种、空间扩散等因素。这种由简入繁的递进能让学习者逐步体会到如何根据实际问题复杂度来调整模型结构这是建模思维的核心训练。2.2 微分方程拟合连接理论与现实的“标尺”建立了微分方程模型里面往往包含一些未知参数比如传染率、恢复率。这些参数无法直接从理论推导必须通过实际观测数据来“反推”这个过程就是参数估计或拟合。微分方程拟合是数学建模中技术含量很高的一环。它不同于简单的曲线拟合如多项式拟合你的拟合对象不是一个显式函数而是一个微分方程组的解。这意味着你需要利用数值方法如龙格-库塔法先求解微分方程组得到一个“数值解”再将这个数值解与真实数据进行比较通过优化算法如最小二乘法调整参数使得数值解尽可能贴近真实数据。这个过程深刻地体现了“用数学模型描述世界并用数据校准模型”的完整科学范式。掌握它你就掌握了量化分析大多数动态过程的能力。3. 从SIR模型出发基础框架与方程建立3.1 SIR模型的基本假设与仓室划分我们以最经典的SIR模型为例来具体看看它是如何构建的。首先我们必须做出明确的假设这是所有建模的起点总人口恒定不考虑出生、死亡和迁移总人口数N S(t) I(t) R(t) 为常数。均匀混合人群充分混合任何一个易感者与任何一个感染者接触的机会均等。传染率与恢复率恒定单位时间内一个感染者能使β个易感者被感染感染者以固定速率γ恢复并具有永久免疫力。疾病传播瞬时完成没有潜伏期。基于这些假设三个仓室间的转移关系就非常清晰了易感者(S)通过接触感染者(I)而变成感染者感染者(I)以一定速率恢复变成移除者(R)。移除者不再参与疾病传播过程。3.2 微分方程组的推导与解释根据上述转移关系我们可以用微分方程来描述每个仓室人数随时间的变化率易感者(S)的变化率 dS/dt易感者只会减少减少的速度取决于当前感染者数量I和易感者数量S以及传染接触的有效性β。并且由于总人口均匀混合一个易感者被感染的概率与感染者占总人口的比例(I/N)成正比。因此单位时间内新增的感染人数为 β * S * (I/N)。由于S在减少所以方程是负的dS/dt -β * S * I / N。感染者(I)的变化率 dI/dt感染者一方面从易感者中补充进来即β * S * I / N另一方面以速率γ恢复移出。所以dI/dt β * S * I / N - γ * I。移除者(R)的变化率 dR/dt移除者只从感染者中恢复而来所以dR/dt γ * I。注意这里β是有效接触率它包含了接触频率和传染概率的综合影响。有时也会看到形式dS/dt -β S I这里的β含义不同是β / N的关系使用时务必明确参数定义这是初期容易混淆的地方。这三个方程就构成了封闭的SIR模型微分方程组。给定初始时刻的S(0), I(0), R(0)和参数β, γ我们就可以通过数值计算预测未来任意时刻各仓室的人数变化。4. 模型求解数值方法入门与编程实现4.1 为什么需要数值求解我们得到了SIR模型的微分方程组但除了极特殊情况这类非线性方程组是找不到解析解即用初等函数公式表达的解的。因此我们必须依赖数值方法来获得方程在离散时间点上的近似解。这就好比我们无法写出一个复杂函数曲线具体的公式但可以用计算机算出这条曲线上成千上万个点的坐标然后用线连起来就能无限逼近真实曲线。4.2 龙格-库塔法最常用的“发动机”在科学计算中四阶龙格-库塔法RK4是求解常微分方程初值问题最经典、最常用的方法之一。它精度和稳定性兼顾对于像SIR这样的非刚性方程非常有效。其核心思想是利用当前点的导数信息通过一个巧妙的加权平均来预测下一个点的函数值。虽然其推导涉及泰勒展开但对于使用者而言我们可以把它理解为一个高度精确的“递推公式”。以SIR模型为例我们的状态变量是一个向量Y [S, I, R]微分方程右边定义了导数dY/dt f(t, Y) [-β*S*I/N, β*S*I/N - γ*I, γ*I]。RK4的每一步计算如下k1 f(t, Y) k2 f(t dt/2, Y dt*k1/2) k3 f(t dt/2, Y dt*k2/2) k4 f(t dt, Y dt*k3) Y_new Y (dt/6) * (k1 2*k2 2*k3 k4)其中dt是我们选择的时间步长。通过循环迭代我们就可以从初始时刻t0算到结束时刻t_end得到一系列[S, I, R]的数值。4.3 编程实现示例Python在实际操作中我们通常不需要自己编写RK4的循环。利用Python的SciPy库可以非常便捷地实现。下面是一个完整的示例import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 1. 定义SIR模型的微分方程 def sir_model(t, y, beta, gamma, N): S, I, R y dSdt -beta * S * I / N dIdt beta * S * I / N - gamma * I dRdt gamma * I return [dSdt, dIdt, dRdt] # 2. 设置参数和初始条件 N 1000 # 总人口 I0, R0 10, 0 # 初始感染者和移除者 S0 N - I0 - R0 # 初始易感者 beta 0.3 # 传染率 gamma 0.1 # 恢复率平均感染期 1/gamma 10天 y0 [S0, I0, R0] # 初始条件向量 # 3. 定义时间范围0到200天 t_span [0, 200] t_eval np.linspace(0, 200, 200) # 希望输出的时间点 # 4. 数值求解 solution solve_ivp(sir_model, t_span, y0, args(beta, gamma, N), t_evalt_eval, methodRK45) # 5. 提取结果并绘图 S solution.y[0] I solution.y[1] R solution.y[2] t solution.t plt.figure(figsize(10,6)) plt.plot(t, S, labelSusceptible, linewidth2) plt.plot(t, I, labelInfected, linewidth2) plt.plot(t, R, labelRecovered, linewidth2) plt.xlabel(Time (days)) plt.ylabel(Number of people) plt.title(SIR Model Simulation (β0.3, γ0.1)) plt.legend() plt.grid(True) plt.show()这段代码清晰地展示了从定义模型到求解可视化的全流程。solve_ivp函数内部默认使用的就是类似RK45自适应步长的龙格-库塔法的算法我们无需操心具体实现细节。实操心得在初次编程时最容易出错的地方是参数单位的一致性。比如如果我们的数据是“每日新增感染数”那么beta和gamma的单位必须是“每天”。如果时间步长dt或求解器输出的时间序列单位是“周”那么参数也需要相应调整为“每周”。务必在代码注释中明确所有参数的单位这是保证结果有意义的前提。5. 核心环节基于真实数据的参数拟合5.1 问题定义与损失函数构建模拟运行只是第一步真正的挑战在于当一场真实的传染病发生时我们如何确定模型中的参数β和γ这就是拟合要解决的问题。假设我们手头有一组时间序列数据例如每天报告的累计感染人数通常对应I(t)R(t)或每日新增感染人数。 我们的目标是找到一组参数(beta, gamma)使得SIR模型在这些参数下运行得到的数值解I_model(t) R_model(t)与真实的累计感染数据C_data(t)之间的差距最小。这个“差距”需要用数学语言量化即构建一个损失函数。最常用的是最小二乘法损失函数定义为所有时间点上模型值与观测值之差的平方和Loss(β, γ) Σ [ (I_model(t_i) R_model(t_i) - C_data(t_i))^2 ]我们的任务就是寻找使Loss最小的(β, γ)组合。这是一个典型的无约束非线性优化问题。5.2 拟合流程与优化算法选择整个拟合流程可以概括为以下步骤数据预处理清洗数据确保时间序列连续处理可能的异常值或报告延迟。将数据整理成(时间, 累计感染数)的数组。定义带参数的模型函数编写一个函数输入参数(beta, gamma)和初始条件输出模型模拟的I(t)R(t)时间序列。定义损失函数如上所述计算模型输出与真实数据之间的误差平方和。选择优化器并求解调用优化算法自动搜索最小化损失函数的参数。对于SIR模型这类规模较小、参数较少2-3个的问题SciPy中的curve_fit或minimize函数就足够强大。curve_fit专为曲线拟合设计接口简单minimize功能更通用可以尝试不同的优化算法如Nelder-Mead, BFGS, L-BFGS-B等。5.3 实战代码使用curve_fit进行拟合假设我们有一份模拟的“真实数据”它是由参数(β0.25, γ0.05)生成的并加上了一些随机噪声。我们现在假装不知道这些参数尝试从数据中反推出来。from scipy.optimize import curve_fit import pandas as pd # 生成带噪声的模拟“真实数据”假设我们不知道真实参数 true_beta, true_gamma 0.25, 0.05 solution_true solve_ivp(sir_model, [0, 150], y0, args(true_beta, true_gamma, N), t_evalnp.linspace(0, 150, 151), dense_outputTrue) C_true solution_true.y[1] solution_true.y[2] # IR即累计感染 np.random.seed(42) C_data C_true * (1 0.05 * np.random.randn(len(C_true))) # 添加5%的高斯噪声 C_data np.maximum(C_data, 0) # 确保非负 # 定义一个用于拟合的包装函数 def fit_function(t, beta, gamma): 对于给定的参数beta, gamma返回模型预测的累计感染数 IR sol solve_ivp(sir_model, [t[0], t[-1]], y0, args(beta, gamma, N), t_evalt, methodRK45) # 返回累计感染数 return sol.y[1] sol.y[2] # 执行拟合 time_points np.linspace(0, 150, 151) # 提供参数的初始猜测值这对收敛很重要 initial_guess [0.2, 0.1] # 设置参数边界例如β和γ都应大于0 bounds ([0, 0], [np.inf, np.inf]) try: popt, pcov curve_fit(fit_function, time_points, C_data, p0initial_guess, boundsbounds, maxfev5000) fitted_beta, fitted_gamma popt print(f真实参数: β{true_beta:.4f}, γ{true_gamma:.4f}) print(f拟合参数: β{fitted_beta:.4f}, γ{fitted_gamma:.4f}) print(f基本再生数 R0 (拟合) {fitted_beta/fitted_gamma:.4f}) except Exception as e: print(f拟合失败: {e}) # 可视化拟合效果 C_fitted fit_function(time_points, fitted_beta, fitted_gamma) plt.figure(figsize(10,6)) plt.scatter(time_points, C_data, alpha0.5, labelNoisy Data (Simulated), s10) plt.plot(time_points, C_true, k--, labelTrue Model (Noise-Free), linewidth2) plt.plot(time_points, C_fitted, r-, labelFitted Model, linewidth2) plt.xlabel(Time (days)) plt.ylabel(Cumulative Infections) plt.title(SIR Model Parameter Fitting) plt.legend() plt.grid(True) plt.show()这段代码演示了完整的拟合过程。curve_fit会调用fit_function多次不断调整beta和gamma直到模型输出的曲线与散点数据C_data的差距最小。最终输出的popt就是最优参数估计值。6. 关键细节与进阶考量6.1 初始条件的敏感性拟合结果不仅依赖于参数(β, γ)也强烈依赖于初始条件[S0, I0, R0]。在实际疫情早期I0和R0往往难以准确估计。一个常见的处理方法是将I0也作为一个待拟合的参数与β, γ一同优化。但这会增加优化问题的复杂度从2维搜索变为3维搜索并且可能使得解的不确定性增大。更稳健的做法是利用疫情早期指数增长阶段的数据通过线性回归粗略估计初始增长率和I0将其作为优化起点。6.2 模型选择与复杂度权衡SIR模型是基础但现实往往更复杂。当拟合效果不佳时可能需要考虑更复杂的模型SEIR模型在S和I之间增加一个潜伏期仓室(E)。这需要多拟合一个参数潜伏期的倒数σ。考虑时变参数在疫情中干预措施如封控、戴口罩会导致传染率β随时间下降。此时可以假设β是一个随时间变化的函数如分段常数或指数衰减这大大增加了拟合难度通常需要更专业的贝叶斯推断方法。考虑人口动力学对于长周期传染病需加入出生和自然死亡项。注意事项切勿盲目追求复杂模型。奥卡姆剃刀原理告诉我们在能同等解释数据的情况下应选择最简单的模型。增加参数虽然可能让曲线拟合得更好损失函数更小但极易导致“过拟合”——即模型过分迎合数据中的噪声反而失去了预测未来或解释现象的能力。通常先用简单模型尝试如果系统性地偏离数据如峰值时间、峰值高度、曲线对称性再考虑增加复杂度。6.3 基本再生数R0的计算与意义在拟合得到β和γ后一个至关重要的衍生指标是基本再生数R0。在SIR模型中R0 β / γ。它的流行病学含义是在一个全部是易感者的人群中一个感染者在其整个传染期内平均能传染的人数。R0 1意味着疾病会蔓延R0 1则疾病会逐渐消失。通过拟合估计R0可以为公共卫生决策如需要将传染率降低多少才能控制疫情提供关键定量依据。7. 常见问题与排查技巧实录在实际操作中从模型建立到成功拟合你会遇到各种“坑”。下面是我总结的一些典型问题及解决方法。7.1 拟合失败或结果不合理问题表现优化算法不收敛或收敛到明显错误的参数值如负值或数量级离谱。排查思路检查参数初始猜测值优化算法像爬山初始点如果离山顶太远可能困在局部洼地或根本找不到路。尝试多个不同的初始猜测值观察结果是否稳定。可以参考文献中类似疾病的参数范围进行设置。检查参数边界像传染率β、恢复率γ这些参数从物理意义上必须大于0。在curve_fit或minimize中设置合理的下界如bounds(0, [np.inf, np.inf])可以防止算法搜索到无意义的区域。缩放问题如果状态变量S,I,R的数量级是百万而参数β、γ在0.1左右数值计算可能会产生精度问题。考虑将人口数归一化例如令N1S,I,R代表比例这样参数和变量的数量级更接近有助于优化器稳定工作。数据本身的问题检查数据是否存在大量零值、突变或平台期。早期数据太少或噪声太大也可能导致拟合失败。有时需要对数据进行平滑处理或选取疫情发展较为明显的阶段进行拟合。7.2 模型模拟曲线与数据形状不匹配问题表现即使拟合出参数模拟曲线的上升速度、峰值位置或下降趋势也与数据有肉眼可见的系统性偏差。可能原因与对策现象可能原因对策建议模型上升比数据慢低估了传染率β或高估了初始易感者S0尝试增大β的初始猜测值检查并修正初始条件I0早期数据可能漏报严重。模型峰值过早到达传染率β估计过高或恢复率γ估计过低尝试引入时变β考虑防控措施的影响。或检查模型是否忽略了潜伏期考虑SEIR。模型下降比数据慢恢复率γ估计过低增大γ的初始猜测值。或考虑是否有部分感染者未被报告即实际I比数据大导致模型中的“感染者池”清空得慢。曲线无法重现数据的“长尾”或“平台”SIR模型假设永久免疫疫情后会结束。长尾可能源于1. 不断有新的输入病例。2. 免疫力非永久如SIRS模型。3. 数据报告延迟或积累效应。根据实际背景判断。如果是输入病例可在模型中加入一个很小的外部输入项。如果怀疑免疫力减弱可改用SIRS模型。7.3 数值求解不稳定问题表现在求解微分方程时解出现剧烈振荡或溢出变成NaN。解决方法减小时间步长对于solve_ivp可以尝试更严格的最大步长max_step或使用更适合刚性问题的求解器如methodRadau或methodBDF。检查方程实现仔细核对微分方程代码确保正负号、分母不为零等。一个常见的错误是在分母中直接使用了I当I0时会导致除零错误。在实际编程中可以给分母加一个极小值epsilon或确保初始I0不为0。归一化变量如前所述将变量归一化到[0,1]区间能有效提升数值稳定性。7.4 关于“过拟合”的再提醒当你使用复杂模型如SEIR或带时变参数的SIR并拟合了大量参数后即使得到了与历史数据完美贴合的曲线也切勿过于乐观。一定要进行交叉验证或后验预测检查。简单来说就是用一部分数据如前80%的时间点来拟合参数然后用拟合的模型去预测剩余20%的数据看预测效果如何。如果预测效果很差说明模型可能只是记住了历史数据的噪声泛化能力不足。在数学建模竞赛中清晰地展示模型验证过程是获得高分的关键。最后我想分享的一点个人体会是传染病模型拟合是一个需要耐心和反复调试的过程。它不像解数学题有标准答案更像是一个与数据对话、不断修正对疾病传播认知的科学探索。第一次尝试可能结果不理想但这正是学习的价值所在——通过分析失败的原因你对模型假设、参数意义和数据特性的理解会深刻得多。不妨从这次直播介绍的经典SIR模型开始亲手运行一遍代码尝试拟合一份公开的疫情数据如某次流感的早期数据你会对整个数学建模的流程有一个非常扎实的感性认识。