MATLAB 数学应用 微分方程 一维偏微分方程 求解单个PDE

本文说明了单个 PDE 的解的构成以及如何对解进行计算和绘图。

以如下偏微分方程为例

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

该方程的定义区间为 0≤x≤1，时间 t≥0。在 t=0 时，解满足初始条件

$u(x,0)=sin(\pi x).$

此外，在 x=0 和 x=1 时，解满足边界条件

$u(0,t)=0,$

$\pi e^{-t}+\frac{\partial u}{\partial x}(1,t)=0.$

要在 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 变为

${\pi }^2\frac{\partial u}{\partial t}=x^0\frac{\partial }{\partial x}(x^0\frac{\partial u}{\partial x})+0。$

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

$m=0$

$c(x,t,u,\frac{\partial u}{\partial x})={\pi }^2$

$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] = pdex1pde(x,t,u,dudx)：

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

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

function [c,f,s] = pdex1pde(x,t,u,dudx)
c = pi^2;
f = dudx;
s = 0;
end

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

代码初始条件

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

对应的函数是

function u0 = pdex1ic(x)
u0 = sin(pi*x);
end

编写边界条件代码

现在，编写一个计算边界条件的函数。对于在区间 a≤x≤b 上提出的问题，边界条件应用于所有 t 以及 x=a 或 x=b。求解器所需的边界条件的标准形式是

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

以这种标准形式重写边界条件，并读取系数值。

对于 x=0，方程为

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

系数是：

• p(0,t,u)=u
• q(0,t)=0

对于 x=1，方程为

$\pi e^{-t}+\frac{\partial u}{\partial x}(1,t)=0\to \pi e^{-t}+1\cdot \frac{\partial u}{\partial x}(1,t)=0.$

系数是：

• $p(1,t,u)=\pi e^{-t}$
• $q(1,t)=1$

由于边界条件函数以 $f(x,t,u,\frac{\partial u}{\partial x})$ 形式表示，并且此项已在主 PDE 函数中定义，因此不需要在边界条件函数中指定方程的此部分。您只需在每个边界处指定 $p(x,t,u)$ 和 $q(x,t)$ 的值。

边界函数应使用函数签名 [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)：

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

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

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

选择解网格

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

对于此问题，请使用一个网格，该网格具有 20 个位于空间区间 [0,1] 中的等距点、5 个位于时间区间 [0,2] 中的 t 值。

x = linspace(0,1,20);
t = linspace(0,2,5);

求解方程

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

m = 0;
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,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 是 5×20 矩阵，但通常您可以使用命令 u = sol(:,:,k) 提取第 k 个解分量。

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

u = sol(:,:,1);

对解进行绘图

创建解的曲面图。

surf(x,t,u)
title('Numerical solution computed with 20 mesh points')
xlabel('Distance x')
ylabel('Time t')

选择此问题的初始条件和边界条件，以得到解析解，由下式给出

$u(x,t)=e^{-t}sin(\pi x).$

绘制采用相同网格点的解析解。

surf(x,t,exp(-t)'*sin(pi*x))
title('True solution plotted with 20 mesh points')
xlabel('Distance x')
ylabel('Time t')

现在，比较在 $t_f$（即 t 的最终值）处的数值解和解析解。在此示例中，$t_f=2。$

plot(x,u(end,:),'o',x,exp(-t(end))*sin(pi*x))
title('Solution at t = 2')
legend('Numerical, 20 mesh points','Analytical','Location','South')
xlabel('Distance x')
ylabel('u(x,2)')

局部函数

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

function [c,f,s] = pdex1pde(x,t,u,dudx) % Equation to solve
c = pi^2;
f = dudx;
s = 0;
end
%----------------------------------------------
function u0 = pdex1ic(x) % Initial conditions
u0 = sin(pi*x);
end
%----------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % Boundary conditions
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;
end
%----------------------------------------------