偏微分方程边值问题的数值解法

求解偏微分方程的边值问题

本实验学习使用MATLAB的图形用户命令pdetool来求解偏微分方程的边值问题。这个工具是用有限元方法来求解的，而且采用三角元。我们用内个例题来说明它的用法。

一、MATLAB支持的偏微分方程类型

考虑平面有界区域D上的二阶椭圆型PDE边值问题：

(cu)uf

(1.1)

其中

(1) , (2) a,f是D上的已知函数 (3) c是标量或22的函数方阵 xy

未知函数为u(x,y) (x,y)D。它的边界条件分为三类：

（1）Direchlet条件：

huf (1.2)

（2）Neumann条件：

n(cu)qug (1.3)

（3）混合边界条件：在边界D上部分为Direchlet条件，另外部分为Neumann条件。

1

其中h,r,q,g,c是定义在边界D的已知函数，另外c也可以是一个2*2的函数矩阵，n是沿边界的外法线的单位向量。

在使用pdetool时要向它提供这些已知参数。

二、例题

例题1 用pdetool求解

u1 D: x2y21uD0

(1.4)

解 ：首先在MATLAB 的工作命令行中键入pdetool 口，选Genenic Scalar模式．

2

按回牟键确定，于是出现PDE Toolbox 窗，

( l ）画区域圆

单击椭圆工具按钮，大致在（0，0）位置单击鼠标右键，拖拉鼠标到适当位置松开。为了保证所绘制的圆是标准的单位园，在所绘园上双击，打开 Object Dialog 对话框，精确地输入圆心坐标X-center 为0 、Y-center 为0 及半径Radius 为l ，然后单击OK 按钮，这样单位画已画好．

3

( 2 ）设置边界条件

单击工具边界模式按钮 ，图形边界变红，逐段双击边界，打开Boundary condition 对话框．输入边界条件．对于同一类型的边界，可以按Shift键，将多个边界同时选择，统一设边界条件．本题选择Dirichlet 条件，输入h 为1 , r 为0。，然后单击OK 按钮．也可以单击Boundary菜单中Spocify Boundary Condition …选项，打开Boundary Condition 对话框输入边界条件．

4

( 3 ）设置方程

单击偏微分方程按钮 ，打开PDE Specification 对话框，选择方程类型· 本题选Ellintic（椭圆型）,输入c为1 , a 为O , f 为1 ，然后单击OK 按钮．

5

( 4 ）网格剖分

单击网格工具，或者单击Mesh 菜单中Initialize Mesh项，可进行初始网格剖分．这时在PDE Toolbox 窗口下方的状态栏内显示出初始网格的节点数和三角形单元数．本题节点数为144 个，三角形单元数为254 个(图?? )。如果要细化网格，单击细化工具，或者单击Mesh 菜单中Refine Mesh 选项，节点数成为541 个，三角形单元数为1016 个。

6

7

( 5 ）解方程

单击解方程工具，或者单击S olve菜单中Solve PDE 选项，可求得方程数值解并用彩色图形显示。单击作图工具，或者单击Plot 菜单中Parameter…选项，出现Plot selection 对话框．从中选择于Height ( 3-D plot) ，然后单击Plot 按钮，方程的图形解如图?? 所示。除了作定解问题解u的图形外，也可以作

grad u等图形·

8

9

10

（6）输出网格节点的编号、单元编号以及节点坐标

单击Mesh 菜单中Show Node Labels选项，再单击网格工具，即可显示节点编号（图？？ ) 。若要输出节点坐标，只需单击Mesh 菜单中Export Mesh … 选项，这时打开的Export对话框中的默认值为p e t,这里p、e、t 分别表示point （点）、edges(边）、triangles(三角形）数据变量，单击OK按钮，然后在MATLAB 命令行键入p,即可以显示按节点编号排列的坐标；键入e再回车则显示边界数据矩阵（7维数组）；键入t按回车则显示三角形单元数据矩阵（4维数组）。

11

点、边、单元的部分输出为：

p =

Columns 1 through 11

-1.0000 0.0000 1.0000 0.0000 -0.7071 0.7071 0.7071 -0.7071 -0.9808 -0.9239 -0.8315

12

-0.0000 -1.0000 0 1.0000 -0.7071 -0.7071 0.7071 0.7071 -0.1951 -0.3827 -0.5556

e =

Columns 1 through 11

1.0000 9.0000 10.0000 2.0000 15.0000 16.0000

9.0000 10.0000 11.0000 15.0000 16.0000 17.0000

0 0.1250 0.2500 0 0.1250 0.2500

0.1250 0.2500 0.3750 0.1250 0.2500 0.3750

1.0000 1.0000 1.0000 2.0000 2.0000 2.0000

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

0 0 0 0 0 0

11.0000 5.0000 0.3750 0.5000 1.0000 1.0000 0 13

5.0000 12.0000 0.5000 0.6250 1.0000 1.0000 0 12.0000 13.0000 0.6250 0.7500 1.0000 1.0000 0 13.0000 0.7500 0.8750 1.0000 1.0000 0 14.0000 2.0000 0.8750 1.0000 1.0000 1.0000 0 14.0000

t =

Columns 1 through 18

32 14 20 26 29 17 23 100 11 66 89 1 9 94 5 12 13 97

1 2 3 4 8 6 7 28 5 32 1 9 10 5 12 13 14 2

89 97 81 98 84 92 99 127 94 89 119 119 95 118 118 90 70 126

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

（7）输出近似解数值

单击Solve菜单中Export Solution。。。选项，在Export对话框中输入u再单击OK按钮，再在MATLAB命令行中输入u并回车，就会显示按节点编号排列的解u的数值。

（8）近似解和准确解的比较

方程(1.4)的准确解为：

1x2y2u(x,y)4

(1.5)

14

为了与准确解比较，单击Plot菜单中Parameters…选项，打开Plot Selection对话框，在Height（3-D plot）行的Property下拉框中选User Entry,并且输入

u-（1-x.^2-y.^2）/4 ,单击Plot按钮，就可以看到误差曲面，其数量级为10。

4

15

练习1 练习2

用pdetool工具求解课本P128的第2、3、4这几题的解，并作出图形。

用pdetool工具求解课本P418的第2、4这两题的解（三角单单元形状不限），并作出16

图形。