偏微分方程边值问题的数值解法论文.docx

偏微分方程边值问题的数值解法论文.docx

《偏微分方程边值问题的数值解法论文.docx》由会员分享，可在线阅读，更多相关《偏微分方程边值问题的数值解法论文.docx（11页珍藏版）》请在冰点文库上搜索。

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

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

本实验学习使用MATLAB的图形用户命令pdetool来求解偏微分方程的边值问题。

这个工具是用有限元方法来求解的，而且采用三角元。

我们用个例题来说明它的用法。

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

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

其中

未知函数为

。

它的边界条件分为三类：

（1）Direchlet条件：

（2）Neumann条件：

（3）混合边界条件：

在边界

上部分为Direchlet条件，另外部分为Neumann条件。

其中

是定义在边界

的已知函数，另外

也可以是一个2*2的函数矩阵，

是沿边界的外法线的单位向量。

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

二、例题

例题1用pdetool求解

解：

首先在MATLAB的工作命令行中键入pdetool，按回牟键确定，于是出现PDEToolbox窗口，选GenenicScalar模式．

（l）画区域圆

单击椭圆工具按钮，大致在（0，0）位置单击鼠标右键，拖拉鼠标到适当位置松开。

为了保证所绘制的圆是标准的单位园，在所绘园上双击，打开ObjectDialog对话框，精确地输入圆心坐标X-center为0、Y-center为0及半径Radius为l，然后单击OK按钮，这样单位画已画好．

（2）设置边界条件

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

，然后单击OK按钮．也可以单击Boundary菜单中SpocifyBoundaryCondition…选项，打开BoundaryCondition对话框输入边界条件．

（3）设置方程

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

（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（三角形）数据变量，单击OK按钮，然后在MATLAB命令行键入p,即可以显示按节点编号排列的坐标；键入e再回车则显示边界数据矩阵（7维数组）；键入t按回车则显示三角形单元数据矩阵（4维数组）。

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

p=

Columns1through11

-1.00000.00001.00000.0000-0.70710.70710.7071-0.7071-0.9808-0.9239-0.8315

-0.0000-1.000001.0000-0.7071-0.70710.70710.7071-0.1951-0.3827-0.5556

e=

Columns1through11

1.00009.000010.000011.00005.000012.000013.000014.00002.000015.000016.0000

9.000010.000011.00005.000012.000013.000014.00002.000015.000016.000017.0000

00.12500.25000.37500.50000.62500.75000.875000.12500.2500

0.12500.25000.37500.50000.62500.75000.87501.00000.12500.25000.3750

1.00001.00001.00001.00001.00001.00001.00001.00002.00002.00002.0000

1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000

00000000000

t=

Columns1through18

3214202629172310011668919945121397

123486728532191051213142

899781988492991279489119119951181189070126

111111111111111111

（7）输出近似解数值

单击Solve菜单中ExportSolution。

。

。

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

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

方程的准确解为：

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

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

。

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

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