MATLAB偏微分方程数值解实战：热传导方程差分格式与稳定性

偏微分方程数值解是很多 MATLAB 学习者从“会算数”走向“会做仿真”的一道分水岭。教材里往往先给推导，再给一个差分格式，但真正打开 MATLAB 后，你面对的是网格怎么取、时间步长怎么定、边界条件如何写进代码、为什么程序跑了一会儿数值就变成 NaN、为什么理论解明明在衰减你的曲线却在发散这一连串问题。这篇免费教程不绕弯，直接用一维热传导方程作为主线，讲清楚偏微分方程数值解在 MATLAB 里从建模、编码、验证到排错的完整链路。内容按概念、环境、求解、差分实现、问题排查和工程实践展开，适合正在学数值方法、计算物理、偏微分方程课程，或者刚接触 MATLAB 仿真但没有系统学过稳定性和边界处理的读者。

1. 先理解偏微分方程数值解在 MATLAB 里到底解决什么问题

1.1 为什么很多偏微分方程没有解析解

偏微分方程描述的是物理量在空间和时间上的变化规律。常见写法是：热传导方程描述温度随时间扩散，波动方程描述振动如何传播，Laplace 方程描述稳态势场分布。解析解的思路是找到一个函数表达式，让它同时满足方程、初始条件和边界条件。问题是，这种方法只对少数规则几何、简单系数、理想边界条件有效。实际工程里的求解域可能是飞机翼型截面，材料系数可能随位置变化，边界形状也不规则，这种情况下几乎不可能写出解析表达式。

数值解的思路完全不同：不找表达式，而是把连续的空间和时间离散成网格，把偏导数改成差分近似，再用计算机推进求解。MATLAB 能做这件事，不是因为它是某种神秘工具，而是因为它天然支持矩阵运算、数组索引、快速绘图和大量数值计算库函数。你不需要自己编写每一个线性代数底层细节，但必须理解离散化后引入了什么误差，时间步长为什么不能随意放大，边界条件由谁来强制施加。

1.2 三类典型偏微分方程决定了不同的求解策略

按照数学上的经典分类，偏微分方程大致分成三种，分类会直接决定你需要初值还是边值、时间上能不能推进、空间上如何迭代。

热传导方程是学习数值方法最合适的切入对象。它写成一维形式是u_t = alpha * u_xx，这里u(x,t)代表温度，alpha是热扩散系数。这个方程既有时间推进，又有空间二阶导数，能让你同时理解初值、边值、网格间距、时间步长和稳定性，而这些概念几乎适用于所有偏微分方程数值解问题。

1.3 在 MATLAB 里有三条主要实现路线

初学时常把工具关系搞混。MATLAB 中解决偏微分方程数值解不是只有一种模式，至少有三条路线可以选，它们定位不同：

• 直接调用pdepe函数。这是 MATLAB 基础环境提供的偏微分方程求解器，适合一维抛物型或椭圆型方程，函数会自动处理空间离散和时间推进，但需要你把方程整理成指定的标准形式。
• 使用 Partial Differential Equation Toolbox。这个工具箱主要处理二维、三维区域上的偏微分方程，底层基于有限元方法。你可以在图形界面或代码里建几何模型、设定边界条件、划分网格并求解。
• 自己手写差分格式或有限元代码。这条路看似老派，却是理解数值稳定性、边界处理、误差来源的最好方式。工程上要修改格式或检查异常时，往往还需要回到这一层。

很多人觉得 MATLAB 里“解偏微分方程”就是找某个神秘函数，其实正确做法是先判断问题是一维还是二维三维，是规则区域还是复杂几何，是需要快速验证还是需要研究算法细节，然后才选择路线。前面的选择错了，后面越改越混乱。

2. 环境与工具箱边界：不要把 pdepe 和 PDE Toolbox 混为一谈

2.1 最低环境准备与版本确认方式

学习偏微分方程数值解并不需要安装很多东西。MATLAB 基础环境本身就可以运行pdepe、差分格式脚本和绘图代码。如果你想做二维三维几何区域的有限元分析，需要额外安装 Partial Differential Equation Toolbox。这里版本不是越新越好，而是要看你的接口习惯和课程要求，但有一点必须检查：安装后确实具备对应的函数能力。

版本信息可以通过命令行直接确认：

% 查看 MATLAB 版本基本信息 version % 查看是否安装了 PDE Toolbox % 如果 createpde 返回路径，说明工具箱可用 which createpde % 查看 pdepe 的帮助 doc pdepe

运行后，如果显示类似C:\Program Files\MATLAB\...\toolbox\pde\...的路径，说明工具箱存在；如果没有路径，只显示“未找到”，说明你的环境中没有启用对应工具链。这个检查在写代码前做一次，能避免你复制了网上的二维求解代码却一直报“未定义函数”的错误。

2.2 pdepe 与 PDE Toolbox 的定位差异

pdepe是 MATLAB 基础自带的求解器，不是额外工具箱，这一点常被误解。它求解的问题限制在一维空间，也就是只有一个空间坐标，但可以同时解多个偏微分方程组成的方程组。它的好处是接口稳定、学习成本低，坏处是无法处理复杂的二维三维形状。

PDE Toolbox 是付费工具箱，核心价值在于支持二维和三维区域。它内部基于有限元方法，可以通过createpde创建模型对象，用geometryFromEdges或geometryFromMesh导入几何，用generateMesh划分网格，再用solvepde求解。适合做热结构耦合、电磁场分布、结构力学这类需要复杂几何的问题。

这里要提醒：不要在学习一维热传导时过早引入 PDE Toolbox。工具箱能生成很漂亮的色图，但它的多层封装会让你很难判断误差来自网格、边界还是物理参数。先跑通pdepe和手写差分，再进入工具箱，效率会高很多。

2.3 学习环境和生产环境的处理差异

如果只是完成课程作业，写脚本即可，不需要考虑部署。你只需要保证 MATLAB 路径切到当前文件夹，脚本和函数文件命名规范，运行时不出现未保存中间变量的问题。建议养成如下习惯：

clear; close all; clc;

这个三连清理了工作区、关闭所有图窗并清空命令行，避免上一次运行的残留变量干扰本次结果。

生产环境或者真正的仿真项目要复杂得多。工程计算不会允许你反复点脚本运行，而应该有可重复流程：

• 模型参数集中写在配置结构体里，而不是散落在多处。
• 每个仿真配置文件都要有版本记录，方便回溯。
• 结果保存成.mat文件，并把图导出版本可控的图片。
• 数值解结果要配套输出诊断信息，比如最大误差、是否满足稳定性条件。
• 大批量参数扫描时，考虑使用parfor并行，但前提是明确结果之间的独立性。

这些不是偏微分方程本身的内容，却是从“作业能跑”走向“结果能用”的必要环节。

3. 第一次实际求解：用 pdepe 计算热传导方程并与解析解对照

3.1 pdepe 要求的方程标准形式

pdepe不能直接接收任意写法，它要求方程整理成如下标准形式：

c(x,t,u,du/dx) * du/dt = x^(-m) * d/dx [ x^m * f(x,t,u,du/dx) ] + s(x,t,u,du/dx)

初看很抽象，先解释各项含义。c是时间导数项的系数，f是通量项，s是源项，m是几何对称系数。对于普通一维热传导方程：

du/dt = alpha * d2u/dx2

可以写成：

m = 0 c = 1 f = alpha * du/dx s = 0

也就是说，热传导方程中的二阶导数项通过通量项来表达。为什么要转换成这样的形式？因为 PDE 求解器可以统一处理后，无论是线性还是非线性方程，都能交给同一套空间离散框架处理。

3.2 pdepe 求解一维热传导的完整代码

考虑一个最简单的热传导问题。求解区间0 <= x <= 1，热扩散系数alpha = 0.02，时间从0到T = 2。

初始条件为：

u(x, 0) = sin(pi * x)

边界条件为两端恒零：

u(0, t) = 0 u(1, t) = 0

这个问题的解析解是：

u(x, t) = exp(-alpha * pi^2 * t) * sin(pi * x)

解析解存在的意义非常大，后面验证数值解是否准确时，不需要靠“看一眼曲线像不像”，而是能计算出一个明确的误差值。

pdepe 的边界条件格式是：

p(x,t,u) + q(x,t) * f(x,t,u,du/dx) = 0

对于要求u = 0的边界，令：

p = u q = 0

因此在左右两端分别写：

pl = ul; ql = 0; pr = ur; qr = 0;

这里ul和ur分别代表左端点和右端点的解值，f是通量。很多人习惯直接把u = 0写在边界条件里，但若写成u = 0而不是p = u, q = 0，会不符合 pdepe 的函数签名，导致无法运行或边界条件不生效。

完整代码可以这样组织：

%% 参数设置 alpha = 0.02; L = 1; T = 2; %% 空间网格和时间向量 xmesh = linspace(0, L, 101); tspan = linspace(0, T, 200); %% 求解 sol = pdepe(0, @(x,t,u,dudx) heatPDE(x,t,u,dudx,alpha), ... @(x) heatIC(x), ... @(xl,ul,xr,ur,t) heatBC(xl,ul,xr,ur,t), ... xmesh, tspan); %% 解析解，用于对照 uA = @(x,t) exp(-alpha * pi^2 * t) .* sin(pi * x); %% 绘制终态对比 figure; plot(xmesh, uA(xmesh, T), 'k-', 'LineWidth', 1.5); hold on; plot(xmesh, sol(end,:), 'ro', 'MarkerSize', 4); xlabel('x'); ylabel('u'); title('pdepe 数值解与解析解对比'); legend('解析解', 'pdepe 数值解', 'Location', 'best'); grid on;

对应的三个子函数：

function [c,f,s] = heatPDE(x,t,u,dudx,alpha) c = 1; f = alpha * dudx; s = 0; end function u0 = heatIC(x) u0 = sin(pi * x); end function [pl,ql,pr,qr] = heatBC(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur; qr = 0; end

代码运行后，图上会看到红色圆点几乎贴合黑色解析解曲线。如果不是这个效果，优先检查子函数的参数顺序，尤其是边界条件里pl要和左端点值ul对应，pr要和右端点值ur对应，这一点最容易写反。

3.3 如何验证结果而不只是“画出图”

画出一条光滑曲线不等于计算正确。需要进一步检查误差。下面这段代码计算每个时间层的最大误差和均方根误差：

U = sol; ErrMax = zeros(length(tspan), 1); ErrRMS = zeros(length(tspan), 1); for n = 1:length(tspan) tNow = tspan(n); uA_now = uA(xmesh, tNow); ErrMax(n) = max(abs(U(n,:) - uA_now)); ErrRMS(n) = sqrt(mean((U(n,:) - uA_now).^2)); end fprintf('最大误差: %.3e\n', max(ErrMax)); fprintf('均方根误差: %.3e\n', max(ErrRMS));

空间网格越密，时间层越多，误差通常会下降。这里也会暴露一个问题：误差太小不一定说明算法好，还可能是解析解本身在长时间后衰减到接近 0，绝对误差自然很小。如果要做严格验证，应该同时观察相对误差或者选择解本身不容易衰减到接近 0 的区间。

绘图时可以顺手使用 MATLAB 的colormap或坐标轴设置，让三维结果更清楚。例如把时间层作为第二维画曲面图：

figure; surf(xmesh, tspan, sol, 'EdgeColor', 'none'); xlabel('x'); ylabel('t'); zlabel('u'); colormap(jet); view(135, 30); colorbar;

曲面图适合观察温度随时间和空间的变化趋势。如果曲线出现锯齿状波动，尤其是波峰附近有明显抖动，说明空间网格过粗或者时间离散步长过大，需要通过后续显式格式实验进一步理解。

4. 亲手实现显式差分：真正理解稳定性门槛在哪里

4.1 从二阶导数到三点差分格式

pdepe能替你完成很多工作，但它不会告诉你为什么时间步长不能随便挑。为了搞懂这个关键问题，很值得手写一次经典显式格式。

把空间区间[0,1]分成Nx - 1段，每个节点编号为i = 1, 2, ..., Nx。二阶导数u_xx可以由相邻三点的值近似：

d2u/dx2 ≈ (u_{i-1} - 2*u_i + u_{i+1}) / dx^2

把热传导方程的时间导数用向前差商近似：

du/dt ≈ (u_i^{n+1} - u_i^n) / dt

于是得到显式递推格式：

u_i^{n+1} = u_i^n + r * (u_{i-1}^n - 2*u_i^n + u_{i+1}^n)

其中：

r = alpha * dt / dx^2

这个格式叫“显式”，因为第n+1层的值可以直接通过第n层的已知值求出，不需要解方程组。但这带来一个重要限制：r不能太大。对于一维扩散方程，稳定性条件通常是：

0 < r <= 0.5

如果超过 0.5，数值解往往会随时间振荡并发散，即使物理上温度本来应该在下降。

4.2 为什么会有稳定性条件而不是任意减小 dt 就行

直观理解：显式格式相当于把当前每一格的变化传给相邻格。如果时间步长太大，一个时间步内热量的传播距离超出了空间网格能承载的范围，数值上就会出现过冲。稳定条件的本质是保证误差不会被逐层放大。

所以这里有一个反直觉的结论：减小空间步长dx会让dx^2变小，从而让r变大；为了让r不超过 0.5，必须同步大幅减小时间步长dt。这导致显式格式在高分辨率网格下计算量增长很快。理解这一点后，你才能解释为什么工程中常会转向隐式格式或pdepe这种自动选择时间步的求解器。

4.3 完整可运行的显式差分脚本

下面脚本求解同一个热传导问题，并同时输出与解析解的误差。

%% 参数设置 alpha = 0.02; L = 1; T = 2; %% 空间网格 Nx = 51; dx = L / (Nx - 1); x = linspace(0, L, Nx)'; %% 时间步选择：取 r = 0.4，保证满足稳定性条件 r = 0.4; dt = r * dx^2 / alpha; Nt = round(T / dt); dt = T / Nt; % 重新规整 dt r = alpha * dt / dx^2; % 重新计算实际 r %% 初始条件 u0 = sin(pi * x); uAnalytic = @(x,t) exp(-alpha * pi^2 * t) .* sin(pi * x); u = u0; uNew = u; errMax = zeros(Nt, 1); %% 显式差分推进 for n = 1:Nt for i = 2:Nx-1 uNew(i) = u(i) + r * (u(i-1) - 2*u(i) + u(i+1)); end % 强制边界条件 uNew(1) = 0; uNew(Nx) = 0; u = uNew; errMax(n) = max(abs(u - uAnalytic(x, n*dt))); end %% 输出 fprintf('空间节点数: %d, dx = %.4f\n', Nx, dx); fprintf('时间步数: %d, dt = %.6f\n', Nt, dt); fprintf('稳定性系数 r = %.4f\n', r); fprintf('最大误差: %.3e\n', max(errMax)); %% 绘图对比 figure; plot(x, uAnalytic(x, T), 'k-', 'LineWidth', 1.5); hold on; plot(x, u, 'ro', 'MarkerSize', 5); xlabel('x'); ylabel('u'); title('显式差分结果与解析解'); legend('解析解', '差分数值解', 'Location', 'best'); grid on;

这段代码中，边界条件在每一时间步都被重新强制为 0。这里的顺序很重要：先更新内部节点，再覆盖边界节点。如果反过来，边界更新后又被内部节点的公式冲掉，边界条件就始终不会生效。

运行后，你会看到r = 0.4时数值解和解析解基本重合。如果只修改一个数字，把r = 0.4改成r = 0.6，然后重跑，最终结果会出现明显的振荡甚至 NaN，这就是稳定条件被破坏的直接证据。能亲手复现这个现象，比背十遍“CFL 条件”更有效。

4.4 差分脚本出现发散时的检查路径

当你修改参数后看见红色圆圈不再贴合理论曲线，而是到处乱跳，不要急着改回原参数。先做如下检查：

1. 打印当前的dx和dt，手动计算r = alpha * dt / dx^2。
2. 确认r <= 0.5。如果大于 0.5，优先减小dt，而不是增大dx。
3. 检查初始条件是否包含过大的高频分量，比如方波或阶跃，这类间断会让数值误差初期偏大。
4. 检查是否在每一层循环结束后强制更新了边界。
5. 如果结果在第一层就 NaN，检查alpha是否为正数，是否出现了除以 0 的网格。

显式格式的优势是代码透明，任何一步出错都能在变量编辑器里逐个节点检查。第一次遇到发散时，建议把Nx设成 11 或 21 这样的小网格，手工算前两三个时间步，这能很快定位是公式、边界还是参数的问题。

5. 数值解异常排查：从“能运行”到“结果可信”

5.1 排查顺序：先查输入，再查格式，别一上来就改算法

很多同学在看到错误曲线时习惯马上去搜索新算法，但实际项目中绝大多数异常来自更基础的问题。按照下面顺序排查通常最有效：

1. 物理量单位是否正确。热扩散系数可能是cm^2/s，而空间长度用的是米，两者相差 1e4 倍，这会直接让时间尺度错乱。
2. 边界条件写的是否是代码实际施加的条件。尤其是pdepe，边界都是以p + q*f = 0形式组合的，直接写u = 0是无效的。
3. 稳定性系数是否满足格式要求。
4. 空间网格分辨率和时间步长是否匹配。细网格配大步长，显式格式极易发散。
5. 看错误日志出现的位置。MATLAB 报错一般会告诉你哪一个函数、哪一行，不要只看最后的红色提示。
6. 检查绘图时的索引。sol是二维矩阵，第一维通常是时间，第二维是空间，索引反了会画出完全混淆的结果。

5.2 高频问题现象与对应处理方案

这张表没有覆盖所有情况，但如果你遇到的是“结果不像噪声却也不像物理”的中间状态，多半要从表格前几行找原因，而不是怀疑方程本身写错。

5.3 MATLAB 语法细节导致的隐蔽错误

偏微分方程代码中，数组运算和矩阵运算的差异经常造成难查的问题。比如你在向量循环中希望逐元素相乘，写成A * u时，MATLAB 会把它当作矩阵乘法，一旦矩阵维度刚好匹配，代码不会报错，但结果完全不是逐元素的物理含义。需要逐元素运算时用.*、./、.^。

下面是一个容易出错的写法：

% 错误示例：矩阵乘，可能报错或得到错误结果 u = alpha * dt / dx^2 * (u(1:end-2) - 2*u(2:end-1) + u(3:end));

问题在于，直接用整段子数组做右端项时，u(1:end-2)的维度已经和你待更新的内部节点个数对齐，这条代码如果写在循环里会导致每次把整体向量赋值的逻辑和索引搞混。建议初学时坚持使用节点循环方式：

for i = 2:Nx-1 uNew(i) = u(i) + r * (u(i-1) - 2*u(i) + u(i+1)); end

这种方式在网格不大时速度可以接受，而且可读性好。当你确认逻辑正确后，再改成下面这种向量化写法以提升速度：

u(2:end-1) = u(2:end-1) + r * (u(1:end-2) - 2*u(2:end-1) + u(3:end));

但注意向量化写法返回的是新数组右端项，而左端仍然引用原数组，实际执行时可能出现先后依赖。推荐先写好节点循环版本，再逐行验证向量化版本的结果，两者一致后才用于正式计算。

6. 工程实践清单：从演示代码变成可信仿真结果

6.1 最小验证清单：发布结果前逐项确认

如果这段代码最终要用于实验报告、课程论文或项目数据，应该执行一遍下面的检查清单：

• [ ] 方程形式与物理问题一致：是一维还是二维，是稳态还是瞬态。
• [ ]alpha和单位正确，空间长度与时间单位匹配。
• [ ] 初始条件在边界处与边界条件相容。若初始条件在端点不为 0，就会立刻形成间断，影响早期解。
• [ ] 显式格式的r不超过 0.5，或明确说明使用其他稳定条件。
• [ ] 数值解与解析解或已知物理规律对照过，误差在可接受范围。
• [ ] 网格加密后误差方向正确。通常二阶差分在网格加密一倍时误差会显著下降。
• [ ] 边界条件在每一步都被强制更新。
• [ ] 关键结果保存为.mat和图片，脚本可以重新运行得到同一结果。

其中“初始条件在边界处是否相容”特别容易被忽略。比如初始温度是sin(pi*x)，在x=0和x=1处天然为 0，这才和两端恒零边界条件相容。如果你初始化为u0 = 1 + 0.1*x，在两端都不等于 0，数值推进早期会产生明显的边界调整层。解决方法是选择更贴近物理的初始分布，或者显式处理边界值上的突变。

6.2 参数选择与结果记录：不要把所有参数散落在脚本各处

偏微分方程数值解项目里，真正容易犯的错误不是代码逻辑，而是参数混乱。一个可维护的项目可以这样组织参数：

%% config.m 或者脚本头部 alpha = 0.02; L = 1; T = 2; Nx = 101; r = 0.4; % 稳定性设计值 scheme = 'explicit'; % explicit / pepe saveResult = true;

把参数集中在文件头部，以后做参数扫描时只需要修改一处，就不需要在不同函数里翻找你曾经硬编码的数值。时间步长尽量用格式推出来，而不是随便给一个值：

dx = L / (Nx - 1); dt = r * dx^2 / alpha; % 由稳定性条件推导

这样如果Nx变了，dt会自动跟着变，不会出现“改了网格但忘记改时间步长”导致发散的经典问题。

如果要做多组参数扫描，把结果写入结构体并按参数命名：

result.config = struct('alpha', alpha, 'Nx', Nx, 'r', r); result.errMax = max(errMax); result.uFinal = u; save(['result_alpha_', num2str(alpha), '_Nx_', num2str(Nx), '.mat'], 'result');

这么做的意义在于：等你第二天回来处理下一个仿真，或者要把结果交给别人复现时，不再需要手动回忆当时的参数。

6.3 从一维课程作业扩展到二维几何模型

当你把一维热传导方程和显式差分完全弄懂后，再走向二维问题就不会手足无措。二维热传导方程多了y方向的二阶导数，离散时每个时间步需要同时计算u_{i-1,j} - 2u_{i,j} + u_{i+1,j}与u_{i,j-1} - 2u_{i,j} + u_{i,j+1}，稳定性条件会变为类似alpha * dt * (1/dx^2 + 1/dy^2) <= 0.5的形式。你会发现，显式格式在高维问题中时间步长越来越受限制，这自然引向隐式格式、交替方向隐式格式 ADI 和预处理迭代法的学习。

如果问题几何比较复杂，就不再适合手写差分，建议转用 PDE Toolbox。它的核心流程是：

model = createpde('thermal', 'transient'); % 导入几何，例如从 STL 文件或几何函数 % 指定材料属性 % 设置初始条件与边界条件 % generateMesh 生成有限元网格 % solvepde 求解 % pdeplot 可视化

这里面每一步都有对应的验证方法，尤其要注意网格划分质量对结果的影响。二维三维仿真中，网格畸形单元往往比算法误差更早击垮结果。做这类问题前，最好先补上有限元方法的节点编号、单元刚度矩阵和边界加载原理，否则遇到工具箱内部的报错信息时会很难定位。

6.4 值得继续深入的方向

一维热传导的显式格式只是入口。继续学习时，可以按以下顺序扩展：

• 隐式格式与 Crank-Nicolson 格式。隐式格式没有严格的时间步稳定性限制，适合长时间推进。
• 二维热传导与 ADI 格式。
• 波动方程的差分格式，你会遇到完全不同的稳定性条件和边界人为反射问题。
• 对流占优问题中的迎风格式和数值耗散。
• 有限元方法如何处理不规则几何。
• 谱方法在高光滑问题中的高精度优势。

每次扩展都建议保留“解析解对照”这一步。解析解不是每个问题都有，但只要存在，它就是最好的验证工具。工程上没有解析解时，则要设计守恒量、单调性或者网格收敛性测试来侧面确认计算结果可信。

如果只挑一个练习来巩固本文内容，建议亲手改一改显式格式中的稳定性系数r，在图上对比r = 0.4和r = 0.6两条完全不同的命运。能解释清楚为什么一个步长之差会让热传导方程从“物理衰减”变成“数值爆炸”，偏微分方程数值解的很多经验就算真正摸到门了。