MATLAB微分方程数值求解实战：从SIR模型到热传导的建模应用

1. 项目概述从数学建模到微分方程求解的核心跨越每年暑假都是数学建模竞赛备战的黄金时期。无论是国赛、美赛还是各类地区性赛事微分方程模型都是工具箱里不可或缺的“重型武器”。从描述传染病传播的SIR模型到模拟热量扩散的热传导方程再到刻画种群竞争的Lotka-Volterra模型其背后都是常微分方程ODE或偏微分方程PDE在支撑。然而很多同学在集训时都会遇到一个尴尬的局面模型方程列出来了理论解却求不出来或者根本不存在解析解。这时候数值求解就成了连接抽象模型与具体结果的唯一桥梁而MATLAB正是搭建这座桥梁最得力的工具之一。我参加过也指导过多次数学建模集训发现大家在学习MATLAB解微分方程时最容易陷入两个极端要么对着几个内置函数死记硬背遇到复杂点的问题就束手无策要么被各种数值算法的理论吓退觉得深不可测。其实对于数学建模而言我们不需要成为数值分析专家但必须成为一个“会调参、懂诊断、能解决问题”的实战派。本次集训的核心目标就是带大家跨越从“知道函数”到“能用函数解决实际问题”这道鸿沟。我们将聚焦MATLAB中求解ODE和PDE的两大核心工具箱通过具体的建模案例拆解每一步操作背后的意图并分享那些只有踩过坑才知道的调试技巧和效率法门。无论你是刚刚接触MATLAB的新手还是想提升求解效率和稳定性的老手相信这些从实战中提炼出的经验都能让你在接下来的建模比赛中更加从容。2. 核心思路与工具箱选型为何是ODE与PDE套件在MATLAB的广阔天地里解决微分方程的函数不止一个。面对具体问题选对工具是成功的第一步。很多初学者会直接搜索“matlab 解微分方程”然后被dsolve,ode45,pdepe等一堆函数搞得眼花缭乱。我们的思路很明确根据方程类型和边界条件快速锁定最合适的求解器并理解其适用场景和局限性。2.1 常微分方程ODE求解器选型从ode45说起对于常微分方程MATLAB提供了一整套以ode为前缀的求解器如ode45,ode23,ode113,ode15s等。数字编号并非随意它暗示了算法的阶数和类型。对于数学建模中的绝大多数初值问题IVPode45是当之无愧的“首发选择”。为什么首选ode45它基于显式Runge-Kutta (4,5)公式即Dormand-Prince算法。这是一种单步法意味着计算下一步只需要前一步的信息编程实现简单。它在精度4阶和计算量之间取得了很好的平衡对于非刚性non-stiff或中等刚性问题表现优异。所谓“刚性”简单类比就是系统里存在变化速度差异巨大的多个过程比如一个化学反应中既有瞬间完成的快速反应又有缓慢进行的慢速反应。刚性方程用普通方法如ode45求解会异常缓慢甚至失败。何时考虑其他求解器对精度要求不高追求速度时可以尝试ode23它使用Bogacki-Shampine公式2,3阶步长更大计算更快适合快速预览解的大致形态。遇到刚性Stiff问题时这是建模中的一个常见坎。如果你的模型用ode45求解时步长变得极小计算时间长得离谱或者直接报错很可能遇到了刚性系统。这时应切换到刚性求解器如ode15s基于数值微分公式适用于中度刚性问题或ode23s基于修正的Rosenbrock公式适用于高度刚性问题。一个典型的刚性系统例子是包含快速衰减瞬态过程的电路模型或化学反应动力学模型。需要更高精度或处理特殊问题时ode113是多步Adams-Bashforth-Moulton算法在允许误差非常严格时可能比ode45更高效。注意不要死记硬背所有求解器。掌握ode45和ode15s这两个最具代表性的非刚性和刚性求解器就能解决95%的ODE建模问题。关键在于学会诊断问题是否为刚性。2.2 偏微分方程PDE求解策略pdepe与有限差分法偏微分方程的世界更复杂MATLAB没有像ODE那样提供“一键通吃”的函数但针对最常见的一类问题——一维空间上的抛物型和椭圆型方程或方程组提供了非常强大的内置求解器pdepe。pdepe的定位与优势pdepe专门用于求解一维空间可以是直线、球体或柱体对称情况上的抛物-椭圆型偏微分方程组。这意味着它非常适合处理诸如一维热传导、物质扩散、反应扩散方程等问题。它的优势在于封装了复杂的空间离散化和时间积分过程用户只需要按照固定格式提供方程系数、初始条件和边界条件函数大大降低了入门门槛。pdepe的局限性它仅限于一维空间问题。对于二维或三维问题或者双曲型PDE如波动方程pdepe就无能为力了。高维或复杂PDE的出路当问题超出pdepe的能力范围时我们通常需要自己实现数值方法。最常用、最直观的就是有限差分法FDM。其核心思想是用网格点上的函数值近似连续空间用差商近似偏导数从而将PDE转化为一个大型的代数方程组对于稳态问题或常微分方程组对于瞬态问题进行求解。虽然实现起来代码量更大但灵活度极高是解决复杂PDE模型的终极手段。MATLAB强大的矩阵运算能力为实现有限差分法提供了极大便利。2.3 整体求解流程设计无论是ODE还是PDE一个稳健的数值求解流程都遵循以下步骤这也是我们后续实操的蓝图方程标准化将你的模型方程整理成MATLAB求解器要求的标准形式。这是最关键的一步形式不对一切白费。编写函数文件根据标准形式编写定义方程、初始条件、边界条件PDE需要的MATLAB函数。调用求解器选择合适的求解器如ode45,pdepe并正确设置时间/空间网格、初始值等参数。结果可视化与验证绘制解随时间/空间的变化图。通过改变网格密度、容差参数等验证解的收敛性和可靠性。模型分析与应用基于数值解进行参数敏感性分析、稳定性分析等为建模结论提供支撑。3. 常微分方程ODE求解实战以传染病SIR模型为例让我们从一个经典的数学建模案例——传染病SIR模型入手完整走一遍ODE的求解流程。SIR模型将人群分为易感者S、感染者I、康复者R三类其微分方程组为 dS/dt -β * S * I / N dI/dt β * S * I / N - γ * I dR/dt γ * I 其中N S I R 为总人口常数β为感染率γ为康复率。3.1 第一步方程标准化与函数编写ode45等求解器要求方程必须写成dy/dt f(t, y)的向量形式。对于SIR模型我们令状态向量 y [S; I; R]。那么右端函数 f(t, y) 就需要计算三个导数。% 文件保存为 sir_ode.m function dydt sir_ode(t, y, beta, gamma, N) % t: 时间未显式使用但格式要求 % y: 状态向量 [S; I; R] % beta, gamma, N: 模型参数 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这里有一个关键技巧我们将参数beta,gamma,N作为函数的额外输入参数而不是在函数内部写死。这样在调用求解器时可以通过匿名函数灵活地传入参数值便于后续进行参数敏感性分析。3.2 第二步调用求解器与参数设置接下来我们在脚本或命令行中设置初始条件、时间区间和参数并调用ode45。% 模型参数 N 1000; % 总人口 I0 1; % 初始感染者 R0 0; % 初始康复者 S0 N - I0 - R0; % 初始易感者 y0 [S0; I0; R0]; % 初始状态向量 beta 0.3; % 感染率 gamma 0.1; % 康复率 (平均感染期 1/gamma 10天) % 时间区间 [0, 150] 天 tspan [0, 150]; % 调用ode45求解 % 使用匿名函数将参数传递给sir_ode [t, y] ode45((t,y) sir_ode(t, y, beta, gamma, N), tspan, y0); % 提取结果 S y(:, 1); I y(:, 2); R y(:, 3);参数设置的讲究时间区间tspan的选取很重要。如果只关心疫情峰值和时间可以设一个较长的区间让系统达到稳定即I趋于0。如果想研究短期爆发区间可以设短一些。初始感染者I0不能为0否则系统不会演化。3.3 第三步结果可视化与初步分析画出三类人群随时间的变化曲线是分析模型的基础。figure(Position, [100, 100, 800, 400]) % 设置图形位置和大小 plot(t, S, b-, LineWidth, 1.5); hold on; plot(t, I, r-, LineWidth, 1.5); plot(t, R, g-, LineWidth, 1.5); hold off; grid on; xlabel(时间 (天)); ylabel(人口数); legend(易感者 S, 感染者 I, 康复者 R, Location, best); title(sprintf(SIR模型动态 (\\beta%.2f, \\gamma%.2f, R0%.2f), beta, gamma, beta/gamma));这里我们在标题中计算并显示了基本再生数 R0 β / γ。R0 1 意味着疫情会扩散这是我们能从图中直观看到的核心结论。3.4 进阶刚性问题的识别与切换求解器假设我们研究一个化学反应模型其中某个中间产物的浓度变化极快。用ode45求解时MATLAB可能会警告Warning: Failure at t... Unable to meet integration tolerances without reducing the step size below the smallest value allowed...或者求解时间异常漫长。这时我们就需要怀疑遇到了刚性系统。一个简单的测试方法是尝试使用刚性求解器ode15s并对比求解时间和结果。% 假设 stiff_ode 是一个刚性ODE的函数 options_ode45 odeset(Stats, on); % 打开统计信息 tic; [t1, y1] ode45(stiff_ode, tspan, y0, options_ode45); time_ode45 toc; fprintf(ode45 求解时间: %.4f 秒\n, time_ode45); options_ode15s odeset(Stats, on); tic; [t2, y2] ode15s(stiff_ode, tspan, y0, options_ode15s); time_ode15s toc; fprintf(ode15s 求解时间: %.4f 秒\n, time_ode15s); % 比较最终结果是否接近 diff norm(y1(end,:) - y2(end,:)); fprintf(最终状态差异范数: %e\n, diff);如果ode15s的求解时间远短于ode45且两者最终结果一致那么就证实了刚性问题的存在后续建模就应选用ode15s。4. 偏微分方程PDE求解实战一维热传导问题我们以一维杆的热传导问题为例展示如何使用pdepe求解。方程是经典的抛物型PDE ∂u/∂t α * ∂²u/∂x², (0 x L, t 0) 其中u(x,t)是温度α是热扩散系数。边界条件设为两端绝热Neumann边界条件∂u/∂x |(x0) 0, ∂u/∂x |(xL) 0。初始条件设为在杆中心有一个高斯分布的高温u(x,0) exp(-(x-L/2)² / (2*σ²))。4.1pdepe的标准形式与函数编写pdepe要求PDE写成如下标准形式 c(x, t, u, ∂u/∂x) * ∂u/∂t x^(-m) * ∂/∂x [ x^m * f(x, t, u, ∂u/∂x) ] s(x, t, u, ∂u/∂x) 其中m0,1,2 分别对应平板、柱对称、球对称几何。f是通量项s是源项。对于我们的热传导方程m 0 平板几何c 1f α * ∂u/∂x 根据傅里叶定律热通量与温度梯度成正比s 0我们需要编写三个函数PDE函数、初始条件函数、边界条件函数。% 1. PDE函数 (保存为 heat_pde.m) function [c, f, s] heat_pde(x, t, u, DuDx, alpha) c 1; % 方程系数 c f alpha * DuDx; % 通量项 f s 0; % 源项 s end % 2. 初始条件函数 (保存为 heat_ic.m) function u0 heat_ic(x, L, sigma) % 在杆中心xL/2处设置一个高斯峰作为初始温度 u0 exp(-(x - L/2).^2 / (2 * sigma^2)); end % 3. 边界条件函数 (保存为 heat_bc.m) function [pl, ql, pr, qr] heat_bc(xl, ul, xr, ur, t, alpha) % 左边界 (x0): 绝热温度梯度为0 pl 0, ql 1 pl 0; ql 1; % 右边界 (xL): 绝热温度梯度为0 pr 0, qr 1 pr 0; qr 1; % p q * f 0 是边界条件形式。对于绝热f alpha * DuDx 0, 所以设置 p0, q1。 end边界条件设置的难点pdepe的边界条件形式为p(x, t, u) q(x, t) * f(x, t, u, ∂u/∂x) 0。对于Dirichlet条件固定温度u常数设 p u - constant, q 0。对于Neumann条件固定热流如绝热时梯度为0设 p 0, q 1因为此时要求 f α * ∂u/∂x 0。这是最容易出错的地方务必理解透彻。4.2 空间与时间网格设置及求解调用空间网格xmesh和时间向量tspan的选取直接影响求解的精度和速度。% 参数设置 L 10; % 杆的长度 alpha 0.1; % 热扩散系数 sigma 0.5; % 初始高斯分布的宽度 % 空间网格在边界附近和初始热点附近可以加密 xmesh linspace(0, L, 101); % 101个空间点通常是个不错的起点 % 时间向量关心初始扩散和最终平衡可以在初期设置密一些 tspan [0:0.1:1, 1.5:0.5:10, 15:5:50]; % 非均匀时间点 % 调用 pdepe sol pdepe(0, ... % 几何参数 m (0平板) (x,t,u,DuDx) heat_pde(x,t,u,DuDx,alpha), ... % PDE函数句柄 (x) heat_ic(x, L, sigma), ... % 初始条件函数句柄 (xl,ul,xr,ur,t) heat_bc(xl,ul,xr,ur,t,alpha), ... % 边界条件函数句柄 xmesh, tspan); % 网格 % 提取结果sol 是一个 3D 数组 (length(tspan) x length(xmesh)) u sol(:,:,1); % 我们只有一个因变量 u网格设置心得空间网格点数不宜过少否则会丢失细节特别是初始温度尖峰也不宜过多否则计算量剧增。可以从50-100点开始尝试。时间点tspan决定了输出解的时间切片。pdepe内部会使用自适应步长积分tspan只是指定了我们需要输出解的那些时刻。为了画出平滑的动画或曲线tspan可以设得密一些。4.3 结果可视化温度时空分布我们可以用多种方式可视化PDE的解。% 方式1时空分布图 (伪彩色图) figure; surf(xmesh, tspan, u, EdgeColor, none); xlabel(位置 x); ylabel(时间 t); zlabel(温度 u); title(一维热传导温度时空演化); colormap(jet); colorbar; view(2); % 俯视图可以看到等高线 % 方式2不同时刻的温度剖面图 figure; hold on; plot_indices [1, find(tspan1), find(tspan5), find(tspan20), length(tspan)]; % 选取几个时刻 colors lines(length(plot_indices)); % 获取不同颜色 for i 1:length(plot_indices) idx plot_indices(i); plot(xmesh, u(idx, :), Color, colors(i,:), LineWidth, 1.5, ... DisplayName, sprintf(t %.1f, tspan(idx))); end hold off; xlabel(位置 x); ylabel(温度 u); legend(show, Location, best); title(不同时刻的温度分布剖面); grid on;时空分布图能全局展示热量如何从中心向两端扩散并最终趋于均匀。剖面图则能更清晰地比较不同时刻分布形态的差异。5. 有限差分法FDM解PDE入门以二维泊松方程为例当问题维度升高或方程形式特殊时pdepe不再适用。例如求解一个二维矩形区域上的稳态泊松方程 ∂²u/∂x² ∂²u/∂y² f(x, y), (0 x a, 0 y b) 边界条件为Dirichlet条件u(0,y)u(a,y)u(x,0)u(x,b)0。 这是一个椭圆型方程我们可以用有限差分法将其离散化求解。5.1 差分格式推导与离散化首先在x方向将区间[0,a]分为M份步长Δx a/My方向将[0,b]分为N份步长Δy b/N。网格点坐标为 (x_i, y_j)其中 x_i iΔx, y_j jΔy, i0,...,M, j0,...,N。 在内部网格点(i,j)处用中心差分近似二阶导数 ∂²u/∂x² ≈ (u_{i-1,j} - 2u_{i,j} u_{i1,j}) / (Δx)² ∂²u/∂y² ≈ (u_{i,j-1} - 2u_{i,j} u_{i,j1}) / (Δy)² 代入泊松方程得到离散方程 (u_{i-1,j} - 2u_{i,j} u_{i1,j})/(Δx)² (u_{i,j-1} - 2u_{i,j} u_{i,j1})/(Δy)² f_{i,j} 对于所有内部点(i1,...,M-1; j1,...,N-1)我们都有这样一个方程。边界点上的u值由边界条件给出此处全为0。5.2 构建线性方程组与MATLAB求解将未知数所有内部点的u值按“行优先”或“列优先”排成一个长向量U。上面的每个差分方程都可以写成一个线性方程。最终整个离散系统可以写成一个大型的稀疏线性方程组A * U F其中A是一个(M-1)*(N-1) 阶的方阵其结构非常有规律带状、对称正定F是由源项f和边界条件贡献构成的右端向量。在MATLAB中我们不需要手动组装巨大的矩阵A。对于这种规则区域上的泊松方程可以使用poisolv针对矩形区域或更通用的pdepe的稳态求解模式但为了理解FDM我们演示一种基于矩阵运算的直观方法适用于较小网格。% 参数设置 a 1; b 1; % 区域大小 [0,1]x[0,1] M 50; N 50; % 网格划分数 dx a / M; dy b / N; x linspace(0, a, M1); y linspace(0, b, N1); % 源项函数 f(x,y) 2*pi^2 * sin(pi*x) * sin(pi*y) 其精确解为 usin(pi*x)*sin(pi*y) [X, Y] meshgrid(x(2:end-1), y(2:end-1)); % 内部点 F 2 * pi^2 * sin(pi*X) .* sin(pi*Y); F_vec F(:); % 将源项矩阵按列展开成向量 % 构建系数矩阵 A (使用稀疏矩阵存储以节省内存和计算量) % 每个内部点(i,j)对应方程涉及自身和上下左右四个邻居 % 我们使用五点差分格式 nx M-1; ny N-1; % 内部点数量 e ones(nx*ny, 1); % 主对角线元素 -2*(1/dx^2 1/dy^2) main_diag -2 * (1/dx^2 1/dy^2) * e; % 次对角线元素对应x方向的邻居 1/dx^2 % 注意在行优先排列下点(i,j)的左边邻居是向量索引 k-1右边邻居是 k1 % 但在矩阵A中这些非零元素的位置需要仔细计算。这里我们使用更简洁的方法 % 利用拉普拉斯算子的离散矩阵具有张量积结构 A kron(Iy, Dxx) kron(Dyy, Ix) % 其中 Dxx 和 Dyy 是一维二阶差分矩阵I是单位矩阵。 Dxx (1/dx^2) * spdiags([ones(nx,1), -2*ones(nx,1), ones(nx,1)], -1:1, nx, nx); Dyy (1/dy^2) * spdiags([ones(ny,1), -2*ones(ny,1), ones(ny,1)], -1:1, ny, ny); Ix speye(nx); Iy speye(ny); A kron(Iy, Dxx) kron(Dyy, Ix); % 这就是离散拉普拉斯算子的矩阵 % 求解线性方程组 A * U F_vec U_vec A \ F_vec; % 将解向量重塑回网格矩阵 U_inner reshape(U_vec, [ny, nx]); % 注意维度对应 % 将内部解嵌入到包含边界零值的完整网格中 U_full zeros(N1, M1); U_full(2:end-1, 2:end-1) U_inner; % 可视化 figure; surf(x, y, U_full, EdgeColor, none); xlabel(x); ylabel(y); zlabel(u(x,y)); title(有限差分法求解二维泊松方程);实操心得对于大规模网格如200x200以上直接使用反斜杠\求解可能内存不足或速度慢。此时应利用A是稀疏、对称正定的特性使用迭代法如共轭梯度法pcg或专门的PDE工具箱。上述代码中构建矩阵A的方法使用kron张量积是处理规则区域标准问题的优雅且高效的方式值得掌握。6. 调试技巧、常见问题与性能优化数值求解微分方程很少能一次成功总会遇到各种报错或不合理的结果。下面分享一些关键的调试经验和优化策略。6.1 ODE求解常见问题与排查错误“矩阵维度必须一致”或“索引超出范围”原因最可能是在定义ODE方程的函数f(t,y)中输出dydt不是列向量。务必检查dydt [dSdt; dIdt; dRdt]用的是分号列向量而非逗号或空格行向量。检查在函数末尾加一行size(dydt)确保输出是[n, 1]而不是[1, n]。错误“在时间t处失败无法满足积分容差”原因这是刚性问题的典型征兆或者方程在某个时间点出现了奇点如除以零。排查检查模型回顾方程是否存在当某个变量为0时分母为零的情况例如在SIR模型中如果总人口N设置为0。可以在函数中加入保护语句if N 0; dSdt0; ...; end。尝试刚性求解器用ode15s替换ode45看是否顺利求解。调整容差使用odeset放宽相对容差RelTol默认1e-3和绝对容差AbsTol默认1e-6。例如options odeset(RelTol, 1e-4, AbsTol, 1e-7);。注意放宽容差会降低精度。检查时间区间是否时间跨度太长导致解的变化尺度跨越多个数量级可以考虑分段求解。解的行为异常如出现负值、爆炸式增长原因可能是模型本身的不稳定性或者数值误差积累导致。排查验证模型检查方程和参数的单位、量纲是否合理。例如人口不应为负可以在ODE函数中对状态变量施加非负约束y(y0)0但这会改变方程需谨慎。减小时间步长通过设置odeset中的InitialStep和MaxStep来限制求解器的步长。例如options odeset(MaxStep, 0.1);。尝试不同求解器换用ode23或ode113看看结果是否一致。6.2 PDE求解常见问题与排查pdepe报错“尝试访问 xx(2)索引超出范围”原因几乎总是因为边界条件函数pdex1bc的输入输出变量数量不匹配。仔细检查函数定义行function [pl, ql, pr, qr] pdex1bc(xl, ul, xr, ur, t)确保输入是5个参数输出是4个参数且顺序正确。解出现非物理振荡或不稳定原因空间网格太粗无法分辨解的空间变化。特别是初始条件或源项有剧烈变化时。解决加密空间网格xmesh。同时对于对流占优的问题中心差分格式可能不稳定需要考虑迎风差分等格式但这已超出pdepe内置能力需要自己实现FDM。计算速度慢原因网格点太多或时间区间太长。优化减少输出点tspan中不要设置过于密集的输出时间点。求解器内部步长是自适应的tspan只控制输出。使用稀疏矩阵如果自己实现FDM矩阵A一定要用sparse或spdiags创建稀疏矩阵。利用对称性如果问题和边界条件是对称的可以只计算一半区域。6.3 性能与精度优化策略向量化编程在定义ODE/PDE的函数中尽量避免使用循环。MATLAB对矩阵和向量运算做了深度优化。例如在计算空间差分时使用矩阵运算代替逐点循环速度可提升数十倍。匿名函数与参数传递如前所述使用匿名函数(t,y) myode(t,y, param1, param2)来传递参数比使用全局变量更清晰、安全。预分配数组在需要存储时间序列结果时比如自己写时间推进的FDM循环务必预先分配好存储数组如U zeros(length(t), length(x))而不是在循环中动态增长数组。精度验证网格收敛性测试将空间网格点数加倍如从50到100时间容差减半比较两次求解结果在关心点上的差异。如果差异很小说明解已收敛。与已知解对比如果问题有解析解或高精度参考解务必进行对比这是检验代码正确性的黄金标准。守恒律检查对于某些物理问题总质量、总能量应该守恒。计算这些量的数值积分看其随时间的变化是否在可接受范围内。7. 在数学建模中的应用拓展与案例点睛掌握了ODE/PDE的求解技术最终要服务于数学建模。在比赛中这不仅仅是“求出解”那么简单。7.1 参数敏感性分析模型的结果往往依赖于参数。以SIR模型为例基本再生数R0 β/γ是关键参数。我们可以通过循环改变β或γ观察疫情峰值、达到时间、最终感染规模等指标如何变化。beta_range 0.1:0.05:0.5; gamma 0.1; peak_infected zeros(size(beta_range)); for i 1:length(beta_range) beta beta_range(i); [t, y] ode45((t,y) sir_ode(t,y,beta,gamma,N), tspan, y0); I y(:,2); peak_infected(i) max(I); end plot(beta_range, peak_infected, o-); xlabel(感染率 \beta); ylabel(疫情峰值感染人数); grid on;这种分析能告诉我们哪个参数对结果影响最大为干预措施如降低β提供定量依据。7.2 模型校准与参数估计当模型需要拟合实际数据时就变成了一个优化问题。例如我们有某地区每日新增感染数据I_data想要估计SIR模型中的β和γ。% 定义误差函数例如最小二乘 error_func (params) sum((simulate_sir(params) - I_data).^2); % params [beta, gamma] initial_guess [0.3, 0.1]; estimated_params fminsearch(error_func, initial_guess);其中simulate_sir(params)是一个封装好的函数用给定的params运行SIR模型并输出与I_data时间点对应的模拟感染人数。fminsearch是MATLAB的无导数优化函数可以用来寻找使误差最小的参数。7.3 耦合模型与多物理场问题真实的建模问题往往是多个过程耦合的。例如一个生态模型可能同时包含种群动力学ODE和空间扩散PDE即反应-扩散系统。这类问题通常需要自己构造数值方法如将PDE空间离散后与ODE部分结合成一个更大的ODE系统再用ode15s等求解或者使用更专业的工具箱如PDE Toolbox。这是数学建模的高阶挑战也是区分队伍水平的关键。7.4 结果的可视化与论文呈现一张好的图胜过千言万语。除了基本的二维线图、三维曲面图可以考虑动画用for循环和getframe制作PDE解随时间演化的动画在论文中提供动画截图或链接。热图用imagesc或pcolor展示二维场比surf图更简洁。参数空间扫描图用contourf或scatter展示不同参数组合下的结果分布。最后在论文中描述数值方法时不必赘述ode45或pdepe的内部算法但必须说明使用了什么求解器、为什么选择它如非刚性/刚性、设置了怎样的容差或网格、并进行了网格无关性验证以确保结果的可靠性。这体现了建模过程的严谨性。