40 偏微分方程工具箱.

练习 40 偏微分方程工具箱
数学知识背景
与解常微分方程一样，求解偏微分方程只有在一定条件下，才能得出其解析解。

在工程数学计算中，经常要求解偏微分方程。

我们可以在一定初始条件和特殊情况下得到偏微分的数值解，对于一定形式的偏微分方程，MA TLAB 提供了求解对应偏微分方程的工具包和相应函数。

主要内容
本练习讲述知识点
本练习主要考查在特定条件下求解偏微分方程的函数用法，而主要运用PDE （Partial Differential Equation ）工具箱中的函数来求对偏微分方程进行求解。

而PDE 工具箱中，则主要提供了求解抛物线型、双曲线型及本征型偏微分方程的函数。

练习过程
（1） 求解抛物线型的函数主要是parabolic,此函数主要是用有限元法求解抛物线型微分方
程以及微分方程组，其用法为：
u1= parabolic(u0,tlist,b,p,e,t,c,a,f,d)
网格参数主要为p 、e 和t,边界条件b 可以用矩阵形式也可以用m 文件格式，主要依赖于时间t,方程参数c 、a 、d 、f 也可是时间的函数。

Tlist 是时间序列，可以在函数末加入误差限来控制。

对于标量形式的偏微分方程，函数返回值为一个矩阵。

u1= parabolic(u0,tlist,K,F,B,ud,M)
这个函数主要用于求解ODE 问题，其中初值为u0。

例：求热传导方程
u t
u ∆=∂∂ 求解的范围为正方形区域： 1,1≤≤-y x ，
初值条件：当12
2≤+y x 时，u(0)=1,在其他条件下，u(0)=0。

在命令区中键入下命令得：
[p,e, t]=initmesh('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);
求得的结果为：
Time: 0.00526316
Time: 0.0105263
Time: 0.0157895
Time: 0.0210526
Time: 0.0263158
Time: 0.0315789
Time: 0.0368421
Time: 0.0421053
Time: 0.0473684
Time: 0.0526316
Time: 0.0578947
Time: 0.0631579
Time: 0.0684211
Time: 0.0736842
Time: 0.0789474
Time: 0.0842105
Time: 0.0894737
Time: 0.0947368
Time: 0.1
96 successful steps
0 failed attempts
194 function evaluations
1 partial derivatives
20 LU decompositions
193 solutions of linear systems
即经过了96步计算，失败次数为0，194次的函数赋值，1次求偏导，20次求LU 分解运算，193组线形系统的解。

（2） 求解双曲线型偏微分方程的主要函数是hyperbolic ，它也是用有限元的方法来求解偏
微分方程或者方程组，其用法为：
u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d)
u0是初始值，网格参数为b 、e 、t,边界条件b 可用矩阵形式也可以用m 文件格式，它依赖于时间，而方程系数c 、a 、d 、f 也可以是时间的函数，可以对函数设置误差限。

u1=hyperbolic(u0,ut0,tlist,K ，F ，B ，ud ，M)
可以用于求解ODE 问题，初值为u0。

（3） 求解特征值主要用函数pdeeig ，函数主要用有限元法来求解特征方程，其用法为：
[v,1]=pdeeig(b,p,e,t,c,a,d,r)
参数p 、e 、t 描述区域，参数b 描述边界条件，r 两个元素的数组。

所得结果中，v 是特征向量矩阵，对于标量形式的特征值方程，v 对应于网格节点的特征值。

求解一般稀疏特征值问题的函数主要是sptarn ，其主要用法是：
[xv,lmb,iresult]=sptarn(a,b,lb,ub,spd,tolconv,jmax,maxmul)
函数主要对特征多项式（0)=-x B A λ在区间[lb ，ub]上的特征值，A 和B 是稀疏矩阵，lb 和ub 分别是要求解的特征值的上界与下界。

若要求解的是ub 左边的所有特征值，则可以令
lb=-inf,若要求解的是lb 右边的所有特征值，则令ub=inf.对于窄区间，可以较快得到结果。

在复数情况下，比较lmb 、lb 、ub 的实数部分。

Xv 是特征向量，它的值使得判断式（a*xv-b*xv*diag(lmb)）最小。

Lmb 是对应的特征值，当iresult 大于或者等于0时，可以进行求解并找到所有特征值，而当iresult 小于0 时，则求解可能是不完全的，有可能有特征值未求出。

如果特征值均为正值，则spd=1，tolconv 是期望的相对精度，缺省值为100*eps, 这里dps 为机器精度。

Jmax 是基向量的最大数目，求解中需要jmax*n 的工作空间。

Maxmul 是Arnoldi 运行次数，应是所有特征值的倍数。

例：求解方程：u u λ=∇-
在L 型区域上的小于100的特征值及其相应的特征模态，命令窗口中演示有：
[p,e,t]=initmesh('lshapeg');
[p,e,t]=refinemesh('lshapeg',p,e,t);
[v,l]=pdeeig('lshapeb',p,e,t,1,0,1,[-Inf100]);
pdesurf(p,t,v(:,16))
运行结果为：
Basic=10, Time=0.65, New cov eig=0
Basic=19, Time=1.09, New cov eig=3
Basic=28 , Time=1.59 , New cov eig=6
Basic=37 , Time=2.25, New cov eig=9
Basic=46 , Time= 3.07, New cov eig=13
Basic= 55, Time=4.01 , New cov eig=27
End of sweep: Basic= 55, Time= 4.01, New cov eig=27
Basic=37, Time=4.77, New cov eig=0
Basic= 46 , Time= 5.54, New cov eig=0
End of sweep: Basic=46, Time= 5.54, New cov eig=0
(4)求解非线性偏微分方程可以用函数pdenonlin ，其用法为：
[u ，res]=pdenonlin （b,p,e,t,c,a,f ）
这个函数主要用来求解非线性标量形式的偏微分方程，其中，c ，a ，f 是依赖于u 的函数，该函数主要用牛顿迭代法求解。

在参数中，主要用于设置方程的迭代次数，迭代中止误差或者初解等。

例：求解最小表面积的问题，在命令框中输入：
g=’circleg ’; b=’circleb 2’; c=’1./sqrt(1+ux.^2)’;
a=0; f=0; rtol=le-3; [p,e,t]=initmesh(g); [p,e,t]=refinemesh(g,p,e,t); u=pdenonlin(b,p,e,t,c,a,f,’tol ’,rtol);
pdesurf(p,t,u);
练习小结
本练习介绍了求解偏微分方程特例的求解，对在一定条件下的偏微分方程MATLAB 有对应的函数来求解，读者在看完本练习后应复习掌握各个函数的用法。