pdepe

求解一维抛物型和椭圆型 PDE

语法

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

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

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

说明

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) 使用一个空间变量 x 和时间 t 求解抛物型和椭圆型 PDE 方程组。至少一个方程必须为抛物型方程。标量 m 表示问题的对称性（平板、柱状或球面）。所求解的方程在 pdefun 中编码，初始值在 icfun 中编码，边界条件在 bcfun 中编码。对空间离散化所得的常微分方程 (ODE) 求积分，以求得 tspan 指定的时间处的近似解。pdepe 函数基于 xmesh 指定的网格返回解值。

示例

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) 还使用由 options（它是使用 odeset 函数创建的）定义的积分设置。例如，使用 AbsTol 和 RelTol 选项指定绝对误差容限和相对误差容限。

示例

[sol,tsol,sole,te,ie] = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) 还求 (t,u(x,t)) 的函数（称为事件函数）在何处为零。在输出中，te 是事件的时间，sole 是事件发生时的解，ie 是触发的事件的索引。tsol 是 tspan 中指定的、第一次终止事件之前的时间的列向量。

对于每个事件函数，应指定积分是否在零点处终止以及过零方向是否重要。为此，请将 odeset 的 Events 选项设置为函数（例如 @myEventFcn），并创建一个对应的函数：[value,isterminal,direction] = myEventFcn(m,t,xmesh,umesh)。xmesh 输入包含空间网格，umesh 是网格点上的解。

示例

全部折叠

求解柱坐标下的热方程

打开实时脚本

使用 pdepe 求解柱坐标下的热方程，并对解绘图。

在角对称的柱坐标中，热方程是

$\frac{\partial u}{\partial t} = \frac{1}{x} \frac{\partial}{\partial x} (x \frac{\partial u}{\partial x}) .$

该方程定义中，时间 $t \geq 0$ 且 $0 \leq x \leq 1$ 。初始条件根据 bessel 函数 $J_{0} (x)$ 及其第一个零点 $n = 2.404825557695773$ 定义为

$u (x, 0) = J_{0} (nx) .$

由于此问题采用柱坐标 (m = 1)，pdepe 会自动在 $x = 0$ 处强制应用对称条件。右边界条件为

$u (1, t) = J_{0} (n) e^{- n^{2} t} .$

初始条件和边界条件已选定为与问题的解析解一致，

$u (x, t) = J_{0} (nx) e^{- n^{2} t} .$

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

编写方程代码

在编写方程代码之前，您需要按照 pdepe 求解器所需的形式对其进行重写。pdepe 所需的标准形式是

$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 变为

$\frac{\partial u}{\partial t} = x^{- 1} \frac{\partial}{\partial x} (x^{1} \frac{\partial u}{\partial x})$

将方程写作适当形式后，可知相关各项为：

$m = 1$
$c (x, t, u, \frac{\partial u}{\partial x}) = 1$
$f (x, t, u, \frac{\partial u}{\partial x}) = \frac{\partial u}{\partial x}$
$s (x, t, u, \frac{\partial u}{\partial x}) = 0$

现在，您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = heatcyl(x,t,u,dudx)：

x 是独立的空间变量。
t 是独立的时间变量。
u 是关于 x 和 t 微分的因变量。
dudx 是偏空间导数 $\partial u / \partial x$ 。
输出 c、f 和 s 对应于 pdepe 所需的标准 PDE 形式中的系数。根据输入变量 x、t、u 和 dudx 对这些系数编写代码。

因此，此示例中的方程可以由以下函数表示：

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

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

编写初始条件代码

接下来，编写一个返回初始条件的函数。初始条件应用于第一个时间值 tspan(1)。该函数应具有签名 u0 = heatic(x)。

$u (x, 0) = J_{0} (nx)$ 的对应函数是

function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end

编写边界条件代码

现在，编写计算以下边界条件的函数

$u (1, t) = J_{0} (n) e^{- n^{2} t} .$

由于此问题采用柱坐标 (m = 1)，pdepe 会自动在 $x = 0$ 处强制应用对称条件，因此您不需要指定左边界条件。

对于在区间 $a \leq x \leq b$ 上提出的问题，边界条件应用于所有 $t$ 以及 $x = a$ 或 $x = b$ 。求解器所需的边界条件的标准形式是

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

以这种形式编写， $u$ 的偏导数的边界条件需要用通量 $f (x, t, u, \frac{\partial u}{\partial x})$ 来表示。因此，此问题的右边界条件是

$u (1, t) - J_{0} (n) e^{- n^{2} t} + 0 \cdot f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

边界函数应具有函数签名 [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)：

对于左边界，输入 xl 和 ul 对应于 $u$ 和 $x$ 。
对于右边界，输入 xr 和 ur 对应于 $u$ 和 $x$ 。
t 是独立的时间变量。
对于左边界，输出 pl 和 ql 对应于 $p_{L} (x, t, u)$ 和 $q_{L} (x, t)$ （对于此问题， $x = 0$ ）。
对于右边界，输出 pr 和 qr 对应于 $p_{R} (x, t, u)$ 和 $q_{R} (x, t)$ （对于此问题， $x = 1$ ）。

此示例中的边界条件由以下函数表示：

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end

选择解网格

在求解方程之前，需要指定希望用 pdepe 计算解的网格点 $(t, x)$ 。将点指定为向量 t 和 x。向量 t 和 x 在求解器中的作用不同。尤其是解的成本和精确度很大程度上依赖于向量 x 的长度。然而，计算对向量 t 中的值并不敏感。

对于此问题，请使用空间区间 [0,1] 内具有 25 个等间距点的网格来指定 x 和 t。

x = linspace(0,1,25);
t = linspace(0,1,25);

求解方程

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

m = 1;
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

pdepe 以三维数组 sol 形式返回解，其中 sol(i,j,k) 是在 t(i) 和 x(j) 处计算的解 $u_{k}$ 的第 k 个分量的逼近值。sol 的大小是 length(t)×length(x)×length(u0)，因为 u0 为每个解分量指定初始条件。对于此问题，u 只有一个分量，因此 sol 是 25×25 矩阵，但通常您可以使用命令 u = sol(:,:,k) 提取第 k 个解分量。

从 sol 中提取第一个解分量。

u = sol(:,:,1);

对解进行绘图

创建解的曲面图。由于问题表现为柱坐标中的一个圆盘，x 值反映圆盘上距中心一定距离处的温度，t 值反映特定位置的温度随时间变化的情况。

surf(x,t,u)
xlabel("x")
ylabel("t")
zlabel("u(x,t)")
view([150 25])

Figure contains an axes object. The axes object with xlabel x, ylabel t contains an object of type surface.

绘制圆盘中心 ( $x = 0$ ) 的温度变化。

plot(t,sol(:,1))
xlabel("Time")
ylabel("Temperature u(0,t)")
title("Temperature change at center of disc")

Figure contains an axes object. The axes object with title Temperature change at center of disc, xlabel Time, ylabel Temperature u(0,t) contains an object of type line.

局部函数

此处列出 PDE 求解器 pdepe 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。

function [c,f,s] = heatcyl(x,t,u,dudx)
c = 1;
f = dudx;
s = 0;
end
%----------------------------------------------
function u0 = heatic(x)
n = 2.404825557695773;
u0 = besselj(0,n*x);
end
%----------------------------------------------
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
n = 2.404825557695773;
pl = 0; %ignored by solver since m=1
ql = 0; %ignored by solver since m=1
pr = ur-besselj(0,n)*exp(-n^2*t);
qr = 0; 
end
%----------------------------------------------

求解振荡 PDE 并记录事件

打开实时脚本

求解偏微分方程，并使用事件函数记录振荡解中的过零点。

假设使用以下方程

$\frac{1}{x} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} (\frac{1}{t} u) .$

该方程定义中， $0 \leq x \leq 1$ 且 $t \geq 0.1$ 。初始条件为

$u (x, 0.1) = 1 .$

边界条件为

$u (0, t) = 1,$

$u (1, t) = \cos (π t) .$

此外，需要关注解的过零点。

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

编写方程代码

在编写方程代码之前，您需要按照 pdepe 求解器所需的形式对其进行重写。pdepe 所需的标准形式是

PDE 方程已采用如下形式：

$\frac{1}{x} \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} (\frac{1}{t} u) .$

由此可知相关各项为：

$m = 0$
$c (x, t, u, \frac{\partial u}{\partial x}) = \frac{1}{x}$
$f (x, t, u, \frac{\partial u}{\partial x}) = \frac{1}{t} u$
$s (x, t, u, \frac{\partial u}{\partial x}) = 0$

现在，您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = oscpde(x,t,u,dudx)：

x 是独立的空间变量。
t 是独立的时间变量。
u 是关于 x 和 t 微分的因变量。
dudx 是偏空间导数 $\partial u / \partial x$ 。
输出 c、f 和 s 对应于 pdepe 所需的标准 PDE 形式中的系数。根据输入变量 x、t、u 和 dudx 对这些系数编写代码。

因此，此示例中的方程可以由以下函数表示：

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

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

编写初始条件代码

接下来，编写一个返回初始条件的函数。初始条件应用于第一个时间值 tspan(1)。该函数应具有签名 u0 = oscic(x)。

$u (x, 0.1) = 1$ 的对应函数是

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

编写边界条件代码

现在，编写计算以下边界条件的函数

$u (0, t) = 1,$

$u (1, t) = \cos (π t) .$

对于在区间 $a \leq x \leq b$ 上提出的问题，边界条件应用于所有 $t$ 以及 $x = a$ 或 $x = b$ 。求解器所需的边界条件的标准形式是

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

采用此形式编写时，此问题的边界条件为

$u (0, t) - 1 + 0 \cdot f (x, t, u, \frac{\partial u}{\partial x}) = 0,$

$u (1, t) - \cos (π t) + 0 \cdot f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

边界函数应具有函数签名 [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)：

对于左边界，输入 xl 和 ul 对应于 $x$ 和 $u$ 。
对于右边界，输入 xr 和 ur 对应于 $x$ 和 $u$ 。
t 是独立的时间变量。
对于左边界，输出 pl 和 ql 对应于 $p_{L} (x, t, u)$ 和 $q_{L} (x, t)$ （对于此问题， $x = 0$ ）。
对于右边界，输出 pr 和 qr 对应于 $p_{R} (x, t, u)$ 和 $q_{R} (x, t)$ （对于此问题， $x = 1$ ）。

此示例中的边界条件由以下函数表示：

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

编写事件函数代码

使用事件函数记录积分中解的过零点。事件函数具有函数签名 [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)：

m 是坐标对称性，指定为 pdepe 的第一个输入。
t 是当前时间（标量）。
xmesh 是空间网格。
umesh 包含网格点上的解。
value 是需要关注的方程，通常用解 umesh 来表示。当 value 等于 0 时，事件发生。
isterminal 指定事件是否导致积分停止。如果 isterminal 为 0，则会记录事件，但积分不会停止。如果 isterminal 为 1，则当事件发生时，积分停止。
direction 指定过零的方向。如果为 1，则只有具有正斜率的过零点才会触发事件。如果为 -1，过零点必须具有负斜率。如果为 0，则任何过零点都会触发事件。

每次积分时，求解器都会调用事件函数来检查过零点。要记录所有过零点，value 应关注解向量 umesh 中的符号变化。将 isterminal 和 direction 指定为大小相同的零向量，因为此示例中的事件不是终止事件，且过零可以以任何斜率出现。

此问题的事件函数是

function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end

选择解网格

在求解方程之前，需要指定希望用 pdepe 计算解的网格点 $(t, x)$ 。对于此问题，请使用区间 $0 \leq x \leq 1$ 和 $0.1 \leq t \leq 1$ 内具有 50 个点的精细网格。精细网格能够清晰的呈现振荡解。

x = linspace(0,1,50);
t = linspace(0.1,pi,50);

求解方程

最后，使用对称性 m、PDE 方程、初始条件、边界条件、事件函数以及 x 和 t 的网格来求解方程。使用 odeset 创建引用事件函数的 options 结构体，并将该结构体作为 pdepe 的最后一个输入参量传入。指定五个输出参量以返回来自事件函数和求解器的信息：

sol 是用 pdepe 计算的解。
tsol 是终止事件前的时间向量。当没有发生终止事件时，tsol 等于 t。
sole 是在每个事件时间点处的解。
te 是每个事件的时间。
ie 是每个事件的索引。由于在事件函数中 values = umesh，ie 也给出了每个时间步上触发事件的 umesh 的索引。

m = 0;
options = odeset(Events=@pdevents);
[sol,tsol,sole,te,ie] = pdepe(m,@oscpde,@oscic,@oscbc,x,t,options);

将解提取为矩阵 u。

u = sol(:,:,1);

对解进行绘图

创建解的曲面图，并采用俯视绘图。

surf(x,t,u)
view(2)

Figure contains an axes object. The axes object contains an object of type surface.

绘制发生事件的点，以曲面 $f (x, t) = 0$ 作为参考。输出索引向量 ie 有助于找出事件位置。表达式 x(ie)' 给出事件发生处的 x 值，表达式 sole(x==x(ie)') 给出对应的解值。

view([39 30])
xlabel("x")
ylabel("t")
zlabel("u(x,t)")
hold on
plot3(x(ie)',te,sole(x==x(ie)'),'r*')
surf(x,t,zeros(size(u)),EdgeColor="flat")
hold off

Figure contains an axes object. The axes object with xlabel x, ylabel t contains 3 objects of type surface, line. One or more of the lines displays its values using only markers

局部函数

此处列出 PDE 求解器 pdepe 为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。

function [c,f,s] = oscpde(x,t,u,dudx)
c = 1/x;
f = u/t;
s = 0;
end
%----------------------------------------------
function u0 = oscic(x)
u0 = 1;
end
%----------------------------------------------
function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)
pl = ul - 1;
ql = 0;
pr = ur - cos(pi*t);
qr = 0;
end
%----------------------------------------------
function [value, isterminal, direction] = pdevents(m,t,xmesh,umesh)
value = umesh;
isterminal = zeros(size(umesh));
direction = zeros(size(umesh));
end
%----------------------------------------------

输入参数

全部折叠

`m` — 对称常量
`0` | `1` | `2`

对称常量，指定为下表中的值之一。m 在 pdefun 指定的方程中指定问题类型。在按照 pdepe 所需的形式重写方程后，您可以从方程中得到 m 的值。

值描述通量的拉普拉斯算子

值	描述	通量的拉普拉斯算子
`0`	没有对称性的一维笛卡尔坐标	$Δ f = \frac{\partial^{2} f}{\partial x^{2}}$
`1`	方位角对称的二维柱坐标	$Δ f = \frac{1}{ρ} \frac{\partial}{\partial ρ} (ρ \frac{\partial f}{\partial ρ}) + \frac{1}{ρ^{2}} \frac{\partial^{2} f}{\partial φ^{2}}$ , 其中 $\frac{\partial f}{\partial φ} = 0$ （方位角对称）。
`2`	方位角和天顶角对称的三维球面坐标	$Δ f = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial f}{\partial r}) + \frac{1}{r^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial f}{\partial θ}) + \frac{1}{r^{2} \sin^{2} θ} \frac{\partial^{2} f}{\partial φ^{2}}$ , 其中 $\frac{\partial f}{\partial θ} = \frac{\partial f}{\partial φ} = 0$ （方位角和天顶角对称）。

0

没有对称性的一维笛卡尔坐标

Δ f = \frac{\partial^{2} f}{\partial x^{2}}

1

方位角对称的二维柱坐标

$Δ f = \frac{1}{ρ} \frac{\partial}{\partial ρ} (ρ \frac{\partial f}{\partial ρ}) + \frac{1}{ρ^{2}} \frac{\partial^{2} f}{\partial φ^{2}}$ ,

其中 $\frac{\partial f}{\partial φ} = 0$ （方位角对称）。

2

方位角和天顶角对称的三维球面坐标

$Δ f = \frac{1}{r^{2}} \frac{\partial}{\partial r} (r^{2} \frac{\partial f}{\partial r}) + \frac{1}{r^{2} \sin θ} \frac{\partial}{\partial θ} (\sin θ \frac{\partial f}{\partial θ}) + \frac{1}{r^{2} \sin^{2} θ} \frac{\partial^{2} f}{\partial φ^{2}}$ ,

其中 $\frac{\partial f}{\partial θ} = \frac{\partial f}{\partial φ} = 0$ （方位角和天顶角对称）。

当 m > 0 时，pdepe 要求左空间边界 a ≥ 0，并且它会自动应用左边界条件以解释原点处的坐标不连续性。在此情况下，pdepe 忽略 bcfun 中指定的任何左边界条件。

`pdefun` — 方程的 PDE 函数
函数句柄

方程的 PDE 函数，指定为定义 PDE 系数的函数句柄。

pdefun 对以下形式的 PDE 的系数进行编码

方程中的项是：

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

PDE 在 t₀ ≤ t ≤ t_f 和 a ≤ x ≤ b 时成立。值 tspan(1) 和 tspan(end) 对应于 t₀ 和 t_f，而 xmesh(1) 和 xmesh(end) 对应于 a 和 b。[a,b] 必须是有限区间。m 可以是 0、1 或 2，分别对应平板、柱状或球面对称性。如果 m > 0，则必须有 a ≥ 0。

编写方程代码

编写函数 pdefun，它计算系数 c、f 和 s 的值。使用函数签名

[c,f,s] = pdefun(x,t,u,dudx)

pdefun 的输入参量是标量 x 和 t 以及向量 u 和 dudx。向量 u 逼近解 u，dudx 逼近它关于 x 的偏导数。输出参量 c、f 和 s 是列向量，其元素数等于方程数。c 存储矩阵 c 的对角线元素。

示例

热方程 $\frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}}$ 可以按照求解器所需的形式重写为

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

写作此形式时，可知系数值 m = 0、c = 1、f = ∂u/∂x 和 s = 0。此方程的函数是：

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

数据类型: function_handle

`icfun` — 初始条件函数
函数句柄

初始条件函数，指定为定义初始条件的函数句柄。

对于 t = t₀ 和所有 x，解分量均满足以下形式的初始条件：

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

编写初始条件代码

编写定义初始条件的函数 icfun。使用函数签名：

function u0 = icfun(x)

当与参量 x 一起调用时，icfun 函数会计算并返回列向量 u0 中 x 处的解分量的初始值。u0 中元素的数量等于方程的数量。

示例

恒定初始条件 $u (x, 0) = 1 / 2$ 对应于以下函数：

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

数据类型: function_handle

`bcfun` — 边界条件函数
函数句柄

边界条件函数，指定为定义边界条件的函数句柄。

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

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

q 的元素全部为零或都不为零。请注意，边界条件以通量 f 的方式而不是 ∂u/∂x 表示。同时，在这两个系数之间，只有 p 可以依赖于 u。

编写边界条件代码

编写函数 bcfun，它定义边界条件的项 p 和 q。使用函数签名

function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

ul 是左边界 xl = a 处的近似解，ur 是右边界 xr = b 处的近似解。
pl 和 ql 是列向量，对应于在 xl 处计算的 p 和 q。同样，pr 和 qr 则与 xr 相对应。
输出 pl、ql、pr 和 qr 中的元素数量等于方程的数量。

当 m > 0 且 a = 0 时，由于解在 x = 0 附近具备有界性，通量 f 需在 a = 0 处消失。pdepe 会自动应用此边界条件，并且忽略 pl 和 ql 中返回的值。

示例

以区间 $0 \leq x \leq 1$ 上的边界条件 $u (0, t) = 0, u (1, t) = 1$ 为例。边界条件以求解器所需形式表示为

$\begin{array}{l} u (0, t) + 0 \cdot f (0, t, u, \frac{\partial u}{\partial x}) = 0, \\ u (1, t) - 1 + 0 \cdot f (0, t, u, \frac{\partial u}{\partial x}) = 0. \end{array}$

在这种形式下，很明显，两个 q 项都为零。p 项非零，以反映 u 的值。对应的函数是

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

数据类型: function_handle

`xmesh` — 空间网格
向量

空间网格，指定为向量 [x0 x1 ... xn]，其中指定了要基于 tspan 中的每个值求数值解的各个点。在 xmesh 中指定的网格上进行解的二阶逼近。pdepe 不会自动选择 x 中的网格。必须在 xmesh 中提供近似的固定网格。

xmesh 的元素必须满足 x0 < x1 < ... < xn。通常，最好使用解快速变化的近距网格点。xmesh 的长度必须为 >=3，解的计算成本在很大程度上依赖于 xmesh 的长度。

对于具有不连续点的问题，应在不连续点处设置网格点，以使问题在每个子区间内均呈平滑。但是，当 m > 0 时，无需使用 x = 0 附近的精细网格来说明坐标奇异性。

示例: xmesh = linspace(0,5,25) 使用一个在 0 到 5 之间包含 25 个点的空间网格。

数据类型: single | double

`tspan` — 积分的时间跨度
向量

积分的时间跨度，指定为向量 [t0 t1 ... tf]，其中指定了要基于 xmesh 中的每个值求解的各个点。函数 pdepe 使用可动态选择时间步长和公式的 ODE 求解器来计算时间积分。tspan 的元素仅指定您需要解的位置，因此 tspan 的长度对解的计算成本影响甚微。

tspan 的元素必须满足 t0 < t1 < ... < tf。tspan 的长度必须 >=3。

示例: tspan = linspace(0,5,5) 使用 0 到 5 之间的五个时间点。

数据类型: single | double

`options` — Options 结构体
结构体数组

Options 结构体，指定为结构体数组。使用 odeset 函数创建或修改 options 结构体。您可以设置 odeset 中的任何选项，但在求解 PDE 时 pdepe 会覆盖以下选项：

Mass
MStateDependence
MvPattern
MassSingular
Jacobian
JPattern
Vectorized

示例: options = odeset(RelTol=1e-5) 指定 1e-5 的相对误差容限。

数据类型: struct

输出参量

全部折叠

`sol` — 解数组
三维数组

解数组，以三维数组形式返回。

pdepe 以多维数组形式返回解。u_i = ui = sol(:,:,i) 是解向量 u 的第 i 个分量的逼近值。元素 ui(j,k) = sol(j,k,i) 在 (t,x) = (tspan(j),xmesh(k)) 处近似于 u_i。

ui = sol(j,:,i) 在 tspan(j) 时和网格点 xmesh(:) 上近似于分量 i。使用 pdeval 计算 xmesh 中未包含的点的逼近值及其偏导数 ∂u_i/∂x。有关详细信息，请参阅 pdeval。

`tsol` — 解时间
列向量

解时间，以 tspan 中指定的、第一次终止事件之前的时间的列向量形式返回。

`sole` — 事件时间处的解
数组

事件时间的解，以数组形式返回。te 中的事件时间对应于 sole 中返回的解，而 ie 指定发生了哪个事件。

`te` — 事件的时间
列向量

事件的时间，以列向量形式返回。te 中的事件时间对应于 sole 中返回的解，而 ie 指定发生了哪个事件。

`ie` — 触发的事件函数的索引
列向量

触发的事件函数的索引，以列向量形式返回。te 中的事件时间对应于 sole 中返回的解，而 ie 指定发生了哪个事件。

提示

如果 uji = sol(j,:,i) 在时间 tspan(j) 和网格点 xmesh 处近似于解的分量 i，则 pdeval 计算在点数组 xout 处的逼近值及其偏导数 ∂u_i/∂x，并将其在 uout 和 duoutdx 中返回：[uout,duoutdx] = pdeval(m,xmesh,uji,xout)。pdeval 函数计算偏导数 ∂u_i/∂x，而不是通量。通量是持续的，但在物质界面处偏导数可能有跳跃。

算法

时间积分使用 ode15s 求解器完成。pdepe 利用 ode15s 的功能解算当 PDE 包含椭圆型方程时产生的微分代数方程，并用来处理具有指定稀疏类型的雅可比矩阵。

经离散化处理后，椭圆型方程会生成代数方程。如果与椭圆型方程对应的初始条件向量的元素与离散化的结果不一致，pdepe 会在开始时间积分之前尝试调整这些元素。因此，为初始时间返回的解可能具有与其他时间类似的离散化误差。如果网格足够精细，pdepe 可找到接近于给定条件的一致初始条件。如果 pdepe 显示消息称找不到一致初始条件，请尝试细化网格。对于与抛物型方程对应的初始条件向量元素，则无需调整。

参考

[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.

扩展功能

全部展开

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

版本历史记录

在 R2006a 之前推出

全部展开

R2025a: 指定输出函数和其他 `odeset` 选项

现在，您可以在使用 pdepe 函数时指定所有 odeset 选项，包括输出函数。但是，pdepe 在定义要求解的 PDE 时会覆盖一些选项。有关详细信息，请参阅 options。

以前，pdepe 仅支持 RelTol、AbsTol、NormControl、InitialStep、MaxStep 和 Events 选项。

另请参阅

pdeval | ode15s | odeset | odeget

pdepe

语法

说明

示例

求解柱坐标下的热方程

求解振荡 PDE 并记录事件

输入参数

m — 对称常量 0 | 1 | 2

pdefun — 方程的 PDE 函数 函数句柄

编写方程代码

示例

icfun — 初始条件函数 函数句柄

编写初始条件代码

示例

bcfun — 边界条件函数 函数句柄

编写边界条件代码

示例

xmesh — 空间网格 向量

tspan — 积分的时间跨度 向量

options — Options 结构体 结构体数组

输出参量

sol — 解数组 三维数组

tsol — 解时间 列向量

sole — 事件时间处的解 数组

te — 事件的时间 列向量

ie — 触发的事件函数的索引 列向量

提示

算法

参考

扩展功能

基于线程的环境 使用 MATLAB® backgroundPool 在后台运行代码或使用 Parallel Computing Toolbox™ ThreadPool 加快代码运行速度。

版本历史记录

R2025a: 指定输出函数和其他 odeset 选项

另请参阅

主题

`m` — 对称常量
`0` | `1` | `2`

`pdefun` — 方程的 PDE 函数
函数句柄

`icfun` — 初始条件函数
函数句柄

`bcfun` — 边界条件函数
函数句柄

`xmesh` — 空间网格
向量

`tspan` — 积分的时间跨度
向量

`options` — Options 结构体
结构体数组

`sol` — 解数组
三维数组

`tsol` — 解时间
列向量

`sole` — 事件时间处的解
数组

`te` — 事件的时间
列向量

`ie` — 触发的事件函数的索引
列向量

基于线程的环境
使用 MATLAB® `backgroundPool` 在后台运行代码或使用 Parallel Computing Toolbox™ `ThreadPool` 加快代码运行速度。

R2025a: 指定输出函数和其他 `odeset` 选项