偏微分方程数值解.pdf

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

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

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

3、型方程：(),c uaufin，椭圆型方程：(),c uauf in 其中是平面有界区域，c,a,f 以及未知数 u 是定义在上的实（或复）函数。抛物型方程：(),.udc uaufint 双曲型方程：22(),uc uaufint.特征值方程：(),c uauduin 其中 d 是定义在上的复函数，是待求特征值。在抛物型方程和双曲型方程中，系数 c,a,f 和 d 可以依赖于时间 t。可以求解非线性椭圆型方程：()()(),c uua uf uin 其中 c,a,f 可以是未知函数 u 的函数。还可以求解如下 PDE 方程组；11112211 1122121122221 12221()(),

4、()()cuucuua ua ufcuucuua ua uf 利用命令行可以求解高阶方程组。对于椭圆型方程，可以用自适应网格算法，还能与非线性解结合起来使用。另外，对于 Poission 方程还有一个矩形网格的快速求解器。122 边界条件（1）Dirichlet 条件：hur (2)Neumann 条件：()nc uqug 其中n是的边界上的单位外法向量，,g q h和r是定义在上的函数。对于特征值问题仅限于齐次条件：0,g 和0r。对于非线性情形系数,g q h和r可以依赖于 u；对于抛物型方程和双曲型方程，系数可以依赖于时间 t。对于方程组情形，边界条件为(1)Dirichlet 条件：1

5、1 11221h uh ur 21 12222h uh ur(2)Neumann 条件：11112211 11221()()ncuncuq uq ug 21122221 12222()()ncuncuq uq ug(3)混合边界条件为:11 11221h uh ur 11112211 1122111()()ncuncuq uq ugh 21122221 1222212()()ncuncuq uq ugh 其中的计算要使得 Dirichlet 条件满足。在有限元法中，Dirichlet 条件也称为本质边界条件，Neumann 条件称为自然边界条件。1.3 如何使用 FDE Toolbox 1.3

6、.1 定解问题的设置 员简单的办法是在 PDE Tool 上直接使用图形用户界面(GUl)。设置定解问题包括三个步骤：(1)Draw 模式：使用 CSG(几何结构实体模型)对话框画几何区域，包括矩形、圆、椭圆和多边形，也可以将它们组合使用。(2)Boundary 模式：在各个边界段上给出边界条件，(3)PDE 模式：确定方程的类型、系数 c，a，f 和 d c。也能够在不同子区域上设置不同的系数(反映材料的性质)。1.3.2 解 PDE 问题 用 GUI 解 PDE 问题主要经过下面两个过程（模式）(1)Mesh 模式；生成网格自动控制网格参数。(2)Solve 模式：对于椭圆型方程还能求非线

7、性和自适应解。对于抛物型和双曲型力程设置初始边值条件后能求出给定 t 时刻的解。对于特征值问题，能求出给定区间内的特征值；求解后可以加密网格再求解。1.3.3 使用 Toolbox 求解非标准的问题 对于非标准的问题。可以用 PDE Too1box 的函数。或者用 FEM(有限元法)求解更为复杂的问题。1.3.4 计算结果的可视化 从 GUI 能够使用 Plot 模式实现可视化。可以使用 Color,Height和 Vector 等作图。对于抛物型和双曲型方程，还可以生成解的动画。这些操作通过命令行都很容易实现。1.3.5 应用领域 在应用界面提供了丁如下应用领域 结构力学平面应力问题 结构力

8、学平面应变问题 静电场问题 静磁场问题 交流电磁场问题 直流导体介质问题 热传导问题 9散问题 这些界面都有对话框，它包括 PDE 的系数、边界条件、解的性质等。1.4 解偏微分方程的一个例子 解 Poisson 方程uf，边界条件为齐次 Dirichlet 类型。第一步：启动 MATLABl,键入 pdetool，按回车键确定便可启动GUI，然后在 Options 菜单下选择 Grid 命令，打开栅格,栅格的使用，能使用户容易确定所绘图形的大小，如图 11 1-1 第二步：分步完成平面几何造型：R1-C1-E1+R2+C2。用菜单或快捷工具，分别画矩形 R1、矩形 R2、椭圆 E1、圆 C1

9、、圆 C2。画圆时，首先选中椭圆工具，按鼠标右键并拖动即可、或者在按 ctrI 的同时，拖动鼠标也可绘制圆。然后在 Set formula 栏，进行编辑并用算术运算将将图形对象名称连接起来，删除默认的表达式键入R1-C1-E1+R2+C2，按等号健得到所需图形。若需要，还可进行储存 形成 M 文件。选择 Boundary 菜单中 Boundary Mode 命令，进入边界模式。单击 Boundary 菜单中 Remove A11 Subdomain Borders 选项，去除子域边界。如果想将几何信息和边界信息进行存储，应选择 Boundary 菜单中的 ExPort Decomposed G

10、eometryBoundary Conds命令，将它们分别储存于 g,b 变量中,通过 MATLAB 形成 M 文件。第三步：选取边界单击 Boundary 菜单中 Specify Bounddy Conditions选项，打开 Boundary conditlons 对话框，输入边界条件，如图 14。本例取缺省条件。即将全部边界设为齐次 Dirichlet 条件，边界颜色显示为红色。第四步：选择 PDE 菜单中 PDE Mode 命令，进入 PDE 模式。单击 PDE 菜单中 PDE Specification选项，打开 PDE 对话框，设置方程类型。本例取缺省设置，类型为椭圆型，参数 c,

11、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 按钮，显示三维图形解。第九步：如果要画等值线图和矢量场图，单击 Plot 菜单中Parameters选项，打开 Plot Sele

12、ction 对话框 选中 Contour 和 Arrows两项。然后单击 P1ot 按钮，可显示解的等值线图和矢量场图。第二章 PDE 图形用户界面 2.1 PDE Toolbox 菜单 File 菜单（如图 1-1）图 1-1 New 新建一个几何结构实体模型（Constructive Solid Geomery,简记为 CSG），默认文件名为“Untitled”。Open 从硬盘装载 M 文件 Save 将在 GUI 内完成的成果储存到一个 M 文件中。Save As 将在 GUI 内完成的成果储存到另外一个 M 文件中。Print 将 PDE 工具箱完成的图形送到打印机内进行硬拷贝。Ex

13、it 退出 PDE 工具图形用户界面。2 Edit 菜单（如图 1-2）图 1-2 Undo 在绘制多边形时退回到上一步操作。Cut 将已选实体剪切到剪贴板上。Copy 将已选实体拷贝到剪贴板上。Paste 将剪贴板上的实体粘贴到当前几何结构实体模型中。Clear 删除已选的实体。Select All 选择当前几何结构实体造型 CSG 中的所有实体及其边界和字域。3 Options 菜单（如图 1-3）图 1-3 Grid 绘图时打开或关闭栅格。Grid Spacing 调整栅格的大小。Snap 打开或关闭捕捉栅格功能。Axes Limits 设置绘图轴的坐标范围。Axes Equal 打开或

14、关闭绘图方轴。Turn off Toolbar Help 关闭工具栏按钮的帮助信息。Zoom 打开或关闭图形缩放功能。Application 选择应用的模式。Refresh 重新显示 PDE 工具箱中的图形实体。4 Draw 菜单（如图 1-4）图 1-4 Draw Mode 进入绘图模式。Rectangle/square 以角点方式画矩形/方行（Ctrl+鼠标）。Rectangle/square（centered）以中心方式画矩形/方行（Ctrl+鼠标）。Ellipse/circle 以矩阵角点方式画椭圆/圆（Ctrl+鼠标）。Ellipse/circle（centered）以中心方式画椭圆

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

16、ary Conditions 对于已选的边界输入条件，如果没有选择边界，则边界条件适用于所有的边界。Show Edge Labels 显示边界区域标识开关，其数据是分解几何矩阵的列数。Show Subdomain Labels 显示子区域标识开关，其数据是分解几何矩阵中的子域数值。Remove Subdomain Border 当图形进行布尔运算时，删除已选取的子域边界。Remove All Subdomain Borders 当图形进行布尔运算时，删除所有的子域边界。Export Decomposed Geometry，Boundary Conds 将分解几何矩阵 g、边界条件矩阵 b 输出

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

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

19、网格。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 Solv

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

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

22、型）对话框。初始化三角形网格。加密三角形网格。解偏微分方程。打开 Plot Selection 对话框，确定后给出解的三维图形。为显示缩放切换按钮。第三章 典型方程及其应用实例 求解 PDE 问题主要有两种方法，一种是使用图形用户界面，另一种是采用命令行编程。前者直观简便，而后者更为灵活。21 求解椭圆方程的例子 例：单位圆上的 Poisson 方程边值问题：221,|1|0ux yxyu 这一问题的精确解为：221,.4xyu x y 若使用图形用户界面（Graphical User Interface,简记为 GUI），则首先在 MATLAB 的工作窗口中键入 pdetool，按回车键确定

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

24、，可以按 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。图 2-1 图 2-2

25、（iv）网格剖分 单击工具，或者单击 Mesh 菜单中 Initialize Mesh 选项，可进行初始网格剖分，这时在 PDE Toolbox 窗口下方的状态栏内显示初始问网格的节点数和三角形单元数。本题节点数为144 个，三角形单元数为 254 个。如果需要网格加密，再单击，或者单击 Mesh 菜单中 Refine Mesh 选项，这时节点数变为 541 个，三角形单元数为 1016 个，如此还可继续加密。（v）解方程 单击工具，或者单击 Solve 菜单中 Solve 菜单中 Solve PDE 选项，可显示方程色彩解。如果单击 Plot 菜单中Parameters选项，出现 Plot

26、Selection 对话框，如图 2-3，从中可以选择 Color，Contour，Arrows,Deformed 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）与精确解作比较 单击 Plot 菜单中 Parameters选项，打开 Plot Selection 对话框，在 Height（3-D pl

27、ot）行 Property 下拉框中选 user entry，且在该行的 User entry 输入框中键入 u-（1-x.2-y.2）/4,单击lot 按钮就可以看到解的绝对误差图形，如图-可见在边界处误差为。图 2-5（vii）输出网格节点的编号、单元编号以及节点坐标 单击 Mesh菜单中 Show Node Labels 选项，再单击工具或，即可显示节点编号。若要输出节点坐标，只需单击 Mesh 菜单中 Export Mesh选项，这时打开的 Export 对话框中默认值为 p，e，t，这里 p，e，t 分别表示 points（点）、edges（边）、triangles（三角形）数据的变

28、量，单击 OK 按钮。然后在 MATLAB 命令窗口键入 p，按回车键确定，即可显示出节点按编号排列的坐标（二维数组）；键入 e，按回车键，则显示边界线段数据矩阵（7 维数组）；输入 t，按回车键，则显示三角形单元数据矩阵（4 维数组）。（viii）输出解的数值 单击 Solve 菜单中 Export Solution选项，在打开的 Export 对话框中输入 u，单击 OK 按钮确定。再在 MATLAB命令窗口中输入 u，按回车确定，即显示按节点编号排列的解的数值。我们也可以用 MATLAB 程序求解 PDE 问题，同时显示解的图形：p,e,t=initmesh(circleg,hmax,1

29、);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;End Pdemesh(p,e,t)Pdesurf(p,t,u)Pdesurf(p,t,u-exact)通过命令行键入 help+命令函数，如 help pdemesh，按回车键，可以调入有关命令函数的定义、参数格式等帮助信息。22 求解抛物型方程的例子 例：考虑一个带有矩

30、形孔的金属板上的热传导问题。板的左边保持在100c，板的右边热量从板向环境空气定常流动，其他边及内孔边界保持绝缘。初始0tt时板的温度为 0，于是概括为如下定解问题：00,100,1,0,|0.t tudutuununu 域的外边界顶点坐标为（-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 窗口的工具栏中选择Generic Scalar 模式。（i）区域设置 单击工具，在窗口拖拉出一个矩形

31、，双击矩形区域，在 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 对话框，设置边界条件。在左边界条件。在左边界上，选择 Dirichlet 条件，输入 h 为 1，r 为 100；

32、右边界上，选择 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 to

33、lerance 为 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-6 MATLAB 程序：p,e,t=initmesh(crackg);u=para

34、bolic(0,0:0.5:5,crackb,p,e,t,1,0,0,1);pdeplot(p,e,t,xydata,u(:,11),mesh,off,colormap,hot)23 求解双曲型方程的例子 例：考虑如下二维波动方程的定解问题：).2sin(1122)sin(3),2(arctan(cos,0,0|,0|,11,11|),(),(,0 xyxextuxutnuuyxyxyxutu 用 GUI 求解。类似前面的例子，首先作正方形区域：设置断点坐标为（-1，1），（-1，-1），（1，-1），（1，1）。在 Object Dialog 对话框中输入 Left 为-1，Bottom 为

35、-1，Width 为 2，Height2，单击 OK 按钮。设置边界条件：左、右边界用 Dirichlet 条件，输入 h 为 1，r 为 0；上、下边界用 Neumann 条件，输入 q 为 0，g 为 0，作网格剖分。设置方程类型为 Hyperbolic（双曲型），键入 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

36、/2*y),Relative tolerance为0.01，Absolute tolerance 为 0.001。最后输出图形解：单击 Plot 菜单中 Parameters选项，打开 Plot Selection 对话框，选 Color，Contour，Arrows，Height（3-D Plot）和Show mesh 五项，单击 Plot 按钮，三维彩色图形的解如图 2-7 所示。图 2-7 MATLAB 的动画程序：p,e,t=initmesh(squareg);x=p(1,:);y=p(2,:);u0=atan(cos(pi/2*x);ut0=3*sin(pi*x).*exp(sin(

37、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;a2;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,

38、colorbar,off,zstyle,continuous);axis(-1 1 1 1 umin umax);caxis(umin umax);M(:,I)=getframe;End Movie(M,10);24 求解特征值问题的例子 例：混合边界条件的特征值问题：,0|43,0|,0|,1,1|),(),(,111xyxununuuyxyxyxuu 首先作区域：方形的 4 个顶点为（-1，1），（-1，-1），（1，-1），（1，1）。在 Object Dialog 对话框中输入 Left 为-1，Bottom 为-1，Width为 2，Height2。单击，逐段双击边界设置边界条件：左

39、边界选 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 Parameters 对话框中键入-inf 10。作网格剖分，然后单击，得到第一个特征值对应的特征函数图形，如图 2-8。图 2-8 再在 Plot Se

40、lection 对话框的 Eigenvalue 项中选取不同的特征值，比如第二个特征值“2.06，并选 Color，Contour，Height（3-D Plot）和 Show mesh 四项，如图 2-9，可作出第二个特征值对应的特征函数的三维图形，如图 2-10。图 2-9 图 2-10 MATLAB 程序：p,e,t=initmesh(squareg);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)25 偏微分方程在一维空间的应用 已知偏微分方程初始条件

41、时间 t 和一维空间变量 x，MATLAB 自身存在对偏微分方程的解函数 pdepe。单一方程 假设我们要求解热传导方程 uuxxt 00,）（tu，11,）（tu xxxu212,0）（1.1）MATLAB 指定抛物型的偏微分方程形式如下：),(),(),(uuxxuuxxmxmtxutxsutxbutxc）（边界条件为：),(),(),(uxxxxlllutbtqutp0 0),(),(),(uxxxxrrrutbtqutp 其中xl为边界左端点，xr为边界右端点，初始条件为)()0(xfxu，。（我们注意到同一函数 b 在方程和边界条件中出现。）通常情况下，每个函数都会被指定到不同的 M

42、ATLAB 文档中。也就是说，与方程有关的函数 c,b 和 s 都会被指定保存到一个 MATLAB 文件中，与边界条件有关的函数 p 和 q 保存到令一个文档中（此外，我们需要注意到函数 b 是相同的只需保存一次即可），最后将初始函数)(xf保存在第三个文档中。执行命令函数 pdepe 将把 m 个文档结合起来进行运算，返回方程的解。在我们的例子中，有：1),(uxutxc uuxxutxb),(0),(uxutxs 我们将它们保存到 MATLAB 文档中，记作 eqn1.m(m=0 以后在加详述)。Functionc,b,s=eqn1(x,t,u,DuDx)%EQN1:MATLAB func

43、tion 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 将它们保存到 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 dimens

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

45、图（如图 2-11）%PDE1：MATLAB script M-file that solves and plots%solutions to the PDE stored in eqn1.m M=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 mesh x=linspace(0,1,20);t=linspace(0,2,10);%Solve

46、the PDE U=pdepe(m,epn1,initial1,bc1,x,t);%Plot solution Surf(x,t,u);Title(Surface plot of solution);Xlabel(Distance x);Ylabel(Time t);图 2-11 通常，我们会发现绘制解的剖面图是十分有用的。此时，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）。其实，我们有

47、一种比较快的方法去绘制剖面图，如下的MATLAB 序列：Fig=plot(x,u(1,:),erase,xor)For k=2:length(t)Set(fig,xdata,x,ydata,u(k,:)Pause(.5)end 试过之后，你会发现它们是如何快速地接近热传导方程的平衡位置的。（平衡位置就是及时停止并要发生改变的位置）。图 2-12 3 PDE Toolboxz 中的命令简介 本节将介绍 PDE Toolbox 中的命令函数。利用这些命令函数，可通过命令行而不是用 PDE Tool 图形用户界面来解偏微分方程。31 PDE Toolbox 中的函数及其分类 表 3-1 PDE 数值

48、计算函数 函 数 目 的 Adaptmesh 生成自适应网格和解 PDE 问题 Assema 组装积分区域贡献 Assemb 组装边界条件贡献 assempde 组装刚度矩阵和方程右边的函数 hyperbolic 解双曲型 PDE 问题 parabolic 解抛物型 PDE 问题 Pdeeig 解特征值 PDE 问题 pdenonlin 解非线性 PDE 问题 poisolv 求矩形网格上泊松方程的快速解 表 3-2 用户界面算法 函 数 目 的 Pdecirc 画圆 Pdeellip 画椭圆 Pdemdlcv 转化为 1.0 版本的 M 文件 Pdepoly 画多边形 Pderect 画矩形

49、 Pdetool 进入 PDE 工具箱的用户界面 表 3-3 几 何 算 法 函 数 目 的 Csgchk 检查几何描述矩阵的有效性 Csgdel 删除最小区域之间的分界线 Decsg 分解结构矩阵成最小区域 Initmesh 创建初始网格 Jigglemesh 调整三角形网格内的点 Pdearcl 在表达式参数和弧长之间插值 Poimesh 在矩形区域上创建规则网格 Refinemesh 加密三角形网格 Wbound 写边界条件说明文件 Wgeom 写几何说明文件 表 3-4 函 数 目 的 Pdecont 绘制等值线图的快速命令 Pdegplot 绘制 PDE 几何图 Pdemesh 绘制

50、 PDE 三角形网格图 Pdeplot 一般 PDE 工具箱的绘制函数 pdesurf 绘制表面图的速写命令 表 3-5 通 用 算 法 函 数 目 的 Dst 离散正弦变换 Pdeadgsc 用相对误差判别准则选择最低劣的三角形 Pdeadworst 相对最差值选择低劣三角形 Pdecgrad 计算 PDE 的通量 Pdeent 与已知三角形单元相邻的三角形单元的索引 Pdegrad 计算 PDE 解的梯度 Pdeintrp 从单元节点数据到三角形单元中点的插值 Pdejmps 为自适应进行误差估计 Pdeprtni 从三角形中点数据到节点数据的插值 Pdesde 子区域集中边缘的索引 Pd