偏微分方程数值解

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

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

3、型方程。椭圆型方程： -'(c、u) a u =,fin椭圆型方程:A © u) au二f,inj其中是平面有界区域，c, a, f以及未知数u是定义在门上的实(或复)函数。抛物型方程：d、u(c，u) au = f, in双曲型方程：-；u) au = f, in门特征值方程：-1 (ci u) au V.du, in j其中d是定义在门上的复函数，是待求特征值。在抛物型方程和双曲型方程中，系数c, a, f和d可以依赖于时间to可以求解非线性椭圆型方程：一7 (c(u)u )+a(u) = f (u), in Q,其中c, a,f可以是未知函数u的函数。还可以求解如下PD

4、E方程组；可 «Cii(u)Vui )可 fci2(u)Vu2 )+aiiui +ai2u2 = fi, V(C2i(u)$ui )可 W22(u)J2 )+a2M +玄22出=匚利用命令行可以求解高阶方程组。对于椭圆型方程，可以用自适应网格算法，还能与非线性解结合起来使用。另外，对于Poission方程还有一个矩形网格的快速求解器。1. 2. 2边界条件(1) Dirichlet 条件:h u= r(2 ) Neumann 条件： n © u) qu =g其中n是门的边界上的单位外法向量，g,q,h和r是定义在"上的函数。对于特征值问题仅限于齐次条件：g =0,

5、和r = 0。对于非线性情形.系数g,q,h和r可以依赖于U；对于抛物型方程和双曲型方 程，系数可以依赖于时间t。对于方程组情形，边界条件为(1 ) Dirichlet 条件：hi 1u < h 1 u 2 r 加*h22u2 = r2(2 ) Neumann 条件：n (qi uj n (s' U2) qu q2 二 g1 n (八 U1) n (c?、出)72鸡 722U2 二 g?(3 )混合边界条件为：hnU1 人2氏=r1n (GN U1) n g 上)w q12U gj IVn(C21I U1) n (c U2) 721U1 722U2 =g2 h12J其中的计算要使

6、得Dirichlet条件满足。在有限元法中，Dirichlet条 件也称为本质边界条件，Neuma nn条件称为自然边界条件。1.3 如何使用 FDE Toolbox1.3.1定解问题的设置员简单的办法是在PDE Tool上直接使用图形用户界面(GUI)。设 置定解问题包括三个步骤：(1) Draw模式：使用CSG(几何结构实体模型)对话框画几何区域， 包括矩形、圆、椭圆和多边形，也可以将它们组合使用。(2) Boundary模式：在各个边界段上给出边界条件，(3) PDE模式：确定方程的类型、系数 c, a, f和d c。也能够在不同子区域上设置不同的系数(反映材料的性质)1.3.2 解 P

7、DE 问题用 GUI 解 PDE 问题主要经过下面两个过程(模式)(1)Mesh 模式；生成网格自动控制网格参数。(2)Solve 模式：对于椭圆型方程还能求非线性和自适应解。对于 抛物型和双曲型力程设置初始边值条件后能求出给定 t 时刻的解。 对于特征值问题， 能求出给定区间内的特征值； 求解后可以加密网格 再求解。1.3.3 使用 Toolbox 求解非标准的问题对于非标准的问题。可以用PDE Toolbox的函数。或者用FEM(有 限元法 )求解更为复杂的问题。1.3.4 计算结果的可视化从 GUI 能够使用 Plot 模式实现可视化。可以使用 Color, Height 和 Vecto

8、r 等作图。对于抛物型和双曲型方程，还可以生成解的动画。 这些操作通过命令行都很容易实现。1.3.5 应用领域在应用界面提供了丁如下应用领域结构力学平面应力问题结构力学平面应变问题静电场问题静磁场问题交流电磁场问题直流导体介质问题.热传导问题.9 散问题这些界面都有对话框，它包括PDE的系数、边界条件、解的性质等。1. 4解偏微分方程的一个例子解Poisson方程二f，边界条件为齐次 Dirichlet类型。第一步：启动 MATLABI,键入pdetool,按回车键确定便可启动 GUI，然后在Options菜单下选择Grid命令，打开栅格,栅格的使用, 能使用户容易确定所绘图形的大小，如图 1

9、 11-1第二步：分步完成平面几何造型： R1-C1-E1+R2+C2。用菜单或 快捷工具，分别画矩形 R1、矩形R2、椭圆E1、圆C1、圆C2。画圆 时，首先选中椭圆工具，按鼠标右键并拖动即可、或者在按 ctrI的同 时，拖动鼠标也可绘制圆。然后在 Set formula栏，进行编辑并用算 术运算将将图形对象名称连接起来，删除默认的表达式键入 R1-C1-E1+R2+C2,按等号健得到所需图形。若需要，还可进行储存.形成 M 文件。选择Boundary菜单中Boundary Mode命令，进入边界模式。单 击 Boundary 菜单中 Remove A11 Subdomain Borders

10、选项，去除子域 边界。如果想将几何信息和边界信息进行存储，应选择 Boundary 菜 单中的 ExPort Decomposed Geometry Bou ndary Co nd'命令，将它 们分别储存于 g,b 变量中 , 通过 MATLAB 形成 M 文件。第三步：选取边界单击 Boundary 菜单中 Specify Bounddy Conditions选项，打开Boundary conditlons对话框，输入边界条件， 如图 14。本例取缺省条件。即将全部边界设为齐次 Dirichlet 条件， 边界颜色显示为红色。第四步：选择PDE菜单中PDE Mode命令，进入PDE模

11、式。单 击PDE菜单中PDE Specification选项，打开PDE对话框，设置方 程类型。本例取缺省设置，类型为椭圆型，参数 c,a,f 分别为 1 ， 0， 10。第五步：选择 Mesh菜单中Initialize Mesh命令，进行网格剖分。 第六步：选择 Mesh菜单中Refine Mesh命令，对网格加密。第七步：选择Solve菜单中So1ve PDE命令，解偏微分方程并显 示图形解。第八步：单击 Plot 菜单中 Parameters, 选项，打开 Plot selection 对话框，选中Color, Height (3D Plot)和Show mesh三项。然后单击 Plot

12、 按钮，显示三维图形解。第九步：如果要画等值线图和矢量场图，单击 Plot 菜单中Parameters 选项，打开 Plot Selection 对话框.选中 Con tour 和 Arrows两项。然后单击Plot按钮，可显示解的等值线图和矢量场图。第二章 PDE图形用户界面2.1 PDE Toolbox 菜单File菜单（如图1-1）File EditOptions DraNewCtrl+N51+0SaveCtrl+SSave M.Print,.Exit匚 trl+wl图1-1New新建一个几何结构实体模型（Constructive Solid Geomery简记为CSG），默认文件名为“

13、 Un titled ”。Open,从硬盘装载M文件Save将在GUI内完成的成果储存到一个 M文件中。Save As将在GUI内完成的成果储存到另外一个 M文件中。Print,将PDE工具箱完成的图形送到打印机内进行硬拷贝Exit退出PDE工具图形用户界面。2 Edit菜单（如图1-2）CutOtrl+X3*Qrl+CPaste.Orl+VClearQrl+RSelect AllQrl+A图1-2Undo在绘制多边形时退回到上一步操作Cut将已选实体剪切到剪贴板上。Paste将剪贴板上的实体粘贴到当前几何结构实体模型中Clear删除已选的实体。Select All选择当前几何结构实体造型 C

14、SG中的所有实体及其边界和字域。Refresh3 Optio ns 菜单（如图 1-3）Grid绘图时打开或关闭栅格。Grid Spac ing,调整栅格的大小。Snap打开或关闭捕捉栅格功能。Axes Limits,设置绘图轴的坐标范围。Axes Equal打开或关闭绘图方轴。图1-3Turn off Toolbar Help关闭工具栏按钮的帮助信息。Zoom打开或关闭图形缩放功能。Applicati on选择应用的模式。Refresh重新显示PDE工具箱中的图形实体。4 Draw菜单(如图1-4)DrakV Boundary PDE 処Hi 创 Plot Window HOrAw Mcdf

15、iiR&ctangle/squareRKtangle/iquare (centered)EJhpse/ 口 ixleElipsc/orcls (centered)敢冷gon图1-4iRotate.,Export Geometn1 Descrptcn. Set Formula., Latwls.Draw Mode进入绘图模式。Rectangle/square以角点方式画矩形/方行（Ctrl+鼠标）。Recta ngle/square （cen tered）以中心方式画矩形/方行（Ctrl+鼠标）。Ellipse/circle以矩阵角点方式画椭圆/圆（Ctrl+鼠标）。Ellipse/ci

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

17、olve Plot Window HelpBaundary MtideCtrl+BSpecify EtaundarY Conditions.Shaw Edge Labelsemoe 5ubdornain Border图1-5RornDve All Sdbdomain 日口rdersExport Decomposed 匚日口meitryboundary Condfs：.Bou ndary Mode进入边界模式。Specify Bou ndary Con ditio ns,对于已选的边界输入条件，如果没有 选择边界，则边界条件适用于所有的边界Show Edge Labels 显示边界区域标识开关，

18、 其数据是分解几何矩阵 的列数。Show Subdomain Labels 显示子区域标识开关， 其数据是分解几何 矩阵中的子域数值。Remove Subdomain Border 当图形进行布尔运算时，删除已选取的子 域边界。Remove All Subdomain Borders 当图形进行布尔运算时， 删除所有的 子域边界。Export Decomposed Geometry， Boundary Cond's,将分解几何矩阵g、边界条件矩阵b输 出到主工作空间。选择 Boundary 菜单中 Specify Boundary Conditions.命令可定义边 界条件。在打开的

19、Boundary condition 对话框，可对已选的边界输入 边界条件。共有如下三种不同的条件类型：NeMmann条件这里边界条件是由方程系数q和g确定的，在方程组的情况下（换成方程组模式），q是2 x 2矩阵，g是2x1矢量。Dirichlet条件 u定义在边界上，边界条件方程是价 h*u=r，这 里 h 是可以选样的权因子 （通常为 1）。在方程组情况下， h 是 2x2 矩 阵， r 是 2x l 矢量，混合边界条件 （仅适合于方程组情形 ） 它是 Dirichlet 和 Neumann 的混合边界条件， q 是 2x 2 矩阵， g 是 2x1 矢量， h 是 1x 2 矢量， r

20、是一个标量。6 PDE菜单（如图1-6）PDE Mesh Solve Plot WindcPDE ModeShow Sub domain LabelsPDE Specification,图1-6PDE ModeShow Subdoma in LabelsPDE Specificati on,Export PDE Coefficie nts,Export PDE Coefficients.,进入偏微分方程模式。显示子区域标识开关。调整PDE参数和类型。将当前PDE参数c,a,f,d输出到主工作空间，其参数变量为字符类型7 Mesh菜单（如图1-7）社esh SoNe Plot 理indc Hel

21、pMesh ModeInitiaize Mbs hiCtrl+1Refine MteshCtrl+MJiggle MeshUndo Mesh ChangeByplay Triangle QualityShow 也ode UbebShow Iningle LabelsMesh Mode输入网格模式。In itialize Mesh建立和显示初始化三角形网格。Refine Mesh加密当前三角型网格。Jiggle Mesh优化网格。Undo Mesh Change退回上一次网格操作。Par-ameters.Export 图1-7Display Triangle Quality 用01之间数字化的颜

22、色显示三角形网格的质量，大于0.6的网格可接受的Show Node Labels 显示网格节点标识开关，节点标识数据是点矩阵p的列Show Tria ngle Labels显示三角形网格标识开关，三角形网格标识图1-8Solve PDE数据是三角形矩阵t的列。Parameters修改网格生成参数。Export Mesh输出节点矩阵p、边界矩阵e和三角形矩阵t到主工作空间。8 Solve菜单（如图1-8）Solve Plot Window HelpSolve RDECtrl+EParameterSiiiExport Solution.对当前的几何结构实体 CSG、三角形网格和图形解偏微分方程。P

23、arameters调整PDE的参数。Export Solutio n,输出PDE的解矢量u。如果可行，将计算的特征值1输出到主工作空间1.1.9 Plot 菜单（如图 1-9）Plot Solution显示图形解。Parameters打开绘图方式对话框。Export Movie,如果动画被录制了，则动画矩阵 M将输出到主工作空间。Plot Window HelpPlot Sulution 匚trl+F1Parameters.Expiort Movie,.图1-910 Window 菜单从Window菜单项，可选择当前打开的所有的 MATLAB图形窗 口，被选择的窗台移至前台。11 Help菜单

24、Help,显示帮助信息About,显示版本信息1. 2 PDE工具栏主菜单下是工具栏，工具栏中喊有许多工具图标按钮，可提 供快速、便捷的操作方式。从左到右5个按钮为绘图模式按钮，紧接 着的6个为边界、网格、解方程和图形显示控制功能按钮，最右边的 为图形缩放功能键。（如图1-10）图 1-10-口"以角点方式画矩形/方行（Ctrl+鼠标）。丨一 I以中心方式画矩形/方行（Ctrl+鼠标）。二以矩形角点长轴方式画椭圆/圆（Ctrl+鼠标）。I旦以中心方式画椭圆/圆（Ctrl+鼠标）。1- 画多边形，按右键可封闭多边形。EI进入边界模式。1- 打开PDE Specification （偏微

25、分方程类型）对话框。匚初始化三角形网格。互J加密三角形网格。解偏微分方程。打开Plot Selection对话框，确定后给出解的三维图形 d为显示缩放切换按钮。第三章典型方程及其应用实例求解PDE问题主要有两种方法，一种是使用图形用户界面，另 一种是采用命令行编程。前者直观简便，而后者更为灵活。2. 1求解椭圆方程的例子例：单位圆上的Poisson方程边值问题：I .''.u =1, I - x, y |x2 y2 : 1 :.uk 尸0这一问题的精确解为：2 2(1_x_y )u x,y 二4若使用图形用户界面(Graphical User In terface简记为GUI)

26、,则 首先在MATLAB的工作窗口中键入pdetool,按回车键确定，于是出 现PDE Toolbox窗口。如果需要坐标网格，单击Options菜单下的Grid 选项即可。下面分步进行操作。(i) 画区域圆 单击工具；三二，大致在(0, 0)位置单击鼠标 右键同时拖拉鼠标到适当位置松开， 绘制圆。为了保证所绘制的圆是 标准的单位圆，在所绘圆上双击，打开Object Dialog对话框，精确地 输入圆心坐标X-center为0、Y-cebter为0及半径Radius为1,然后 单击OK按钮，这样单位远已画好。(ii) 设置边界条件 单击工具还，图形边界变红，逐段双击 边界，打开Boundary

27、Condition对话框，输入边界条件。对于同一类 型的边界，可以按Shift键，将多个边界同时选择，统一设置边界条 件。本题选择Dirichlet条件，输入h为1, r为0,然后单击OK按钮。 也可以单击 Boundary 菜单中 Specify Boundary Conditions,选项，打 开Boundary Condition对话框输入边界条件，如图 2-1。(iii) 设置方程 单击PDE菜单中PDE Specification,选项，打 开PDE Specification对话框，选项方程类型。本题单击Elliptic，输入 c为1，a为0，f为1,然后单击OK按钮，如图2-2&

28、#163; llounidaiy ComdltronBoundiiy concftkrt广 Neurrjmq|oFh'ir厂OL|Cmd图2-1图2-2(iv )网格剖分 单击工具一，或者单击 Mesh菜单中Initialize Mesh选项，可进行初始网格剖分，这时在 PDE Toolbox窗口下方的 状态栏内显示初始问网格的节点数和三角形单元数。本题节点数为 144个，三角形单元数为254个。如果需要网格加密，再单击商 或者单击 Mesh菜单中Refine Mesh选项，这时节点数变为 541个， 三角形单元数为1016个，如此还可继续加密。(v) 解方程单击工具=，或者单击Sol

29、ve菜单中Solve菜单中Solve PDE选项，可显示方程色彩解。如果单击Plot菜单中Parameters选项，出现Plot Selection对话框，如图2-3,从中可以选 择 Color, Con tour, ArrowsQeformed mesh, Height (3-D polt), 还可 以设置等值线的数目等。本例中选择Color,Contour, Height(3-D polt) 和Show mesh四项，然后单击Plot按钮，方程的图形解如图2-4所示。 除了作定解问题解u的图形外，也可以作|grad u|, |cgrad u等图形。图2-3图2-4(vi) 与精确解作比较单

30、击Plot菜单中Parameters选项，打开 Plot Selection 对话框，在 Height ( 3-D plot)行 Property下拉框中选user entry,且在该行的User entry输入框中键入u- (1-x.A2-y.A2)/4,单击P lot按钮就可以看到解的绝对误差图形，如图2- 5.可见在边界处误差为0图2-5（vii）输出网格节点的编号、单元编号以及节点坐标单击Mesh菜单中Show Node Labels选项，再单击工具亠或亠：即可显示节点 编号。若要输出节点坐标，只需单击Mesh菜单中Export Mesh,选项，这时打开的Export对话框中默认值为p

31、, e, t，这里p, e, t分 别表示points （点）、edges （边）、triangles （三角形）数据的变量, 单击0K按钮。然后在MATLAB命令窗口键入p,按回车键确定， 即可显示出节点按编号排列的坐标（二维数组）；键入e,按回车键， 则显示边界线段数据矩阵（7维数组）;输入t,按回车键，则显示三 角形单元数据矩阵（4维数组）。（viii ）输出解的数值 单击Solve菜单中Export Solution,选项， 在打开的Export对话框中输入u,单击0K按钮确定。再在MATLAB 命令窗口中输入u,按回车确定，即显示按节点编号排列的解的数值。我们也可以用MATLAB程序

32、求解PDE问题，同时显示解的图形：p,e,t=initmesh( circleg ' , ' hmax' ,1);Error=;err=1;While err>0.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)

33、Pdesurf(p,t,u-exact)通过命令行键入help+命令函数，女口 help pdemesh按回车键，可 以调入有关命令函数的定义、参数格式等帮助信息。2. 2求解抛物型方程的例子例：考虑一个带有矩形孔的金属板上的热传导问题。板的左边保持在100 c，板的右边热量从板向环境空气定常流动，其他边及内孔边界 保持绝缘。初始t"。时板的温度为0,于是概括为如下定解问题：d 色 4=0, 盘u =100,岂，cnS,nu =°.域门的外边界顶点坐标为(-0.5, -0.8), (0.5, -0.8), (0.5, 0.8), (-0.5,0.8)o 内边界顶点坐标为(-

34、0.005, -0.4), (0.05, -0.4), (0.05, 0.4),(-0.05, 0.4)使用GUI求解这一问题。在PDE Toolbox窗口的工具栏中选择Gen eric 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栏

35、中键入R1-R2。(ii) 设置边界条件单击瓯，使边界变红色，然后分别双击每 段边界，打开Boundary Cvondition对话框，设置边界条件。在左边 界条件。在左边界上，选择Dirichlet条件，输入h为1, r为100；右 边界上，选择Neumann条件，输入g为-1, q为0;其他边界上，选 择Neumann条件，输入g为0, q为0。(iii) 设置方程类型单击Z,打开PDE Specification对话框， 设置方程类型为Parabolic (抛物型)，d=1, c=1, a=0, f=0,单击OK按钮（iv ）网格剖分单击葩，或者加密网格，单击亠。（v） 初值和误差的设置

36、单击Solve菜单中Parameters选项，打开 Solve Parameters对话框，输入 Time 为 0:5, u（t0）为 0，Relative toleranee为 0.01，Absolute toleranee为 0.001，然后单击 OK 按钮。（vi） 数值解的输出单击Solve菜单中Export Solution,选项，在Export对话框中输入u，单击OK按钮。再在 MATLAB命令窗口 中键入u，按回车键，这时显示出按节点编号的数值解。（vii ）解的图形 单击Plot菜单中Parameters选项，打开PlotSelection 对话框，选 Color，Con to

37、ur，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,0,1);pdeplot(p,e,t, Xydata；u(:,11), mesh", off', colormap' hot'2. 3求解双曲型方程的例子例：考虑如下二维波动方程的定解问题:.：t2-.：u =0,(x, y) i - "(x, y) | 一仁：x : 1, 一1 ： y : 1

38、 ；Ux = 1=0,.:u.:ny = 1 = 0,nt =0,u =arcta n( cos(x),2亠"sin(二x)eSin(3.:t用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 (

39、双曲型)，键入 c=1, a=0, f=0, d=1。单击Solve菜单中Parameters选项，打开 Solve Parameters寸话 框，在Time栏中键入linspace (0, 5, 31),设置u的初始值u (t0) 为ata n(cos(pi/2*x),u' 的 初 始 值 u'(t0) 为3*sin(pi*x).*exp(sin(pi/2*y),Relative toleranee 为 0.01 , Absolute toleranee为 0.001。最后输出图形解：单击 Plot菜单中Parameters选项，打开 PlotSelection 对话框，选

40、Color, Contour, Arrows , Height （3-D Plot） 和Show mesh五项，单击Plot按钮，三维彩色图形的解如图 2-7所示。图2-7MATLAB的动画程序：p,e,t=initmesh( Squareg);x=p(1,:) ;y=p(2,:)'u0=ata n(cos(pi/2*x);ut0=3*si n(pi*x).*exp(si n(pi/2*y);n=3l;tlist=li nspace(0,5, n);uu二hyperbolic(uO,utO,tlist, squareb3,p,e,t,1,0,0,1); delta=-1:0.1:1;u

41、xy,t n,a2,a3=tri2grid(p,t,uu(:,1),delta,delta);gp=tn;a2;a3;umax=max(max(uu);umi n=min(min (uu);n ewplotM=moviei n(n);For 1=1: n,pdeplot(p,e,t,ydata：uu(:,l), Zdata；uu(:,l),mesh", off', ygrid ',On' gridparam',gp, Colorbar； Off' zstyle', Continuous); axis(-1 1 1 umin umax);

42、caxis(umin umax); M(:,I)=getframe;EndMovie(M,10);2. 4求解特征值问题的例子例：混合边界条件的特征值问题:-u = u,(x, y) i(x, y)| 一1 ：x,y ： 1u |x二 0,3 u4,；:n -首先作区域：方形的4个顶点为(-1, 1), (-1 , -1) , (1,-1), (1,1)。在 Object Dialog 对话框中输入 Left 为-1，Bottom 为-1，Width 为 2, Height2。单击还，逐段双击边界设置边界条件：左边界选 Dirichlet条件, 输入h为1, r为0;右边界选Neumann条件

43、，输入q为-3/4, g为0； 上、下边界选Neumann条件，输入q为0, g为0。在PDE Specification对话框中设置方程类型为Eigenmodes(特征 值模式)，PDE的系数设置为c=1, a=0, d=1。由于右边界上Neumann条件，特征值可以出现负值，因此在求小于10的特征值时应在 Solve Parameters对话框中键入-inf 10。作网格剖分，然后单击 二，得到第一个特征值对应的特征函数图形，如图2-8图2-8再在Plot Selection对话框的Eigenvalue项中选取不同的特征值， 比如第二个特征值“ 2.06,并选 Color, Con tou

44、r, Height (3-D Plot) 和Show mesh四项，如图2-9，可作出第二个特征值对应的特征函数 的三维图形，如图2-10。图2-9图 2-10MATLAB 程序：p,e,t=initmesh( s'quareg'p,e,t二refinemesh squareg,p,e,t);v,l=pdeeig( squareb2,p,e,t,1,0,1,-inf 10);pdesurf(p,t,v(:,4)2. 5偏微分方程在一维空间的应用已知偏微分方程初始条件时间t和一维空间变量x, MATLAB自身存在对偏微分方程的解函数pdepeb单一方程假设我们要求解热传导方程Ut

45、 二 Uxxu (t,0) = 0 u( t,1) = 1u (0,x)二2x21 X(1.1)MATLAB指定抛物型的偏微分方程形式如下：_m mc（x,t,u，u）ut 二 X '（X b（x,t,u，ux）s（x,t,u，ux）一 x边界条件为：P（Xi，t，u）q（Xi，t）b（Xi，t,u，Ux）=0p（Xr，t,u） q（Xr，t） b（Xr，t,u，Ux）= 0其中X为边界左端点，Xr为边界右端点，初始条件为u（0, x） = f（x）。（我们注意到同一函数b在方程和边界条件中出现。）通常情况下，每个函数都会被指定到不同的MATLAB文档中。也就是说，与方程有关的函数c,

46、b和s都会被指定保存到一个 MATLAB文件中，与边 界条件有关的函数p和q保存到令一个文档中（此外，我们需要注意 到函数b是相同的只需保存一次即可），最后将初始函数f（x）保存在 第三个文档中。执行命令函数 pdepe将把m个文档结合起来进行运 算，返回方程的解。在我们的例子中，有：c（x,t,u，Ux）=1 b（X，t，u，Ux）二 Ux s（x,t,u，Ux）=0我们将它们保存到MATLAB文档中，记作eqn 1.m （m=0以后在加详 述）。Fun cti on c,b,s=e qn 1（x,t,u,DuDx）%EQN1:MATLAB fun ctio n M-file that sp

47、ecifies %a PDE in time and one space dime nsion.C=1;B=DuDx;S=0;对于边界条件，我们有:P(O,t,u) = u ; q(O,t) = 0P(1,t,u) = u -1 ;q(1,t) = 0将它们保存到MATLAB文档中，记作bcl.m.Fun cti on p1,q1,pr,qr二bc1(x1,u1,xr,ur,t) %BC1:MATLAB fun ctio n M-file that specifies bou ndary con diti ons %for a PDE in time and one space dimensi

48、on.P1=u1;Q仁0;Pr=ur-1;Qr=0;对初始条件，有:f（ X）二2x将它保存到MATLAB文档中，记作initial1.m.Fun cti on value二 in itial1(x)%INITIAL1:MATLAB function M-file that specifies the initial condition %for a PDE in time and one space dime nsion.Value=2*x/(1+xA2);最终我们准备好了用pdepe函数处理偏微分方程。在如下的MATLAB 文档中，我们选择坐标值x和t,求解偏微分方程并绘制解的表面图(如图

49、 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=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

50、PDEU=pdepe(m,epn1,initial1,bc1,x,t);%Plot solutionSurf(x,t,u);Title( Surface plot of solution');Xlabel( Distance x');Ylabel( Time t');Sw£blZ-« fka： cif1 1P-fi- 31 plh 1 £ IT Hf'III. Il i-iTI 1 (>11 1" -HI. |' 11 | )图 2-11通常，我们会发现绘制解的剖面图是十分有用的。此时，t是固定的，u是x的函

51、数。解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,:), 'rasW,'or'For k=2:le ngth(t)Set(fig, 'data',x, ydata；u(k,:)Pause(.5) end试过之后，你会发现它们是如何快速地接近热传导方程的平衡位置 的。(平衡位置就是及时停止并要发生改变的

52、位置)。Viulldlb Fl H-O- l.l frn>图 2-123 PDE Toolboxz中的命令简介本节将介绍PDE Toolbox中的命令函数。利用这些命令函数，可通过 命令行而不是用PDE Tool图形用户界面来解偏微分方程。3. 1 PDE Toolbox中的函数及其分类表3-1PDE数值计算函数函数目的Adaptmesh生成自适应网格和解PDE问题Assema组装积分区域贡献Assemb组装边界条件贡献assempde组装刚度矩阵和方程右边的函数hyperbolic解双曲型PDE问题parabolic解抛物型PDE问题Pdeeig解特征值PDE问题pde non li

53、n解非线性PDE问题poisolv求矩形网格上泊松方程的快速解表3-2用户界面算法函数目的Pdecirc画圆Pdeellip画椭圆Pdemdlcv转化为1.0版本的M文件Pdepoly画多边形Pderect画矩形Pdetool进入PDE工具箱的用户界面表3-3几何算法函数目的Csgchk检查几何描述矩阵的有效性Csgdel删除最小区域之间的分界线Decsg分解结构矩阵成最小区域In itmesh创建初始网格Jigglemesh调整三角形网格内的点Pdearcl在表达式参数和弧长之间插值Poimesh在矩形区域上创建规则网格Refin emesh加密三角形网格Wbou nd写边界条件说明文件Wg

54、eom写几何说明文件表3-4函数目的Pdec ont绘制等值线图的快速命令Pdegplot绘制PDE几何图Pdemesh绘制PDE三角形网格图Pdeplot一般PDE工具箱的绘制函数pdesurf绘制表面图的速写命令表3-5通用算法函数目的Dst离散正弦变换Pdeadgsc用相对误差判别准则选择最低劣的三角形Pdeadworst相对最差值选择低劣三角形Pdecgrad计算PDE的通量Pdee nt与已知三角形单元相邻的三角形单元的索引Pdegrad计算PDE解的梯度Pdei ntrp从单元节点数据到三角形单元中点的插值Pdejmps为自适应进行误差估计Pdeprt ni从三角形中点数据到节点数据的插值Pdesde子区域集中边缘的索引Pdesdp子区域集中点的索引Pdesdt子区域集中三角形的索引Pdesmesh计算结构力学张量函数Pdetrg三角形几何数据Pdetriq测量三角形网格质量Poiasma泊松方程快速解的边界点矩阵的贡献Poicalc矩形网格上泊松方程的快速解Poii ndex矩形网格正规次序的索引Sptar n解一般稀疏矩阵特征值问题Tri2grid由PDE三角形网格转换成矩形网格表3-6用户定义算法函数目的Pdebound创建边界M文件Pdegeom创建几何M文件3. 2 PDE数值计算函数简介表3-7属性名属性值缺省值说明Maxt