第九章偏微分方程差分方法.doc

|第 9 章 偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。9.1 椭圆型方程边值问题的差分方法9.1.1 差分方程的建立最典型的椭圆型方程是 Poisson（泊松）方程（9.1）Gyxfyuxu),(,)(2G 是 x,y 平面上的有界区域，其边界 为分段光滑的闭曲线。当 f(x,y)0 时，方程（9.1）称为 Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件第一边值条件 （9.2）),(yxu第二边值条件 （9.3）,n第三边值条件 （9.4）),()(yxku这里，n 表示 上单位外法向， (x,y), (x,y), (x,y)和 k(x,y)都是已知的函数，k(x,y)0 。满足方程（9.1）和上述三种边值条件之一的光滑函数 u(x,y)称为椭圆型方程边值问题的解。用差分方法求解偏微分方程，就是要求出精确解 u(x,y)在区域 G 的一些离散节点（x i，y i）上的近似值 ui,j( xi，y i） 。差分方法的基本思想是，对求解区域 G 做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。设 G=00, B(x,y) B min >0, E(x,y) 0。引进半节点,121hii利用一阶中心差商公式，在节点（i,j ）处可有221yii| )(2,1(),(), )(),1(),),(,12121)(),( 21 21, 2hOhjiujijixu hOhjiujiAjijiAxuxuhjixu jiji 对 类似处理，就可推得求解方程（9.9）的差分方程yB),(（9.10）h jijijijijijiiGfuauaua ),(, ,1,1,1,1其中（9.11） jijijijijiji jijiji jijiji jijiji jijiji EBhAhaDBhCAa ,21,2,21,21, ,2, ,1,21, ,1, ,21, )()()()(显然，当系数函数 A(x,y)=B(x,y)=1, C(x,y)=D(x,y)=E(x,y)=0 时，椭圆型方程（9.9）就成为 Poisson 方程（9.1） ，而差分方程（9.10）就成为差分方程（9.6） 。容易看出，差分方程（9.10）的截断误差为 阶。)21hO9.1.2 一般区域的边界条件处理前面已假设 G 为矩形区域，现在考虑 G 为一般区域情形，这里主要涉及边界条件的处理。考虑 Poisson 方程第一边值问题（9.12）),(,yxuf其中 G 可为平面上一般区域，例如为曲边区域。仍然用两组平行直线：x=x0+ih1,y=y0+jh2,i,j=0,±1,对区域 G 进行矩形网格剖分，见图 9-3。|如果一个内节点（i,j）的四个相邻节点（i+1,j ） ， （i-1, j） ， （i,j+1）和（i,j -1）属于 ，则称其为正则内点，见图 9-3 中打“。 ”号者；如果一个节点G（i,j）属于 且不为正则内点，则称其为非正则内点，见图 9-3 中打“.”号者。记正则内点集合为 ，非正则内点集合为 。显然，当 G 为矩形区域时，h h成立。h,在正则内点（i,j）处，完全同矩形区域情形，可建立五点差分格式hjijijijijijiji Gfuuhuh ),(2121 ,1,2,2（9.13）在方程（9.13）中，当（i,j）点临近边界时，将出现非正则内点上的未知量，因此必须补充非正则内点处的方程。若非正则内点恰好是边界点，如图 9-4 中D 点，则利用边界条件可取 uD= (D)对于不是边界点的非正则内点，如图 9-4 中B 点，一般可采用如下两种处理方法。a.直接转移法.取与点 B 距离最近的边界点（如图 9-4 中 E 点）上的 u 的值作为 u(B)的近似值 uB，即 uB=u(E)= (E)直接转移法的优点是简单易行，但精度较低，只为一阶近似。b.线性插值法.取 B 点的两个相邻点（如图 9-4 中边界点 A 和正则内点 C 作为插值节点对 u(B)进行线性插值| )()()( 21hOCuxAuxBuABCB则得到点 B 处的方程ABCxhh,)(11线性插值法精度较高，为二阶近似。对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般区域上边值问题（9.12）的差分近似解。对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。第二、三边值条件的处理较为复杂，这里不再讨论。9.2 抛物型方程的差分方法本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。9.2.1 一维问题作为模型，考虑一维热传导的初边值问题（9.14）Ttlxtfxuat 0,),(2（9.15）l0,),(（9.16）ttgtutgtu0),(021其中 a 是正常数， 都是已知的连续的函数。),(1xf和现在讨论求解问题（9.14）-(9.18)的差分方法。首先对求解区域 G=0xl, 0t T进行网格剖分。取空间步长 h=l/N,时间步长 = T/M,其中 N,M 是正整数，作两族平行直线Mktjxkj ,10,将区域 G 剖分成矩形网格，见图 9-5，网格交点（x j，t k）称为节点。|用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解 u(x,t)在每个节点（x j， tk）处的近似值 。为简化记号，简记节点（x j，t k）),(kjkjtu=u(j,k)。 利用一元函数的 Taylor 展开公式，可推出下列差商表达式（9.17）)(,()1,(),( Ojjjt （9.18）kuk（9.19）)(2,(),(),( 2jjjtu（9.20）)(,1),1, 22 hOhkjukjjkjx 1.古典显格式在区域 G 的内节点(j,k)处，利用公式（9.17）和（9.20） ，可将偏微分方程（9.14）离散为 )(),1(),2),1(1, 2hfhkjujkjuau j 其中 。舍去高阶小项 ，就得到节点近似值（差分解）),(kikjtxfO所满足的差分方程ju（9.21）kjjkjkjkjj fhuau2111显然，在节点(j,k )处，差分方程（9.21）逼近偏微分方程（9.14）的误差为，这个误差称为 截断误差 ，它反映了差分方程逼近偏微分方程的精度。2hO|现将（9.21）式改写为便于计算的形式，并利用初边值条件（9.15）与（9.16）补充上初始值和边界点方程，则得到（9.22） MktgutguNjxj frurkkj kjjkjkjk ,10),(),(2,1,)(10 11其中 称为网比。2har与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程（9.22） ，当第 k 层节点值 已知时，可直接计算出第 k+1 层节点值 。这kju 1kju样，从第 0 层已知值 开始，就可逐层求出各时间层的节点值。差分方程)(0ijx（9.22）的求解计算是显式的，无须求解方程组，故称为古典显格式。此外，在式（9.22）中，每个内节点处方程仅涉及 k 和 k+1 两层节点值，称这样的差分格式为双层格式。差分方程（9.22）可表示为矩阵形式（9.23）01 1,0,uMFAkk 其中rrrA2121 TNTkkNkkkxtgfftrgfFu)(,)( )(,11 212211 2. 古典隐格式在区域 G 的内节点(j,k)处，利用公式（9.18）和（9.20） ，可将偏微分方程（9.14）离散为 )(),1(),2),1( 2hOfhkjujkjuau j |舍去高阶小项 ，则得到如下差分方程)(2hO（9.24）kjjkjkjkjj fuau2111它的截断误差为 ，逼近精度与古典显格式相同。改写（9.24）式为便于)(h计算的形式，并补充上初始值与边界点方程，则得到（9.25）MktgutguNjxj furrkkj kjjkjjk ,10),(),(2,1,)(10 11与古典显格式不同，在差分方程（9.25）的求解中，当第 k-1 层值 已知时，1kju必须通过求解一个线性方程组才能求出第 k 层值 ，所以称（9.25）式为古典ju隐格式，它也是双层格式。差分方程（9.25）的矩阵形式为（9.26）01,2,uMkFBk其中rrrB2121 向量 同（9.23）式中定义。从（9.26）式看到，古典隐格式在每一层计,kFu算时，都需求解一个三对角形线性方程组，可采用追赶法求解。3.Crank-Nicolson 格式（六点对称格式）利用一元函数 Taylor 展开公式可得到如下等式|)()1,(),(21),(, 222 2Okjxujkjxujkjt 使用这两个公式，在 点离散偏微分方程（9.14） ，然后利用（9.20）式进,j一步离散二阶偏导数，则可导出差分方程（9.27）21211211 kjkjjkjjkjkjkjj fhuhuuau其截断误差为 ，在时间方向的逼近阶较显格式和隐格式高出一阶。这)(2O个差分格式称为 Crank-Nicolson 格式，有时也称为六点对称格式，它显然是双层隐式格式。改写（9.27）式，并补充初始值和边界点方程得到（9.28）MktgutguNjxj frruukkj kjkjjkjjj ,10),(),(2,1,)(10 2111它的矩阵形式为 （9.29）1,0, ,2,)()(0 21kuFAIBIk 在每层计算时，仍需求解一个三对角形方程组。4. Richardson 格式利用公式（9.19）和（9.20） ，可导出另一个截断误差为 阶的差分)(2hO方程 kjjkjkjkjj fhuau21112称之为 Richardson 格式。可改写为(9.32)kjkjjkjkjj fr)(111