偏微分方程边值问题的数值解法.docx
《偏微分方程边值问题的数值解法.docx》由会员分享，可在线阅读，更多相关《偏微分方程边值问题的数值解法.docx（11页珍藏版）》请在冰豆网上搜索。
偏微分方程边值问题的数值解法
求解偏微分方程的边值问题
本实验学习使用MATLAB的图形用户命令pdetool来求解偏微分方程的边值问题。
这个工具是用有限元方法来求解的，而且采用三角元。
我们用内个例题来说明它的用法。
一、MATLAB支持的偏微分方程类型
考虑平面有界区域D上的•二阶椭圆型PDE边值问题：
-VHcVw）+aw=f
其中
（1）V=
（2）叮退D上的己知函数（3）c是标赞或2x2的函数方阵
未知函数为
“3）（x^）e£>
O它的边界条件分为三类：
（1）Direchlet条件：
加=/
（2）Neumann条件：
wC（cVw）+^w=g
（3）混合边界条件：
在边界
dD
上部分为Direchlet条件，另外部分为Neumann条件。
其中是定义在边界
dD的已知函数，另外
C也可以是一个2*2的函数矩阵，
n
是沿边界的外法线的单位向量。
在使用pdetool时要向它提供这些已知参数。
二、例题
例题1用pdetool求解
-Am=1D:
x?
+y2<1
wLd=0
解：
首先在MATLAB的工作命令行中键入pdetool,按回牟键确定，于是出现PDEToolbox窗口，选GenenicScalar模式.
（1）画区域圆
单击椭圆工具按钮，大致在（0,0）位置单击鼠标右键，拖拉鼠标到适当位置松开。
为了保证所绘制的圆是标准的单位园，在所绘园上双击，打开ObjectDialog对话框，精确地输入圆心坐标X-center为0、Y-center为0及半径Radius为1,然后单击0K按钮，这样单位画已画好.
（2）设置边界条件
单击工具边界模式按钮，图形边界变红，逐段双击边界，打开Boundarycondition对话框.输入边界条件.对于同一类型的边界，可以按Sh辻t键，将多个边界同时选择，统一设边界条件.本题选择Dirichlet条件，输入h为1,r为0。
，然后单击0K按钮.也可以单击Boundary菜单中Spoc辻yBoundaryCondition…选项，打开BoundaryCondition对话框输入边界条件.
（3）设冒方程
单击偏微分方程按钮，打开PDESpec辻ication对话框，选择方程类型・本题选Ellintic（椭圆型），输入c为1,&为0,f为1，然后单击0K按钮.
（4）网格剖分
单击网格工具，或者单击Mesh菜单中InitializeMesh项，可进行初始网格剖分.这时在PDEToolbox窗口下方的状态栏内显示出初始网格的节点数和三角形单元数.本题节点数为144个，三角形单元数为254个（图?
？
）。
如果要细化
网格，单击细化工具，或者单击Mesh菜单中RefineMesh选项，节点数成为541个，三角形单元数为1016个。
（5）解方程
单击解方程工具，或者单击Solve菜单中SolvePDE选项，可求得方程数值解并用彩色图形显示。
单击作图工具，或者单击Plot菜单中Parameter…选项，出现Plotselection对话框.从中选择于Height（3~Dplot）,然后单击Plot按钮，方程的图形解如图?
？
所示。
除了作定解问题解U的图形外，也可以作
等图形•
（6）输出网格节点的编号、单元编号以及节点坐标
单击Mesh菜单中ShowNodeLabels选项，再单击网格匸具，即可显示节点编号（图？
？
）。
若要输出节点坐标，只需单击Mesh菜单中ExportMesh•••选项，这时打开的Export对话框中的默认值为pet,这里p、e、t分别表示point（点）、edges（边）、triangles（三角形）数据变量，单击0K按钮，然后在MATLAB命令行键入p,即可以显示按节点编号排列的坐标；键入e再回车则显示边界数据矩阵（7维数组）；键入t按回车则显示三角形单元数据矩阵（4维数组）。
点、边、单元的部分输出为:
P二
Columns1through11
-1.00000.00001.00000.0000一
0.70710.7071
0.7071
-0.7071
-0.9808
-0.9239-0・8315
-0.0000-1.0000
0
1.0000-o.7071-
0.70710.7071
0.7071
-0・1951
-0.3827
-0.5556
e二
Columns1through
11
1.0000
9.000010.0000
11.0000
5.0000
12.0000
13.0
00014.0000
2.000015.0000
16.0000
9.0000
10.000011.0000
5.0000
12.0000
13.0000
14.000
0
2.000015.000016.000017.0000
00.12500.25000.37500.5000
0.6250
0.75000.8750
0
0.1250
0.2500
0.1250
0.2500
0.3750
0.5000
0.6250
0.7500
0.8750
1.0000
0.1250
0.2500
0.3750
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
2.0000
2.0000
2.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.
0000
0
0
0
0
0
0
0
0
0
0
0
t-
Columns1
through18
3214
20
26
29
17
23
100
11
6689
19
94
5
12
13
97
12
34
8
67
28
5
32
19
10
512
13
14
2
8997
81
98
84
92
99
127
94
8911911995
118118
90
70
126
11
11
1
11
1
1
1
11
11
1
1
1
1
（7）输出近似解数值
单击Solve菜单中ExportSolution。
。
。
选项，在Export对话框中输入u再单击OK按钮，再在MATLAB命令行中输入u并回车，就会显示按节点编号排列的解u的数值。
（8）近似解和准确解的比较
方程的准确解为：
1—x2—V2
飒x,y）=—厂・
为了与准确解比较，单击Plot菜单中Parameters…选项，打开PlotSelection对话框，在Height（3~Dplot）行的Property下拉框中选UserEntry,并且输入
u-（1-x."2-y.⑵/4,单击Plot按钮，就可以看到误差曲面，其数量级为
104
练习1用pdetool11具求解课本P128的第2、3、4这儿题的解，并作岀图形。
练习2用pdetool工具求解课本P418的第2、4这两题的解（三角单单元形状不限），并作出图形。
升级会员 