Main Content

求解单个 PDE

此示例说明单个 PDE 的解的构成以及如何对解进行计算和绘图。

以如下偏微分方程为例

π2ut=2ux2.

该方程的定义区间为 0x1,时间 t0。在 t=0 时,解满足初始条件

u(x,0)=sin(πx).

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

u(0,t)=0,

πe-t+ux(1,t)=0.

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

编写方程代码

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

c(x,t,u,ux)ut=x-mx(xmf(x,t,u,ux))+s(x,t,u,ux).

以这种形式编写的 PDE 变为

π2ut=x0x(x0ux)+0.

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

m=0

c(x,t,u,ux)=π2

f(x,t,u,ux)=ux

s(x,t,u,ux)=0

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

  • x 是独立的空间变量。

  • t 是独立的时间变量。

  • u 是关于 xt 微分的因变量。

  • dudx 是偏空间导数 u/x

  • 输出 cfs 对应于 pdepe 所需的标准 PDE 形式中的系数。根据输入变量 xtududx 对这些系数编写代码。

因此,此示例中的方程可以由以下函数表示:

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

编写边界条件代码

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

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

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

对于 x=0,方程为

u(0,t)=0u+0ux=0.

系数是:

  • p(0,t,u)=u

  • q(0,t)=0

对于 x=1,方程为

πe-t+ux(1,t)=0πe-t+1ux(1,t)=0.

系数是:

  • p(1,t,u)=πe-t

  • q(1,t)=1

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

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

  • 对于左边界,输入 xlul 对应于 ux

  • 对于右边界,输入 xrur 对应于 ux

  • t 是独立的时间变量。

  • 对于左边界,输出 plql 对应于 p(x,t,u)q(x,t)(对于此问题,x=0)。

  • 对于右边界,输出 prqr 对应于 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)。将点指定为向量 tx。向量 tx 在求解器中的作用不同。尤其是解的成本和精确度很大程度上依赖于向量 x 的长度。然而,计算对向量 t 中的值并不敏感。

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

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

求解方程

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

m = 0;
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

pdepe 以三维数组 sol 形式返回解,其中 sol(i,j,k) 是在 t(i)x(j) 处计算的解 uk 的第 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')

Figure contains an axes object. The axes object with title Numerical solution computed with 20 mesh points, xlabel Distance x, ylabel Time t contains an object of type surface.

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

u(x,t)=e-tsin(πx).

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

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

Figure contains an axes object. The axes object with title True solution plotted with 20 mesh points, xlabel Distance x, ylabel Time t contains an object of type surface.

现在,比较在 tf(即 t 的最终值)处的数值解和解析解。在此示例中,tf=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)')

Figure contains an axes object. The axes object with title Solution at t = 2, xlabel Distance x, ylabel u(x,2) contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Numerical, 20 mesh points, Analytical.

局部函数

此处列出 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
%----------------------------------------------

另请参阅

相关主题