求解偏微分方程

在偏微分方程 (PDE) 中，要求解的函数取决于几个变量，微分方程可以包括关于每个变量的偏导数。偏微分方程可用于对波浪、热流、流体扩散和其他空间行为随时间变化的现象建模。

使用 MATLAB 可求解哪些类型的 PDE？

MATLAB^® PDE 求解器 pdepe 使用一个空间变量 x 和时间 t 对 PDE 方程组的初始边界值问题求解。您可以将这些看作一个变量的 ODE，它们也会随着时间而变化。

pdepe 要求解的一维方程大概可分为以下两类：

带时间导数的方程是抛物型方程。例如热方程 $\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}}$ 。
不带时间导数的方程是椭圆型方程。例如拉普拉斯方程 $\frac{\partial^{2} u}{\partial x^{2}} = 0$ 。

pdepe 要求方程组中存在至少一个抛物型方程。换句话说，方程组中至少一个方程必须包含时间导数。

pdepe 还可求解某些二维和三维问题，这些问题由于角对称而简化为一维问题（有关详细信息，请参阅对称常量 m 的参数描述）。

Partial Differential Equation Toolbox™ 将此功能扩展到狄利克雷和诺伊曼边界条件下的二维和三维广义问题。

求解一维 PDE

一维 PDE 包含函数 u(x,t)，该函数依赖于时间 t 和一个空间变量 x。MATLAB PDE 求解器 pdepe 求解以下形式的一维抛物型和椭圆型 PDE 的方程组

$c (x, t, u, \frac{\partial u}{\partial x}) \frac{\partial u}{\partial t} = x^{- m} \frac{\partial}{\partial x} (x^{m} f (x, t, u, \frac{\partial u}{\partial x})) + s (x, t, u, \frac{\partial u}{\partial x}) .$

方程具有以下属性：

PDE 在 t₀ ≤ t ≤ t_f 和 a ≤ x ≤ b 时成立。
空间区间 [a, b] 必须为有限值。
m 可以是 0、1 或 2，分别对应平板、柱状或球面对称性。如果 m > 0，则 a ≥ 0 也必须成立。
系数 $f (x, t, u, \frac{\partial u}{\partial x})$ 是通量项， $s (x, t, u, \frac{\partial u}{\partial x})$ 是源项。
通量项必须取决于偏导数 ∂u/∂x。

关于时间的偏导数耦合只限于与对角矩阵 $c (x, t, u, \frac{\partial u}{\partial x})$ 相乘。此矩阵的对角线元素为零或正数。为零的元素对应于椭圆型方程，任何其他元素对应于抛物型方程。必须至少存在一个抛物型方程。如果 x 的某些孤立值是网格点（即计算解的位置），那么在这些值处，抛物型方程对应的 c 元素可能消失。当物质界面上有网格点时，允许 c 和 s 中出现界面导致的不连续点。

求解过程

要使用 pdepe 求解 PDE，您必须定义 c、f 和 s 的方程系数、初始条件、解在边界处的行为以及在其上计算解的点网格。函数调用 sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) 使用以下信息计算指定网格上的一个解：

m 是对称常量。
pdefun 定义要求解的方程。
icfun 定义初始条件。
bcfun 定义边界条件。
xmesh 是 x 的空间值向量。
tspan 是 t 的时间值向量。

xmesh 和 tspan 向量共同构成一个二维网格，pdepe 在该网格上计算解。

方程

您必须按照 pdepe 所需的标准形式表示 PDE。以这种形式编写，您可以读取系数 c、f、s 的值。

在 MATLAB 中，您可以用以下形式的函数编写方程代码

function [c,f,s] = pdefun(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end

在本例中，pdefun 定义方程

\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}}

。如果有多个方程，则 c、f 和 s 均为向量，其中每个元素对应于一个方程。

初始条件

在初始时间 t = t₀ 时，针对所有 x，解分量均满足以下格式的初始条件

$u (x, t_{0}) = u_{0} (x) .$

在 MATLAB 中，您可以用以下形式的函数对初始条件进行编码

function u0 = icfun(x)
u0 = 1;
end

在本例中，u0 = 1 定义 u₀(x,t₀) = 1 的初始条件。如果有多个方程，则 u0 是一个向量，其中每个元素定义一个方程的初始条件。

边界条件

在边界 x = a 或 x = b 时，针对所有 t，解分量满足以下形式的边界条件

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0.$

q(x,t) 是对角线矩阵，其元素全部是零或全部是非零。请注意，边界条件以通量 f（而非关于 x 的 u 的偏导数）形式表示。同时，在 p(x,t,u) 和 q(x,t) 这两个系数之间，只有 p 可以依赖于 u。

在 MATLAB 中，您可以用以下形式的函数对边界条件进行编码

function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
pL = uL;
qL = 0;
pR = uR - 1;
qR = 0;
end

pL 和 qL 是左边界的系数，pR 和 qR 是右边界的系数。在本例中，bcfun 定义边界条件

$\begin{array}{l} u_{L} (x_{L}, t) = 0 \\ u_{R} (x_{R}, t) = 1 \end{array}$

如果有多个方程，则输出 pL、qL、pR 和 qR 是向量，其中每个元素定义一个方程的边界条件。

积分选项

可以选择 MATLAB PDE 求解器中的默认积分属性来处理常见问题。在某些情况下，可以通过覆盖这些默认值来提高求解器的性能。为此，请使用 odeset 创建一个 options 结构体。然后，将该结构体作为最后一个输入参数传递给 pdepe：

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

在基础 ODE 求解器 ode15s 的选项中，只有下表中所示的选项可用于 pdepe。

类别	选项名称
误差控制	`RelTol`, `AbsTol`, `NormControl`
步长	`InitialStep`, `MaxStep`
事件日志记录	`Events`

解的计算

在您用 pdepe 求解方程后，MATLAB 将以三维数组 sol 返回解，其中 sol(i,j,k) 包含在 t(i) 和 x(j) 处计算的解的第 k 个分量。通常，您可以使用命令 u = sol(:,:,k) 提取第 k 个解分量。

您指定的时间网格仅用于输出目的，不影响求解器采用的内部时间步。但是，您指定的空间网格会影响解的质量和速度。求解方程后，您可以使用 pdeval 计算 pdepe 采用不同空间网格返回的解结构体。

示例：热方程

打开实时脚本

抛物型 PDE 的一个示例是一维热方程：

$\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}} .$

此方程描述 $0 \leq x \leq L$ 和 $t \geq 0$ 的散热情况。目标是求解 $u (x, t)$ 温度问题。温度最初是一个非零常量，因此初始条件是

$u (x, 0) = T_{0} .$

此外，左边界的温度为零，右边界的温度不为零，因此边界条件为

$u (0, t) = 0,$

$u (L, t) = 1 .$

要在 MATLAB® 中求解该方程，您需要对方程、初始条件和边界条件编写代码，然后在调用求解器 pdepe 之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾（如本示例所示），或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。

编写方程代码

在编写方程代码之前，您需要确保它的形式符合 pdepe 求解器的要求：

在此形式中，热方程是

$1 \cdot \frac{\partial u}{\partial t} = x^{0} \frac{\partial}{\partial x} (x^{0} \frac{\partial u}{\partial x}) + 0 .$

因此，系数的值如下：

$m = 0$
$c = 1$
$f = \frac{\partial u}{\partial x}$
$s = 0$

$m$ 的值作为参数传递给 pdepe，而其他系数编写为方程的一个函数，即

function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end

（注意：所有函数都作为局部函数包含在示例的末尾。）

代码初始条件

热方程的初始条件函数对 $u_{0}$ 赋给一个常量值。此函数必须接受 $x$ 的输入，即使它未使用。

function u0 = heatic(x)
u0 = 0.5;
end

编写边界条件代码

pdepe 求解器所需的边界条件的标准形式是

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

以这种形式编写的此问题的边界条件是

$u (0, t) + (0 \cdot f) = 0,$

$(u (L, t) - 1) + (0 \cdot f) = 0 .$

因此， $p$ 和 $q$ 的值是

$p_{L} = u_{L}, q_{L} = 0 .$
$p_{R} = u_{R} - 1, q_{R} = 0$ .

则对应的函数是

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end

选择解网格

使用包含 20 个点的空间网格和包含 30 个点的时间网格。由于解快速达到稳态， $t = 0$ 附近的时间点间隔更近以将此行为捕获到输出中。

L = 1;
x = linspace(0,L,20);
t = [linspace(0,0.05,20), linspace(0.5,5,10)];

求解方程

最后，使用对称性值 $m$ 、PDE 方程、初始条件、边界条件以及 $x$ 和 $t$ 的网格来求解方程。

m = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);

对解进行绘图

使用 pcolor 可视化解矩阵。

colormap hot
pcolor(x,t,sol)
colorbar
xlabel('Distance x','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')

Figure contains an axes object. The axes object with title Heat Equation for 0 less equals x less equals 1 and 0 less equals t less equals 5, xlabel Distance x, ylabel Time t contains an object of type surface.

局部函数

function [c,f,s] = heatpde(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
function u0 = heatic(x)
u0 = 0.5;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur - 1;
qr = 0;
end

PDE 示例和文件

我们提供几个示例文件，它们可作为求解最常见的一维 PDE 问题的绝佳起点。要查看和运行示例，请使用 Differential Equations Examples App。要运行此 App，请键入

odeexamples

要打开单独的文件进行编辑，请键入

edit exampleFileName.m

要运行示例，请键入

exampleFileName

下表包含可用 PDE 示例文件的列表。

示例文件	描述	示例链接
`pdex1`	简单的 PDE，用于说明解的公式、计算和绘图。	求解单个 PDE
`pdex2`	涉及不连续性的问题。	求解具有不连续性的 PDE
`pdex3`	需要计算偏导数的值的问题。	求解 PDE 并计算偏导数
`pdex4`	含两个 PDE 的方程组，其解在区间两端具有边界层，并且对于较小的 t 值，解的变化很快。	求解 PDE 方程组
`pdex5`	阶跃函数为初始条件的 PDE 方程组。	使用初始条件阶跃函数求解 PDE 方程组

参考

[1] Skeel, R. D. and M. Berzins, "A Method for the Spatial Discretization of Parabolic Equations in One Space Variable," SIAM Journal on Scientific and Statistical Computing, Vol. 11, 1990, pp. 1–32.