偏微分方程数值解

以中心方式画椭圆/圆（Ctrl+鼠标）。 画多边形，按右键可封闭多边形。 进入边界模式。

打开PDE Specification（偏微分方程类型）对话框。 初始化三角形网格。 加密三角形网格。 解偏微分方程。

打开Plot Selection对话框，确定后给出解的三维图形。 为显示缩放切换按钮。

第三章 典型方程及其应用实例

求解PDE问题主要有两种方法，一种是使用图形用户界面，另一种是采用命令行编程。前者直观简便，而后者更为灵活。 2．1 求解椭圆方程的例子

例：单位圆上的Poisson方程边值问题：

22????u?1,????x,y?|x?y?1???u|???0 ?

这一问题的精确解为：

u?x,y?1?x??2?y2?4.

若使用图形用户界面（Graphical User Interface,简记为GUI），则首先在MATLAB的工作窗口中键入pdetool，按回车键确定，于是出现PDE Toolbox窗口。如果需要坐标网格，单击Options菜单下的Grid选项即可。下面分步进行操作。

（i）画区域圆 单击工具

，大致在（0，0）位置单击鼠标右键同时拖拉鼠标到适当位置松开，

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

（ii）设置边界条件 单击工具，图形边界变红，逐段双击边界，打开Boundary 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。

图2-1

图2-2

（iv）网格剖分 单击工具

，或者单击Mesh菜单中Initialize Mesh选项，可进行初始网格剖分，

这时在PDE Toolbox窗口下方的状态栏内显示初始问网格的节点数和三角形单元数。本题节点数为144个，三角形单元数为254个。如果需要网格加密，再单击

，或者单击Mesh菜单中Refine Mesh选项，这

时节点数变为541个，三角形单元数为1016个，如此还可继续加密。

（v）解方程 单击工具

，或者单击Solve菜单中Solve菜单中Solve PDE选项，可显示方程色

彩解。如果单击Plot菜单中Parameters?选项，出现Plot 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 plot）行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（三角形）数据的变量，单击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); 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]; End

Pdemesh(p,e,t) Pdesurf(p,t,u) Pdesurf(p,t,u-exact)

通过命令行键入help+命令函数，如help pdemesh，按回车键，可以调入有关命令函数的定义、参数格式等帮助信息。

2．2 求解抛物型方程的例子

例：考虑一个带有矩形孔的金属板上的热传导问题。板的左边保持在100c，板的右边热量从板向环境空气定常流动，其他边及内孔边界保持绝缘。初始

?t?t0时板的温度为0，于是概括为如下定解问题：