很多初学者学习 MATLAB 时，往往把大量时间花在矩阵运算、绘图和 Simulink 仿真上，真正遇到偏微分方程数值解的时候，会突然发现不知道从哪里入手：方程形式五花八门，教材里全是差分格式和收敛性证明，直接照抄代码又经常报错。这篇文章想帮你打通“物理问题 → 数学模型 → MATLAB 数值实现”的完整链路。
文中不追求严格的数学证明，而是把最常用的几类偏微分方程数值求解方法讲清楚：先用pdepe求解一维抛物型方程，再用手写有限差分法处理二维稳态问题，然后介绍 PDE Toolbox 的有限元流程，最后用一个线方法案例说明时间依赖问题如何借助 ODE 求解器高效推进。如果你正准备做 MATLAB 偏微分方程相关的作业、课程设计或工程验证，这篇文章可以作为一份可以直接对照操作的入门资料。
1. 偏微分方程数值解与 MATLAB：先说清楚背景
1.1 什么是偏微分方程数值解
偏微分方程，英文全称是 Partial Differential Equation，简称 PDE。它描述的是包含多个自变量的函数及其偏导数之间的关系。最常见的自变量是空间坐标x, y, z和时间t。
例如一维热传导方程：
∂u/∂t = α · ∂²u/∂x²其中u(x,t)表示温度场，α是热扩散系数。这个方程同时包含时间一阶导数和空间二阶导数，因此属于典型的抛物型偏微分方程。
与常微分方程 ODE 不同，PDE 的解是一个场，而不是一条简单的曲线。绝大多数实际 PDE 都无法求出解析解，即使能求出，也通常只在极简单的几何区域和边界条件下成立。工程中更常用的做法是对区域进行离散，把连续方程转化为有限维代数方程或常微分方程组，再借助计算机求解。这就是偏微分方程数值解的基本思想。
按方程性质，PDE 大致可以分成三类：
| 类型 | 典型方程 | 物理背景 | 常用处理思路 |
|---|---|---|---|
| 椭圆型 | 拉普拉斯方程Δu = 0 | 稳态热传导、静电场 | 有限差分法、有限元法 |
| 抛物型 | 热传导方程u_t = αΔu | 扩散、瞬态导热 | pdepe、线方法、有限元 |
| 双曲型 | 波动方程u_tt = c²Δu | 声波、弹性波 | 特征线法、CFL 条件控制 |
不同方程对时间推进格式、空间离散阶数和稳定性条件的要求不同，所以数值求解第一件事是判断方程类型。
1.2 MATLAB 在 PDE 数值解中的优势
很多编程语言都能做 PDE 数值模拟，但 MATLAB 的优势在于交互效率高、矩阵运算方便、画图工具成熟。对于教学和科研验证阶段，我们通常先在小规模网格上验证算法，再逐步扩大计算规模。MATLAB 可以让你用较短的代码完成“建模 → 求解 → 可视化 → 对比解析解”的闭环。
MATLAB 自带的pdepe函数专门用于求解一维抛物型和椭圆型 PDE，不需要额外安装工具箱。二维和三维有限元问题可以借助 PDE Toolbox 完成。如果不使用工具箱，也可以自己写有限差分法，MATLAB 的稀疏矩阵和向量化操作能明显降低代码难度。
1.3 不是所有偏微分方程都适合硬推公式
数值方法的本质是用离散近似替代微分运算。我们不需要掌握所有差分格式的详细推导，但至少要理解：网格越细，逼近精度通常越高；时间步长和空间步长之间的比例会影响稳定性；边界条件处理错误会导致结果完全不可信。
下面从最简单的环境准备开始，然后逐步进入代码实战。
2. 环境准备与 MATLAB 基础检查
2.1 软件版本与工具箱
本文的示例主要使用 MATLAB 基础模块和 PDE Toolbox。pdepe、ode15s、spdiags、surf等函数在基础模块中已经包含。
如果你使用的是近几年的 MATLAB 版本，以下代码基本可以直接运行。如果版本较老，例如 R2014a 之前，部分函数在参数解析上会略有不同。建议先运行：
ver查看当前环境的版本信息和已安装工具箱。
PDE Toolbox 并非 MATLAB 基础安装自带，需要通过 Add-Ons 或学校授权获取。如果运行createpde时提示：
Undefined function or variable 'createpde'说明当前环境没有安装 PDE Toolbox。此时可以跳过本文第 5 章，先使用pdepe和手写有限差分法完成练习。
如果还没有可用的 MATLAB 环境，建议优先使用学校或单位提供的正版授权，也可以申请 MathWorks 官方试用版。不要使用来路不明的破解包，一方面存在安全风险，另一方面也会影响后续系统稳定性。
2.2 建议的项目文件结构
为了让代码更清晰，建议在本地新建一个独立目录，例如matlab_pde_demo。后续所有脚本都放在这个目录中：
matlab_pde_demo/ main_heat_pdepe.m laplace2d_jacobi.m pde_toolbox_poisson.m line_method_diffusion.mMATLAB 对中文路径的支持虽然还可以，但某些函数在读取文件或导出数据时可能因中文路径出现问题。工程推荐使用英文目录名，代码内部的注释可以使用中文。
2.3 入门必须注意的 MATLAB 语法习惯
PDE 数值解代码通常会涉及大量数组运算。先检查下面几个基础点。
第一，矩阵乘法*和数组乘法.*完全不同。差分矩阵作用到解向量时使用*，边界条件和系数场逐点计算时使用.*。
第二，linspace生成的默认是行向量。如果函数需要的是列向量，例如后续ode15s的初值向量，建议显式转置。
x = linspace(0,1,101); % 行向量 xcol = x(:); % 列向量第三，匿名函数是 PDE 求解中非常方便的工具。
d = 0.01; rhs = @(t,u) d * u;匿名函数可以捕获当前工作区中的变量，例如这里的d，从而避免频繁使用全局变量。
3. 用 pdepe 求解一维热传导方程
3.1 pdepe 的标准形式
pdepe求解的是如下形式的一维偏微分方程：
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)这里：
m表示坐标系类型，0 代表直角坐标，1 代表柱坐标，2 代表球坐标；c是时间项系数；f是通量项；s是源项。
对应的边界条件需要写成：
p(x,t,u) + q(x,t) · f(x,t,u,∂u/∂x) = 0初值形式为：
u(x,0) = u0(x)pdepe的调用方式为：
sol = pdepe(m, pdefun, icfun, bcfun, xmesh, tspan);其中pdefun定义偏微分方程和源项，icfun定义初始条件，bcfun定义边界条件。
3.2 热传导方程的标准转换
现在考虑一个具体问题。长度为 1 的细杆，两端温度始终为 0，初始温度分布为：
u(x,0) = sin(πx)热传导方程为：
∂u/∂t = ∂²u/∂x²该方程有解析解：
u(x,t) = exp(-π²t) · sin(πx)因此我们可以用解析解来验证数值方法的准确性。
对应pdepe标准形式，令：
m = 0; c = 1; f = ∂u/∂x; s = 0;边界条件在x=0和x=1处都取温度为零，即：
p = u; q = 0;3.3 完整代码与函数拆解
新建脚本文件main_heat_pdepe.m：
function main_heat_pdepe() % 一维热传导方程：u_t = u_xx % 区间：0 < x < 1，0 < t < 1 % 边界：u(0,t)=0，u(1,t)=0 % 初始：u(x,0)=sin(pi*x) clear; close all; clc; x = linspace(0, 1, 101); t = linspace(0, 1, 21); sol = pdepe(0, @heatPDE, @heatIC, @heatBC, x, t); % sol 的大小为 length(t) x length(x) x 方程个数 u = sol(:,:,1); % 绘制三维曲面图 figure; surf(x, t, u); xlabel('x'); ylabel('t'); zlabel('u'); title('pdepe 求解一维热传导方程'); % 对比数值解和解析解 uExact = exp(-pi^2 * t(end)) * sin(pi * x); figure; plot(x, u(1,:), 'k--', 'LineWidth', 1.5); hold on; plot(x, u(end,:), 'r-', 'LineWidth', 1.5); plot(x, uExact, 'bo', 'LineWidth', 1); legend('t=0 数值解', 't=1 数值解', 't=1 解析解', 'Location', 'North'); xlabel('x'); ylabel('u'); title('数值解与解析解对比'); end function [c, f, s] = heatPDE(x, t, u, DuDx) % 方程主体函数 c = 1; f = 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在 MATLAB 中运行该脚本后，会得到两个图像。第一张图是u(x,t)随时间变化的曲面，第二张图对比了t=1时刻的数值解和解析解。
3.4 结果解读
从曲面上可以看出，初始温度分布sin(πx)随时间推移逐渐衰减，最终趋于 0。由于两端温度固定为 0，热量不断从边界导出。
这里要注意两个关键点：
pdepe自动选择时间步长和空间离散方案，但对非常尖锐的初值或复杂源项，你需要增加网格点数量，并仔细检查结果是否出现非物理振荡；sol(:,:,1)表示第一个因变量。如果方程组包含多个未知函数，第三个维度会对应多个变量。
pdepe适合一维问题，但如果是二维或三维区域，就需要换用其他方法。
4. 手写有限差分法求解二维稳态问题
4.1 有限差分法的核心思想
有限差分法的基本思路是用差商近似导数。对于二维拉普拉斯方程：
∂²u/∂x² + ∂²u/∂y² = 0假设空间步长都为h，那么在内部节点(i,j)处可以写成：
(u(i+1,j) - 2u(i,j) + u(i-1,j)) / h² + (u(i,j+1) - 2u(i,j) + u(i,j-1)) / h² ≈ 0整理后得到经典的五点差分格式：
u(i,j) ≈ 0.25 * (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1))这个公式说明，在拉普拉斯方程中，任意内部节点的值等于它上下左右四个邻居的平均值。这就是 Jacobi 迭代求解 Laplace 方程的基础。
4.2 具体问题与边界条件
考虑单位正方形区域：
0 ≤ x ≤ 1，0 ≤ y ≤ 1边界条件设置为：
u(0,y) = 0 u(1,y) = 1 u(x,0) = x u(x,1) = x该问题的解析解恰好是：
u(x,y) = x因此非常适合检验差分格式是否正确。
4.3 MATLAB 完整代码
新建laplace2d_jacobi.m：
clear; close all; clc; % 内部网格剖分数 n = 50; % 包含边界点，因此坐标点数为 n+2 x = linspace(0, 1, n + 2); y = x; % 解矩阵：行对应 y，列对应 x u = zeros(n + 2, n + 2); % 边界条件 u(1, :) = x; % y = 0 u(end, :) = x; % y = 1 u(:, 1) = 0; % x = 0 u(:, end) = 1; % x = 1 % Jacobi 迭代 maxIter = 30000; tol = 1e-6; for iter = 1:maxIter uNew = u; uNew(2:end-1, 2:end-1) = 0.25 * ( ... u(1:end-2, 2:end-1) + ... u(3:end, 2:end-1) + ... u(2:end-1, 1:end-2) + ... u(2:end-1, 3:end)); err = max(abs(uNew(:) - u(:))); u = uNew; if err < tol fprintf('Jacobi 迭代在第 %d 次收敛\n', iter); break; end end % 解析解 uAnalytic = repmat(x, n + 2, 1); % 查看最大误差 maxErr = max(abs(u(:) - uAnalytic(:))); fprintf('最大误差：%.3e\n', maxErr); % 可视化 figure; surf(x, y, u); xlabel('x'); ylabel('y'); zlabel('u'); title('二维拉普拉斯方程有限差分解'); shading interp;运行后可以看到，迭代解和线性分布u = x几乎完全一致，最大误差通常在1e-6量级。
4.4 为什么使用 Jacobi 迭代
这里没有直接求解大规模线性方程组，而是使用 Jacobi 迭代。这么做是为了让代码更容易理解：每次更新时，内部节点只依赖上一层迭代的周围节点值。
Jacobi 迭代的缺点是收敛速度较慢。网格越细，迭代次数越多。实际项目中更推荐组装稀疏矩阵后用：
uVec = A \ b;求解。但上一段代码展示的思路是有限差分法的起点，适合教学和验证二维 Laplace 方程。
4.5 扩展到 Poisson 方程
如果把拉普拉斯方程变为 Poisson 方程：
-Δu = f(x,y)五点差分格式会变为：
4u(i,j) - u(i-1,j) - u(i+1,j) - u(i,j-1) - u(i,j+1) = h² · f(i
