求解 PDE 方程组
此示例说明由两个偏微分方程构成的方程组的解的构成,以及如何对解进行计算和绘图。
以如下 PDE 方程组为例
(函数 用作速记形式。)
该公式在区间 上对于时间 成立。初始条件为
边界条件为
要在 MATLAB® 中求解该方程,您需要对方程、初始条件和边界条件编写代码,然后在调用求解器 pdepe
之前选择合适的解网格。您可以将所需的函数作为局部函数包含在文件末尾(如本处所示),或者将它们作为单独的命名文件保存在 MATLAB 路径上的目录中。
编写方程代码
在编写方程代码之前,您需要确保它的形式符合 pdepe
求解器的要求:
在此形式中,PDE 系数是矩阵值,方程变为
因此,方程中的系数的值是
(仅对角线值)
现在,您可以创建一个函数以编写方程代码。该函数应具有签名 [c,f,s] = pdefun(x,t,u,dudx)
:
x
是独立的空间变量。t
是独立的时间变量。u
是关于x
和t
微分的因变量。它是二元素向量,其中u(1)
是 ,u(2)
是 。dudx
是偏空间导数 。它是二元素向量,其中dudx(1)
是 ,dudx(2)
是 。输出
c
、f
和s
对应于pdepe
所需的标准 PDE 形式中的系数。
因此,此示例中的方程可由以下函数表示:
function [c,f,s] = pdefun(x,t,u,dudx) c = [1; 1]; f = [0.024; 0.17] .* dudx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; end
(注意:所有函数都作为局部函数包含在示例的末尾。)
编写初始条件代码
接下来,编写一个返回初始条件的函数。初始条件应用在第一个时间值处,并为 x 的任何值提供 的值。初始条件的数量必须等于方程的数量,因此对于此问题,有两个初始条件。使用函数签名 u0 = pdeic(x)
编写函数。
初始条件为
对应的函数是
function u0 = pdeic(x) u0 = [1; 0]; end
编写边界条件代码
现在,编写计算以下边界条件的函数
对于在区间 上提出的问题,边界条件应用于所有 以及 或 。求解器所需的边界条件的标准形式是
以这种形式编写, 的偏导数的边界条件需要用通量 来表示。因此,此问题的边界条件是
对于 ,方程为
系数是:
同样,对于 ,方程是
系数是:
边界函数应使用函数签名 [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
:
对于左边界,输入
xl
和ul
对应于 和 。对于右边界,输入
xr
和ur
对应于 和 。t
是独立的时间变量。对于左边界,输出
pl
和ql
对应于 和 (对于此问题,)。对于右边界,输出
pr
和qr
对应于 和 (对于此问题,)。
此示例中的边界条件由以下函数表示:
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; end
选择解网格
当 较小时,此问题的解会快速变化。虽然 pdepe
选择了适合解析急剧变化的时间步,但要在输出绘图中显示该行为,您需要选择适当的输出时间。对于空间网格,在 两端的解中都存在边界层,因此您需要在那里指定网格点来解析急剧变化。
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
求解方程
最后,使用对称性值 、PDE 方程、初始条件、边界条件以及 和 的网格来求解方程。
m = 0; sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
pdepe
以三维数组 sol
形式返回解,其中 sol(i,j,k)
是在 t(i)
和 x(j)
处计算的解 的第 k
个分量的逼近值。将每个解分量提取到一个单独变量中。
u1 = sol(:,:,1); u2 = sol(:,:,2);
对解进行绘图
创建在 和 的所选网格点上绘制的 和 的解的曲面图。
surf(x,t,u1) title('u_1(x,t)') xlabel('Distance x') ylabel('Time t')
surf(x,t,u2) title('u_2(x,t)') xlabel('Distance x') ylabel('Time t')
局部函数
此处列出 PDE 求解器 pdepe
为计算解而调用的局部辅助函数。您也可以将这些函数作为它们自己的文件保存在 MATLAB 路径上的目录中。
function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve c = [1; 1]; f = [0.024; 0.17] .* dudx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; end % --------------------------------------------- function u0 = pdeic(x) % Initial Conditions u0 = [1; 0]; end % --------------------------------------------- function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; end % ---------------------------------------------