matlab使用教程(30)—求解偏微分方程(1)

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

1.简介

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

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

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

        pdepe 还可求解某些可以简化成一维问题的二维和三维问题 。

2.求解一维 PDE

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

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

2.1 求解过程

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

• m 是对称常量。

• pdefun 定义要求解的方程。

• icfun 定义初始条件。

• bcfun 定义边界条件。

• xmesh 是 x 的空间值向量。

• tspan 是 t 的时间值向量。

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

2.2 方程的表示

        您必须按照 pdepe 所需的标准形式表示 PDE。以这种形式编写，您可以读取系数 c、f、s 的值。 在 MATLAB 中，您可以用以下形式的函数编写方程代码

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

2.3 初始条件

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

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

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

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

2.4 边界条件

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

        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 定义边界条件

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

2.5 积分选项

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

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

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

2.6 解的计算

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

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

3. 求解一维热方程

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

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

3.1 编写方程代码

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

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

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

3.2 代码初始条件

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

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

3.3 编写边界条件代码

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

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)];

3.3 求解方程

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

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

3.4 可视化运行结果

        使用 imagesc 可视化解矩阵。

colormap hot
imagesc(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')

        

3.5 其他所需函数

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

4.其他示例文件

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

odeexamples

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

edit exampleFileName.m

        要运行示例，请键入

exampleFileName

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