尧图网站建设 
大家在科研和工程计算里遇到偏微分方程Partial Differential EquationPDE时常常会陷入两难数学推导太复杂、手写求解器太耗时而网上找到的 MATLAB 代码又往往只给片段没有完整思路。本文围绕“偏微分方程数值解”这个主题整理一份偏理论、偏代码、又偏实用的 MATLAB 教程。对于正在学习“数学实验”“数值分析”“计算物理”等课程的同学或者需要用热传导、波动、扩散等模型做仿真的工程师这篇文章都能直接提供可运行的 MATLAB 脚本和拆解思路。文章覆盖三种主流求解路线有限差分法手写实现、MATLAB 内置 pdepe 函数求解、PDE Toolbox 有限元求解。每一条路线都配了最小例子目的不是“炫技”而是让你在拿到一个偏微分方程后知道第一步写什么、第二步看什么、报错时查什么。先说一下本文的约定所有代码都按常见 MATLAB 版本编写如果你用的是 R2020a 之前的版本个别函数名可能需要轻微调整。示例中没有酷炫的 3D 动画但会把结果可视化和后续扩展的方向讲清楚。1. 偏微分方程数值解的基础概念1.1 什么是偏微分方程偏微分方程是包含未知多元函数及其偏导数的方程。相比常微分方程中未知函数只依赖一个自变量偏微分方程中的未知函数通常依赖时间 (t) 和空间坐标 (x, y, z) 中的若干个。常见形式如下热传导方程抛物型 [ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2} ]波动方程双曲型 [ \frac{\partial^2 u}{\partial t^2} c^2 \frac{\partial^2 u}{\partial x^2} ]拉普拉斯方程或泊松方程椭圆型 [ \frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2} f(x, y) ]这三个方程是很多 MATLAB 偏微分方程数值解教程的经典案例因为它们分别对应不同的物理背景扩散、振动、稳态场分布。1.2 为什么需要数值解偏微分方程的解析解只在极少数规则边界、简单初值条件下才能求出。实际工程中的几何形状、材料参数、边界载荷通常非常复杂继续“手推解析解”并不现实。所以数值解成为主流将连续问题离散成有限个未知量。用代数方程组近似替代微分方程。借助 MATLAB 矩阵运算能力快速求解。数值解并不等于近似很差。只要网格足够密、格式足够稳定计算精度可以满足工程需求。1.3 三大类主要解法在实际学习中接触最多的方法有三个方向方法核心思想MATLAB 实现方式优缺点有限差分法用差商代替导数手写矩阵和循环简单直观适合规则网格谱方法用全局基函数逼近解手写 FFT 或 Chebyshev 变换精度高适合光滑解有限元法将区域剖分并构造插值函数MATLAB PDE Toolbox适应复杂区域和边界条件pdepe 函数MATLAB 内置的 PDE 求解器直接调用函数适合一维抛物型/椭圆型问题下面从最容易理解的有限差分法开始。2. 环境准备与 MATLAB 版本说明2.1 基础环境要求编写偏微分方程数值解示例不需要额外购买工具箱只要安装 MATLAB 基础环境即可。如果后续要用 PDE Toolbox才需要确认许可证中已经包含该工具箱。版本方面可以这样理解本文示例基于较新的 MATLAB 版本编写R2020a 及以上版本基本都能直接运行。 函数如果存在版本差异我会在代码注释中说明替代写法。建议在学习时使用clear; close all; clc;开启干净的工作区。做数值实验时尽量减少变量污染。2.2 常用函数和脚本组织建议一个规范的 MATLAB PDE 数值实验项目建议这样组织目录heat_pde_demo/ ├── main_heat_explicit.m % 主脚本显式有限差分法 ├── pde_heat_ic.m % 初始条件函数必要时单独放 ├── pde_heat_bc.m % 边界条件函数 ├── plot_result.m % 绘图功能 └── readme.md % 说明文档初学者不必一开始就拆成多个文件。第一个完整示例可以只写一个脚本方便运行和调试。2.3 确认工具箱是否可用如果你要使用pdepe不需要安装额外工具箱它是 MATLAB 基础函数。如果你要调用solvepde、geometryFromEdges这些函数则需要 PDE Toolbox。可以使用ver命令查看ver % 或者用 which 检查函数是否存在 which solvepde如果输出显示路径正常就说明 PDE Toolbox 可用。如果返回built-in或实际安装路径也可以确认基础函数已经存在。3. 有限差分法求解热传导方程的 MATLAB 实现3.1 有限差分法的核心思路以一维热传导方程为例[ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2} ]设空间网格步长为 ( \Delta x )时间步长为 ( \Delta t )。常用的差分格式包括前向差分( \frac{u_i^{n1} - u_i^n}{\Delta t} \approx \frac{\partial u}{\partial t} )二阶中心差分( \frac{u_{i1}^n - 2u_i^n u_{i-1}^n}{\Delta x^2} \approx \frac{\partial^2 u}{\partial x^2} )显式格式可以写成[ u_i^{n1} u_i^n \alpha \frac{\Delta t}{\Delta x^2} \left( u_{i1}^n - 2u_i^n u_{i-1}^n \right) ]这个格式简单但必须满足稳定性条件[ r \frac{\alpha \Delta t}{\Delta x^2} \le 0.5 ]如果不满足计算结果会发散出现数值振荡。3.2 完整显式有限差分 MATLAB 脚本下面给出一个一维热传导方程的完整求解示例。物理场景是一根长度为 (L1) 的细杆两端温度恒为 0初始时刻中间区域温度为 1材料热扩散系数 ( \alpha 0.02 )。%% main_heat_explicit.m % 显式有限差分法求解一维热传导方程 % du/dt alpha * d2u/dx2 % 边界条件u(0,t)0, u(L,t)0 % 初始条件中间矩形分布 clear; close all; clc; % 参数设置 L 1.0; % 杆长 alpha 0.02; % 热扩散系数 nx 101; % 空间网格数网格点编号 1..101 dx L/(nx-1); % 空间步长 x linspace(0, L, nx); % 时间设置 T 1.0; % 总计算时间 dt 0.001; % 时间步长 nt round(T/dt); % 时间步数 r alpha * dt / dx^2; % 傅里叶数用于稳定性判断 fprintf(傅里叶数 r %.4f\n, r); if r 0.5 error(稳定性条件不满足r 0.5请减小 dt 或增大 dx); end % 初始条件 u zeros(nx,1); % 中间 0.4 到 0.6 区域初始温度为 1 u(x0.4 x0.6) 1; % 保存用于绘图 u_all zeros(nx, nt1); u_all(:,1) u; % 显式时间迭代 for n 1:nt % 内点更新 u_new u; for i 2:nx-1 u_new(i) u(i) r * (u(i1) - 2*u(i) u(i-1)); end % 边界条件直接置 0 u_new(1) 0; u_new(nx) 0; u u_new; u_all(:, n1) u; end % 绘制结果 figure; [X, Tgrid] meshgrid(0:nt, x); % 这里只需绘制三个时刻 plot_times [1, round(nt/5), round(nt/2), nt1]; colors lines(length(plot_times)); figure; hold on; for k 1:length(plot_times) idx plot_times(k); plot(x, u_all(:, idx), LineWidth, 1.5, Color, colors(k,:), ... DisplayName, sprintf(t %.3f, (idx-1)*dt)); end xlabel(x); ylabel(u(x,t)); title(显式有限差分法求解热传导方程); legend(show); grid on; hold off;对这个代码做几点解释稳定性参数r在这里非常重要。当r大于 0.5 时程序会直接报错防止算出无意义结果。时间循环中我们遍历每个内点使用上一时刻的u(i-1)、u(i)、u(i1)计算下一时刻这种格式称为显式格式。边界条件通过强制赋值实现。Dirichlet 边界条件下直接修改端点值是简单的处理方式。运行以上脚本后可以看到温度从中间矩形区域向两侧扩散最终整体温度趋于 0与无热源且两端恒温的物理直觉一致。3.3 隐式格式与 Crank-Nicolson 格式简介显式格式虽然简单但时间步长受限制。实际计算中更常用隐式格式因为它在参数上更稳定。对同一方程后向欧拉格式可以写为[ u_i^{n1} - r \left( u_{i1}^{n1} - 2u_i^{n1} u_{i-1}^{n1} \right) u_i^n ]每步需要解一个三对角方程组。MATLAB 可以高效用\\求解% 核心片段隐式后向欧拉 A eye(nx); for i 2:nx-1 A(i, i-1) -r; A(i, i) 1 2*r; A(i, i1) -r; end % 边界行不变对应 u(1)0, u(nx)0 A(1,:) 0; A(1,1) 1; A(nx,:) 0; A(nx,nx) 1; for n 1:nt b u; b(1) 0; b(nx) 0; u A \ b; end隐式格式的优势是不再受 (r \le 0.5) 的限制可以取更大时间步长适合长时间演化问题。缺点是每步需要解线性方程组不过对一维问题而言纯 MATLAB 矩阵求解速度仍然很快。4. 使用 MATLAB 内置 pdepe 函数求解偏微分方程4.1 pdepe 可以解什么MATLAB 提供的pdepe函数可以求解如下形式的一维偏微分方程组[ c\left(x,t,u,\frac{\partial u}{\partial x}\right) \frac{\partial u}{\partial t} x^{-m} \frac{\partial}{\partial x} \left( x^m f\left(x,t,u,\frac{\partial u}{\partial x}\right) \right) s\left(x,t,u,\frac{\partial u}{\partial x}\right) ]参数含义如下(m)问题的几何对称类型。(m0) 代表平板/直角坐标(m1) 代表柱对称(m2) 代表球对称。(c)时间导数项系数。(f)通量项。(s)源项。要使用pdepe一共需要提供三个函数pdefun定义方程中的系数函数。icfun定义初始条件。bcfun定义边界条件。调用形式如下sol pdepe(m, pdefun, icfun, bcfun, xmesh, tspan);sol是一个三维数组[ \text{sol}(i,j,k) ] 表示第 (k) 个因变量在时间 (t_i)、空间位置 (x_j) 处的值。如果只求一个变量会把所有时间层和空间层的解都存进来。4.2 pdepe 求解热传导方程完整案例仍然使用一维热传导方程[ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2} ]把它改写成pdepe能识别的标准形式[ 1 \cdot \frac{\partial u}{\partial t} \frac{\partial}{\partial x} \left( \alpha \frac{\partial u}{\partial x} \right) 0 ]因此得到(c1)(f\alpha \cdot \partial u / \partial x)(s0)完整代码如下%% main_pdepe_heat.m % 使用 pdepe 求解一维热传导方程 clear; close all; clc; % 参数 alpha 0.02; L 1.0; % 定义方程 function [c, f, s] heat_pdefun(x, t, u, dudx) c 1; f alpha * dudx; s 0; end % 初始条件 function u0 heat_icfun(x) if x 0.4 x 0.6 u0 1; else u0 0; end end % 边界条件 function [pl, ql, pr, qr] heat_bcfun(xl, ul, xr, ur, t) % 边界条件p q*f 0 % 对两端 u0相当于 pul-0、q0 pl ul; ql 0; pr ur; qr 0; end % 生成网格 xmesh linspace(0, L, 101); tspan linspace(0, 1, 101); % 求解 m 0; sol pdepe(m, heat_pdefun, heat_icfun, heat_bcfun, xmesh, tspan); % sol 维度为 length(tspan) x length(xmesh) x 1 u sol(:,:,1); % 绘图选取几个时间层 figure; hold on; plot_times [1, 21, 51, 101]; for idx plot_times plot(xmesh, u(idx,:), LineWidth, 1.5, ... DisplayName, sprintf(t %.2f, tspan(idx))); end xlabel(x); ylabel(u(x,t)); title(pdepe 求解热传导方程); legend(show); grid on; hold off;注意到一个小技术点在函数内部访问alpha时MATLAB 嵌套函数可以直接引用主函数工作区的变量。这里我把主函数和子函数放在同一个main_pdepe_heat.m文件中利用function ... end的局部函数特性比较适合初学者阅读。如果写成单独的.m文件则需要把alpha、L作为全局变量或通过参数传递。pdepe 会自动处理时间和空间离散用户不需要手动考虑稳定性条件。但一定要保证边界条件写法正确否则结果会出现奇怪振荡。4.3 pdepe 求解结果的输出技巧sol的维度容易让人困惑。这里用代码验证disp(size(sol)); % 如果要得到所有时间层在 x0.5 处的温度曲线 u_mid sol(:, 51); figure; plot(tspan, u_mid); xlabel(t); ylabel(u(t, x0.5)); title(x0.5 处温度随时间变化);如果换了一套网格x0.5 不一定正好落在节点上。更稳妥的做法是先找最近索引[~, idx] min(abs(xmesh - 0.5)); u_mid sol(:, idx);在处理数值解时这种“先找索引再取值”的写法比硬编码节点编号更安全。5. 使用 PDE Toolbox 求解二维稳态问题5.1 PDE Toolbox 是什么当问题推广到二维或三维尤其边界形状不是矩形时有限差分法的手写成本快速上升。MATLAB 的 PDE Toolbox偏微分方程工具箱采用有限元方法可以处理更复杂的几何区域。PDE Toolbox 的基本流程可以概括为创建几何模型createpde。建立几何形状geometryFromEdges或从 STL 导入三维几何。指定方程系数和边界条件。生成网格generateMesh。求解并可视化。这里用一个经典二维拉普拉斯方程案例来说明。问题为在一个边长为 1 的正方形区域内求解泊松方程[ -\nabla^2 u 1 ]边界条件设为 (u0)。这是一个有源场问题相当于在均匀边界接地、内部均匀激励的物理模型。5.2 PDE Toolbox 完整示例%% main_pdetoolbox_poisson.m clear; close all; clc; % 1. 创建 PDE 模型 model createpde(); % 2. 创建几何对象单位正方形 g geometryFromEdges(model, squareg); % 3. 指定方程类型系数形式 % solvepde 默认求解 -div(c*grad(u)) a*u f specifyCoefficients(model, m, 0, d, 0, c, 1, a, 0, f, 1); % 4. 边界条件u 0 applyBoundaryCondition(model, dirichlet, Edge, 1:model.Geometry.NumEdges, u, 0); % 5. 生成网格 generateMesh(model, Hmax, 0.05); % 6. 求解 results solvepde(model); % 7. 查看结果 u results.NodalSolution; figure; pdeplot(model, XYData, u, ZData, u, ColorBar, on); title(PDE Toolbox 求解 -nabla^2 u 1); xlabel(x); ylabel(y);说明一点squareg是 PDE Toolbox 自带的单位正方形几何函数。如果你导入自定义多边形可以改用更灵活的方式例如gd [2; 4; 0; 1; 1; 0; 0; 0; 1; 1]; % 四边形坐标描述 g decsg(gd, R1, char(R1)); geometryFromEdges(model, g);这种底层几何描述对新手不太友好。更推荐使用pdegplot或pdeModeler工具先绘图确认边界结构。5.3 PDE Toolbox 求解三维问题扩展PDE Toolbox 也支持三维几何例如从 STL 文件导入模型并做热传导仿真。核心代码结构类似只是把geometryFromEdges替换为importGeometry(model, your_model.stl);然后使用generateMesh和solvepde。三维问题真正难点通常不是求解器而是几何模型是否闭合、材料参数是否合理、网格数量是否在内存允许范围内。初学者建议先使用内置几何函数比如multicuboid创建长方体组合模型。6. 偏微分方程数值解的常用边界条件处理6.1 Dirichlet、Neumann、Robin 边界条件在 MATLAB 中书写偏微分方程时边界条件的处理直接决定结果是否正确。边界条件类型数学表达式物理含义pdepe 写法Dirichlet(u g)边界值固定(p u - g), (q0)Neumann(\partial u / \partial n g)边界通量给定(p -g), (q1) 或按对应形式Robin(a u b \partial u / \partial n g)边界导热等写成 (pq f0) 的组合形式在 pdepe 边界条件书写中要求满足[ p q \cdot f 0 ]其中 (f) 就是方程中的通量函数。以热传导方程为例如果左端绝热意味着 (\partial u / \partial x 0)那么取 (p0)、(q1)这样[ 0 1 \cdot f f \alpha \frac{\partial u}{\partial x} 0 ]正好符合绝热条件。如果右端是恒温 (u1)则取 (pu-1)、(q0)。初学者最容易犯的错误是把 Neumann 边界理解成直接使用导数赋值但在 pdepe 格式中要写成通量条件必须知道 (f) 的表达式。不同类型方程之间这种写法可能略有差异。6.2 初始条件与边界条件的协调设定初始条件时还有个很容易忽略的问题初始条件最好和边界条件协调一致否则在 (t0) 附近可能出现非物理的突变。比如在前面的案例中初始时刻 (x\in[0.4,0.6]) 的温度为 1而边界温度恒为 0。虽然边界与初始高温区域没有直接重叠但离散节点在第一时刻就会经历快速变化这会导致初始几个时间层斜率较大。如果要测数值精度可以改成更光滑的初始条件例如高斯分布sigma 0.05; u0 exp(-((x-0.5).^2)/(2*sigma^2));这样在物理上更自然时间演化也更平滑。7. 结果可视化与精度对比7.1 一维结果作图MATLAB 中常用的 PDE 结果作图有三种plot(x, u)某时刻曲线。surf(x, t, u)或pcolor全部时空分布。contourf二维等值线图。对于热传导方程时空分布图最直观figure; [Tgrid, Xgrid] meshgrid(tspan, xmesh); surf(Xgrid, Tgrid, u); xlabel(x); ylabel(t); zlabel(u(x,t)); shading interp; title(热传导方程时空分布);注意surf要求u的行列与坐标矩阵匹配。上面的代码中用转置来处理。7.2 验证程序是否写对的基本方法对偏微分方程数值解有一个通用的“验算三件套”检查守恒量或单调趋势是否合理。检查稳态解是否和解析解一致。检查减小步长后结果是否变化很小。以热传导方程为例当 (t) 足够大时解会趋于边界条件决定的稳态。如果边界条件为零则最终全场接近 0。如果做的是绝热边界则总热量应保持不变可以计算total_heat trapz(x, u);不同时刻的总热量变化不应太大。如果热量明显流失说明边界条件可能写错。8. 常见报错与排查思路偏微分方程数值解的 MATLAB 程序中报错信息和你预想的可能不太一样。下面列几个频率较高的现象。问题现象常见原因解决思路显式差分结果直接变成 NaN 或 Inf不满足稳定性条件即 (r0.5)减小 (dt) 或增大 (dx)检查代码中r是否计算正确pdepe求解报错“空间离散化失败”方程系数函数格式写错或边界条件返回了错误数组维度打印c、f、s的维度确认返回列向量检查p和q是否与输入点一一对应pdepe结果剧烈振荡时间网格太粗或初边界条件不协调加密tspan或改用更平滑的初始条件Undefined function solvepde当前 MATLAB 没有安装 PDE Toolbox使用ver检查工具箱或者改用pdepe和有限差分法geometryFromEdges报错几何描述矩阵格式不对先用pdegplot画几何确认 Edge 编号正确结果看起来“一直不变”总时间太短或边界与初始状态已经平衡增加 T或设置tspan跨度更大8.1 pdepe 隐式边界条件错误示例错误代码function [pl, ql, pr, qr] heat_bcfun(xl, ul, xr, ur, t) pl 0; ql 1; pr ul - 1; % 类型错误右端条件被写成了左端变量 qr 0; end这里右端温度应为常数 1但pr却使用了ul导致边界条件与空间位置无关时仍可能造成奇怪的耦合。正确写法是function [pl, ql, pr, qr] heat_bcfun(xl, ul, xr, ur, t) pl ul; ql 0; pr ur - 1; qr 0; end所以排查边界条件时一个技巧是如果方程只有一个因变量先检查pl是否只依赖xl和ulpr是否只依赖xr和ur。8.2 显式差分中索引越界很多新手的循环会写成for i 2:nx u_new(i) u(i) r * (u(i1) - 2*u(i) u(i-1)); end当i nx时会访问u(nx1)越界报错。这类问题并不难修关键是写循环前先画一条网格索引轴线明确边界点是 1 和 nx内点是 2 到 nx-1。9. 最佳实践与工程建议9.1 尽量使用无量纲化模型实际工程问题中长度、时间、温度等物理量可能相差多个数量级。直接带入 SI 单位可能导致矩阵条件数很大收敛变慢。建议先做无量纲化比如用特征长度去除所有空间坐标用特征时间去除所有时间项。整理后再用 MATLAB 处理你会更容易判断网格和时间步长是否合理。9.2 把“数值实验”和“正式仿真”分开写代码时可以把脚本设计成两块第一部分手算/已知解析解做验证。第二部分在自己关心的参数范围内做正式计算。这种分离能帮你尽早发现离散实现中的 bug。9.3 涉及几何模型时先可视化再求解使用 PDE Toolbox 时先执行pdegplot(model, EdgeLabels, on)查看边界编号。如果边界条件应用在错误的边上计算不会报错但结果完全不正确。边界条件不是“报错友好”的错误所以一定要先可视化。9.4 大批量参数扫描时避免重复生成网格在做一个参数的扫描实验时网格可以保持不变只需要重复修改系数或边界条件。PDE Toolbox 中generateMesh只需要执行一次后续用solvepde(model)反复求解即可。如果每次都重新生成网格耗时可能翻倍而且如果网格变化比较不同参数的结果时也会引入额外误差。9.5 保存关键状态与结果数值实验经常需要对比多组参数。建议把关键结果保存为.mat文件save(result_alpha002.mat, x, tspan, u);后期绘图时单独写一个脚本加载数据而不是重新运行整个求解。这在大规模仿真中能节省大量时间。9.6 真实项目中的边界条件安全边界如果程序涉及真实物理设备的仿真必须遵循两条原则单元测试先行用已知解析解验证代码。参数合法性检查材料参数如热扩散系数、导热系数不能为负边界温度在合理范围内。仿真程序不报错不代表结果有意义。10. 下一步学习路线与总结数学上看偏微分方程数值解是计算数学的核心主题工程上看它又是从物理模型到软件仿真的必经环节。MATLAB 把很多底层逻辑封装得非常方便但这并不意味着可以忽略基本概念。学完本文的显式差分和pdepe后建议按下面顺序继续深入。第一掌握稳定性分析。你会慢慢发现显式格式的r 0.5并不是数值计算的经验法则而是来源于 Von Neumann 稳定性分析。理解这个概念后你会更从容地应对各种非线性问题。第二学习如何将一维问题推广到二维。有限差分在二维矩形区域仍然容易实现但要注意交替方向隐式格式ADI等更高效的解法。如果区域不规则再去学习 PDE Toolbox 的网格生成和有限元组装原理。第三重视误差分析和收敛阶验证。数值解是否可信不只是“图看起来像”就行。通过比较不同网格步长下的数值结果估计收敛阶是数值方法实验中非常有价值的技能。第四把问题上升到“计算流程”。真实项目不只是运行一段代码而是包含几何建模、网格生成、求解、后处理、优化迭代整个流程。MATLAB Live Script 也是很好的记录工具可以把文字、公式、代码放在同一个文档里方便梳理思路和写实验报告。希望这篇偏微分方程数值解 MATLAB 教程帮到正在学习或准备做数值仿真的你。学 PDE 数值解关键不在于一次性理解所有理论而在于把最简单的热传导方程真正调通再去挑战更复杂的耦合方程组。 偏微分方程数值解MATLAB教程：差分法、pdepe与PDE Toolbox指南
大家在科研和工程计算里遇到偏微分方程Partial Differential EquationPDE时常常会陷入两难数学推导太复杂、手写求解器太耗时而网上找到的 MATLAB 代码又往往只给片段没有完整思路。本文围绕“偏微分方程数值解”这个主题整理一份偏理论、偏代码、又偏实用的 MATLAB 教程。对于正在学习“数学实验”“数值分析”“计算物理”等课程的同学或者需要用热传导、波动、扩散等模型做仿真的工程师这篇文章都能直接提供可运行的 MATLAB 脚本和拆解思路。文章覆盖三种主流求解路线有限差分法手写实现、MATLAB 内置 pdepe 函数求解、PDE Toolbox 有限元求解。每一条路线都配了最小例子目的不是“炫技”而是让你在拿到一个偏微分方程后知道第一步写什么、第二步看什么、报错时查什么。先说一下本文的约定所有代码都按常见 MATLAB 版本编写如果你用的是 R2020a 之前的版本个别函数名可能需要轻微调整。示例中没有酷炫的 3D 动画但会把结果可视化和后续扩展的方向讲清楚。1. 偏微分方程数值解的基础概念1.1 什么是偏微分方程偏微分方程是包含未知多元函数及其偏导数的方程。相比常微分方程中未知函数只依赖一个自变量偏微分方程中的未知函数通常依赖时间 (t) 和空间坐标 (x, y, z) 中的若干个。常见形式如下热传导方程抛物型 [ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2} ]波动方程双曲型 [ \frac{\partial^2 u}{\partial t^2} c^2 \frac{\partial^2 u}{\partial x^2} ]拉普拉斯方程或泊松方程椭圆型 [ \frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2} f(x, y) ]这三个方程是很多 MATLAB 偏微分方程数值解教程的经典案例因为它们分别对应不同的物理背景扩散、振动、稳态场分布。1.2 为什么需要数值解偏微分方程的解析解只在极少数规则边界、简单初值条件下才能求出。实际工程中的几何形状、材料参数、边界载荷通常非常复杂继续“手推解析解”并不现实。所以数值解成为主流将连续问题离散成有限个未知量。用代数方程组近似替代微分方程。借助 MATLAB 矩阵运算能力快速求解。数值解并不等于近似很差。只要网格足够密、格式足够稳定计算精度可以满足工程需求。1.3 三大类主要解法在实际学习中接触最多的方法有三个方向方法核心思想MATLAB 实现方式优缺点有限差分法用差商代替导数手写矩阵和循环简单直观适合规则网格谱方法用全局基函数逼近解手写 FFT 或 Chebyshev 变换精度高适合光滑解有限元法将区域剖分并构造插值函数MATLAB PDE Toolbox适应复杂区域和边界条件pdepe 函数MATLAB 内置的 PDE 求解器直接调用函数适合一维抛物型/椭圆型问题下面从最容易理解的有限差分法开始。2. 环境准备与 MATLAB 版本说明2.1 基础环境要求编写偏微分方程数值解示例不需要额外购买工具箱只要安装 MATLAB 基础环境即可。如果后续要用 PDE Toolbox才需要确认许可证中已经包含该工具箱。版本方面可以这样理解本文示例基于较新的 MATLAB 版本编写R2020a 及以上版本基本都能直接运行。 函数如果存在版本差异我会在代码注释中说明替代写法。建议在学习时使用clear; close all; clc;开启干净的工作区。做数值实验时尽量减少变量污染。2.2 常用函数和脚本组织建议一个规范的 MATLAB PDE 数值实验项目建议这样组织目录heat_pde_demo/ ├── main_heat_explicit.m % 主脚本显式有限差分法 ├── pde_heat_ic.m % 初始条件函数必要时单独放 ├── pde_heat_bc.m % 边界条件函数 ├── plot_result.m % 绘图功能 └── readme.md % 说明文档初学者不必一开始就拆成多个文件。第一个完整示例可以只写一个脚本方便运行和调试。2.3 确认工具箱是否可用如果你要使用pdepe不需要安装额外工具箱它是 MATLAB 基础函数。如果你要调用solvepde、geometryFromEdges这些函数则需要 PDE Toolbox。可以使用ver命令查看ver % 或者用 which 检查函数是否存在 which solvepde如果输出显示路径正常就说明 PDE Toolbox 可用。如果返回built-in或实际安装路径也可以确认基础函数已经存在。3. 有限差分法求解热传导方程的 MATLAB 实现3.1 有限差分法的核心思路以一维热传导方程为例[ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2} ]设空间网格步长为 ( \Delta x )时间步长为 ( \Delta t )。常用的差分格式包括前向差分( \frac{u_i^{n1} - u_i^n}{\Delta t} \approx \frac{\partial u}{\partial t} )二阶中心差分( \frac{u_{i1}^n - 2u_i^n u_{i-1}^n}{\Delta x^2} \approx \frac{\partial^2 u}{\partial x^2} )显式格式可以写成[ u_i^{n1} u_i^n \alpha \frac{\Delta t}{\Delta x^2} \left( u_{i1}^n - 2u_i^n u_{i-1}^n \right) ]这个格式简单但必须满足稳定性条件[ r \frac{\alpha \Delta t}{\Delta x^2} \le 0.5 ]如果不满足计算结果会发散出现数值振荡。3.2 完整显式有限差分 MATLAB 脚本下面给出一个一维热传导方程的完整求解示例。物理场景是一根长度为 (L1) 的细杆两端温度恒为 0初始时刻中间区域温度为 1材料热扩散系数 ( \alpha 0.02 )。%% main_heat_explicit.m % 显式有限差分法求解一维热传导方程 % du/dt alpha * d2u/dx2 % 边界条件u(0,t)0, u(L,t)0 % 初始条件中间矩形分布 clear; close all; clc; % 参数设置 L 1.0; % 杆长 alpha 0.02; % 热扩散系数 nx 101; % 空间网格数网格点编号 1..101 dx L/(nx-1); % 空间步长 x linspace(0, L, nx); % 时间设置 T 1.0; % 总计算时间 dt 0.001; % 时间步长 nt round(T/dt); % 时间步数 r alpha * dt / dx^2; % 傅里叶数用于稳定性判断 fprintf(傅里叶数 r %.4f\n, r); if r 0.5 error(稳定性条件不满足r 0.5请减小 dt 或增大 dx); end % 初始条件 u zeros(nx,1); % 中间 0.4 到 0.6 区域初始温度为 1 u(x0.4 x0.6) 1; % 保存用于绘图 u_all zeros(nx, nt1); u_all(:,1) u; % 显式时间迭代 for n 1:nt % 内点更新 u_new u; for i 2:nx-1 u_new(i) u(i) r * (u(i1) - 2*u(i) u(i-1)); end % 边界条件直接置 0 u_new(1) 0; u_new(nx) 0; u u_new; u_all(:, n1) u; end % 绘制结果 figure; [X, Tgrid] meshgrid(0:nt, x); % 这里只需绘制三个时刻 plot_times [1, round(nt/5), round(nt/2), nt1]; colors lines(length(plot_times)); figure; hold on; for k 1:length(plot_times) idx plot_times(k); plot(x, u_all(:, idx), LineWidth, 1.5, Color, colors(k,:), ... DisplayName, sprintf(t %.3f, (idx-1)*dt)); end xlabel(x); ylabel(u(x,t)); title(显式有限差分法求解热传导方程); legend(show); grid on; hold off;对这个代码做几点解释稳定性参数r在这里非常重要。当r大于 0.5 时程序会直接报错防止算出无意义结果。时间循环中我们遍历每个内点使用上一时刻的u(i-1)、u(i)、u(i1)计算下一时刻这种格式称为显式格式。边界条件通过强制赋值实现。Dirichlet 边界条件下直接修改端点值是简单的处理方式。运行以上脚本后可以看到温度从中间矩形区域向两侧扩散最终整体温度趋于 0与无热源且两端恒温的物理直觉一致。3.3 隐式格式与 Crank-Nicolson 格式简介显式格式虽然简单但时间步长受限制。实际计算中更常用隐式格式因为它在参数上更稳定。对同一方程后向欧拉格式可以写为[ u_i^{n1} - r \left( u_{i1}^{n1} - 2u_i^{n1} u_{i-1}^{n1} \right) u_i^n ]每步需要解一个三对角方程组。MATLAB 可以高效用\\求解% 核心片段隐式后向欧拉 A eye(nx); for i 2:nx-1 A(i, i-1) -r; A(i, i) 1 2*r; A(i, i1) -r; end % 边界行不变对应 u(1)0, u(nx)0 A(1,:) 0; A(1,1) 1; A(nx,:) 0; A(nx,nx) 1; for n 1:nt b u; b(1) 0; b(nx) 0; u A \ b; end隐式格式的优势是不再受 (r \le 0.5) 的限制可以取更大时间步长适合长时间演化问题。缺点是每步需要解线性方程组不过对一维问题而言纯 MATLAB 矩阵求解速度仍然很快。4. 使用 MATLAB 内置 pdepe 函数求解偏微分方程4.1 pdepe 可以解什么MATLAB 提供的pdepe函数可以求解如下形式的一维偏微分方程组[ c\left(x,t,u,\frac{\partial u}{\partial x}\right) \frac{\partial u}{\partial t} x^{-m} \frac{\partial}{\partial x} \left( x^m f\left(x,t,u,\frac{\partial u}{\partial x}\right) \right) s\left(x,t,u,\frac{\partial u}{\partial x}\right) ]参数含义如下(m)问题的几何对称类型。(m0) 代表平板/直角坐标(m1) 代表柱对称(m2) 代表球对称。(c)时间导数项系数。(f)通量项。(s)源项。要使用pdepe一共需要提供三个函数pdefun定义方程中的系数函数。icfun定义初始条件。bcfun定义边界条件。调用形式如下sol pdepe(m, pdefun, icfun, bcfun, xmesh, tspan);sol是一个三维数组[ \text{sol}(i,j,k) ] 表示第 (k) 个因变量在时间 (t_i)、空间位置 (x_j) 处的值。如果只求一个变量会把所有时间层和空间层的解都存进来。4.2 pdepe 求解热传导方程完整案例仍然使用一维热传导方程[ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2} ]把它改写成pdepe能识别的标准形式[ 1 \cdot \frac{\partial u}{\partial t} \frac{\partial}{\partial x} \left( \alpha \frac{\partial u}{\partial x} \right) 0 ]因此得到(c1)(f\alpha \cdot \partial u / \partial x)(s0)完整代码如下%% main_pdepe_heat.m % 使用 pdepe 求解一维热传导方程 clear; close all; clc; % 参数 alpha 0.02; L 1.0; % 定义方程 function [c, f, s] heat_pdefun(x, t, u, dudx) c 1; f alpha * dudx; s 0; end % 初始条件 function u0 heat_icfun(x) if x 0.4 x 0.6 u0 1; else u0 0; end end % 边界条件 function [pl, ql, pr, qr] heat_bcfun(xl, ul, xr, ur, t) % 边界条件p q*f 0 % 对两端 u0相当于 pul-0、q0 pl ul; ql 0; pr ur; qr 0; end % 生成网格 xmesh linspace(0, L, 101); tspan linspace(0, 1, 101); % 求解 m 0; sol pdepe(m, heat_pdefun, heat_icfun, heat_bcfun, xmesh, tspan); % sol 维度为 length(tspan) x length(xmesh) x 1 u sol(:,:,1); % 绘图选取几个时间层 figure; hold on; plot_times [1, 21, 51, 101]; for idx plot_times plot(xmesh, u(idx,:), LineWidth, 1.5, ... DisplayName, sprintf(t %.2f, tspan(idx))); end xlabel(x); ylabel(u(x,t)); title(pdepe 求解热传导方程); legend(show); grid on; hold off;注意到一个小技术点在函数内部访问alpha时MATLAB 嵌套函数可以直接引用主函数工作区的变量。这里我把主函数和子函数放在同一个main_pdepe_heat.m文件中利用function ... end的局部函数特性比较适合初学者阅读。如果写成单独的.m文件则需要把alpha、L作为全局变量或通过参数传递。pdepe 会自动处理时间和空间离散用户不需要手动考虑稳定性条件。但一定要保证边界条件写法正确否则结果会出现奇怪振荡。4.3 pdepe 求解结果的输出技巧sol的维度容易让人困惑。这里用代码验证disp(size(sol)); % 如果要得到所有时间层在 x0.5 处的温度曲线 u_mid sol(:, 51); figure; plot(tspan, u_mid); xlabel(t); ylabel(u(t, x0.5)); title(x0.5 处温度随时间变化);如果换了一套网格x0.5 不一定正好落在节点上。更稳妥的做法是先找最近索引[~, idx] min(abs(xmesh - 0.5)); u_mid sol(:, idx);在处理数值解时这种“先找索引再取值”的写法比硬编码节点编号更安全。5. 使用 PDE Toolbox 求解二维稳态问题5.1 PDE Toolbox 是什么当问题推广到二维或三维尤其边界形状不是矩形时有限差分法的手写成本快速上升。MATLAB 的 PDE Toolbox偏微分方程工具箱采用有限元方法可以处理更复杂的几何区域。PDE Toolbox 的基本流程可以概括为创建几何模型createpde。建立几何形状geometryFromEdges或从 STL 导入三维几何。指定方程系数和边界条件。生成网格generateMesh。求解并可视化。这里用一个经典二维拉普拉斯方程案例来说明。问题为在一个边长为 1 的正方形区域内求解泊松方程[ -\nabla^2 u 1 ]边界条件设为 (u0)。这是一个有源场问题相当于在均匀边界接地、内部均匀激励的物理模型。5.2 PDE Toolbox 完整示例%% main_pdetoolbox_poisson.m clear; close all; clc; % 1. 创建 PDE 模型 model createpde(); % 2. 创建几何对象单位正方形 g geometryFromEdges(model, squareg); % 3. 指定方程类型系数形式 % solvepde 默认求解 -div(c*grad(u)) a*u f specifyCoefficients(model, m, 0, d, 0, c, 1, a, 0, f, 1); % 4. 边界条件u 0 applyBoundaryCondition(model, dirichlet, Edge, 1:model.Geometry.NumEdges, u, 0); % 5. 生成网格 generateMesh(model, Hmax, 0.05); % 6. 求解 results solvepde(model); % 7. 查看结果 u results.NodalSolution; figure; pdeplot(model, XYData, u, ZData, u, ColorBar, on); title(PDE Toolbox 求解 -nabla^2 u 1); xlabel(x); ylabel(y);说明一点squareg是 PDE Toolbox 自带的单位正方形几何函数。如果你导入自定义多边形可以改用更灵活的方式例如gd [2; 4; 0; 1; 1; 0; 0; 0; 1; 1]; % 四边形坐标描述 g decsg(gd, R1, char(R1)); geometryFromEdges(model, g);这种底层几何描述对新手不太友好。更推荐使用pdegplot或pdeModeler工具先绘图确认边界结构。5.3 PDE Toolbox 求解三维问题扩展PDE Toolbox 也支持三维几何例如从 STL 文件导入模型并做热传导仿真。核心代码结构类似只是把geometryFromEdges替换为importGeometry(model, your_model.stl);然后使用generateMesh和solvepde。三维问题真正难点通常不是求解器而是几何模型是否闭合、材料参数是否合理、网格数量是否在内存允许范围内。初学者建议先使用内置几何函数比如multicuboid创建长方体组合模型。6. 偏微分方程数值解的常用边界条件处理6.1 Dirichlet、Neumann、Robin 边界条件在 MATLAB 中书写偏微分方程时边界条件的处理直接决定结果是否正确。边界条件类型数学表达式物理含义pdepe 写法Dirichlet(u g)边界值固定(p u - g), (q0)Neumann(\partial u / \partial n g)边界通量给定(p -g), (q1) 或按对应形式Robin(a u b \partial u / \partial n g)边界导热等写成 (pq f0) 的组合形式在 pdepe 边界条件书写中要求满足[ p q \cdot f 0 ]其中 (f) 就是方程中的通量函数。以热传导方程为例如果左端绝热意味着 (\partial u / \partial x 0)那么取 (p0)、(q1)这样[ 0 1 \cdot f f \alpha \frac{\partial u}{\partial x} 0 ]正好符合绝热条件。如果右端是恒温 (u1)则取 (pu-1)、(q0)。初学者最容易犯的错误是把 Neumann 边界理解成直接使用导数赋值但在 pdepe 格式中要写成通量条件必须知道 (f) 的表达式。不同类型方程之间这种写法可能略有差异。6.2 初始条件与边界条件的协调设定初始条件时还有个很容易忽略的问题初始条件最好和边界条件协调一致否则在 (t0) 附近可能出现非物理的突变。比如在前面的案例中初始时刻 (x\in[0.4,0.6]) 的温度为 1而边界温度恒为 0。虽然边界与初始高温区域没有直接重叠但离散节点在第一时刻就会经历快速变化这会导致初始几个时间层斜率较大。如果要测数值精度可以改成更光滑的初始条件例如高斯分布sigma 0.05; u0 exp(-((x-0.5).^2)/(2*sigma^2));这样在物理上更自然时间演化也更平滑。7. 结果可视化与精度对比7.1 一维结果作图MATLAB 中常用的 PDE 结果作图有三种plot(x, u)某时刻曲线。surf(x, t, u)或pcolor全部时空分布。contourf二维等值线图。对于热传导方程时空分布图最直观figure; [Tgrid, Xgrid] meshgrid(tspan, xmesh); surf(Xgrid, Tgrid, u); xlabel(x); ylabel(t); zlabel(u(x,t)); shading interp; title(热传导方程时空分布);注意surf要求u的行列与坐标矩阵匹配。上面的代码中用转置来处理。7.2 验证程序是否写对的基本方法对偏微分方程数值解有一个通用的“验算三件套”检查守恒量或单调趋势是否合理。检查稳态解是否和解析解一致。检查减小步长后结果是否变化很小。以热传导方程为例当 (t) 足够大时解会趋于边界条件决定的稳态。如果边界条件为零则最终全场接近 0。如果做的是绝热边界则总热量应保持不变可以计算total_heat trapz(x, u);不同时刻的总热量变化不应太大。如果热量明显流失说明边界条件可能写错。8. 常见报错与排查思路偏微分方程数值解的 MATLAB 程序中报错信息和你预想的可能不太一样。下面列几个频率较高的现象。问题现象常见原因解决思路显式差分结果直接变成 NaN 或 Inf不满足稳定性条件即 (r0.5)减小 (dt) 或增大 (dx)检查代码中r是否计算正确pdepe求解报错“空间离散化失败”方程系数函数格式写错或边界条件返回了错误数组维度打印c、f、s的维度确认返回列向量检查p和q是否与输入点一一对应pdepe结果剧烈振荡时间网格太粗或初边界条件不协调加密tspan或改用更平滑的初始条件Undefined function solvepde当前 MATLAB 没有安装 PDE Toolbox使用ver检查工具箱或者改用pdepe和有限差分法geometryFromEdges报错几何描述矩阵格式不对先用pdegplot画几何确认 Edge 编号正确结果看起来“一直不变”总时间太短或边界与初始状态已经平衡增加 T或设置tspan跨度更大8.1 pdepe 隐式边界条件错误示例错误代码function [pl, ql, pr, qr] heat_bcfun(xl, ul, xr, ur, t) pl 0; ql 1; pr ul - 1; % 类型错误右端条件被写成了左端变量 qr 0; end这里右端温度应为常数 1但pr却使用了ul导致边界条件与空间位置无关时仍可能造成奇怪的耦合。正确写法是function [pl, ql, pr, qr] heat_bcfun(xl, ul, xr, ur, t) pl ul; ql 0; pr ur - 1; qr 0; end所以排查边界条件时一个技巧是如果方程只有一个因变量先检查pl是否只依赖xl和ulpr是否只依赖xr和ur。8.2 显式差分中索引越界很多新手的循环会写成for i 2:nx u_new(i) u(i) r * (u(i1) - 2*u(i) u(i-1)); end当i nx时会访问u(nx1)越界报错。这类问题并不难修关键是写循环前先画一条网格索引轴线明确边界点是 1 和 nx内点是 2 到 nx-1。9. 最佳实践与工程建议9.1 尽量使用无量纲化模型实际工程问题中长度、时间、温度等物理量可能相差多个数量级。直接带入 SI 单位可能导致矩阵条件数很大收敛变慢。建议先做无量纲化比如用特征长度去除所有空间坐标用特征时间去除所有时间项。整理后再用 MATLAB 处理你会更容易判断网格和时间步长是否合理。9.2 把“数值实验”和“正式仿真”分开写代码时可以把脚本设计成两块第一部分手算/已知解析解做验证。第二部分在自己关心的参数范围内做正式计算。这种分离能帮你尽早发现离散实现中的 bug。9.3 涉及几何模型时先可视化再求解使用 PDE Toolbox 时先执行pdegplot(model, EdgeLabels, on)查看边界编号。如果边界条件应用在错误的边上计算不会报错但结果完全不正确。边界条件不是“报错友好”的错误所以一定要先可视化。9.4 大批量参数扫描时避免重复生成网格在做一个参数的扫描实验时网格可以保持不变只需要重复修改系数或边界条件。PDE Toolbox 中generateMesh只需要执行一次后续用solvepde(model)反复求解即可。如果每次都重新生成网格耗时可能翻倍而且如果网格变化比较不同参数的结果时也会引入额外误差。9.5 保存关键状态与结果数值实验经常需要对比多组参数。建议把关键结果保存为.mat文件save(result_alpha002.mat, x, tspan, u);后期绘图时单独写一个脚本加载数据而不是重新运行整个求解。这在大规模仿真中能节省大量时间。9.6 真实项目中的边界条件安全边界如果程序涉及真实物理设备的仿真必须遵循两条原则单元测试先行用已知解析解验证代码。参数合法性检查材料参数如热扩散系数、导热系数不能为负边界温度在合理范围内。仿真程序不报错不代表结果有意义。10. 下一步学习路线与总结数学上看偏微分方程数值解是计算数学的核心主题工程上看它又是从物理模型到软件仿真的必经环节。MATLAB 把很多底层逻辑封装得非常方便但这并不意味着可以忽略基本概念。学完本文的显式差分和pdepe后建议按下面顺序继续深入。第一掌握稳定性分析。你会慢慢发现显式格式的r 0.5并不是数值计算的经验法则而是来源于 Von Neumann 稳定性分析。理解这个概念后你会更从容地应对各种非线性问题。第二学习如何将一维问题推广到二维。有限差分在二维矩形区域仍然容易实现但要注意交替方向隐式格式ADI等更高效的解法。如果区域不规则再去学习 PDE Toolbox 的网格生成和有限元组装原理。第三重视误差分析和收敛阶验证。数值解是否可信不只是“图看起来像”就行。通过比较不同网格步长下的数值结果估计收敛阶是数值方法实验中非常有价值的技能。第四把问题上升到“计算流程”。真实项目不只是运行一段代码而是包含几何建模、网格生成、求解、后处理、优化迭代整个流程。MATLAB Live Script 也是很好的记录工具可以把文字、公式、代码放在同一个文档里方便梳理思路和写实验报告。希望这篇偏微分方程数值解 MATLAB 教程帮到正在学习或准备做数值仿真的你。学 PDE 数值解关键不在于一次性理解所有理论而在于把最简单的热传导方程真正调通再去挑战更复杂的耦合方程组。