4.8.2  偏微分方程

在自然科学的很多领域内，都会遇到微分方程初值问题，特别是偏微分方程，它的定解问题是描述自然界及科学现象的最重要的工具。可以说，几乎自然界和各种现象都可以通过微分方程(特别是偏微分方程)来描述。

MATLAB提供了一个专门用于求解偏微分方程的工具箱PDE Toolbox。本小节仅介绍一些最简单、经典的偏微分方程，如椭圆型、双曲型、抛物型等偏微分方程，并给出求解方法。用户可以从中了解其解题的基本方法，从而解决类似的问题。

1．椭圆型问题

assempde函数是PDE工具箱中的一个基本函数，它使用有限元法组合PDE问题。该函数用来有选择地生成PDE问题的解。可以用assempde函数求解下面的标量椭圆型问题：

或系统椭圆型问题：

对于标量的情况，用解的列向量代表解矢量u，列矢量中的值对应于p的对应节点处的解。对于具有np个节点的N维系统，u1的前np行描述u的第1个元素，接下来的np行描述u的第2个元素，依次类推。这样，u的元素就作为N块节点行放到u中。assempde函数的调用语法如下。

(1)u=assempde(b,p,e,t,c,a,f)： 通过在线性方程组中剔除Dirichlet边界条件来组合和求解PDE问题。

(2)[K,F]=assempde(b,p,e,t,c,a,f)： 通过刚性位移近似Dirichlet边界条件来组合和求解PDE问题。K和F分别为刚性矩阵和右边项。PDE问题的有限元解为u1=K\F。

(3)[K,F,B,ud]=assempde(b,p,e,t,c,a,f)：通过从线性方程组剔除Dirichlet边界条件来组合PDE问题。u1=K\F，则返回非Dirichlet点上的解。完整的PDE问题可以通过MATLAB中的表达式u=B*u1+ud求解。

(4)[K,M,F,Q,G,H,R]=assempde(b,p,e,t,c,a,f)：给出PDE问题的分离表示。

(5)u=assempde(K,M,F,Q,G,H,R)：将PDE问题的分离表示转换为单一矩阵或矢量的形式，然后通过从线性方程组中剔除Dirichlet边界条件来求解。

(6)[K1,F1]=assempde(K,M,F,Q,G,H,R)：用很大的常数修改Dirichlet边界条件，从而将PDE问题的分离表示转换为单一矩阵或矢量的形式。

(7)[K1,F1,B,ud]=assempde(K,M,F,Q,G,H,R)：从线性方程组中剔除Dirichlet边界条件，从而将PDE问题的分离表示转换为单一矩阵或矢量的形式。

b描述PDE问题的边界条件。b也可以是边界条件矩阵或边界M文件的文件名。PDE问题几何模型由网络数据p,e,t描述。网格数据的生成可以查询help文档中的initmesh函数。

系数c，a，d，f可以通过多种方式给定。这些系数也可以与时间t相关，在assempde函数中可以看到所有选项的列表。

【例4-47】  求解L型薄膜的方程，为Dirichlet边界条件。最后绘图显示结果。

>> [p,e,t]=initmesh('lshapeg','hmax',0.2);    %生成初始三角形网格，hmax指网格大小

>> [p,e,t]=refinemesh('lshapeg',p,e,t);       %  将初始的三角形网格细化

>> u=assempde('lshapeb',p,e,t,1,0,1);         %  求解方程

>> pdesurf(p,t,u)                                %  绘制结果图形

lshapeg和lshapeb分别为表示对象几何模型和边界条件的M文件，为工具箱自带文件。initmesh函数和 refinemesh函数分别对网格模型进行初始化和细化。pdesurf函数绘制解的表面图。L型薄膜泊松方程的解如图4-15所示。

2．抛物型问题

MATLAB提供了parabolic函数求解标量抛物型问题：

图4-15  L型薄膜泊松方程的解

或系统PDE问题：

该函数的调用语法如下。

(1)u1=parabolic(u0,tlist,g,b,p,e,t,c,a,f,d)：p,e,t为网格数据，b为边界条件，初值为u0。

(2)u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d,rtol,atol)：atol和rtol为绝对和相对容限。

(3)u1=parabolic(u0,tlist,K,F,B,ud,M)：生成下面PDE问题的解。

的初始值为u0。

【例4-48】  求解热传导方程：。其中。在的区域内令，在其他区域内令。使用Dirichlet边界条件。要求计算时间linspace(0,0.1,20)处的解。

%  生成初始三角形网格数据

>> [p,e,t]=initmesh('squareg');                %  参数squareg指计算区域是方形的

>> [p,e,t]=refinemesh('squareg',p,e,t);       %  将初始的三角形网格细化

>> u0=zeros(size(p,2),1);

>> ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);    %  查找区域内的元素

>> u0(ix)=ones(size(ix));                       %  令

>> tlist=linspace(0,0.1,20);                   %  时间列表

>> u1=parabolic(u0,tlist,'squareb1',p,e,t,1,0,1,1);   %  求解偏微分方程

本例首先调用initmesh函数生成偏微分方程的初始网格，然后调用parabolic函数求解偏微分方程，运行结束将显示如下信息：

96 successful steps

0 failed attempts

194 function evaluations

1 partial derivatives

20 LU decompositions

193 solutions of linear systems

具体的结果u1是一个665*20的矩阵，这里就略去不显示了。

3．双曲型问题

MATLAB提供了hyperbolic函数来求解标量双曲型问题：

或系统双曲型问题：

hyperbolic函数的调用语法如下。

u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d,rtol,atol)：p,e,t为网格数据，b为边界条件，u0为初值，初始导数为ut0。atol和rtol为绝对和相对容限。

u1=hyperbolic(u0,ut0,tlist,K,F,B,ud,M)生成下面PDE问题的解。

u的初始值为u0和ut0。

在MATLAB命令窗口输入：

>> [p,e,t]=initmesh('squareg');                            %   生成初始三角形网格

>> x=p(1,:)';

>> y=p(2,:)';

>> u0=atan(cos(pi/2*x));

>> ut0=3*sin(pi*x).*exp(cos(pi*y));

>> tlist=linspace(0,5,31);                                  %  时间列表

>> uu=hyperbolic(u0,ut0,tlist,'squareb3',p,e,t,1,0,0,1);   %  求解方程

本例首先调用initmesh函数生成偏微分方程的初始网格，然后调用hyperbolic函数求解偏微分方程，运行结束将显示如下信息：

462 successful steps

70 failed attempts

1066 function evaluations

1 partial derivatives

156 LU decompositions

1065 solutions of linear systems

具体的结果uu是一个177*31的矩阵，这里就略去不显示了。