偏微分方程数值解

第一章 概 述

大家采用下面的方法求解Terzaghi一维固结方程。 1．1 偏微分方程工具箱的功能

偏微分方程工具箱(PDE Toolbox)提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境。PDE Toolbox的功能包括：

(1) 设置PDE (偏微分方程)定解问题，即设置二维定解区域、边界条件以及方程的形式和系数； (2) 用有限元法 (FEM) 求解PDE数值解； (3) 解的可视化。

无论是高级研究人员还是初学者，在使用PDE Too1box时都会感到非常方便。只要PDE定解问题的提法正确，那么，启动MATLAB后，在MATLAB工作空间的命令行中键人pdetool，系统立即产生偏微分方程工具箱(PDE Toolbox)的图形用户界面(Graphical User Interface，简记为GUI)，即PDE解的图形环境，这时就可以在它上面画出定解区域、设置方程和边界条件、作网格剖分、求解、作图等工作，详见1．4节中的例子。我们将在第二章详细介绍GUI的使用，在第二章给出大量典型例子和应用实例。除了用GUI求解PDE外，也可以用M文件的编程计算更为复杂的问题，详见第三章和第四章 的内容。

1．2 PDE Toolbox求解的问题及其背景 1．2．1 方程类型

PDE Toolbox求解的基本方程有椭圆型方程、抛物型方程、双曲型方程、特征值方程、椭圆型方程组以及非线性椭圆型方程。

椭圆型方程： ???(c?u)?au?f, in ?,， 椭圆型方程：???(c?u)?au?f,in?,

其中?是平面有界区域，c,a,f以及未知数u是定义在?上的实（或复）函数。 抛物型方程：d?u???(c?u)?au?f, in ?. ?t?2u双曲型方程：?2???(c?u)?au?f, in ?.

?t特征值方程：??(c?u)?au??du, in ?,

其中d是定义在?上的复函数，?是待求特征值。在抛物型方程和双曲型方程中，系数c,a,f和d可以依赖于时间t。

可以求解非线性椭圆型方程：

????c(u)?u??a(u)?f(u), in ?, 其中c,a,f可以是未知函数u的函数。

还可以求解如下PDE方程组；

??????c11(u)?u1?????c12(u)?u2??a11u1?a12u2?f1,? ??????c21(u)?u1?????c22(u)?u2??a21u1?a22u2?f1

利用命令行可以求解高阶方程组。对于椭圆型方程，可以用自适应网格算法，还能与非线性解结合起

来使用。

另外，对于Poission方程还有一个矩形网格的快速求解器。 1．2．2 边界条件

（1）Dirichlet条件 ： hu?r

( 2 ) Neumann 条件： n?(c?u)?qu?g

其中n是?的边界??上的单位外法向量，g,q,h和r是定义在??上的函数。对于特征值问题仅限于齐次条件：g?0,和r?0。对于非线性情形．系数g,q,h和r可以依赖于u；对于抛物型方程和双曲型方程，系数可以依赖于时间t。

对于方程组情形，边界条件为

( 1 ) Dirichlet 条件： h11u1?h12u?2 r h21u1?h22u?2 r( 2 ) Neumann 条件： n?(c11?u1)?n?(c12?u2)?q11u1?q12u2?g1

??????n?(c21?u1)?n?(c22?u2)?q21u1?q22u2?g2

( 3 ) 混合边界条件为: h11u1?h12u2?r1

n?(c11?u1)?n?(c12?u2)?q11u1?q12u2?g1?h11? n?(c21?u1)?n?(c22?u2)?q21u1?q22u2?g2?h12?

其中?的计算要使得Dirichlet条件满足。在有限元法中，Dirichlet条件也称为本质边界条件，Neumann条件称为自然边界条件。 1.3 如何使用FDE Toolbox 1.3.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模式：对于椭圆型方程还能求非线性和自适应解。对于抛物型和双曲型力程．设置初始边值条件后能求出给定t时刻的解。对于特征值问题，能求出给定区间内的特征值；求解后可以加密网格再求解。

1.3.3 使用Toolbox求解非标准的问题

对于非标准的问题。可以用PDE Too1box的函数。或者用FEM(有限元法)求解更为复杂的问题。 1.3.4 计算结果的可视化

从GUI能够使用Plot模式实现可视化。可以使用Color, Height和Vector等作图。对于抛物型和双曲型方程，还可以生成解的动画。这些操作通过命令行都很容易实现。

????

1.3.5 应用领域

在应用界面提供了丁如下应用领域 ．结构力学——平面应力问题 ．结构力学——平面应变问题 ．静电场问题 ．静磁场问题 ．交流电磁场问题 ．直流导体介质问题 ．热传导问题 ．9‘散问题

这些界面都有对话框，它包括PDE的系数、边界条件、解的性质等。 1. 4 解偏微分方程的一个例子

解Poisson方程??u?f，边界条件为齐次Dirichlet类型。

第一步：启动MATLABl, 键入pdetool，按回车键确定便可启动GUI，然后在Options菜单下选择Grid命令，打开栅格, 栅格的使用，能使用户容易确定所绘图形的大小，如图1—1

1--1

第二步：分步完成平面几何造型：R1-C1-E1+R2+C2。用菜单或快捷工具，分别画矩形R1、矩形R2、椭圆E1、圆C1、圆C2。画圆时，首先选中椭圆工具，按鼠标右键并拖动即可、或者在按ctrI的同时，拖动鼠标也可绘制圆。然后在Set formula栏，进行编辑并用算术运算将将图形对象名称连接起来，删除默认的表达式键入R1-C1-E1+R2+C2，按等号健得到所需图形。若需要，还可进行储存．

形成M文件。

选择Boundary菜单中Boundary Mode命令，进入边界模式。单击Boundary菜单中Remove A11 Subdomain Borders选项，去除子域边界。如果想将几何信息和边界信息进行存储，应选择Boundary菜单中的ExPort Decomposed Geometry．Boundary Cond’s…命令，将它们分别储存于g,b变量中, 通过MATLAB形成M文件。

第三步：选取边界．单击Boundary菜单中Specify Bounddy Conditions…选项，打开Boundary conditlons对话框，输入边界条件，如图1—4。本例取缺省条件。即将全部边界设为齐次Dirichlet条件，边界颜色显示为红色。

第四步：选择PDE菜单中PDE Mode命令，进入PDE模式。单击PDE菜单中PDE Specification…选项，打开PDE对话框，设置方程类型。本例取缺省设置，类型为椭圆型，参数c,a,f分别为1，0，10。

第五步：选择Mesh菜单中Initialize Mesh命令，进行网格剖分。 第六步：选择Mesh菜单中Refine Mesh命令，对网格加密。

第七步：选择Solve菜单中So1ve PDE命令，解偏微分方程并显示图形解。

第八步：单击Plot菜单中Parameters?选项，打开Plot selection对话框，选中Color, Height (3—D Plot)和Show mesh三项。然后单击Plot按钮，显示三维图形解。

第九步：如果要画等值线图和矢量场图，单击Plot菜单中Parameters?选项，打开Plot Selection对话框．选中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文件中。