求解单个 PDE
此示例说明单个 PDE 的解的构成以及如何对解进行计算和绘图。
以如下偏微分方程为例
该方程的定义区间为 ,时间 。在 时,解满足初始条件
此外,在 和 时,解满足边界条件
要在 MATLAB® 中求解该方程,您需要对方程、初始条件和边界条件编写代码,然后在调用求解器 pdepe
之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写方程代码
在编写方程代码之前,您需要按照 pdepe
求解器所需的形式对其进行重写。pdepe
所需的标准形式是
.
以这种形式编写的 PDE 变为
.
将方程写作适当形式后,可知相关各项为:
现在,您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = pdex1pde(x,t,u,dudx)
:
x
是独立的空间变量。t
是独立的时间变量。u
是关于x
和t
微分的因变量。dudx
是偏空间导数 。输出
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
编写边界条件代码
现在,编写一个计算边界条件的函数。对于在区间 上提出的问题,边界条件应用于所有 以及 或 。求解器所需的边界条件的标准形式是
以这种标准形式重写边界条件,并读取系数值。
对于 ,方程为
系数是:
对于 ,方程为
系数是:
由于边界条件函数以 形式表示,并且此项已在主 PDE 函数中定义,因此不需要在边界条件函数中指定方程的此部分。您只需在每个边界处指定 和 的值。
边界函数应使用函数签名 [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
:
对于左边界,输入
xl
和ul
对应于 和 。对于右边界,输入
xr
和ur
对应于 和 。t
是独立的时间变量。对于左边界,输出
pl
和ql
对应于 和 (对于此问题,)。对于右边界,输出
pr
和qr
对应于 和 (对于此问题,)。
此示例中的边界条件由以下函数表示:
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
在求解器中的作用不同。尤其是解的成本和精确度很大程度上依赖于向量 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)
处计算的解 的第 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')
选择此问题的初始条件和边界条件,以得到解析解,由下式给出
绘制采用相同网格点的解析解。
surf(x,t,exp(-t)'*sin(pi*x)) title('True solution plotted with 20 mesh points') xlabel('Distance x') ylabel('Time t')
现在,比较在 (即 的最终值)处的数值解和解析解。在此示例中,。
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 %----------------------------------------------