MATLAB的偏微分方程工具箱PDEToolbox

第二十二章数学物理方程的计算机仿真求解本章将主要讲述如何用MATLAB实现对偏微分方程的仿真求解．MATLAB的偏微分方程工具箱（PDEToolbox）的出现，为偏微分方程的求解以及定性研究提供了捷径．主要步骤为：(1)设置PDE的定解问题．即设置二维定解区域、边界条件以及方程的形式和系数；(2)用有限元法（FEM）求解PDE．即网格的生成、方程的离散以及求出数值解；(3)解的可视化．用PDEToolbox可以求解的基本方程有：椭圆方程、抛物方程、双曲方程、特征值方程、椭圆方程组以及非线性椭圆方程．具体操作上可用两个途径：(1)直接使用图形用户界面（GraphicalUserInterface，简记作GUI）求解．计算机仿真求解的偏微分方程类型分为：椭圆型方程：()cuauf抛物型方程：()udcuauft双曲型方程：22()udcuauft特征值问题：()cuaudu特征值偏微分方程中不含参数f．（2）用M文件编程求解．本章首先对可视化方法（GUI）求解作初步介绍，然后详细介绍用M文件编程解几类基本偏微分方程，并对典型偏微分的解的静态（或动态）显示曲线分布进行了讨论。22.1用偏微分方程工具箱求解微分方程直接使用图形用户界面（GraphicalUserInterface，简记作GUI）求解．【解】第一步：启动MATLAB，键入命令pdetool并回车，就进入GUI．在Options菜单下选择Grid命令，打开栅格．栅格使用户容易确定所绘图形的大小．第二步：选定定解区域本题为自定区域：自拟定解区域如图22.1所示：E1-E2+R1-E3．具体用快捷工具分别画椭圆E1、圆E2、矩形R1、圆E3．然后在Setformula栏中进行编辑并用算术运算符将图形对象名称连接起来.(或删去默认的表达式，直接键入E1-E2+R1-E3)例22.1.1解热传导方程tuuf边界条件是齐次类型，定解区域自定。计算机仿真图22.1所讨论定解问题的区域第三步：选取边界首先选择Boundary菜单中BoundaryMode命令，进入边界模式．然后单击Boundary菜单中RemoveAllSubdomainBorders选项，从而去掉子域边界，如图22.2．单击Boundary菜单中SpecifyBoundaryConditions选项，打开BoundaryConditions对话框，输入边界条件．本例取默认条件，即将全部边界设为齐次Dirichlet条件，边界显示为红色．如果想将几何与边界信息存储，可选择Boundary菜单中的ExportDecomposedGeometry,BoundaryCond’s命令，将它们分别存储在g、b变量中，并通过MATLAB形成M文件．图22.2定解问题的边界第四步：设置方程类型选择PDE菜单中PDEMode命令，进入PDE模式．再单击PDE菜单中PDESecification选项，打开PDESecification对话框，设置方程类型．本例取抛物型方程()udcuauft，故参数,,,abcd分别是1,0,10,1.第五步：选择Mesh菜单中InitializeMesh命令，进行网格剖分．选择Mesh菜单中RefineMesh命令，使网格密集化，如图22.3．图22.3网格密集化第六步：解偏微分方程并显示图形解选择Solve菜单中SolvePDE命令，解偏微分方程并显示图形解，如图22.4所示。图22.4偏微分方程的图解图第七步：单击Plot菜单中Parameter选项，打开PlotSelection对话框，选中Color,Height(3-Dplot)和Showmesh三项．再单击Polt按钮，显示三维图形解，如图22.5所示．图24.15图22.5偏微分方程的三维图形解第八步：若要画等值线图和矢量场图，单击Plot菜单中Parameter选项，在Plotselection对话框中选中Contour和Arrows两项．然后单击Plot按钮，可显示解的等值线图和矢量场图，如图22.6所示。图22.6解的等值线图和矢量场图22.2计算机仿真编程求解偏微分方程求解偏微分方程除了上节介绍的直接使用偏微分方程工具箱外，还可以用编写程序（对于MATLAB仿真，可以编写M文件）的方法求解偏微分方程。22.2.1双曲型：波动方程的求解1.求解双曲型方程下面将讨论标准波动方程的求解问题．波动方程属于双曲型方程，即fauuctud)(22（a、c、d、f是参数）．求解双曲型方程调用形式如下：(1)u1=hyperbolic（u0，ut0，tlist，b，p，e，t，c，a，f，d）其中：参数c、a、f、d决定了方程的类型．b代表求解域的边界条件．b可以是边界条件矩阵，也可以是相应的PDE边界条件M文件名．网格坐标描述矩阵p,e,t是由网格初始化命令得到的：(2)[p,e,t]=initmesh(g)其中：g代表求解区域几何形状，是相应的PDE几何分类函数M文件名．initmesh函数的作用是将求解区域进行三角形网格化，网格大小由区域的几何形状决定．输出的p、e、t都是网格数据．点阵p的第1、2行分别包含了网格中点的x、y坐标．e是边缘矩阵，其各行的意义与求解步骤无直接关系，从略．t是三角矩阵，其中的几行描述了区域的顶点．有时还会用到修整网格（精细化）命令(3)[p,e,t]=refinemesh(g,p,e,t)即是迭代过程，得到更细小的网格，使结果更精确．其中：u0、ut0（ut即是tu/）是初始条件．tlist是t=0时刻以后均匀的时间矩阵．hyperbolic函数返回的是u在tlist中各个时间点、在区域各三角形网格处的值u1．u1的每行是由p中相应列上的坐标值所得的函数值，u1的每一列是矩阵tlist中相应的时间项所对应的函数值．2.动画图形显示为了将所得的解形象地表示出来，还要通过一些动画图形命令．为了加速绘图，首先把三角形网格转化成矩形网格．调用形式如下：(1)uxy=tri2grid（p，t，u1，x，y）p、t是描述三角形网格的矩阵，x、y是求解区域中矩形网格的坐标点（矩阵x、y必须都是递增顺序），u1是各时刻三角形网格中的解．输出矩阵uxy是用线性插值法在矩形网格点上得出的相应u值．(2)[uxy,tn,a2,a3]=tri2grid(p,t,u,x,y)uxy、p、t、u、x、y意义同上，tn是格点的指针矩阵，a2、a3是内插法的系数．(3)uxy=tri2grid(p,t,u,tn,a2,a3)用此命令之前，应先用一个tri2grid命令得出矩阵tn、a2、a3．用此方法可以加快速度．注意:如果矩形网格点在三角形网格之外，则结果中将会出现出错信息‘NAN’．主要的绘图（包括动画）命令函数有：moviein、movie、pedplot、pdesurf等，其具体应用见下面的例题．例24.2.1用MATLAB求解下面波动方程定解问题并动态显示解的分布2222221111()0||0,0π(,,0)atan[sin()],(,,0)2cos(π2xxyytuuutxyuuuuyyuxyxuxyπ)exp[cos()]2xy已知求解域是方形区域，空间坐标的个数由具体问题确定．【解】采用步骤如下（注：在MATLAB环境下运行时，下面的中文不输入。凡是“%”后的语句为解释语句,MATLAB不执行）（1）题目定义g='squareg';%定义单位方形区域b='squareb3';%左右零边界条件，顶底零导数边界条件c=1;a=0;f=0;d=1;（2）初始的粗糙网格化[p,e,t]=initmesh('squareg');（3）初始条件x=p(1,:)';%注意坐标向量都是列向量y=p(2,:)';u0=atan(sin(pi/2*x));ut0=2*cos(pi*x).*exp(cos(pi/2*y));（4）在时间段0～5内的31个点上求解n=31;tlist=linspace(0,5,n);%在0～5之间产生n个均匀的时间点（5）求解此双曲问题u1=hyperbolic(u0,ut0,tlist,b,p,e,t,c,a,f,d);%得到如下结果：%Time:0.166667%Time:0.333333%Time:0.5%…%Time:4.5%Time:4.66667%Time:4.83333%Time:5%428successfulsteps%62failedattempts%982functionevaluations%1partialderivatives%142LUdecompositions%981solutionsoflinearsystems%现在把解u1用动态图形表示出来．（6）矩形网格插值delta=-1:0.1:1;[uxy,tn,a2,a3]=tri2grid(p,t,u1(:,1),delta,delta);gp=[tn;a2;a3];（7）在0～5时间内动画显示newplot;%建立新的坐标系newplot;M=moviein(n);umax=max(max(u1));umin=min(min(u1));fori=1:n,...%注意‘…’符号不可省略ifrem(i,10)==0,...%当n是10的整数倍时，在命令窗口打印出相应的数字fprintf('%d',i);...图22.7某一瞬时的波动方程的解图end,...pdeplot(p,e,t,'xydata',u1(:,i),'zdata',u1(:,i),'zstyle','continuous','mesh','on','xygrid','on','gridparam',gp,'colorbar','off');...axis([-1,1,-1,1uminumax]);caxis([uminumax]);...M(:,i)=getframe;...ifi==n,...fprintf('done\n');...end,...end%运行结果是：102030done%若要显示持续不断的动画，则再加上下面语句：nfps=5;movie(M,10,nfps);动态解图可以直接通过MATLAB仿真程序执行看出，图22.7是动态图的某一瞬间的解的分布。图22.7某一瞬时的波动方程的解图例22.1.2求解下列热传导定解问题．222211()0(,,)||01(0.4)(,,0)0(0.4)xyxyuuutxyuxyturuxyr求解域是方形区域，其中空间坐标的个数由具体问题确定．图22.8某一瞬时的热传导方程解分布22.1.3椭圆型：稳定场方程的求解稳定场方程属于椭圆方程，下面我们来求解含有源项的标准稳定场方程，也即是泊松方程．在MATLAB中是指如下形式：()cuauf求解椭圆方程用命令：u=assempde(b,p,e,t,c,a,f)各输入量的意义同前．由于是稳定场方程，则输出的矩阵u是坐标矩阵p相应的解．例22.1.3用MATLAB求解下列泊松方程，并将计算机仿真解与精确解比较.11|0ruu22.9(b)精确解与仿真解的误差解图图22.9(a)泊松方程的动态解图22.1.4点源泊松方程的适应解本节介绍函数的适应性网格解法．例22.1.4用MATLAB求解含点源的泊松方程，并与精确解比较．方程如下：0),(12ruyxu图22.10(a)点源泊松方程的解图图22.10(b)仿真解与精确解的误差图22.3定解问题的计算机仿真本小节主要对定解问题的求解给出计算机仿真显示，直观明了地给出解动态（或静态）分布