企业官网建设流程全解析

大家在科研和工程计算里遇到偏微分方程（Partial Differential Equation，PDE）时，常常会陷入两难：数学推导太复杂、手写求解器太耗时，而网上找到的 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 三大类主要解法

在实际学习中，接触最多的方法有三个方向：

下面从最容易理解的有限差分法开始。

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^{n+1} - u_i^n}{\Delta t} \approx \frac{\partial u}{\partial t} )
• 二阶中心差分：( \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2} \approx \frac{\partial^2 u}{\partial x^2} )

显式格式可以写成：

[ u_i^{n+1} = u_i^n + \alpha \frac{\Delta t}{\Delta x^2} \left( u_{i+1}^n - 2u_i^n + u_{i-1}^n \right) ]

这个格式简单，但必须满足稳定性条件：

[ r = \frac{\alpha \Delta t}{\Delta x^2} \le 0.5 ]

如果不满足，计算结果会发散，出现数值振荡。

3.2 完整显式有限差分 MATLAB 脚本

下面给出一个一维热传导方程的完整求解示例。物理场景是一根长度为 (L=1) 的细杆，两端温度恒为 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(x>=0.4 & x<=0.6) = 1; % 保存用于绘图 u_all = zeros(nx, nt+1); u_all(:,1) = u; % 显式时间迭代 for n = 1:nt % 内点更新 u_new = u; for i = 2:nx-1 u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1)); end % 边界条件直接置 0 u_new(1) = 0; u_new(nx) = 0; u = u_new; u_all(:, n+1) = u; end % 绘制结果 figure; [X, Tgrid] = meshgrid(0:nt, x); % 这里只需绘制三个时刻 plot_times = [1, round(nt/5), round(nt/2), nt+1]; 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;

对这个代码做几点解释：

1. 稳定性参数r在这里非常重要。当r大于 0.5 时，程序会直接报错，防止算出无意义结果。
2. 时间循环中，我们遍历每个内点，使用上一时刻的u(i-1)、u(i)、u(i+1)计算下一时刻，这种格式称为显式格式。
3. 边界条件通过强制赋值实现。Dirichlet 边界条件下，直接修改端点值是简单的处理方式。

运行以上脚本后，可以看到温度从中间矩形区域向两侧扩散，最终整体温度趋于 0，与无热源且两端恒温的物理直觉一致。

3.3 隐式格式与 Crank-Nicolson 格式简介

显式格式虽然简单，但时间步长受限制。实际计算中更常用隐式格式，因为它在参数上更稳定。对同一方程，后向欧拉格式可以写为：

[ u_i^{n+1} - r \left( u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1} \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, i+1) = -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)：问题的几何对称类型。(m=0) 代表平板/直角坐标，(m=1) 代表柱对称，(m=2) 代表球对称。
• (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 ]

因此得到：

• (c=1)
• (f=\alpha \cdot \partial u / \partial x)
• (s=0)

完整代码如下：

%% 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 % 对两端 u=0，相当于 p=ul-0、q=0 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)); % 如果要得到所有时间层在 x=0.5 处的温度曲线 u_mid = sol(:, 51); figure; plot(tspan, u_mid); xlabel('t'); ylabel('u(t, x=0.5)'); title('x=0.5 处温度随时间变化');

如果换了一套网格，x=0.5 不一定正好落在节点上。更稳妥的做法是先找最近索引：

[~, idx] = min(abs(xmesh - 0.5)); u_mid = sol(:, idx);

在处理数值解时，这种“先找索引再取值”的写法比硬编码节点编号更安全。

5. 使用 PDE Toolbox 求解二维稳态问题

5.1 PDE Toolbox 是什么

当问题推广到二维或三维，尤其边界形状不是矩形时，有限差分法的手写成本快速上升。MATLAB 的 PDE Toolbox（偏微分方程工具箱）采用有限元方法，可以处理更复杂的几何区域。

PDE Toolbox 的基本流程可以概括为：

1. 创建几何模型（createpde）。
2. 建立几何形状（geometryFromEdges或从 STL 导入三维几何）。
3. 指定方程系数和边界条件。
4. 生成网格（generateMesh）。
5. 求解并可视化。

这里用一个经典二维拉普拉斯方程案例来说明。问题为：在一个边长为 1 的正方形区域内，求解泊松方程：

[ -\nabla^2 u = 1 ]

边界条件设为 (u=0)。这是一个有源场问题，相当于在均匀边界接地、内部均匀激励的物理模型。

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 边界条件书写中，要求满足：

[ p + q \cdot f = 0 ]

其中 (f) 就是方程中的通量函数。以热传导方程为例，如果左端绝热，意味着 (\partial u / \partial x = 0)，那么取 (p=0)、(q=1)，这样：

[ 0 + 1 \cdot f = f = \alpha \frac{\partial u}{\partial x} = 0 ]

正好符合绝热条件。如果右端是恒温 (u=1)，则取 (p=u-1)、(q=0)。

初学者最容易犯的错误是：把 Neumann 边界理解成直接使用导数赋值，但在 pdepe 格式中要写成通量条件，必须知道 (f) 的表达式。不同类型方程之间，这种写法可能略有差异。

6.2 初始条件与边界条件的协调

设定初始条件时还有个很容易忽略的问题：初始条件最好和边界条件协调一致，否则在 (t=0) 附近可能出现非物理的突变。

比如在前面的案例中，初始时刻 (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 验证程序是否写对的基本方法

对偏微分方程数值解有一个通用的“验算三件套”：

1. 检查守恒量或单调趋势是否合理。
2. 检查稳态解是否和解析解一致。
3. 检查减小步长后结果是否变化很小。

以热传导方程为例，当 (t) 足够大时，解会趋于边界条件决定的稳态。如果边界条件为零，则最终全场接近 0。如果做的是绝热边界，则总热量应保持不变，可以计算：

total_heat = trapz(x, u);

不同时刻的总热量变化不应太大。如果热量明显流失，说明边界条件可能写错。

8. 常见报错与排查思路

偏微分方程数值解的 MATLAB 程序中，报错信息和你预想的可能不太一样。下面列几个频率较高的现象。

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和ul，pr是否只依赖xr和ur。

8.2 显式差分中索引越界

很多新手的循环会写成：

for i = 2:nx u_new(i) = u(i) + r * (u(i+1) - 2*u(i) + u(i-1)); end

当i = nx时，会访问u(nx+1)，越界报错。这类问题并不难修，关键是写循环前先画一条网格索引轴线，明确边界点是 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 数值解，关键不在于一次性理解所有理论，而在于把最简单的热传导方程真正调通，再去挑战更复杂的耦合方程组。