40%20%20偏微分方程工具箱

1、 练习 40 偏微分方程工具箱 数学知识背景 与解常微分方程一样，求解偏微分方程只有在一定条件下，才能得出其解析解。在工程数学计算中，经常要求解偏微分方程。我们可以在一定初始条件和特殊情况下得到偏微分的数值解，对于一定形式的偏微分方程，MATLAB提供了求解对应偏微分方程的工具包和相应函数。 主要内容本练习讲述知识点 本练习主要考查在特定条件下求解偏微分方程的函数用法，而主要运用PDE（Partial Differential Equation）工具箱中的函数来求对偏微分方程进行求解。而PDE工具箱中，则主要提供了求解抛物线型、双曲线型及本征型偏微分方程的函数。 练习过程（1） 求解抛物线型的

2、函数主要是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(0)=1,

3、在其他条件下，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.00526316Time: 0.0

4、105263Time: 0.0157895Time: 0.0210526Time: 0.0263158Time: 0.0315789Time: 0.0368421Time: 0.0421053Time: 0.0473684Time: 0.0526316Time: 0.0578947Time: 0.0631579Time: 0.0684211Time: 0.0736842Time: 0.0789474Time: 0.0842105Time: 0.0894737Time: 0.0947368Time: 0.196 successful steps0 failed attempts194 funct

5、ion evaluations1 partial derivatives20 LU decompositions193 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

6、、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,tolco

7、nv,jmax,maxmul)函数主要对特征多项式（在区间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=

8、1，tolconv是期望的相对精度，缺省值为100*eps,这里dps为机器精度。Jmax是基向量的最大数目，求解中需要jmax*n的工作空间。Maxmul是Arnoldi运行次数，应是所有特征值的倍数。例：求解方程：在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)运行结果为： Basi

9、c=10, Time=0.65, New cov eig=0 Basic=19, Time=1.09, New cov eig=3Basic=28 , Time=1.59 , New cov eig=6Basic=37 , Time=2.25, New cov eig=9Basic=46 , Time= 3.07, New cov eig=13Basic= 55, Time=4.01 , New cov eig=27End of sweep: Basic= 55, Time= 4.01, New cov eig=27 Basic=37, Time=4.77, New cov eig=0 Bas

10、ic= 46 , Time= 5.54, New cov eig=0End 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=circleb2; c=1./sqrt(1+ux.2);a=0; f=0; rtol=le-3; p,e,t=initmes