企业官网建设流程全解析

很多学习 MATLAB 的朋友都有类似经历：用ode45解常微分方程时还能跟上节奏，一到“偏微分方程数值解”这一章，突然就看不进去了。原因不复杂——常微分方程的未知量是y(t)，本质是沿着时间轴推进；偏微分方程的未知量是u(x,t)，空间和时间两个维度同时出现，光是把一个连续问题变成程序能算的离散问题，就涉及网格、差分格式、边界条件、收敛性和稳定性。

如果你关注过“大谦MATLAB”这类免费 MATLAB 教程，会发现偏微分方程数值解几乎是所有进阶路线的共同分水岭。学到这里，算法本身倒不是最大的障碍，真正的障碍是：你脑子里没有一个清晰的“从方程到代码”的转换流程。市面上的教程要么上来就堆差分格式推导，要么只告诉你怎么调用工具箱，很少有人先把三件事讲清楚：这条方程在数学上属于哪一类？边界条件和初始条件应该怎么给？MATLAB 里到底有哪几条路可以走，分别适合什么场景？

这篇文章想帮你解决的就是这个问题。我会用一维热传导方程作为主线，分别演示 MATLAB 中求解偏微分方程数值解的三种典型路线：自带函数pdepe、PDE Toolbox 的有限元流程、以及自己写有限差分代码。读完以后，你应该能判断手里的方程该走哪条路，并且能自己排查常见的运行错误。

1. 偏微分方程数值解：为什么值得单独花时间学？

先说判断：学习偏微分方程数值解，重点不在把公式背熟，而在建立“连续方程 → 离散系统 → 数值验证”的实验意识。这门技术在工程里的价值，要远大于考试里“会解一个方程”的价值。

实际项目中的热传导、流体扩散、电磁场分布、结构应力，基本都由偏微分方程描述。能用解析方法求解的问题只占极少数：几何必须规则、系数必须常值、边界条件必须简单。真实工程问题通常是不规则几何加非线性系数，这时候唯一可行的办法就是数值解。

MATLAB 在数值计算领域的优势在于，它把很多底层实现封装成了可调用函数。新手不需要从零写一个有限元求解器，也能先跑通二维区域上的泊松方程；同时，它又允许你手动指定网格密度、边界条件和误差容限。换句话说，MATLAB 能把你的精力引导到数学建模和结果分析上，而不是陷入底层矩阵组装和线性方程求解细节。

另外值得强调的是：偏微分方程数值解的很多概念，比如网格收敛性、时间步长稳定性、边界条件处理，换到其他语言和工具依然适用。你在 MATLAB 里理解清楚一次，后面转到 C++、Python、COMSOL、FEniCS，思路都能复用。

2. MATLAB 里求解 PDE 的三条主流路线

很多新手最大的困惑在于，不知道自己该用哪个工具。这里先给你一张总览，后面每一节再展开细节。

2.1pdepe：一维问题的最快入口

pdepe的最大价值是让你不用关心有限元、有限差分具体怎么实现。你只需要把一个一维 PDE 描述成它规定的标准形式，再提供初始条件和边界条件，它就能在时间方向上帮你完成积分。

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)

第一次看这个形式的人会有点头晕，但实际使用中，我们通常只需要定义四个函数：pdefun描述方程里的c、f、s，icfun描述初始条件，bcfun描述边界条件。边界条件一律写成：

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

这个约定是新手最容易出错的地方。很多人直接写“左边界等于 0”，却不理解p和q的取值组合。后面我会用代码演示。

2.2 PDE Toolbox：面向任意几何的有限元流程

如果你的问题在二维区域上，或者几何形状不规则，pdepe就不够用了。PDE Toolbox 提供了一条偏微分方程数值解的经典有限元流程：

创建 PDE 模型 → 绘制或导入几何 → 设置方程系数 → 设置边界条件 → 生成网格 → 求解 → 后处理

这条流程和商用有限元软件高度相似。学会它以后，哪怕你以后换用别的软件，操作思路也是通用的。

2.3 自写差分法：理解问题本质的必修课

自写差分法在生产环境中看起来“效率不高”，因为你要自己处理时间推进、空间网格和稳定性。但从学习角度讲，只有自己写过一次差分代码，你才会真正理解：为什么计算会发散？为什么要限制时间步长？为什么边界条件必须单独处理？

3. 正式动手前必须理清的三个关键概念

不管是哪条路线，开始做偏微分方程数值解之前，有几个关键概念必须搞清楚。否则代码写得再漂亮，结果也解释不了。

3.1 方程类型决定了问题“胖瘦”

偏微分方程按数学性质可以分很多类，但你做数值解时只要关心三种大类：

• 椭圆型方程：例如泊松方程-Δu = f，只包含空间二阶导数，问题只在空间上“定解”。
• 抛物型方程：例如热传导方程ut = uxx，包含时间一阶导，需要初值和边界条件。
• 双曲型方程：例如波动方程utt = c²uxx，包含时间二阶导，既需要初值，也可能涉及波前传播。

不同类型方程的最优数值算法差别很大。椭圆型最终落到一个大规模线性方程组；抛物型通常沿时间推进，每一步可能都要解一个线性系统；双曲型则要额外关注数值耗散和振荡。

3.2 定解条件不是“辅助信息”

没有初始条件和边界条件的偏微分方程，在数值上通常会得到毫无意义的结果，甚至直接报错。

边界条件又分成大家熟悉的几类：

拿到一个模型，第一件事不是写代码，而是确认自己遇到的是哪一类定解问题。缺少边界条件的模型，在数学上可能是不适定的。

3.3 网格和时间步长不是“随便选的”

在数值解里，空间离散的粗细用网格步长h表示，时间推进步长用dt表示。减小h会让空间分辨率更高，但也会导致线性方程组规模变大、计算耗时增长。更麻烦的是，显式时间推进格式通常有一个稳定性限制，比如显式有限差分求解热传导，要求：

dt / h^2 ≤ 0.5

如果不满足这个条件，数值误差会被指数放大，前面还正常的波形，几步之后就会出现高频振荡，最终变成 NaN。因此，网格和时间步长是偏微分方程数值解里第一优先级的问题。

4. 环境准备：跑 PDE 数值解需要 MATLAB 里有什么？

先说结论：做偏微分方程数值解，不建议只用免费在线版或未安装工具箱的精简版。最好安装完整版 MATLAB，并确保需要的工具箱处于可用状态。

你可以在 MATLAB 命令窗口输入：

ver

在输出的列表里找到是否有Partial Differential Equation Toolbox。如果没有，但你的问题只是一维经典热传导方程，那么可以直接使用pdepe，它在 MATLAB 基础环境中就可以使用；如果要做二维三维有限元，就需要安装 PDE Toolbox。

安装或激活工具箱时，注意下面几点：

• 登录 MathWorks 账号，在 License Center 中确认授权范围内是否包含 PDE Toolbox。
• 如果安装后工具箱仍不可用，优先查看ver列表，而不是怀疑代码本身。
• 涉及许可证的服务问题，请通过学校或公司正版渠道联系管理员，不要使用来路不明的激活方式。

使用免费 MATLAB 教程学习时，建议把版本差异放在心上。不同 MATLAB 版本的 PDE Toolbox API 会有变化，比如老版本常用pdetool图形界面，新版本则推荐命令行方式createpde、solvepde。下面的代码我以目前主流的命令行 API 为例。

5. 路线一：用 pdepe 求解一维热传导方程

我们从一个最简单的物理模型开始：一根长度为 1 的细杆，初始温度分布为sin(pi*x)，两端温度始终为 0。方程是：

ut = uxx

它描述的是热量从中间向两端传导、最终温度全部变成 0 的过程。这个问题的解析解是u_exact = exp(-pi^2 * t) * sin(pi*x)，适合用来验证数值结果。

完整的 MATLAB 函数文件如下。

% 文件名：heat_pdepe_demo.m function heat_pdepe_demo() m = 0; xmesh = linspace(0, 1, 50); tspan = linspace(0, 0.2, 100); sol = pdepe(m, @pdefun, @icfun, @bcfun, xmesh, tspan); % pdepe 可能返回多分量方程，这里只取第 1 个分量 u = sol(:,:,1); surf(xmesh, tspan, u); xlabel('x'); ylabel('t'); zlabel('u(x,t)'); title('pdepe 求解一维热传导方程'); end % ------------------------------------------------------------------ % 方程系数定义： % c * ut = x^(-m) * d/dx [ x^m * f ] + s % 对 ut = uxx：c=1，f=dudx，s=0 function [c, f, s] = pdefun(x, t, u, dudx) c = 1; f = dudx; s = 0; end % ------------------------------------------------------------------ % 初始条件：t=0 时，u0 = sin(pi*x) function u0 = icfun(x) u0 = sin(pi * x); end % ------------------------------------------------------------------ % 边界条件： % p + q * f = 0 % 左边界和右边界都要求 u = 0 function [pl, ql, pr, qr] = bcfun(xl, ul, xr, ur, t) pl = ul; % 左边界满足 ul = 0 ql = 0; pr = ur; % 右边界满足 ur = 0 qr = 0; end

运行方式很简单：把文件保存为heat_pdepe_demo.m，在命令窗口输入heat_pdepe_demo。如果一切正常，会弹出一个三维曲面图，初始温度中间高、两端低，随时间向下衰减，最终趋于 0。

这里一定要重点理解bcfun中的写法和教材里的区别。在教材上，两端固定温度 0 会写成u(0,t)=0、u(1,t)=0；在pdepe里，每个边界条件都要求写成p + q*f = 0的形式。要让u本身等于 0，就应该令p = u、q = 0，这样边界条件退化成u = 0。如果写成p=0, q=1，含义就变成了f = 0，也就是dudx = 0，这是绝热边界，不是你要的固定温度边界。

为了检查数值解是否正确，可以继续在函数末尾追加误差分析代码：

% 与解析解对比 u_exact = exp(-pi^2 * tspan)' * sin(pi * xmesh); err = max(max(abs(u - u_exact))); fprintf('pdepe 最大绝对误差: %.3e\n', err);

运行后，误差一般会很小。这个结果也能说明：当你正确设置边界条件和初始条件后，pdepe对这类经典热传导问题的求解精度是相当高的。

6. 路线二：用 PDE Toolbox 处理二维问题

接下来我们把问题的维度从一维提升到二维。考虑这样一个经典问题：在一个单位圆区域内求解泊松方程：

-Δu = 1

边界条件为 Dirichlet 条件u = 0。物理上可以理解为一个圆形薄板上均匀加热，边缘温度保持 0，最终形成的稳态温度分布。

在 PDE Toolbox 中，解决这个问题的标准流程如下。

% 文件名：poisson_disk_demo.m % 在单位圆内求解 -Δu = 1, 边界 u = 0 model = createpde(1); % 使用 MATLAB 自带的单位圆几何函数 geometryFromEdges(model, @circleg); % 查看几何边界编号（调试时取消注释） % pdegplot(model, 'EdgeLabels', 'on'); % 所有边界都使用 Dirichlet 条件 applyBoundaryCondition(model, 'dirichlet', ... 'Edge', 1:model.Geometry.NumEdges, ... 'u', 0); % 定义方程系数： % 对于 -Δu = 1，写成 m*utt + d*ut - div(c*grad(u)) + a*u = f % 所以 m=0, d=0, c=1, a=0, f=1 specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', 1); % 生成网格，Hmax 控制最大网格尺寸 generateMesh(model, 'Hmax', 0.05); % 求解 result = solvepde(model); % 绘制结果 figure; pdeplot(model, 'XYData', result.NodalSolution, ... 'ZData', result.NodalSolution, ... 'ColorBar', 'on'); title('单位圆内 Poisson 方程数值解'); xlabel('x'); ylabel('y');

运行这段代码前，确保 MATLAB 中已经安装了 PDE Toolbox。如果命令窗口报错说createpde未定义，或者找不到geometryFromEdges，基本可以确定是工具箱缺失。

这段代码中值得注意的点有三个。

第一，几何模型里的边界编号不一定从 1 连续排到某个数，但在@circleg这样简单的内置几何中，用1:model.Geometry.NumEdges就能覆盖所有边界。更复杂的几何，建议先用pdegplot(model, 'EdgeLabels', 'on')查看边界编号。

第二，specifyCoefficients的方程系数是带符号约定的。PDE Toolbox 内部统一使用：

m * utt + d * ut - div(c * grad(u)) + a * u = f

所以，要让方程等价于-Δu = 1，需要把c设为 1，f设为 1，而不是把u的符号手动写成负号。

第三，Hmax越小，网格越密，计算结果越接近真实解，但耗时会增加。做正式计算前，可以先从Hmax=0.1开始试跑，确认模型没问题后再加密。

对于这个圆域问题，你可以用解析解u_exact = (1 - r^2)/4来验证中心点的数值结果，其中r是到圆心的距离。中心点解析值应为0.25，网格加密后，有限元解的中心值会越来越接近这个数。

7. 路线三：用有限差分法自己写一遍

工具函数虽然方便，但如果你想真正理解偏微分方程数值解的核心，我建议至少自己写一次显式有限差分。这里仍然以热传导方程为例，但故意不使用pdepe，而是用最原始的向量化写法。

思路很简单：

• 把[0,1]空间区域分成nx-1段；
• 把时间[0,T]分成nt-1步；
• 空间二阶导数用中心差分：uxx ≈ (u_{i+1} - 2u_i + u_{i-1}) / dx^2
• 时间一阶导数用前向欧拉。

代码实现如下：

% 文件名：fd_heat_1d_demo.m function fd_heat_1d_demo() L = 1; T = 0.2; nx = 50; dx = L / (nx - 1); % 显式格式稳定性条件：dt / dx^2 <= 0.5 dt = 0.4 * dx^2; nt = ceil(T / dt); dt = T / nt; x = linspace(0, L, nx)'; % 初始条件 u = sin(pi * x); u(1) = 0; u(end) = 0; figure; plotStep = max(1, round(nt / 5)); for k = 1:nt u_old = u; % 只更新内部节点，边界节点始终为 0 u(2:end-1) = u_old(2:end-1) + ... dt / dx^2 * (u_old(3:end) - 2 * u_old(2:end-1) + u_old(1:end-2)); % 每隔一段时间画一条曲线 if mod(k, plotStep) == 0 plot(x, u, 'LineWidth', 1.2); hold on; drawnow; end end xlabel('x'); ylabel('u'); grid on; title('显式有限差分：u_t = u_{xx}'); legend(compose('t = %.3f', (plotStep:plotStep:nt)' * dt)); end

保存后直接运行fd_heat_1d_demo，你应该能看到一条初始正弦曲线逐渐被压平，最终接近直线u=0。

最后，可以加上一步验证：

% 最终时刻解析解 u_true = exp(-pi^2 * T) * sin(pi * x); fprintf('有限差分最终时刻最大误差: %.3e\n', max(abs(u - u_true)));

显式有限差分虽然简单，但它把偏微分方程数值解的两个核心问题暴露得很清楚。首先是稳定性：如果dt / dx^2超过 0.5，程序很容易发散。其次是边界处理：如果你的差分更新公式不小心把边界节点也更新了，那么固定温度边界就会被破坏，求解出来的温度会不断偏离真实值。

我建议你做一个实验：把dt = 0.4 * dx^2改成dt = 0.6 * dx^2，再运行一遍。你会看到数值解在几个时间步之后出现明显的高频振荡，最终结果变成 NaN。亲眼看到这个发散过程，比背十遍稳定性公式都管用。

8. 数值解可信吗？网格收敛性与可视化验证

代码不报错，不代表数值解正确。偏微分方程数值解里最重要的检查，就是看加密网格后结果是否稳定收敛。

8.1 用网格收敛性检查解

网格越细，理论上数值解应该越接近真实解。如果同一问题在nx=20、40、80、160下计算，解的误差不断下降，说明你的算法实现和边界条件大概率是正确的；如果误差不降反升，说明代码里很可能有隐藏 bug。

在heat_pdepe_demo.m的函数体内，你可以追加这样的循环：

% 网格收敛性测试 nxList = [20 40 80 160]; errList = zeros(size(nxList)); for i = 1:numel(nxList) nx = nxList(i); xmesh = linspace(0, 1, nx); sol = pdepe(m, @pdefun, @icfun, @bcfun, xmesh, tspan); u = sol(:,:,1); u_exact = exp(-pi^2 * tspan)' * sin(pi * xmesh); errList(i) = max(max(abs(u - u_exact))); end table(nxList', errList', 'VariableNames', {'nx', 'maxError'})

运行后你会看到一个两列的表：网格数量翻倍时，最大误差通常会明显下降。如果你看到的是误差反而增大，那就要优先检查边界条件函数和初始条件函数是否写错。

8.2 PDE 数值解的可视化建议

偏微分方程数值解的结果通常是一个二维或三维数组，直接用命令行打印数据并不直观。以下三个命令可以提高你的结果查看效率：

% 二维场图用 surf 或 imagesc surf(xmesh, tspan, u); xlabel('x'); ylabel('t'); zlabel('u'); % 或者用平面俯视图 figure; imagesc(tspan, xmesh, u'); axis xy; colorbar; xlabel('t'); ylabel('x');

对于 PDE Toolbox 的有限元结果，刚才代码里已经用了pdeplot。需要调整颜色风格时，可以在绘图后加：

colormap(parula); colorbar;

如果你需要把结果导出到论文或者报告里，高版本的 MATLAB 推荐使用：

exportgraphics(gcf, 'heat_solution.png', 'Resolution', 300);

这个命令比老式的print更容易控制图片尺寸和清晰度，是偏微分方程数值解后处理中很实用的小技巧。如果你的结果是在圆域或极坐标几何中得到的，需要画半径方向随角度变化时，可以先把某条半径或边界上的场值提取成theta和u(theta)两个向量，再用：

polarplot(theta, u_theta);

不要把整个二维场直接传给polarplot，那样只会得到形状混乱的图。

9. 常见问题、工程建议与学习路线图

9.1 常见问题排查表

9.2 数组运算中的 .* 与 * 到底怎么选

很多人在写有限差分或解析解验证时，会在 MATLAB 的.*和*上栽跟头。这两者的区别很简单：

• *是矩阵乘法；
• .*是对应元素相乘。

比如下面这个例子：

a = [1 2 3]; b = [4 5 6]; c = a .* b; % 结果 [4 10 18]，对应元素相乘 d = a * b'; % 结果 32，矩阵乘法

在偏微分方程数值解里，我们通常操作的是网格向量和矩阵，绝大多数情况下需要的是元素级运算。比如构造解析解矩阵时：

u_exact = exp(-pi^2 * tspan)' * sin(pi * xmesh);

这里我故意用了矩阵乘法*，目的是把时间向量和空间向量做外积，生成一个二维网格矩阵。如果你把它误写成.*，虽然 MATLAB 新版在维度兼容时会做隐式扩展，但你的本意是生成二维网格，不是对应元素相乘，代码的可读性会大大降低。

9.3 给初学者的工程建议

第一条建议：从有解析解的问题入手。偏微分方程数值解最怕的是你算出一个数却不知道对不对。先用ut = uxx、`-