修订版偏微分方程差分方法汇总修订版

1、第9章偏微分方程的差分方法含有偏导数的微分方程称为偏微分方程。由于变量的增多和区域的复杂性，求偏微分方程的精确解一般是不可能的，经常采用数值方法求方程的近似解。偏微分方程的数值方法种类较多，最常用的方法是差分方法。差分方法具有格式简单，程序易于实现，计算量小等优点，特别适合于规则区域上偏微分方程的近似求解。本章将以一些典型的偏微分方程为例，介绍差分方法的基本原理和具体实现方法。9.1 椭圆型方程边值问题的差分方法9.1.1 差分方程的建立最典型的椭圆型方程是Poisson(泊松)方程：2u：2u一一Au三一(r+:)=f(x,y),(x,y)eG(9):xcyG是x,y平面上的有界区域，其边界

2、r为分段光滑的闭曲线。当f(x,y)三0时，方程(9.1)称为Laplace(拉普拉斯)方程。椭圆型方程的定解条件主要有如下三种边界条件第一边值条件ur=a(x,y)(9.2)第二边值条件r=P(x,y)(9.3)二n第三边值条件("+ku)=¥(x,y)(9.4)：:n这里，n表示r上单位外法向，a(x,y),3(x,y),丫(x,y)和k(x,y)都是已知的函数，k(x,y)>0O满足方程(9.1)和上述三种边值条件之一的光滑函数u(x,y)称为椭圆型方程边值问题的解。用差分方法求解偏微分方程，就是要求出精确解u(x,y)在区域G的一些离散节点(xy。上的近似值口

3、一"(为，y)。差分方法的基本思想是，对求解区域G做网格剖分，将偏微分方程在网格节点上离散化，导出精确解在网格节点上近似值所满足的差分方程，最终通过求解差分方程，通常为一个线性方程组，得到精确解在离散节点上的近似值。设G=0<x<a,0<y<b为矩形区域，在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,yO称为剖分节点(区域内节点集合记为Gh=(X,y);(xi,y。CG),网格线与边界r的交点称为边界

4、点，边界点集合记为现在将微分方程（9.1）在每一个内节点（Xi,yi）上进行离散。在节点（Xi,yi）处，方程（9.1）为22-u-U咨（Xi,yi）十日（Xi,yi）=f（Xi,yi）,（Xi,yi）WGh（9.5）:x:y需进一步离散（9.5）中的二阶偏导数。为简化记号，简记节点（Xi,yi）=（i,j）,节点函数值u（Xi,yi）=u（i,j）o利用一元函数的Taylor展开公式，推得二阶偏导数的差商表达式.2/W（i,j）=3u（i1.j）-2u（i.j）u（i-1,j）0（h12）x'lu12F(i,j)Fu(i,j"MJ)u(i,i0(h2)代入（9.5）式中，得

5、到方程（9.1）在节点（i,j）处的离散形式11-2u(i1,j)-2u(i,j)u(i-1,j)-yu(i,j1)-2u(i,j)u(i,j-1)h1h20(h12h2),(i,j)Gh其中=f(为，v)。舍去高阶小项0(h2+h；),就导出了u(i,j)的近似值ui,j所满足的差分方程11._；Tui+,j-2ui,j+ui4,j-ui,-2ui,j+ui,jU=fi,j,(i,j)Gh(9.6)h1h2在节点（i,j）处方程（9.6）逼近偏微分方程（9.1）的误差为O（h2+h；）,它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。在差分方

6、程（9.6）中，每一个节点（i,j）处的方程仅涉及五个节点未知量u，j,ui+1,j,Ui-1,j,Ui,j+1,Ui,j-1,因此通常称（9.6）式为五点差分格式，当h1=h2=h时，它简化为h2ui1,jUi,jui,j1Ui,j-4Ui,j=fi,j,(i,j)Gh差分方程（9.6）中，方程个数等于内节点总数，但未知量除内节点值Ui,j,（i,j）£Gh外，还包括边界点值。例如，点（1,j）处方程就含有边界点未知量U0,j。因此，还要利用给定的边值条件补充上边界点未知量的方程。对于第一边值条件式（9.2）,可直接取Ui,j=a（Xi,yi）,（i,j）Th（9.7）对于第三（k

7、=0时为第二）边值条件式（9.4）,以左边界点（1,j）为例，见图9-2,利用一阶差商公式«肛口cuu（Qj）u（1,j）<9j）=h+OSi）图9-一nh1-则得到边界点（0,j）处的差分方程u0,j-也h1-k0,ju0,j=r0,j(9.8)联立差分方程（9.6）与（9.7）或（9.8）就形成了求解Poisson方程边值问题的差分方程组，它实质上是一个关于未知量u”的线性代数方程组，可采用第2,3章介绍的方法进行求解。这个方程组的解就称为偏微分方程的差分近似解，简称差分解。考虑更一般形式的二阶椭圆型方程uu、c：uu-(A-)+丁(Bt)+Ct+D丁+Eu=f(x,y),

8、(x,y产G(9.9)二x二x二y二y二x二y1.其中A(x,y)>Amn>0,B(x,y)>Bmn>0,E(x,y)>0。引进半节点x1=xi±-h1i二22，1y1=yi士一h2利用一阶中心差商公式，在下点（i,j）处可有V2u1：u1：u1-2一(A一)(i,j)(A一)(i;，j)-(A一)(i：,j)O(h2)x二x%二x20x2J-1u(i1,j)-u(i,j)_A1u(i,j)-u(i-1,j)O(h12)ni2,jh1i-2,jh1卫(u(i-7-1)2),x2h1:u:u.对（B）,类似处理，就可推得求解方程（9.9）的差分方程.y二y

9、二y-aii,jui,j.ai,jui,j-ai,jdui,jiai,jdui,ji-ai,jui,ji(9.iU)=f(i,j),(i,j)Gh其中ai1,j2=h(Aii2ai,jai,jH=hi(AiiJ=h21B.ii,j-2hi:Ci'j)hi-*)ai,j1ai,j=hi(Ai.i一i,j-2i.A2，jih22h22Di,j)Di,j)、_2，一，一、一.i)h(B.iBi)Ei,ji,ji,j-i,j222(9.11)显然，当系数函数A（x,y）=B（x,y）=i,C（x,y）=D（x,y）=E（x,y）=0时，椭圆型方程（9.9）就成为Poisson方程（9.i）,而

10、差分方程（9.i0）就成为差分方程（9.6）。容易看出，差分方程（9.i0）的截断误差为O（h；+h；）阶。9.i.2一般区域的边界条件处理前面已假设G为矩形区域，现在考虑G为一般区域情形，这里主要涉及边界条件的处理。考虑Poisson方程第一边值问题(9.i2)M=f(x,y),(x,y)wGU=s(x,y),(x,y)w其中G可为平面上一般区域，例如为曲边区域。仍然用两组平行直线：x=xo+ihi,y=yo+jh2,i,j=0,±i,，对区域G进行矩形网格剖分,见图9-3。如果一个内节点（i,j）的四个相邻节点（i+i,j）,（i-i,j）,（i,j+i）和（i,j-i）属于G=

11、G，j，则称其为正则内点，见图9-3中打”号者；如果一个节点（i,j）属于G且不为正则内点，则称其为非正则内点，见图9-3中打".”号者。记正则内点集合为G非正则内点集合为"。显然，当G为矩形区域时，Gh=Gh,甲=曙成立。（如图9-4中E点）上的u的值作为Ubh1二(A)6uC,h1、C,"Xb-Xa在正则内点（i,j）处，完全同矩形区域情形，可建立五点差分格式1 .八r1.工卅为-2切/工卅皿Ui”=fi,j'(i'J尸Gh(9.13)在方程（9.13）中，当（i,j）点临近边界时，将出现非正则内点上的未知量，因此必须补充非正则内点处的方程。若

12、非正则内点恰好是边界点，如图9-4中D点，则利用边界条件可取Ud=oc（D）对于不是边界点的非正则内点，如图9-4中点，一般可采用如下两种处理方法。a.直接转移法.取与点B距离最近的边界点u（B）的近似值UB,即Ub=u（E）="（E）直接转移法的优点是简单易行，但精度较低，只为一阶近似。b.线性插值法.取B点的两个相邻点（如图9-4中边界点A和正则内点C作为插值节点对u（B）进行线性插值u(B)Xc-Xbu(A)Xb-Xau(C)O(h12)Xc-XaXc-Xa则得到点B处的方程线性插值法精度较高，为二阶近似。对每一个非正则内点进行上述处理，将所得到的方程与（9.13）式联立，就组

13、成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般区域上边值问题（9.12）的差分近似解。对于一般区域上二阶椭圆型方程（9.9）的第一边值问题，可完全类似处理。第二、三边值条件的处理较为复杂，这里不再讨论。9.2抛物型方程的差分方法本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。9.2.1 一维问题作为模型，考虑一维热传导的初边值问题2*=a-2+f(x,t),0cxMl,0<t<T(9.14).:tjxu(x,0)=*(x),0<x<l(9.15)u（0,t）=g（t）,u（l,t）=g2（t）,0<t<T（9.16

14、）其中a是正常数，f（x,t）,中（x）,g（t）和g2（t）都是已知的连续的函数。现在讨论求解问题（9.14）-（9.18）的差分方法。首先对求解区域G=0WxWl,0<t<T进行网格剖分。取空间步长h=l/N,时间步长t=T/M,其中N,M是正整数，作两族平行直线x=xj=jh,j=0,1,Ntf=k,k=0,1,M将区域G剖分成矩形网格，见图9-5,网格交点（为，t。称为节点。0123_N-1图95用差分方法求解初边值问题（9.14）-（9.16）就是要求出精确解u（x,t）在每k个下点（xj,tk）处的近似值uju（xj,tk）。为简化记号，简记节点（xj,tk）=u（j,

15、k）。利用一元函数的Taylor展开公式，可推出下列差商表达式(9.17)(9.18)(9.19)(9.20)卫u(j,k1)-u(j,k)0()Zt卫此水)川IDo()由卫(j,k)=u(j,k1)-u(j,k-1)0(.2).122_,.、£u(jk)u(j1,k)-2u(j,k)u(j-1,k)0(h2).:x2'h21.古典显格式在区域G的内节点(j,k)处，利用公式(9.17)和(9.20),可将偏微分方程(9.14)离散为u(j,k1)-u(j,k)au(j1,k)-2u(j,k)u(j-1,k)fk岛一h2j其中f；=f(Xi,tk)o舍去高阶小项O(w+h2)

16、,就得到节点近似值(差分解)uk所满足的差分方程k1kkkk(9.21)Uj-ujuj1-2UjUj.k二a：2fjh显然，在节点(j,k)处，差分方程(9.21)逼近偏微分方程(9.14)的误差为0。+h2),这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。现将(9.21)式改写为便于计算的形式，并利用初边值条件(9.15)与(9.16)补充上初始值和边界点方程，则得到k+k4/c、kkJuj=ruj由十(12r)uj+ruj十%j=1,2,，N1,k=0,1，M-1/、,'',(9.22uj=5(Xj),j=1,2,，N-1kkuo=g1(tk),un=g2(t

17、k),k=0,1,m其中r=a=称为网比。h2与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程kk1(9.22),当第k层节点值山已知时，可直接计算出第k+1层节点值uj。这样，从第0层已知值u0=邛(为)开始，就可逐层求出各时间层的节点值。差分方程(9.22)的求解计算是显式的，无须求解方程组，故称为古典显格式。此外，在式(9.22)中，每个内节点处方程仅涉及k和k+1两层节点值，称这样的差分格式为双层格式。差分方程（9.22）可表示为矩阵形式'uk*=Auk+Fk,k=0,1，M-1/、Jc（9.23）U0=巾其中1-2rr-r1-2rr+A=,'+，二

18、rr1-2rkkk.kTU二（U1,UN）Fk=（f1k-rg1（tk）,f2k，-2一rg2（tk）T：=r（X1）J；：（XN）T2.古典隐格式在区域G的内节点(j,k)处，利用公式(9.18)和(9.20),可将偏微分方程(9.14)离散为u(j,k)-u(j,k-1)u(j1,k)-2u(j,k)u(j-1,k)k2a2fjO(h)h舍去高阶小项O。+h2),则得到如下差分方程kkkUj-UjUj1二akk_2ujUjJfkh2fj（9.24）它的截断误差为O(T+h2),逼近精度与古典显格式相同。改写(9.24)式为便于计算的形式，并补充上初始值与边界点方程，则得到kkkk上士卜一r

19、uj由十(1十2r)Ujruj=Uj+f,丁飞小-1,k=0,，M(9.25)uj”(Xj),j=1,2,，N-1k.kU0=g1(tk),un=g2(tk),k=0,1,M与古典显格式不同，在差分方程（9.25）的求解中，当第k-1层彳1u：已知时,必须通过求解一个线性方程组才能求出第k层彳！u：,所以称（9.25）式为古典隐格式，它也是双层格式。差分方程（9.25）的矩阵形式为ckk_1,Ikr4Bu=u+F,k=1,2，M/、c'''（9.26）U°二巾其中12rt-r1+2r-r+B=、丁+''-rI-r1+2r向量uk,FkW同（9.

20、23）式中定义。从（9.26）式看到，古典隐格式在每一层计算时，都需求解一个三对角形线性方程组，可采用追赶法求解。3.Crank-Nicolson格式（六点对称格式）利用一元函数Taylor展开公式可得到如下等式k.1）=u（j，k1）-u（j，k）O（2）:t'2-22-212（j，k）=12（j，k）2（j,k1）O（2）x22exex1使用这两个公式，在（j,k+）点离散偏微分方程（9.14）,然后利用（9.20）式进2一步离散二阶偏导数，则可导出差分方程k"1k"1k1kkk11ui1-2uiui4ui1-2uiuiJk-=a1jjL±+_L1J+

21、32（9.27）2h2h2j其截断误差为O。2+h2）,在时间方向的逼近阶较显格式和隐格式高出一阶。这个差分格式称为Crank-Nicolson格式，有时也称为六点对称格式，它显然是双层隐式格式。改写（9.27）式，并补充初始值和边界点方程得到k1k1-ruj12(1r)Ujkk=ruj12(1-r)ujj=1,2,N-1,ru；,2小k=0,1,M-1(9.28)u0=(Xj),j=1,2,N-1-kkU0=g1(tk),un=g2(tk),k=0,1,M它的矩阵形式为1_k1k_k(IB)uk1=(IA)ukF2,k=1,2,MU0=,k-0,1,M-1(9.29)在每层计算时，仍需求解一

22、个三对角形方程组。4.Richardson格式利用公式（9.19）和（9.20）,可导出另一个截断误差为O。2+h2）阶的差分方程k1k1Uj-Uj二a2h2-fjk称之为Richardson格式。可改写为k1Uj=u2r(u：12u：u：)2fjk(9.32)这是一个三层显式差分格式。在逐层计算时，需用到k1一uj和ukj两层值才能得到k+1层彳tuk*。这样，从第0层已知值u0=（Xj）开始，还须补充上第一层值u1,才能逐层计算下去。可采用前述的双层格式计算u1。9.14）的差分格式,除上述四种差分格式外，还可构造出许多逼近偏微分方程（但并不是每个差分格式都是可用的。一个有实用价值的差分格

23、式应具有如下性质:（1）收敛性。对任意固定的节点（Xj,tk）,当剖分步长NhT0时，差分解uk应收敛到精确解U（Xj,tk）。（2）稳定性。当某一时间层计算产生误差时，在以后各层的计算中，这些误差的传播积累是可控制而不是无限增长的。理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。因此，这里主要研究差分格式的稳定性。作为例子，先考查Richardson格式的稳定性。设U；是当计算过程中带有误差时，按Richardson格式（9.30）得到的实际计算值。u：是理论值，误差e：=U；-U：。假定右端项f：的计算是精确的，网比r=1,则e：满足j2jk1k1kkkej=ej+（0书一20+

24、e）（9.31）设前k-1层计算时精确的，误差只是在第k层j0点发生，即e=0,ek0=名岩=0,j,j。则利用（9.31）式可得到误差名的传播情况，见表9-1。表9-1r=1/2时Richardson格式的误差传播j0-4j0-3j0-2j0-1j0jo+1j0+2j0+3j0+4k0000£0000k+1000£-2££000k+200£-4£7£-4££00k+30£-6s17£-24£17£-6s£0k+4£-8£31s-68s

25、89s-68s31s-8££k+5-10S49s-144£273s-388s273s-144£49s-10Sk+671£-260s641£-1096s1311£-1096s641£-260s71£1从表中看出，误差是逐层无限增长的。表中的计算虽然是就网比r=进行的,2实际上对任何r>0都会产生类似现象，所以Richardson格式是不稳定的。利用误差传播图表方法考查差分格式的稳定性虽然直观明了，但只能就具体取定的r值进行，并且也不适用于隐式差分格式。9.2.2差分格式的稳定性前节构造的几种双层差分格

26、式都可以表示为如下的矩阵方程形式kk4k(9.32)u=HuFu°二其中H称为传播矩阵。对于显格式H=A,隐格式H=B-1,六点对称格式H=（I+B）（I+A）。一般的三层格式也可以转化为双层格式。为了讨论方便，设在初始层产生误差总°,且假定右端项Fk的计算是精确的。用uk表示当初始层存在误差/时，由差分格式（9.32）得到的计算解，则k口kku=Hu+F-0A,0u=牛+名记误差向量J=Uk_u;则Sk满足方程7k=H”，k=1,2,1°为初始误差定义9.1称差分格式（9.32）是稳定的，如果对任意初始误差s0在某种范数下满足6klMc卜0|,k>0,0&

27、lt;T<T0其中C为与h,°无关的常数。这个定义表明，当差分格式稳定时，它的误差传播是可控制的。从（9.34）式递推得到kk0T；=H；,0-k-因此，差分格式稳定的充分必要条件是Hk<C,0<k<-定理9.3（稳定性必要条件）差分格式稳定的必要条件是存在与常数M,使谱半径：（H）-1M、.-.*定理9.4（稳定性充分条件）设H为正规矩阵，即HH=Huk满足方程(9.33)(9.34)，误差向量飞(9.35)(9.36)h,t无关的(9.37)H,则(9.37)式也是差分格式稳定的充分条件。下面讨论几种差分格式的稳定性。为便于讨论，引进-0111 01+S='I101I10N-1阶矩阵这个特殊矩阵的特征值为a-j：以=2cos匚，j=1,2,-,N-1(9.38)jN例9-1古典显格式此时H=A=(1-2r)I+rS。利用(9.38)式和三角函数公式，可求得H的特征值为2 ji=1-4rsin(),j=1,2,N-1j2N1为使稳定性条件式(9.39)成立，必须且只须rM2。由于H=A为实对称矩阵，所1以古典显