《偏微分方程概述及运用matlab求解偏微分方程常见问题》.

北京航空航天大学

偏微分方程概述及运用matlab求解微分方

程求解常见问题

姓名徐敏

学号57000211

班级380911班

2011年6月

偏微分方程概述及运用matlab求解偏微分

方程常见问题

徐敏

摘要偏微分方程简介，matlab偏微分方程工具箱应用简介，用这个工具箱解方程的过程是：确定待解的偏微分方程；确定边界条件；确定方程所在域的几何形状；划分有限元；解方程

关键词MATLAB 偏微分方程程序

如果一个微分方程中出现的未知函数只含有一个自变量，这个方程叫做常微分方程，也简称微分方程：如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。

一，偏微分方程概述

偏微分方程是反映有关的未知变量关于时间的导数和关于空间变量的导数之间制约关系的等式。许多领域中的数学模型都可以用偏微分方程来描述，很多重要的物理、力学等学科的基本方程本身就是偏微分方程。早在微积分理论刚形成后不久，人们就开始用偏微分方程来描述、解释或预见各种自然现象，并将所得到的研究方法和研究成果运用于各门科学和工程技术中，不断地取得了显著的成效，显示了偏微分方程对于人类认识自然界基本规律的重要性。逐渐地，以物

理、力学等各门科学中的实际问题为背景的偏微分方程的研究成为传统应用数学中的一个最主要的内容，它直接联系着众多自然现象和实际问题，不断地提出和产生出需要解决的新课题和新方法，不断地促进着许多相关数学分支(如泛函分析、微分几何、计算数学等)的发展，并从它们之中引进许多有力的解决问题的工具。偏微分方程已经成为当代数学中的一个重要的组成部分，是纯粹数学的许多分支和自然科学及工程技术等领域之间的一座重要的桥梁。

在国外，对偏微分方程的应用发展是相当重视的。很多大学和研究单位都有应用偏微分方程的研究集体，并得到国家工业、科学部门及军方、航空航天等方面的大力资助。比如在国际上有重大影响的美国的Courant研究所、法国的信息与自动化国立研究所等都集中了相当多的偏微分方程的研究人员，并把数学模型、数学方法、应用软件及实际应用融为一体，在解决实际课题、推动学科发展及加速培养人才等方面都起了很大的作用。

在我国，偏微分方程的研究起步较晚。但解放后，在党和国家的大力号召和积极支持下，我国偏微分方程的研究工作发展比较迅速，涌现出一批在这一领域中做出杰出工作的数学家，如谷超豪院士、李大潜院士等，并在一些研究方向上达到了国际先进水平。但总体来说，偏微分方程的研究队伍的组织和水平、研究工作的广度和深度与世界先进水平相比还有很大的差距。因此，我们必须继续努力，大力加强应用偏微分方程的研究，逐步缩小与世界先进水平的差距

二，偏微分方程的内容

偏微分方程是什么样的？它包括哪些内容？这里我们可从一个例子的研究加以介绍。

弦振动是一种机械运动，当然机械运动的基本定律是质点力学的 F=ma，但是弦并不是质点，所以质点力学的定律并不适用在弦振动的研究上。然而，如果我们把弦细细地分成若干个极小极小的小段，每一小段抽象地看作是一个质点，这样我们就可以应用质点力学的基本定律了。

弦是指又细又长的弹性物质，比如弦乐器所用的弦就是细长的、柔软的、带有弹性的。演奏的时候，弦总是绷紧着具有一种张力，这种张力大于弦的重量几万倍。当演奏的人用薄片拨动或者用弓在弦上拉动，虽然只因其所接触的一段弦振动，但是由于张力的作用，传播到使整个弦振动起来。

用微分的方法分析可得到弦上一点的位移是这一点所在的位置和时间为自变量的偏微分方程。偏方程又很多种类型，一般包括椭圆型偏微分方程、抛物型偏微分方程、双曲型偏微分方程。上述的例子是弦振动方程，它属于数学物理方程中的波动方程，也就是双曲型偏微分方程。

偏微分方程的解一般有无穷多个，但是解决具体的物理问题的时候，必须从中选取所需要的解，因此，还必须知道附加条件。因为偏微分方程是同一类现象的共同规律的表示式，仅仅知道这种共同规律还不足以掌握和了解具体问题的特殊性，所以

就物理现象来说，各个具体问题的特殊性就在于研究对象所处的特定条件，就是初始条件和边界条件。

拿上面所举的弦振动的例子来说，对于同样的弦的弦乐器，如果一种是以薄片拨动弦，另一种是以弓在弦上拉动，那么它们发出的声音是不同的。原因就是由于“拨动”或“拉动”的那个“初始”时刻的振动情况不同，因此产生后来的振动情况也就不同。

天文学中也有类似情况，如果要通过计算预言天体的运动，必须要知道这些天体的质量，同时除了牛顿定律的一般公式外，还必须知道我们所研究的天体系统的初始状态，就是在某个起始时间，这些天体的分布以及它们的速度。在解决任何数学物理方程的时候，总会有类似的附加条件。

就弦振动来说，弦振动方程只表示弦的内点的力学规律，对弦的端点就不成立，所以在弦的两端必须给出边界条件，也就是考虑研究对象所处的边界上的物理状况。边界条件也叫做边值问题。

当然，客观实际中也还是有“没有初始条件的问题”，如定场问题（静电场、稳定浓度分布、稳定温度分布等），也有“没有边界条件的问题”，如着重研究不靠近两端的那段弦，就抽象的成为无边界的弦了。

在数学上，初始条件和边界条件叫做定解条件。偏微分方程本身是表达同一类物理现象的共性，是作为解决问题的依据；

定解条件却反映出具体问题的个性，它提出了问题的具体情况。方程和定解条件合而为一体，就叫做定解问题。

求偏微分方程的定解问题可以先求出它的通解，然后再用定解条件确定出函数。但是一般来说，在实际中通解是不容易求出的，用定解条件确定函数更是比较困难的。

偏微分方程的解法还可以用分离系数法，也叫做傅立叶级数；还可以用分离变数法，也叫做傅立叶变换或傅立叶积分。分离系数法可以求解有界空间中的定解问题，分离变数法可以求解无界空间的定解问题；也可以用拉普拉斯变换法去求解一维空间的数学物理方程的定解。对方程实行拉普拉斯变换可以转化成常微分方程，而且初始条件也一并考虑到，解出常微分方程后进行反演就可以了。

应该指出，偏微分方程的定解虽然有以上各种解法，但是我们不能忽视由于某些原因有许多定解问题是不能严格解出的，只可以用近似方法求出满足实际需要的近似程度的近似解。

常用的方法有变分法和有限差分法。变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，在数学上是拉普拉斯方程的边值问题，由于求

解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。

随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。从这个角度说，偏微分方程变成了数学的中心。

三，用matlab解偏微分方程

解偏微分方程不是一件轻松的事，但是偏微分方程在自然科学和工程领域中应用很广，因此，我们可以运用matlab这个软件来解决一些常见的偏微分方程。

（一）Matlab偏微分方程工具箱简介。

1，概述。本文只给出该工具箱的函数列表

2，偏微分方程算法函数列表。

adaptmesh 生成自适应网络及偏微分方程的解

assemb 生成边界质量和刚度矩阵

assema 生成积分区域上质量和刚度矩阵

assempde 组成偏微分方程的刚度矩阵及右边

hyperbolic 求解双曲线型偏微分方程

parabolic 求解抛物线型偏微分方程

pdeeig 求解特征型偏微分方程

pdenonlin 求解非线性型微分方程

poisolv 利用矩阵格式快速求解泊松方程3，图形界面函数。

pdecirc 画圆

pdeellip 画椭圆

pdemdlcv 转化为版本1.0式的*.m文件pdepoly 画多边形

pderect 画矩形

pdetool 偏微分方程工具箱的图形用户界面4，几何处理函数。

csgchk 检查几何矩阵的有效性

csgdel 删除接近边界的小区

decsg 将固定的几何区域分解为最小区域initmesh 产生最初的三角形网络jigglemesh 微调区域内的三角形网络poimesh 在矩形区域上产生规则的网络refinemesh 细化三角形网络

wbound 写一个边界描述文件

wgeom 写一个几何描述文件

pdecont 画轮廓图

pdemesh 画偏微分方程的三角形网络pdeplot 画偏微分方程的三角形网络

pdesurf 画表面图命令

5，通用函数 。

pdetriq 三角形单元的品性度量

poiasma 边界点对快速求解泊松方程的“贡献”矩阵

poicalc 规范化的矩阵格式的点索引

（二）Matlab 偏微分方程工具箱应用。

可以用词工具箱求解如椭圆方程，双曲线方程，特征值方程，抛

物线方程。

椭圆型偏微分方程

椭圆型偏微分方程的一般形式为

()(,)div c u au f x t -?+=

其中：若12(,,,,)(,)n u u x x x t u x t ==，u ?为u 的梯度，则其定义为

12

,,,n u u x x x ??????=??????? 散度()div v 的定义为

12

()n div v v x x x ?????=+++ ??????

这样，()div c u ?可以更明确地表示为 1122()n n u u u div c u c c c x x x x x x ???????????????=++

+?? ? ? ???????????

????

若c 为常数，则进一步化简为 22

222212()n div c u c u c u x x x ??????=+++=? ??????

其中，?又称为Laplace 算子。这样椭圆型偏微分方程可以简单地写

为

22

222212(,)n c u au f x t x x x ?????-++++= ?????? 抛物型偏微分方程

抛物型偏微分方程的一般形式为

()(,)u d div c u au f x t t

?-?+=? 根据上面叙述，若c 为常数，则该方程可以更简单地写为

22

222212(,)n u d c u au f x t t x x x ??????-++++= ???????

双曲型偏微分方程

双曲型偏微分方程的一般形式为

22()(,)u d div c u au f x t t

?-?+=? 若c 为常数，则可以将该方程简化为

222

2222212(,)n u d c u au f x t t x x x ??????-++++= ???????

三类方程的直接的区别在于u 对t 的导数的阶次。

若对t 没有求导，可以理解为其值为常数，故称为椭圆型的。

若取u 对时间t 的一阶导数，则与u 对x 的二阶导数直接构成了抛

物线关系，故称为抛物型偏微分方程。

若取u 对时间t 的二阶导数，称其为双曲型偏微分方程。

特征值型偏微分方程

特征值型偏微分方程为

()div c u au du λ-?+=

对常数c 该方程可以化简为

22

222212

(,)n c u au f x t x x x λ?????-++++= ?????? 该方程是椭圆型偏微分方程的特例。

pdetool 的使用：在matlab 命令窗口中键入pdetool 窗口打开进入

工作状态，pdetool 提供两种解方程的方法，一种是通过函数，利用函数可以编程也可以用命令行的方式解方程，两一种是对pdetool 窗口进行交互操作。一般来说，用函数解方程比较繁琐，但是比较灵活：通过窗口交互操作比较简单。解方程的全部过程以及结果都可以输出保存为文本文件。限于文本的篇幅，我们主要介绍交互操作偏微分方程的方法。

1.确定待解的偏微分方程。利用函数assempde 可以对待解的偏微

分方程加以描述。

在交互操作中，为了方便用户，pdetool 把常见问题归结为及各

类型，可以再pdetool 窗口的工具栏上找到选择类型的弹出菜单，这些类型如下：

通用问题

通用系统（二维的偏微分方程组）

结构力学：平面应力

结构力学：平面应变

静电学

静磁学

交流电电磁学

直流电导电介质

热传导

扩散

确定问题类型后，可以再PDE Specification对话框填入c,a,f,d 等系数，这样就确定了待解的偏微分方程。

2. 确定边界条件

用函数assemb可以描述边界条件。

用pdetool提供的边界条件对话框，在对话框里填入g,h,q,r等边界条件。

3. 确定偏微分方程所在Ω域的几何图形

平面上波得散射问题。

按照上面所说的解方程的过程，首先确定带解的偏微分方程。

散射是介质对入射波的反射。假定介质是均匀的，那么入射波在介质中传播的速度是一个常数c。

可以用函数画出Ω域的几何图形。Pdecirc:画圆：pdeellip:画椭圆：pderect:画矩形：pdepoly：画多边形。

无论哪种画法，图形一经画出，pdetool就为这个图形自动取名，并把代表图形的名字放入Set formula窗口，在这个窗口，可以通过+，-图形的名字现在对图形的拓扑运算，以便构造复杂的Ω域几何图形

4. 划分有限元

对Ω域进行有限元的划分函数有，initmesh ：基本划分；

jigglemesh ：采用jiggle 方法划分；refinemesh ：精细划分。

在pdetool 窗口中直接点击划分有限元的按钮划分有限元，划分

的方法与上面的函数想对应。

5.解方程

经过1.1—1.4步骤后就可以解方程。解方程的函数有，

adaptunesh ：解方程的通用函数；poicalc ：矩形有限元快速解椭圆形方程；poisolv ：矩形有限元解椭圆形方程；parabolic ：解抛物线型方程：hyperbolic ：解双曲线型方程。

在pdetool 窗口中直接点击解方程的按钮即可解方程，解方程所

耗费的时间在于有限元划分的多少。

（三）实例

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

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

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

()c u u f α-??+= (0.1)

其中

(1) , (2) a,f D c x y ?????=? ?????是上的已知函数(3)是标量或22的函数方阵 未知函数为(,) (,)u x y x y D ∈。它的边界条件分为三类：

（1）Direchlet 条件：

hu f = (0.2) （2）Neumann 条件：

()n c u qu g ?+= (0.3)

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

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

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

二、例题

例题1 用pdetool 求解 22D 1 D: 10u x y u ??-?=+≤??=?? (0.4)

解 ：首先在MATLAB 的工作命令行中键入pdetool ，按回牟键确定，于是出现PDE Toolbox 窗口，选Genenic Scalar 模式．

( l ）画区域圆

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

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

开Boundary Condition 对话框输入边界条件．

( 3 ）设置方程

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

( 4 ）网格剖分

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

Mesh 菜单中Refine Mesh 选项，节点数成为541 个，三角形单元数为1016 个。

( 5 ）解方程

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