求解具有不连续性的 PDE
此示例说明如何求解涉及物质界面的 PDE。物质界面使得问题在 处具有不连续点,初始条件在右边界 处具有不连续点。
以如下分段 PDE 为例
初始条件为
边界条件为
要在 MATLAB® 中求解该方程,您需要对方程、初始条件和边界条件编写代码,然后在调用求解器 pdepe
之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写方程代码
在编写方程代码之前,您需要确保它的形式符合 pdepe
求解器的要求。pdepe
所需的标准形式是
.
在本例中,PDE 采用了正确形式,因此您可以读取系数的值。
通量项 和源项 的值根据 的值而变化。系数是:
现在,您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = pdex2pde(x,t,u,dudx)
:
x
是独立的空间变量。t
是独立的时间变量。u
是关于x
和t
微分的因变量。dudx
是偏空间导数 。输出
c
、f
和s
对应于pdepe
所需的标准 PDE 形式中的系数。根据输入变量x
、t
、u
和dudx
对这些系数编写代码。
因此,此示例中的方程可以由以下函数表示:
function [c,f,s] = pdex2pde(x,t,u,dudx) c = 1; if x <= 0.5 f = 5*dudx; s = -1000*exp(u); else f = dudx; s = -exp(u); end end
(注意:所有函数都作为局部函数包含在示例的末尾。)
编写初始条件代码
接下来,编写一个返回初始条件的函数。初始条件应用在第一个时间值处,并为 的任何值提供 的值。使用函数签名 u0 = pdex2ic(x)
编写函数。
初始条件为
对应的函数是
function u0 = pdex2ic(x) if x < 1 u0 = 0; else u0 = 1; end end
编写边界条件代码
现在,编写一个计算边界条件的函数。对于区间 上的问题,边界条件应用于所有 t 以及 或 的情形。求解器所需的边界条件的标准形式是
由于此示例具有球面对称性 (),pdepe
求解器会自动强制执行左边界条件以约束在原点的解,并忽略在边界函数中为左边界指定的任何条件。因此,对于左边界条件,您可以指定 。对于右边界条件,您可以用标准形式重写边界条件,并读取 和 的系数值。
对于 ,方程为 。系数是:
边界函数应使用函数签名 [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t)
:
对于左边界,输入
xl
和ul
对应于 u 和 x。对于右边界,输入
xr
和ur
对应于 u 和 x。t
是独立的时间变量。对于左边界,输出
pl
和ql
对应于 和 (对于此问题,)。对于右边界,输出
pr
和qr
对应于 和 (对于此问题,)。
此示例中的边界条件由以下函数表示:
function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) pl = 0; ql = 0; pr = ur - 1; qr = 0; end
选择解网格
空间网格应包括 附近的几个值以表示不连续界面,并包括 附近的点,因为在该点上具有不一致的初始值 () 和边界值 ()。对于较小的 ,解的变化很快,因此请使用可以解析这种急剧变化的时间步。
x = [0 0.1 0.2 0.3 0.4 0.45 0.475 0.5 0.525 0.55 0.6 0.7 0.8 0.9 0.95 0.975 0.99 1]; t = [0 0.001 0.005 0.01 0.05 0.1 0.5 1];
求解方程
最后,使用对称性值 m
、PDE 方程、初始条件、边界条件以及 x
和 t
的网格来求解方程。
m = 2; sol = pdepe(m,@pdex2pde,@pdex2ic,@pdex2bc,x,t);
pdepe
以三维数组 sol
形式返回解,其中 sol(i,j,k)
是在 t(i)
和 x(j)
处计算的解 的第 k
个分量的逼近值。sol
的大小是 length(t)
×length(x)
×length(u0)
,因为 u0
为每个解分量指定初始条件。对于此问题,u
只有一个分量,因此 sol
是 8×18 矩阵,但通常您可以使用命令 u = sol(:,:,k)
提取第 k
个解分量。
从 sol
中提取第一个解分量。
u = sol(:,:,1);
对解进行绘图
创建在 和 的所选网格点上绘制的解 的曲面图。由于 问题是在具有球面对称性的球面几何中提出的,因此解仅在径向 方向上变化。
surf(x,t,u) title('Numerical solution with nonuniform mesh') xlabel('Distance x') ylabel('Time t') zlabel('Solution u')
现在,只需绘制 和 即可获得曲面图中等高线的侧视图。在 处添加一条线,以突出材料接口的效果。
plot(x,u,x,u,'*') line([0.5 0.5], [-3 1], 'Color', 'k') xlabel('Distance x') ylabel('Solution u') title('Solution profiles at several times')
局部函数
此处列出 PDE 求解器 pdepe
为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function [c,f,s] = pdex2pde(x,t,u,dudx) % Equation to solve c = 1; if x <= 0.5 f = 5*dudx; s = -1000*exp(u); else f = dudx; s = -exp(u); end end %---------------------------------------------- function u0 = pdex2ic(x) %Initial conditions if x < 1 u0 = 0; else u0 = 1; end end %---------------------------------------------- function [pl,ql,pr,qr] = pdex2bc(xl,ul,xr,ur,t) % Boundary conditions pl = 0; ql = 0; pr = ur - 1; qr = 0; end %----------------------------------------------