MATLAB 偏微分方程 (PDE) 求解¶

偏微分方程 (Partial Differential Equation, PDE) 是描述未知多元函数及其偏导数之间关系的方程。它们是模拟热传导、流体力学、波传播、电磁场等几乎所有物理现象的基础。

MATLAB 求解 PDE 的核心是 Partial Differential Equation Toolbox™。该工具箱主要使用有限元法 (Finite Element Method, FEM) 来求解 PDE。

数学原理：有限元法 (FEM)¶

有限元法的核心思想是“分而治之”： 1. 离散化 (Meshing): 将复杂的求解域（几何形状）分割成许多简单、互不重叠的子域，称为“有限元”（如三角形或四面体）。所有有限元的顶点被称为“节点”。 2. 局部近似: 在每个有限元内部，假设未知函数的解可以由一个简单的基函数（通常是低阶多项式）来近似。 3. 组装方程组: 将 PDE 转化为一个等效的积分形式（弱形式），然后通过将每个单元上的近似解“组装”起来，最终形成一个大型的、稀疏的代数方程组 \(KU=F\)。 4. 求解: 求解这个代数方程组，得到在每个节点上的解的近似值。

MATLAB 的 PDE 工具箱将这个复杂的过程封装在了一套直观的工作流程中。

PDE 的一般形式¶

工具箱求解的 PDE 通常可以写成以下一般形式：

\[ m \frac{\partial^2 u}{\partial t^2} + d \frac{\partial u}{\partial t} - \nabla \cdot (c \nabla u) + au = f \]

其中 \(u\) 是待求解的函数，系数 \(m, d, c, a, f\) 分别代表质量、阻尼、扩散、吸收/反应和源项，它们可以是空间 \((x, y, z)\)、时间 \(t\)、解 \(u\) 甚至解的梯度 \(\nabla u\) 的函数。

交互式工作流：PDE Modeler App¶

PDE Modeler App - 图形化PDE求解工具¶

功能对于二维 PDE 问题，PDE Modeler App 提供了一个完整的图形化用户界面 (GUI)，允许用户通过交互式操作完成从几何建模、网格生成到求解和后处理的全过程。
启动方式
1. 在 MATLAB APPS 选项卡中，找到并点击 PDE Modeler 图标。
2. 在命令行窗口中输入 pdeModeler。
工作流程
1. 绘制几何: 使用内置的绘图工具（矩形、圆形、多边形）通过构造实体几何 (CSG) 的方式创建求解域。
2. 设置边界条件: 在边界模式下，为每条边界指定边界条件。
3. 指定 PDE 系数: 在 PDE 模式下，输入方程的系数。
4. 生成网格: 初始化并细化有限元网格。
5. 求解: 运行求解器。
6. 后处理: 以多种方式（如彩色图、等高线、箭头图）可视化结果。

示例

% 打开 PDE Modeler App
pdeModeler

% 在弹出的 App 界面中，可以手动操作：
% 1. 绘制一个单位圆。
% 2. 在 "Boundary Mode" 下，为所有边界设置 Dirichlet 条件 u=0。
% 3. 在 "PDE Specification" 中，设置 c=1, a=0, f=10 (求解泊松方程)。
% 4. 点击 "Solve" 按钮即可看到结果。

编程工作流 (Programming Workflow)¶

对于需要自动化、参数化或求解三维问题的场景，编程工作流是必经之路。

`createpde` - 创建PDE模型容器¶

功能创建一个 PDE 模型对象，作为容纳几何、网格、系数、边界条件等所有问题信息的容器。
语法 matlab model = createpde() model = createpde(N) 语法说明*
- createpde(): 创建一个标量 PDE 模型（求解一个因变量）。
- createpde(N): 创建一个包含 N 个耦合方程的 PDE 系统。

示例

% 创建一个标量 PDE 模型
heatModel = createpde();

几何定义与网格生成¶

`decsg` - 分解构造实体几何 (CSG) 矩阵¶

功能 decsg (Decomposed Solid Geometry) 是一个核心函数，用于处理通过“构造实体几何”(Constructive Solid Geometry, CSG) 方法定义的二维几何。它的主要作用是将多个基本形状（如矩形、圆形、多边形）通过集合运算（并、交、差）组合成的复杂几何，分解成一组互不重叠的最小区域（minimal regions）。这个分解后的几何描述是后续生成网格（generateMesh）和进行有限元分析的基础。

语法

g = decsg(gd)
g = decsg(gd, sf, ns)
[g, bt] = decsg(...)

语法说明
- 基本形状定义 (gd 矩阵) gd (Geometry Description) 是几何描述矩阵，矩阵的每一列代表一个基本形状。
  - 圆形 (Circle): 第1行为 1。第2、3行分别是圆心的x, y坐标。第4行是半径。
```
% 圆心(0,0), 半径0.5
C1 = [1; 0; 0; 0.5];
```
  - 多边形 (Polygon): 第1行为 2。第2行是顶点的数量。后续行成对出现，表示每个顶点的x, y坐标。
```
% 有3个顶点的三角形
% P1 = [2; 3; x1; x2; x3; y1; y2; y3];
```
  - 矩形 (Rectangle): 第1行为 3。更简单的定义方式是 [3, 4, xmin, xmax, ymin, ymax] 的转置形式。
```
% 左下角(-1,-1), 右上角(1,1)的矩形
R1 = [3; 4; -1; 1; 1; -1; -1; -1; 1; 1];
```
  - 椭圆 (Ellipse): 第1行为 4。第2、3行是中心坐标。第4、5行是半轴长度。第6行是旋转角度（弧度）。
```
% 中心(0,0), 半轴(2,1), 旋转pi/4
% E1 = [4; 0; 0; 2; 1; pi/4];
```
    注意: gd 矩阵的所有列必须有相同的行数，如果某个形状的参数较少，则需用零填充使其对齐。
- sf (Set Formula): 一个字符串，定义了如何通过集合运算（+ 代表并集, * 代表交集, - 代表差集）来组合 gd 中定义的基本形状。
- ns (Name Space): 一个字符矩阵，将 gd 中的形状（列）与 sf 中使用的名称对应起来。
- g: 分解后的几何矩阵，可以被 geometryFromEdges 函数使用。
- bt (Boolean Table): 一个布尔表，显示了每个最小区域隶属于哪个原始基本形状。

示例：创建一个带孔的矩形

% 1. 定义基本形状：一个矩形(R1)和一个圆形(C1)
% 矩形 R1: 角点(-1, -1)，长2，宽2
R1 = [3, 4, -1, 1, 1, -1, -1, -1, 1, 1]'; 
% 圆形 C1: 圆心(0, 0)，半径0.5
C1 = [1, 0, 0, 0.5]';

% 2. 将形状名称与几何描述关联
% 注意：gd矩阵的列数要对齐，这里R1有10行，C1只有4行，需要将C1用0填充至10行
gd = [R1, [C1; zeros(6,1)]];
ns = (['R1'; 'C1'])';

% 3. 设置集合公式：矩形减去圆形
sf = 'R1 - C1';

% 4. 分解几何
g = decsg(gd, sf, ns);

% 5. 在 PDE 模型中使用并可视化
model = createpde();
geometryFromEdges(model, g);

figure;
pdegplot(model, 'EdgeLabels', 'on', 'FaceLabels', 'on');
title('带孔矩形的分解几何');

`geometryFromEdges` - 从边创建二维几何¶

功能通过定义几何体的边界（边）来创建二维几何描述。
语法
```
g = geometryFromEdges(model, edgeData)
```

示例

% 创建一个矩形几何
model = createpde();
R1 = [3, 4, 0, 1, 1, 0, 0, 0, 1, 1]'; % 矩形的几何描述矩阵
g = decsg(R1); % 分解 CSG 矩阵
geometryFromEdges(model, g);

% figure; pdegplot(model, 'EdgeLabels', 'on'); % 可视化几何

`generateMesh` - 生成有限元网格¶

功能在已添加到模型中的几何上创建并存储一个二维或三维的有限元网格。

语法

generateMesh(model)
generateMesh(model, 'Name', Value)

语法说明
- 'Hmax': 控制网格单元的最大尺寸，是控制网格密度的最常用参数。

示例

% ... 延续上一个示例 ...
generateMesh(model, 'Hmax', 0.1); % 生成最大单元尺寸为 0.1 的网格
% figure; pdeplot(model); % 可视化网格

PDE 问题定义¶

`specifyCoefficients` - 定义PDE方程系数¶

功能为模型指定 PDE 的系数（\(m, d, c, a, f\))。

语法

specifyCoefficients(model, 'Name', Value, ...)

示例: 求解泊松方程 \(- \nabla \cdot (\nabla u) = 1\)

% 对应于一般形式，c=1, a=0, f=1
specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', 1);

`applyBoundaryCondition` - 施加边界条件¶

功能为几何的指定边界（边或面）施加边界条件。
数学原理
- 狄利克雷 (Dirichlet): 直接指定解在边界上的值。 \(hu = r\)
- 诺伊曼 (Neumann): 指定解在边界法向上的通量或梯度。 \(\bm{n} \cdot (c \nabla u) + qu = g\)

语法

applyBoundaryCondition(model, type, location, value, ...)

示例

% ... 延续上一个示例 ...
% 在所有边界上施加 Dirichlet 条件 u=0
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:4, 'u', 0);

`setInitialConditions` - 设置初始条件¶

功能为瞬态（时域）PDE 问题指定初始条件。

语法

setInitialConditions(model, u0) % t=0 时的解
setInitialConditions(model, u0, ut0) % t=0 时的解和解的导数

PDE 求解器¶

`solve` - 求解静态PDE问题¶

功能求解静态（稳态）PDE 问题，如静电场、稳态热传导等（椭圆型 PDE）。
语法
```
results = solve(model)
```

示例 (完整的泊松方程求解)

% 1. 创建模型
model = createpde();
% 2. 创建几何 (L 形膜)
gd = [2 6 0 1 1 -1 -1 0 0 0 -1 -1 1 1]';
g = decsg(gd);
geometryFromEdges(model, g);
% 3. 生成网格
generateMesh(model, 'Hmax', 0.05);
figure; % 创建新图形窗口
pdeplot(model); % 绘制网格
title('生成的网格');
% 4. 定义系数
specifyCoefficients(model, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', 10);
% 5. 定义边界条件 (所有边界 u=0)
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:model.Geometry.NumEdges, 'u', 0);
% 6. 求解
results = solvepde(model);
% 7. 可视化
figure;
pdeplot(model, 'XYData', results.NodalSolution, 'Contour', 'on');
title('solve - 求解泊松方程');

`solvepde` - 求解时域PDE问题¶

功能求解瞬态（时域）PDE 问题，如热传导过程、波传播等（抛物线型和双曲线型 PDE）。
语法
```
results = solvepde(model, tlist)
```
语法说明
- tlist: 一个时间向量，指定了需要输出解的时间点。

示例 (热传导)

此示例程序求解的热传导方程为：

\[ \frac{\partial u}{\partial t} - \nabla \cdot (c \nabla u) = 0 \]

在下面的代码中，通过 specifyCoefficients 设置 d=1 (对应 \(\partial u/\partial t\) 项) 和 c=0.5 (扩散项) 来实现。

% 1. 创建模型和几何 (同上)
model = createpde();
geometryFromEdges(model, g);
generateMesh(model, 'Hmax', 0.05);
% 2. 定义热传导系数 (抛物线型)
% du/dt - div(c*grad(u)) = 0  =>  d=1, c=0.5, a=0, f=0
specifyCoefficients(model, 'm', 0, 'd', 1, 'c', 0.5, 'a', 0, 'f', 0);
% 3. 边界条件 (u=0)
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:model.Geometry.NumEdges, 'u', 0);
% 4. 初始条件 (初始温度为 1)
setInitialConditions(model, 1);
% 5. 定义求解时间
tlist = 0:0.005:0.2;
% 6. 求解
results = solvepde(model, tlist);
% 7. 可视化并创建动图
figure;
filename = 'temperature_animation.gif'; % 定义GIF文件名
total_duration = 5; % 动图总时长 (秒)
delay_time = total_duration / length(tlist); % 计算每帧的延迟
for i = 1:length(tlist)
    % 绘制当前时间的温度分布
    pdeplot(model, 'XYData', results.NodalSolution(:, i), 'Contour', 'on', 'ColorMap', 'hot');
    title(['t = ', num2str(tlist(i)), ' s']);
    xlabel 'x';
    ylabel 'y';
    clim([0, 1]); % 固定颜色范围以便观察变化

    % 捕获帧
    frame = getframe(gcf);
    im = frame2im(frame);
    [imind, cm] = rgb2ind(im, 256);

    % 写入 GIF
    if i == 1
        imwrite(imind, cm, filename, 'gif', 'Loopcount', inf, 'DelayTime', delay_time);
    else
        imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delay_time);
    end
end
disp(['动图已保存为 ', filename]);

`solvepdeeig` - 求解特征值PDE问题¶

功能求解 PDE 特征值问题，如结构振动的模态、量子力学中的能级等。
语法
```
results = solvepdeeig(model, range)
```
语法说明
- range: 一个二元向量 [rmin, rmax]，指定了搜索特征值的范围。

示例 (一维量子谐振子)

求解定态薛定谔方程，找到一维谐振子的特征能量 (Energy Eigenvalues) 和特征波函数 (Eigenfunctions)。

定态薛定谔方程为：

\[ -\frac{\hbar^2}{2m} \frac{\mr{d}^2\psi}{{\d{x}}^2} + V(x)\psi = E\psi \]

对于谐振子，势能 \(V(x) = \frac{1}{2}m\omega^2x^2\)。为了简化，我们取 \(\hbar=1, m=1, \omega=1\)，方程变为：

\[ -\frac{1}{2} \frac{\mr{d}^2\psi}{{\d{x}}^2} + \frac{1}{2}x^2\psi = E\psi \]

这对应于 PDE 特征值问题的一般形式 \(-\nabla \cdot (c \nabla u) + au = \lambda d u\)，其中：

\(c = 1/2\)
\(a = x^2/2\)
\(d = 1\)
特征值 \(\lambda\) 对应能量 \(E\)

% 1. 创建模型
model = createpde(1);

% 2. 定义一维几何 (区间 [-L, L])
L = 5;
model.Geometry = pde.IntervalGeometry(-L, L);

% 3. 生成网格
generateMesh(model, 'Hmax', 0.05);

% 4. 定义 PDE 系数
% -d/dx(c*du/dx) + a*u = lambda*d*u
% c=1/2, a=x^2/2, d=1
specifyCoefficients(model, 'm', 0, 'd', 1, 'c', 0.5, 'a', @(location, state) 0.5*location.x.^2, 'f', 0);

% 5. 边界条件 (波函数在无穷远处为 0)
applyBoundaryCondition(model, 'dirichlet', 'Edge', 1:2, 'u', 0);

% 6. 求解特征值
% 搜索范围，理论能量为 E_n = n + 0.5
range = [0, 10]; 
results = solvepdeeig(model, range);

% 7. 可视化
figure;
% 提取特征值和特征向量
eigenvalues = results.Eigenvalues;
eigenvectors = results.Eigenvectors;

% 绘制前 5 个波函数
plot(results.Mesh.Nodes, eigenvectors(:, 1:5));
grid on;
title('一维谐振子的前5个波函数');
xlabel('x');
ylabel('\psi(x)');
legend(strcat('E = ', num2str(eigenvalues(1:5), '%.4f')));

% 打印能量值
disp('计算得到的能量特征值:');
disp(eigenvalues);
disp('理论能量特征值 (n+0.5):');
disp((0:length(eigenvalues)-1)' + 0.5);

`pdepe` - 一维PDE专用求解器¶

功能一个专门用于求解一类一维抛物线和椭圆 PDE 的函数。其语法与现代 PDE 工具箱工作流完全不同。

语法

sol = pdepe(m, pdefun, icfun, bcfun, xmesh, tspan)

语法说明
- m: 对称性常数（0-笛卡尔，1-圆柱，2-球）。
- pdefun, icfun, bcfun: 分别定义 PDE 系数、初始条件和边界条件的函数句柄。

示例 (一维热传导)

m = 0;
x = 0:0.1:1;
t = 0:0.1:2;

sol = pdepe(m, @pdex1pde, @pdex1ic, @pdex1bc, x, t);
u = sol(:,:,1);

figure;
surf(x, t, u);
title('pdepe - 一维热传导方程');
xlabel('x'); ylabel('t'); zlabel('u(x,t)');

function [c,f,s] = pdex1pde(x,t,u,dudx) % 定义 PDE
c = 1; f = dudx; s = 0;
end
function u0 = pdex1ic(x) % 定义初始条件
u0 = sin(pi*x);
end
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % 定义边界条件
pl = ul; ql = 0; pr = ur; qr = 0;
end

MATLAB 偏微分方程 (PDE) 求解¶

数学原理：有限元法 (FEM)¶

PDE 的一般形式¶

交互式工作流：PDE Modeler App¶

PDE Modeler App - 图形化PDE求解工具¶

编程工作流 (Programming Workflow)¶

createpde - 创建PDE模型容器¶

几何定义与网格生成¶

decsg - 分解构造实体几何 (CSG) 矩阵¶

geometryFromEdges - 从边创建二维几何¶

generateMesh - 生成有限元网格¶

PDE 问题定义¶

specifyCoefficients - 定义PDE方程系数¶

applyBoundaryCondition - 施加边界条件¶

setInitialConditions - 设置初始条件¶