第九章-偏微分方程差分方法

170第9章偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。9.1椭圆型方程边值问题的差分方法9.1.1差分方程的建立最典型的椭圆型方程是Poisson（泊松）方程Gyxyxfyuxuu),(),,()(2222（9.1）G是x,y平面上的有界区域，其边界Γ为分段光滑的闭曲线。当f(x,y)≡0时，方程（9.1）称为Laplace(拉普拉斯）方程。椭圆型方程的定解条件主要有如下三种边界条件第一边值条件),(yxu（9.2）第二边值条件),(yxnu（9.3）第三边值条件),()(yxkunu（9.4）这里，n表示Γ上单位外法向，α(x,y),β(x,y),γ(x,y)和k(x,y)都是已知的函数，k(x,y)≥0。满足方程（9.1）和上述三种边值条件之一的光滑函数u(x,y)称为椭圆型方程边值问题的解。用差分方法求解偏微分方程，就是要求出精确解u(x,y)在区域G的一些离散节点（xi，yi）上的近似值ui,j≈(xi，yi）。差分方法的基本思想是，对求解区域G做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。设G={0xa,0yb}为矩形区域，在x,y平面上用两组平行直线x=ih1,i=0,1,…,N1,h1=a/N1y=jh2,j=0,1,…,N2,h2=b/N2将G剖分为网格区域，见图9-1。h1，h2分别称为x方向和y方向的剖分步长，网格交点（xi，yi）称为剖分节点（区域内节点集合记为Gh={（xi，yi）;（xi，yi）∈G}）,网格线与边界Γ的交点称为边界点，边界点集合记为Γh。171现在将微分方程（9.1）在每一个内节点（xi，yi）上进行离散。在节点（xi，yi）处，方程（9.1）为hiiiiiiiiGyxyxfyxyuyxxu),(),,()],(),([2222（9.5）需进一步离散（9.5）中的二阶偏导数。为简化记号，简记节点（xi，yi）=(i,j)，节点函数值u（xi，yi）=u(i,j)。利用一元函数的Taylor展开公式，推得二阶偏导数的差商表达式)(0)]1,(),(2)1,([1),()(0)],1(),(2),1([1),(222222212122hjiujiujiuhjiyuhjiujiujiuhjixu代入（9.5）式中，得到方程（9.1）在节点(i,j)处的离散形式hjiGjihhfjiujiujiuhjiujiujiuh),(),(0)]1,(),(2)1,([1)],1(),(2),1([12221,2221其中),(,iijiyxff。舍去高阶小项)(02221hh，就导出了u(i,j)的近似值ui,j所满足的差分方程hjijijijijijijiGjifuuuhuuuh),(,]2[1]2[1,1,,1,22,1,,121（9.6）在节点(i,j)处方程（9.6）逼近偏微分方程（9.1）的误差为)(2221hhO，它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。在差分方程（9.6）中，每一个节点(i,j)处的方程仅涉及五个节点未知量ui,j，ui+1,j，ui-1,j，ui,j+1，ui,j-1，因此通常称（9.6）式为五点差分格式，当h1=h2=h时，它简化为172hjijijijijijiGjifuuuuuh),(,]4[1,,1,1,,1,12差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值ui,j,(i,j)∈Gh外，还包括边界点值。例如，点(1,j)处方程就含有边界点未知量u0,j。因此，还要利用给定的边值条件补充上边界点未知量的方程。对于第一边值条件式（9.2），可直接取ui,j=α（xi，yi），(i,j)∈Γh（9.7）对于第三（k=0时为第二）边值条件式（9.4），以左边界点(1,j)为例，见图9-2,利用一阶差商公式)(),1(),0(),0(11hOhjujujnu则得到边界点(0,j)处的差分方程jjjjjrukhuu,0,0,01,1,0（9.8）联立差分方程（9.6）与（9.7）或（9.8）就形成了求解Poisson方程边值问题的差分方程组，它实质上是一个关于未知量{ui,j}的线性代数方程组，可采用第2，3章介绍的方法进行求解。这个方程组的解就称为偏微分方程的差分近似解，简称差分解。考虑更一般形式的二阶椭圆型方程GyxyxfEuyuDxuCyuByxuAx),(),,(])()([（9.9）其中A(x,y)≥Amin0,B(x,y)≥Bmin0,E(x,y)≥0。引进半节点,12121hxxii,22121hyyii利用一阶中心差商公式，在节点（i,j）处可有)(2),1(),1(),()(]),1(),(),(),1([1)()],21)((),21)([(1),)((211211,211,211211hOhjiujiujixuhOhjiujiuAhjiujiuAhhOjixuAjixuAhjixuAxjiji173对yuyuBy),(类似处理，就可推得求解方程（9.9）的差分方程hjijijijijijijijijijiGjijifuauauauaua),(),,(][,,1,1,1,1,,1,1,1,1（9.10）其中jijijijijijijijijijijijijijijijijijiEBBhAAhaDhBhaDhBhaChAhaChAha,21,21,22,21,2121,,221,221,,221,221,,1,2121,1,1,2121,1)()()2()2()2()2(（9.11）显然，当系数函数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）的截断误差为)(2221hhO阶。9.1.2一般区域的边界条件处理前面已假设G为矩形区域，现在考虑G为一般区域情形，这里主要涉及边界条件的处理。考虑Poisson方程第一边值问题),(),,(),(),,(yxyxuGyxyxfu（9.12）其中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）属于GG，则称其为正则内点，见图9-3中打“。”号者；如果一个节点（i,j）属于G且不为正则内点，则称其为非正则内点，见图9-3中打“.”号者。记正则内点集合为hG，非正则内点集合为h。显然，当G为矩形区域时，174hhhhGG,成立。在正则内点（i,j）处，完全同矩形区域情形，可建立五点差分格式hjijijijijijijiGjifuuuhuuuh),(,]2[1]2[1,1,,1,22,1,,121（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)进行线性插值)()()()(21hOCuxxxxAuxxxxBuACABACBC则得到点B处的方程ABCBxxuhAhhu,)(111线性插值法精度较高，为二阶近似。对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般175区域上边值问题（9.12）的差分近似解。对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。第二、三边值条件的处理较为复杂，这里不再讨论。9.2抛物型方程的差分方法本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。9.2.1一维问题作为模型，考虑一维热传导的初边值问题Ttlxtxfxuatu0,0),,(22（9.14）lxxxu0),()0,(（9.15）Tttgtlutgtu0),(),(),(),0(21（9.16）其中a是正常数，)()(),(),,(21tgtgxtxf和都是已知的连续的函数。现在讨论求解问题（9.14）-(9.18)的差分方法。首先对求解区域G={0≤x≤l,0≤t≤T}进行网格剖分。取空间步长h=l/N,时间步长τ=T/M,其中N,M是正整数，作两族平行直线MkkttNjjhxxkj,1,0,,,1,0,将区域G剖分成矩形网格，见图9-5，网格交点（xj，tk）称为节点。用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解u(x,t)在每个节点（xj，tk）处的近似值),(kjkjtxuu。为简化记号，简记节点（xj，tk）=u(j,k)。利用一元函数的Taylor展开公式，可推出下列差商表达式176)(),()1,(),(Okjukjukjtu（9.17）)()1,(),(),(Okjukjukjtu（9.18）)(2)1,()1,(),(2Okjukjukjtu（9.19）)(),1(),(2),1(),(2222hOhkjukjukjukjxu（9.20）1.古典显格式在区域G的内节点(j,k)处，利用公式（9.17）和（9.20），可将偏微分方程（9.14）离散为)(),1(),(2),1(),()1,(22hOfhkjukjukjuakjukjukj其中),(kikjtxff。舍去高阶小项)(2hO，就得到节点近似值（差分解）kju所满足的差分方程kjkjkjkjkjkjfhuuuauu21112（9.21）显然，在节点(j,k)处，差分方程（9.21）逼近偏微分方程（9.14）的误差为)(2hO，这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。现将（9.21）式改写为便于计算的形式，并利用初边值条件（9.15）与（9.16）补充上初始值和边界点方程，则得到MktgutguNjxuMkNjfruurruukkNkkjjkjkjkjkjkj,1,0),(),(1,,2,1),(1,,1,0,1,,2,1)21(2