偏微分方程数值解.doc

《偏微分方程数值解.doc》由会员分享，可在线阅读，更多相关《偏微分方程数值解.doc（31页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第一章 概 述大家采用下面的方法求解Terzaghi一维固结方程。11 偏微分方程工具箱的功能 偏微分方程工具箱(PDE Toolbox)提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境。PDE Toolbox的功能包括： (1) 设置PDE (偏微分方程)定解问题，即设置二维定解区域、边界条件以及方程的形式和系数； (2) 用有限元法 (FEM) 求解PDE数值解； (3) 解的可视化。 无论是高级研究人员还是初学者，在使用PDE Too1box时都会感到非常方便。只要PDE定解问题的提法正确，那么，启动MATLAB后，在MATLAB工作空间的命令行中键人pdetool，系

2、统立即产生偏微分方程工具箱(PDE Toolbox)的图形用户界面(Graphical User Interface，简记为GUI)，即PDE解的图形环境，这时就可以在它上面画出定解区域、设置方程和边界条件、作网格剖分、求解、作图等工作，详见14节中的例子。我们将在第二章详细介绍GUI的使用，在第二章给出大量典型例子和应用实例。除了用GUI求解PDE外，也可以用M文件的编程计算更为复杂的问题，详见第三章和第四章的内容。 12 PDE Toolbox求解的问题及其背景 121 方程类型 PDE Toolbox求解的基本方程有椭圆型方程、抛物型方程、双曲型方程、特征值方程、椭圆型方程组以及非线性椭

3、圆型方程。椭圆型方程： ，椭圆型方程：其中是平面有界区域，c,a,f以及未知数u是定义在上的实（或复）函数。抛物型方程：双曲型方程：.特征值方程：其中d是定义在上的复函数，是待求特征值。在抛物型方程和双曲型方程中，系数c,a,f和d可以依赖于时间t。可以求解非线性椭圆型方程： 其中c,a,f可以是未知函数u的函数。还可以求解如下PDE方程组； 利用命令行可以求解高阶方程组。对于椭圆型方程，可以用自适应网格算法，还能与非线性解结合起来使用。另外，对于Poission方程还有一个矩形网格的快速求解器。122 边界条件（1）Dirichlet条件 ： ( 2 ) Neumann 条件： 其中是的边界

4、上的单位外法向量，和是定义在上的函数。对于特征值问题仅限于齐次条件：和。对于非线性情形系数和可以依赖于u；对于抛物型方程和双曲型方程，系数可以依赖于时间t。对于方程组情形，边界条件为( 1 ) Dirichlet 条件： ( 2 ) Neumann 条件： ( 3 ) 混合边界条件为: 其中的计算要使得Dirichlet条件满足。在有限元法中，Dirichlet条件也称为本质边界条件，Neumann条件称为自然边界条件。1.3 如何使用FDE Toolbox 1.3.1 定解问题的设置 员简单的办法是在PDE Tool上直接使用图形用户界面(GUl)。设置定解问题包括三个步骤： (1)Draw

5、模式：使用CSG(几何结构实体模型)对话框画几何区域，包括矩形、圆、椭圆和多边形，也可以将它们组合使用。 (2)Boundary模式：在各个边界段上给出边界条件， (3)PDE模式：确定方程的类型、系数c，a，f和d c。也能够在不同子区域上设置不同的系数(反映材料的性质)。1.3.2 解PDE问题用GUI解PDE问题主要经过下面两个过程（模式）(1)Mesh模式；生成网格自动控制网格参数。(2)Solve模式：对于椭圆型方程还能求非线性和自适应解。对于抛物型和双曲型力程设置初始边值条件后能求出给定t时刻的解。对于特征值问题，能求出给定区间内的特征值；求解后可以加密网格再求解。1.3.3 使用

6、Toolbox求解非标准的问题对于非标准的问题。可以用PDE Too1box的函数。或者用FEM(有限元法)求解更为复杂的问题。1.3.4 计算结果的可视化从GUI能够使用Plot模式实现可视化。可以使用Color, Height和Vector等作图。对于抛物型和双曲型方程，还可以生成解的动画。这些操作通过命令行都很容易实现。1.3.5 应用领域在应用界面提供了丁如下应用领域结构力学平面应力问题结构力学平面应变问题静电场问题静磁场问题交流电磁场问题直流导体介质问题热传导问题9散问题这些界面都有对话框，它包括PDE的系数、边界条件、解的性质等。1. 4 解偏微分方程的一个例子解Poisson方程

7、，边界条件为齐次Dirichlet类型。第一步：启动MATLABl, 键入pdetool，按回车键确定便可启动GUI，然后在Options菜单下选择Grid命令，打开栅格, 栅格的使用，能使用户容易确定所绘图形的大小，如图11 1-1第二步：分步完成平面几何造型：R1-C1-E1+R2+C2。用菜单或快捷工具，分别画矩形R1、矩形R2、椭圆E1、圆C1、圆C2。画圆时，首先选中椭圆工具，按鼠标右键并拖动即可、或者在按ctrI的同时，拖动鼠标也可绘制圆。然后在Set formula栏，进行编辑并用算术运算将将图形对象名称连接起来，删除默认的表达式键入R1-C1-E1+R2+C2，按等号健得到所需

8、图形。若需要，还可进行储存形成M文件。选择Boundary菜单中Boundary Mode命令，进入边界模式。单击Boundary菜单中Remove A11 Subdomain Borders选项，去除子域边界。如果想将几何信息和边界信息进行存储，应选择Boundary菜单中的ExPort Decomposed GeometryBoundary Conds命令，将它们分别储存于g,b变量中, 通过MATLAB形成M文件。第三步：选取边界单击Boundary菜单中Specify Bounddy Conditions选项，打开Boundary conditlons对话框，输入边界条件，如图14。本

9、例取缺省条件。即将全部边界设为齐次Dirichlet条件，边界颜色显示为红色。第四步：选择PDE菜单中PDE Mode命令，进入PDE模式。单击PDE菜单中PDE Specification选项，打开PDE对话框，设置方程类型。本例取缺省设置，类型为椭圆型，参数c,a,f分别为1，0，10。第五步：选择Mesh菜单中Initialize Mesh命令，进行网格剖分。第六步：选择Mesh菜单中Refine Mesh命令，对网格加密。第七步：选择Solve菜单中So1ve PDE命令，解偏微分方程并显示图形解。第八步：单击Plot菜单中Parameters选项，打开Plot selection对话

10、框，选中Color, Height (3D Plot)和Show mesh三项。然后单击Plot按钮，显示三维图形解。第九步：如果要画等值线图和矢量场图，单击Plot菜单中Parameters选项，打开Plot Selection对话框选中Contour和Arrows两项。然后单击P1ot按钮，可显示解的等值线图和矢量场图。第二章 PDE图形用户界面2.1 PDE Toolbox菜单 File菜单（如图1-1） 图1-1New 新建一个几何结构实体模型（Constructive Solid Geomery,简记为CSG），默认文件名为“Untitled”。Open 从硬盘装载M文件Save 将

11、在GUI内完成的成果储存到一个M文件中。Save As 将在GUI内完成的成果储存到另外一个M文件中。Print 将PDE工具箱完成的图形送到打印机内进行硬拷贝。Exit 退出PDE工具图形用户界面。2 Edit菜单（如图1-2） 图1-2Undo 在绘制多边形时退回到上一步操作。Cut 将已选实体剪切到剪贴板上。Copy 将已选实体拷贝到剪贴板上。Paste 将剪贴板上的实体粘贴到当前几何结构实体模型中。Clear 删除已选的实体。Select All 选择当前几何结构实体造型CSG中的所有实体及其边界和字域。3 Options菜单（如图1-3） 图1-3Grid 绘图时打开或关闭栅格。Gr

12、id Spacing 调整栅格的大小。Snap 打开或关闭捕捉栅格功能。Axes Limits 设置绘图轴的坐标范围。Axes Equal 打开或关闭绘图方轴。Turn off Toolbar Help 关闭工具栏按钮的帮助信息。Zoom 打开或关闭图形缩放功能。Application 选择应用的模式。Refresh 重新显示PDE工具箱中的图形实体。4 Draw菜单（如图1-4） 图1-4Draw Mode 进入绘图模式。Rectangle/square 以角点方式画矩形/方行（Ctrl+鼠标）。Rectangle/square（centered） 以中心方式画矩形/方行（Ctrl+鼠标）。

13、Ellipse/circle 以矩阵角点方式画椭圆/圆（Ctrl+鼠标）。Ellipse/circle（centered） 以中心方式画椭圆/圆（Ctrl+鼠标）。Polygon 画多边形，单击鼠标右键可封闭多边形。Rotate 旋转已选的图形。Export Geometry Description,Set Formula，Labels 将几何描述矩阵gd、公式设置字符sf和标识空间矩阵ns输出到主工作空间去。 单击Draw菜单中Rotate选项，可打开Rotate比对活框，通过输入旋转的角度，可使选择的物体按输入的角度逆时针旋转。旋转中心的选择如果缺省，则为图形的质心，也可以输入旋转中心坐标

14、。5 Boundary菜单（如图1-5） 图1-5Boundary Mode 进入边界模式。Specify Boundary Conditions 对于已选的边界输入条件，如果没有选择边界，则边界条件适用于所有的边界。Show Edge Labels 显示边界区域标识开关，其数据是分解几何矩阵的列数。Show Subdomain Labels 显示子区域标识开关，其数据是分解几何矩阵中的子域数值。Remove Subdomain Border 当图形进行布尔运算时，删除已选取的子域边界。Remove All Subdomain Borders 当图形进行布尔运算时，删除所有的子域边界。Expo

15、rt Decomposed Geometry，Boundary Conds 将分解几何矩阵g、边界条件矩阵b输出到主工作空间。 选择Boundary菜单中Specify Boundary Conditions命令可定义边界条件。在打开的Boundary condition对话框，可对已选的边界输入边界条件。共有如下三种不同的条件类型：NeMmann条件 这里边界条件是由方程系数q和g确定的，在方程组的情况下(换成方程组模式)，q是22矩阵，g是2x1矢量。Dirichlet条件 u定义在边界上，边界条件方程是价h*u=r，这里h是可以选样的权因子(通常为1)。在方程组情况下，h是2x2矩阵，r

16、是2x l矢量，混合边界条件(仅适合于方程组情形) 它是Dirichlet和Neumann的混合边界条件，q是2x 2矩阵，g是2x1矢量，h是1x 2矢量，r是一个标量。6 PDE菜单（如图1-6） 图1-6PDE Mode 进入偏微分方程模式。Show Subdomain Labels 显示子区域标识开关。PDE Specification 调整PDE参数和类型。Export PDE Coefficients 将当前PDE参数c,a,f,d输出到主工作空间，其参数变量为字符类型。7 Mesh菜单（如图1-7） 图1-7Mesh Mode 输入网格模式。Initialize Mesh 建立和

17、显示初始化三角形网格。Refine Mesh 加密当前三角型网格。Jiggle Mesh 优化网格。Undo Mesh Change 退回上一次网格操作。Display Triangle Quality 用01之间数字化的颜色显示三角形网格的质量，大于0.6的网格可接受的。Show Node Labels 显示网格节点标识开关，节点标识数据是点矩阵p的列。Show Triangle Labels 显示三角形网格标识开关，三角形网格标识数据是三角形矩阵t的列。Parameters 修改网格生成参数。Export Mesh 输出节点矩阵p、边界矩阵e和三角形矩阵t到主工作空间。8 Solve菜单（

18、如图1-8） 图1-8Solve PDE 对当前的几何结构实体CSG、三角形网格和图形解偏微分方程。Parameters 调整PDE的参数。Export Solution 输出PDE的解矢量u。如果可行，将计算的特征值1输出到主工作空间。1.1.9 Plot菜单（如图1-9） 图1-9Plot Solution 显示图形解。Parameters 打开绘图方式对话框。Export Movie 如果动画被录制了，则动画矩阵M将输出到主工作空间。10 Window菜单从Window菜单项，可选择当前打开的所有的MATLAB图形窗口，被选择的窗台移至前台。11 Help菜单Help 显示帮助信息Abo

19、ut 显示版本信息12 PDE工具栏 主菜单下是工具栏，工具栏中喊有许多工具图标按钮，可提供快速、便捷的操作方式。从左到右5个按钮为绘图模式按钮，紧接着的6个为边界、网格、解方程和图形显示控制功能按钮，最右边的为图形缩放功能键。（如图1-10） 图1-10 以角点方式画矩形/方行（Ctrl+鼠标）。 以中心方式画矩形/方行（Ctrl+鼠标）。 以矩形角点长轴方式画椭圆/圆（Ctrl+鼠标）。 以中心方式画椭圆/圆（Ctrl+鼠标）。 画多边形，按右键可封闭多边形。 进入边界模式。 打开PDE Specification（偏微分方程类型）对话框。 初始化三角形网格。 加密三角形网格。 解偏微分方

20、程。 打开Plot Selection对话框，确定后给出解的三维图形。 为显示缩放切换按钮。第三章 典型方程及其应用实例求解PDE问题主要有两种方法，一种是使用图形用户界面，另一种是采用命令行编程。前者直观简便，而后者更为灵活。21 求解椭圆方程的例子例：单位圆上的Poisson方程边值问题： 这一问题的精确解为：若使用图形用户界面（Graphical User Interface,简记为GUI），则首先在MATLAB的工作窗口中键入pdetool，按回车键确定，于是出现PDE Toolbox窗口。如果需要坐标网格，单击Options菜单下的Grid选项即可。下面分步进行操作。（i）画区域圆

21、单击工具，大致在（0，0）位置单击鼠标右键同时拖拉鼠标到适当位置松开，绘制圆。为了保证所绘制的圆是标准的单位圆，在所绘圆上双击，打开Object Dialog对话框，精确地输入圆心坐标X-center为0、Y-cebter为0及半径Radius为1，然后单击OK按钮，这样单位远已画好。（ii）设置边界条件 单击工具，图形边界变红，逐段双击边界，打开Boundary Condition对话框，输入边界条件。对于同一类型的边界，可以按Shift键，将多个边界同时选择，统一设置边界条件。本题选择Dirichlet条件，输入h为1，r为0，然后单击OK按钮。也可以单击Boundary菜单中Specif

22、y Boundary Conditions选项，打开Boundary Condition对话框输入边界条件，如图2-1。（iii）设置方程 单击PDE菜单中PDE Specification选项，打开PDE Specification对话框，选项方程类型。本题单击Elliptic，输入c为1，a为0，f为1，然后单击OK按钮，如图2-2。 图2-1图2-2（iv）网格剖分 单击工具，或者单击Mesh菜单中Initialize Mesh选项，可进行初始网格剖分，这时在PDE Toolbox窗口下方的状态栏内显示初始问网格的节点数和三角形单元数。本题节点数为144个，三角形单元数为254个。如果需

23、要网格加密，再单击，或者单击Mesh菜单中Refine Mesh选项，这时节点数变为541个，三角形单元数为1016个，如此还可继续加密。（v）解方程 单击工具，或者单击Solve菜单中Solve菜单中Solve PDE选项，可显示方程色彩解。如果单击Plot菜单中Parameters选项，出现Plot Selection对话框，如图2-3，从中可以选择Color，Contour，Arrows,Deformed mesh，Height（3-D polt），还可以设置等值线的数目等。本例中选择Color，Contour，Height（3-D polt）和Show mesh四项，然后单击Plot按

24、钮，方程的图形解如图2-4所示。除了作定解问题解u的图形外，也可以作|grad u|，|cgrad u|等图形。图2-3图2-4（vi）与精确解作比较 单击Plot菜单中Parameters选项，打开Plot Selection对话框，在Height（3-D plot）行Property下拉框中选user entry，且在该行的User entry输入框中键入u-（1-x.2-y.2）/4,单击lot按钮就可以看到解的绝对误差图形，如图-可见在边界处误差为。图2-5（vii）输出网格节点的编号、单元编号以及节点坐标 单击Mesh菜单中Show Node Labels选项，再单击工具或，即可显示

25、节点编号。若要输出节点坐标，只需单击Mesh菜单中Export Mesh选项，这时打开的Export对话框中默认值为p，e，t，这里p，e，t分别表示points（点）、edges（边）、triangles（三角形）数据的变量，单击OK按钮。然后在MATLAB命令窗口键入p，按回车键确定，即可显示出节点按编号排列的坐标（二维数组）；键入e，按回车键，则显示边界线段数据矩阵（7维数组）；输入t，按回车键，则显示三角形单元数据矩阵（4维数组）。（viii）输出解的数值 单击Solve菜单中Export Solution选项，在打开的Export对话框中输入u，单击OK按钮确定。再在MATLAB命令

26、窗口中输入u，按回车确定，即显示按节点编号排列的解的数值。我们也可以用MATLAB程序求解PDE问题，同时显示解的图形：p,e,t=initmesh(circleg,hmax,1);Error=;err=1;While err0.001,p,e,t=refinemesh(circleg,p,e,t);U=assempde(circleb1,p,e,t,1,0,1);Exact=(1-p(1,:).2-p(2,:).2)/4;Err=norm(u-exact,inf);Error=error err;EndPdemesh(p,e,t)Pdesurf(p,t,u)Pdesurf(p,t,u-exa

27、ct)通过命令行键入help+命令函数，如help pdemesh，按回车键，可以调入有关命令函数的定义、参数格式等帮助信息。22 求解抛物型方程的例子例：考虑一个带有矩形孔的金属板上的热传导问题。板的左边保持在100，板的右边热量从板向环境空气定常流动，其他边及内孔边界保持绝缘。初始时板的温度为0，于是概括为如下定解问题：域的外边界顶点坐标为（-0.5，-0.8），（0.5，-0.8），（0.5，0.8），（-0.5，0.8）。内边界顶点坐标为（-0.005，-0.4），（0.05，-0.4），（0.05，0.4），（-0.05，0.4）。 使用GUI求解这一问题。在PDE Toolbox窗

28、口的工具栏中选择Generic Scalar模式。 （i）区域设置 单击工具，在窗口拖拉出一个矩形，双击矩形区域，在Object Dialog对话框中输入Left为-0.5，Bottom为-0.8，Width为1，Height为1.6，单击OK按钮，显示矩形区域R1。用同样方法作内孔R2，只要设置Left=-0.05，Bottom=-0.4，Width=0.1，Height=0.8即可。然后在Set formula栏中键入R1-R2。（ii）设置边界条件 单击，使边界变红色，然后分别双击每段边界，打开Boundary Cvondition对话框，设置边界条件。在左边界条件。在左边界上，选择Di

29、richlet条件，输入h为1，r为100；右边界上，选择Neumann条件，输入g为-1，q为0；其他边界上，选择Neumann条件，输入g为0，q为0。（iii）设置方程类型 单击，打开PDE Specification对话框，设置方程类型为Parabolic（抛物型），d=1，c=1，a=0，f=0，单击OK按钮。（iv）网格剖分 单击，或者加密网格，单击。（v）初值和误差的设置 单击Solve菜单中Parameters选项，打开Solve Parameters对话框，输入Time为0:5，u（t0）为0，Relative tolerance为0.01，Absolute toleranc

30、e为0.001，然后单击OK按钮。（vi）数值解的输出 单击Solve菜单中Export Solution选项，在Export对话框中输入u，单击OK按钮。再在MATLAB命令窗口中键入u，按回车键，这时显示出按节点编号的数值解。（vii）解的图形 单击Plot菜单中Parameters选项，打开Plot Selection对话框，选Color，Contour，Arrows，单击Plot按钮，窗口显示出Time=5时解的彩色图形，如图2-6。图2-6MATLAB程序： p,e,t=initmesh(crackg); u=parabolic(0,0:0.5:5,crackb,p,e,t,1,0,

31、0,1); pdeplot(p,e,t,xydata,u(:,11),mesh,off,colormap,hot)23 求解双曲型方程的例子例：考虑如下二维波动方程的定解问题：用GUI求解。类似前面的例子，首先作正方形区域：设置断点坐标为（-1，1），（-1，-1），（1，-1），（1，1）。在Object Dialog对话框中输入Left为-1，Bottom为-1，Width为2，Height2，单击OK按钮。设置边界条件：左、右边界用Dirichlet条件，输入h为1，r为0；上、下边界用Neumann条件，输入q为0，g为0，作网格剖分。设置方程类型为Hyperbolic（双曲型），键入

32、c=1，a=0，f=0，d=1。单击Solve菜单中Parameters选项，打开Solve Parameters对话框，在Time栏中键入linspace（0，5，31），设置u的初始值u（t0）为atan(cos(pi/2*x),u的初始值u(t0)为3*sin(pi*x).*exp(sin(pi/2*y),Relative tolerance为0.01，Absolute tolerance为0.001。最后输出图形解：单击Plot菜单中Parameters选项，打开Plot Selection对话框，选Color，Contour，Arrows，Height（3-D Plot）和Show

33、mesh五项，单击Plot按钮，三维彩色图形的解如图2-7所示。图2-7MATLAB的动画程序： p,e,t=initmesh(squareg); x=p(1,:);y=p(2,:); u0=atan(cos(pi/2*x); ut0=3*sin(pi*x).*exp(sin(pi/2*y); n=3l; tlist=linspace(0,5,n); uu=hyperbolic(u0,ut0,tlist,squareb3,p,e,t,1,0,0,1); delta=-1:0.1:1; uxy,tn,a2,a3=tri2grid(p,t,uu(:,1),delta,delta); gp=tn;a

34、2;a3; umax=max(max(uu); umin=min(min(uu); newplot M=moviein(n); For I=1:n,pdeplot(p,e,t,xydata,uu(:,I),zdata,uu(:,I), mesh,off,xygrid, on, gridparam,gp,colorbar,off,zstyle, continuous);axis(-1 1 1 1 umin umax);caxis(umin umax);M(:,I)=getframe;EndMovie(M,10);24 求解特征值问题的例子例：混合边界条件的特征值问题： 首先作区域：方形的4个顶点

35、为（-1，1），（-1，-1），（1，-1），（1，1）。在Object Dialog对话框中输入Left为-1，Bottom为-1，Width为2，Height2。单击，逐段双击边界设置边界条件：左边界选Dirichlet条件，输入h为1，r为0；右边界选Neumann条件，输入q为-3/4，g为0；上、下边界选Neumann条件，输入q为0，g为0。在PDE Specification对话框中设置方程类型为Eigenmodes（特征值模式），PDE的系数设置为c=1，a=0，d=1。由于右边界上Neumann条件，特征值可以出现负值，因此在求小于10的特征值时应在Solve Paramet

36、ers对话框中键入-inf 10。作网格剖分，然后单击，得到第一个特征值对应的特征函数图形，如图2-8。 图2-8再在Plot Selection对话框的Eigenvalue项中选取不同的特征值，比如第二个特征值“2.06，并选Color，Contour，Height（3-D Plot）和Show mesh四项，如图2-9，可作出第二个特征值对应的特征函数的三维图形，如图2-10。图2-9图2-10MATLAB程序：p,e,t=initmesh(squareg);p,e,t=refinemesh（squareg,p,e,t）;v,l=pdeeig(squareb2,p,e,t,1,0,1,-i

37、nf 10);pdesurf(p,t,v(:,4)25 偏微分方程在一维空间的应用已知偏微分方程初始条件时间t和一维空间变量x，MATLAB自身存在对偏微分方程的解函数pdepe。单一方程假设我们要求解热传导方程 ， （1.1）MATLAB指定抛物型的偏微分方程形式如下：边界条件为：其中为边界左端点，为边界右端点，初始条件为 。（我们注意到同一函数b在方程和边界条件中出现。）通常情况下，每个函数都会被指定到不同的MATLAB文档中。也就是说，与方程有关的函数c,b和s都会被指定保存到一个MATLAB文件中，与边界条件有关的函数p和q保存到令一个文档中（此外，我们需要注意到函数b是相同的只需保存

38、一次即可），最后将初始函数保存在第三个文档中。执行命令函数pdepe将把m个文档结合起来进行运算，返回方程的解。在我们的例子中，有： 我们将它们保存到MATLAB文档中，记作eqn1.m (m=0以后在加详述)。Functionc,b,s=eqn1(x,t,u,DuDx)%EQN1:MATLAB function M-file that specifies%a PDE in time and one space dimension.C=1;B=DuDx;S=0;对于边界条件，我们有： P(0,t,u) = u ； q(0,t) = 0 P(1,t,u) = u 1 ； q(1,t) = 0将它

39、们保存到MATLAB文档中，记作bc1.m.Function p1,q1,pr,qr=bc1(x1,u1,xr,ur,t)%BC1:MATLAB function M-file that specifies boundary conditions%for a PDE in time and one space dimension.P1=u1;Q1=0;Pr=ur-1;Qr=0;对初始条件，有： 将它保存到MATLAB文档中，记作initial1.m.Function value=initial1(x)%INITIAL1:MATLAB function M-file that specifies

40、 the initial condition%for a PDE in time and one space dimension.Value=2*x/(1+x2);最终我们准备好了用pdepe函数处理偏微分方程。在如下的MATLAB文档中，我们选择坐标值x和t，求解偏微分方程并绘制解的表面图（如图2-11）%PDE1：MATLAB script M-file that solves and plots%solutions to the PDE stored in eqn1.mM=0;%NOTE:m=0 specifies no symmetry in the problem.Taking%m=

41、1 specifies cylindrical symmetry,while m=2 specifies%spherical symmetry.%Define the solution meshx=linspace(0,1,20);t=linspace(0,2,10);%Solve the PDEU=pdepe(m,epn1,initial1,bc1,x,t);%Plot solutionSurf(x,t,u);Title(Surface plot of solution);Xlabel(Distance x);Ylabel(Time t); 图2-11通常，我们会发现绘制解的剖面图是十分有用

42、的。此时，t是固定的，u是x的函数。解u（t,x）作为t和x的向量指标被保存为指标矩阵。例如，u(1,5)返回在点（t(1),x(5)）出u的值。我们可以使用命令plot(x,u(1,:)绘制u的最初值的图像（t=0）（见图2-12）。其实，我们有一种比较快的方法去绘制剖面图，如下的MATLAB序列： Fig=plot(x,u(1,:),erase,xor) For k=2:length(t) Set(fig,xdata,x,ydata,u(k,:) Pause(.5) end试过之后，你会发现它们是如何快速地接近热传导方程的平衡位置的。（平衡位置就是及时停止并要发生改变的位置）。 图2-123 PDE Toolboxz中的命令简介本节将介绍PDE Toolbox中的命令函数。利用这些命令函数，可通过命令行而不是用PDE Tool图形用户界面来解偏微分方程。31 PDE Toolbox中的函数及其分类表3-1 PDE数值计算函数函 数 目 的 Adaptmesh 生成自适应网格和解PDE问题 Assema 组装积分区域贡献 Assemb 组装边界条件贡献 assempde 组装刚度矩阵和方程右边的函数 hyperbolic 解双曲型PDE问题 parabolic 解抛物型PDE问题