第九章偏微分方程差分方法汇总

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

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

3、 yj上的近似值Uj,j (xi, yi) o差分方法的基本思想是，对求解区域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,yi)称为剖分节点(区域内节点集合记为Gh=(为,yi);(Xi,y

4、i) G),网格线与边界r的交点称为边界点，边界点集合记为rho*40hla图g-】现在将微分方程（9.1 ）在每一个内节点（Xi, yi）上进行离散。在节点（Xi, yj处，方程（9.1 ）为：2u：2u-芋（Xi,y） U（Xi,yi）二 f （Xi,y） （Xi,yJ Gh（ 9.5 ）:x: y需进一步离散（9.5 ）中的二阶偏导数。为简化记号，简记节点（Xi，yi）=（ i, j），节:：2u2(i, J)X；：2u(i,j) y点函数值u （Xi，y） =u（i,j）。利用一元函数的Taylor展开公式，推得二阶偏导数的 差商表达式=2 u(i 1,j)-2u(i,j) u(i-1

5、,j) 0(h2)=2 u(i,j 1)-2u(i,j) u(i,j-1) 0(hf) h2代入（9.5 ）式中，得到方程（9.1 ）在节点（i,j）处的离散形式11-2u(i 1,j)-2u(i,j) u(i -1,j)-；ju(i,j 1)-2u(i,j) u(i, j -1) h1h2二 fi,j0(h12h|), (i,j) Gh其中=f *, yj。舍去高阶小项0（h2+h；），就导出了 u（i,j）的近似值ui,j所满足的差分方程11 / 、-尹-2叽u7一严十皿 5"心）Gh （ 9'）在节点（i,j）处方程（9.6 ）逼近偏微分方程（9.1 ）的误差为0（叶*

6、 h；），它关于剖分步长是二阶的。这个误差称为差分方程逼近偏微分方程的截断误差，它的大小将影响近似解的精度。在差分方程（9.6 ）中，每一个节点（i,j）处的方程仅涉及五个节点未知量 ui,j， ui+1, j， ui-1, j，ui, j+1，ui, j-1，因此通常称（9.6 ）式为五点差分格式，当h1= h；=h时， 它简化为1出1 un，（i,j）3差分方程（9.6 ）中，方程个数等于内节点总数，但未知量除内节点值u,j ,（ i,j） Gh外，还包括边界点值。例如，点 （1,j）处方程就含有边界点未知量 u°,j。因此，还 要利用给定的边值条件补充上边界点未知量的方程。对于

7、第一边值条件式（9.2），可直接取 对于第三（k=0时为第二）边值条件式（ 以左边界点（1, j）为例，见图9-2, 利用一阶差商公式.:uu（0, j）-u（1, j）（0, j）O（h1）一nh|则得到边界点（0, j）处的差分方程uo,j "j“-k0,ju0,j -r0,jhiUi, j=a (Xi, yi) ,( i,j) r h (9.7 )9.4 ),图42(9.8 )189联立差分方程（9.6 ）与（9.7 ）或（9.8 ）就形成了求解 Poisson方程边值问题 的差分方程组，它实质上是一个关于未知量u,j的线性代数方程组，可采用第2, 3章介绍的方法进行求解。这个

8、方程组的解就称为偏微分方程的差分近似解，简称差分解。考虑更一般形式的二阶椭圆型方程uu u u-(A ) (B ) C D Eu = f(x, y), (x,y) G (9.9 ) :x :x :y :y:y其中 A(x, y) >人听>0, B(x, y) > Bn >0, E( x, y)> 0。引进半节点x 1 =Xj 士丄g电 2 ,1»尹2,利用一阶中心差商公式,在节点i, j）处可有u1；:u1u 12多=)(i,j)飞(n,w)(匕,j) 5)2u(i 1,j)-u(i,j)_A 1 u(i,ju(1,j)O(h2) h1i_2,jdYi,

9、j)3 冒"x2h1O(h12)u - u对-(B' ), '类似处理，就可推得求解方程(9.9 )的差分方程 .y:-y:-y-ai1,juij ' ai j,jUi j,j - ai,j .1Ui,j 1- ai,j jUi,j j - ai,jUi, j二 f(i,j), (i,j) Gh(9.10 )其中2* (A.i 2，j=h-i (A 2 Pji,j七2二 h2，(B.-i,j p_2ai,j = hi (A.+ A. 透 iai 1,jai,j 1ai,jh1叱)2 I丿hi巧 Cij)+ h2 D )2 i,J 丿 -S)2 i,J+B. i

10、) + Ei,jTi，j 2 i，j(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 )的截断误差为 0(h； h；)阶。9.1.2 一般区域的边界条件处理前面已假设 G为矩形区域，现在考虑G为一般区域情形，这里主要涉及边界条 件的处理。考虑Poisson方程第一边值问题(9.12 )”Au = f(x,y),(x,y)w GU=a(x,y),(x,y)EF其中G可为平面上一

11、般区域，例如为曲边区域。仍然用两组平行直线： x=xo+ih1, y=yo+jh2, i, j=0, ± 1,对区域G进行矩形网格剖分，见图9-3。如果一个内节点(i,j)的四个相邻节点(i+1,j), (i-1, j), (i,j+1 )和(i,j-1 ) 属于G二G _丨，则称其为正则内点，见图9-3中打“。”号者；如果一个节点(i,j)属于G且不为正则内点，则称其为非正则内点，见图9-3中打“.”号者。记正则 内点集合为Gh ,非正则内点集合为-h。显然，当G为矩形区域时，Gh 二 Gh,- h 二-h 成立。r91G* 11t)4I19<L T00图43在正则内点（i,

12、j）处，完全同矩形区域情形，可建立五点差分格式Uiijhi1-2Ui,j *丄j 一，2 *,j i -2比)ui,jH fi,j,(i,j) Gh(9.13)在方程（9.13 ）中，当（i, j）点临近边界时，将 出现非正则内点上的未知量，因此必须补充非正 则内点处的方程。若非正则内点恰好是边界点，如图9-4中D点，则利用边界条件可取Uo=a （D）对于不是边界点的非正则内点，如图9-4中B点，一般可采用如下两种处理方法。Dg1'czn图9-4a.直接转移法.取与点B距离最近的边界点（如图9-4中E点）上的u的值作为u（B）的近似值 UB，即 Ub=U（E）= a （E）直接转移法的

13、优点是简单易行，但精度较低，只为一阶近似。b.线性插值法.取B点的两个相邻点（如图 9-4中边界点A和正则内点C作为 插值节点对u（B）进行线性插值u(B) = Xc _Xb u(A)Xb _Xa u(C) - O(h12)Xc _XaXc _Xa则得到点B处的方程UBUc,线性插值法精度较高，为二阶近似。对每一个非正则内点进行上述处理，将所得到的方程与（9.13 ）式联立，就组成了方程个数与未知量个数相一致的线性代数方程组。求解此方程组就可得到一般区域上边值问题(9.12 )的差分近似解。对于一般区域上二阶椭圆型方程(9.9 )的第一边值问题，可完全类似处理。第二、三边值条件的处理较为复杂，

14、这里不再讨论。9.2抛物型方程的差分方法本节介绍抛物型方程的差分方法，重点讨论差分格式的构造和稳定性分析。9.2.1一维问题作为模型，考虑一维热传导的初边值问题f(x,t),(9.14)(9.15)(9.16)u(x,0) C(x),0 ：. x : lu(O,t) =g(t), u(l,t) =g2(t), OEt 汀其中a是正常数,f (x,t),(x), g(t)和g2(t)都是已知的连续的函数。现在讨论求解问题(9.14 ) -(9.18)的差分方法。首先对求解区域G=0 < xw l,0< t< T进行网格剖分。取空间步长h=l/ N,时间步长t =T/M,其中N,

15、M是正整数，作两族平行直线x = Xj = jh, j = 0,1, Nt =tk = k , k = 0,1, M将区域G剖分成矩形网格，见图9-5，网格交点(为,tk)称为节点。图9-5用差分方法求解初边值问题(9.14 ) - ( 9.16 )就是要求出精确解 u(x, t)在每k个节点(Xj,tk)处的近似值Uj： u(Xj,tk)。为简化记号，简记节点(Xj,tk)=u(j, k)。利用一元函数的 Taylor展开公式，可推出下列差商表达式3(j,k)=u(j,k D-uJk) 0(.)Zt-u(j,kHu(j，ku(j，k-1) 0(.) a(j,k)= u(j，k DW1) 0(

16、.2) jt2.g u (j k) _u(j +1,k)_2u(j,k)+u(j _1，k) 土0仆2) _：x2'h2(9.17 )(9.18 )(9.19 )(9.20 )1. 古典显格式在区域G的内节点（j,k）处，利用公式（9.17 ）和（9.20 ），可将偏微分方程（9.14 ） 离散为u（j,k 1） -u（j,k） _au（j 1,k）-2u（j,k） u（j-1,k） . f k . 0（. h2）th2J其中fjk =f（Xi，tk）。舍去高阶小项0（. h2），就得到节点近似值（差分解）uk所k 1 kuj ujkc k 丄 kUj 1 -2Uj Uj4fjk(9.

17、21 )满足的差分方程显然，在节点（j,k）处，差分方程（9.21 ）逼近偏微分方程（9.14 ）的误差为0G "h2），这个误差称为截断误差，它反映了差分方程逼近偏微分方程的精度。现将（9.21 ）式改写为便于计算的形式，并利用初边值条件（9.15 ）与（9.16 ）补充上初始值和边界点方程，则得到'k*k丄"c、k 丄k丄"kuj =ru j+(12r)uj +ruj+f(9.22 )j =1,2,，N -1, k =0,1，,M -1 u；N(Xj), J =1,2，,N -1kkw =g1(tk), uN=g2(tk), k = 0,1Mh2称为网

18、比。与时间相关问题差分方程的求解通常是按时间方向逐层进行的。对于差分方程kk 1（9.22 ），当第k层节点值Uj已知时，可直接计算出第k+1层节点值Uj 。这样，从第0层已知值u0 =（Xi）开始，就可逐层求出各时间层的节点值。差分方程（9.22 ）的求解计算是显式的，无须求解方程组，故称为古典显格式。此外，在式（9.22 ）中，每个内节点处方程仅涉及k和k+1两层节点值，称这样的差分格式为双层格式。差分方程(9.22 )可表示为矩阵形式.k二 AuFk,k =0,1, ,M -1(9.23 )其中1 -2r r_r 1 -2rr+ + +A =- J -+ +rr 1 -2rk / k .

19、 k TU(u1 , ,U N)F k = ( f/ - rg1(tk), f2k，,f,2 代rg 2(tk)T：WxJ,，(Xn)t2. 古典隐格式在区域G的内节点(j, k)处，利用公式(9.18 )和(9.20 )，可将偏微分方程(9.14) 离散为u(j,k)u(j,k1)u(j +1,k)-2u(j,k)+u(j-1,k) k +CF +以、a2fj 0( h )th舍去高阶小项OC，h2)，则得到如下差分方程k k 4 Uj -Ujkku j 1 -2u jh2(9.24)它的截断误差为 0(. h2)，逼近精度与古典显格式相同。改写(9.24 )式为便于计 算的形式，并补充上初

20、始值与边界点方程，则得到k丄"丄kkJ 丄ru j 卅+(1 十2r)uj - ru j = u j 中 fj(9.25 )j =1,2,，N -1, k=0,1，M u0"(Xj), j =1,2,，N1kk卩0 =g1(tk),un =g2(tk), k=0,1，M与古典显格式不同，在差分方程（9.25 ）的求解中，当第 k-1层值u：已知时,必须通过求解一个线性方程组才能求出第k层值u：，所以称（9.25 ）式为古典隐 格式，它也是双层格式。差分方程（9.25 ）的矩阵形式为(9.26)Bu： =u：+F：, k =1,2，M其中1 2r -r -r 1 +2r -r

21、+ +B=J-+J- "-r-r 1+2r向量U：,F：,同（9.23 ）式中定义。从（9.26 ）式看到，古典隐格式在每一层计算 时，都需求解一个三对角形线性方程组，可采用追赶法求解。3. Cra nk-Nicolson 格式（六点对称格式）利用一元函数Taylor展开公式可得到如下等式詐 k»u(j,k"u(j,k)20(.),2:uC 2x(j,k1 12"2;=2u§2(j,k)U(j,k 1) 0(.2)ex使用这两个公式，在1（J,k第）点离散偏微分万程（9.14），然后利用（9.20）式进k 1 kUj UjT1k 1 k 11川

22、_2山Uj 4 山1_2山二a_22h2h2h2k¥fj 2(9.27 )步离散二阶偏导数，则可导出差分方程其截断误差为0（ 2h2），在时间方向的逼近阶较显格式和隐格式高出一阶。这个差分格式称为Crank-Nicolson格式，有时也称为六点对称格式，它显然是双层隐式格式。改写（9.27 ）式，并补充初始值和边界点方程得到心，、k*kd1_ru j书 +2(1 +r)Uj -ru j /k=ru 爲 +2(1 _r)u： + ru 廿 +2计 j 2j =1,2，,N -1, k =0,1，,M -1u0 =®(Xj), j =1,2，N -1kk=g(tk),un =g

23、2(tk), k =0,1，M它的矩阵形式为f1k -1kk 2(I B)uk 1 =(l A)uk F 2, k =1,2 ,M u0 二，k =0,1, ,M -1(9.28)(9.29)在每层计算时，仍需求解一个三对角形方程组。4. Richards。n格式利用公式（9.19 ）和（9.20），可导出另一个截断误差为O（2 h2）阶的差分方程k 1 k 1Uj _uj2kc k 丄 kUj 1 -2Uj uj4fjk称之为 Richardson格式。可改写为k 1Uj二 u： 2r(uk 十-2u： u：)2 f ；(9.32)这是一个三层显式差分格式。在逐层计算时，需用到u：和u：两层

24、值才能得到k+1层值u： 1。这样，从第0层已知值u： h（Xj）开始，还须补充上第一层值 u1，才能逐层计算下去。可采用前述的双层格式计算u1。除上述四种差分格式外，还可构造出许多逼近偏微分方程（9.14 ）的差分格式，但并不是每个差分格式都是可用的。一个有实用价值的差分格式应具有如下性质：（1）收敛性。对任意固定的节点（Xj，tk）,当剖分步长.，h； 0时，差分解u： 应收敛到精确解 u （Xj，tk）。（2）稳定性。当某一时间层计算产生误差时，在以后各层的计算中，这些误差 的传播积累是可控制而不是无限增长的。理论上可以证明，在一定条件下，稳定的差分格式必然是收敛的。因此，这里主要研究差

25、分格式的稳定性。作为例子，先考查 Richards。n格式的稳定性。设 U：是当计算过程中带有误差时，按Richardson格式（9.30 ）得到的实际计算值。u：是理论值，误差e：-U：。1假定右端项fjk的计算是精确的，网比r ，则e：满足j 2 jk 1 k 1kk kej二 ej © i -2qq（9.31）设前k-1层计算时精确的，误差只是在第k层j0点发生，即k 1kkej0, ej0 = ；' ej 0, j j 0。则利用（9.31 ）式可得到误差:的传播情况，见表9-1。表9-1 r=1/2 时Richardson 格式的误差传播Xj0-4j0-3j0-2j

26、0-1j0j0+1j0+2j0+3j0+4k0000£0000k+1000£-2 ££000k+200£-4 £7 £-4 ££00k+30£-6 £17 £-24 £17 £-6 ££0k+4£-8 £31 £-68 £89 £-68 £31 £-8 ££k+5-10 £49 £-144 £273 £-38

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

28、差分格式都可以表示为如下的矩阵方程形式 kU02-Huk4 Fk(9.32)其中H称为传播矩阵。对于显格式 H =A,隐格式H =B-1,六点对称格式 H =（ I +B） （I+ A）。一般的三层格式也可以转化为双层格式。为了讨论方便，设在初始层产生误差；°，且假定右端项Fk的计算是精确的。用uk表示当初始层存在误差；°时，由差分格式（9.32 ）得到的计算解，则k LJ 2ku =Hu + F一° 木丄 °u =伞+名记误差向量；k =Jk -uk，则；*满足方程言=H 4, k =1,2,呂°为初始误差定义9.1称差分格式（9.32 ）是

29、稳定的，如果对任意初始误差；°在某种范数下满足兰k 3 0,° s < 1°其中C为与h, T无关的常数。这个定义表明，当差分格式稳定时，它的误差传播是可控制的。从（9.34 ）式递推得到kk °TH ； ,0 _ k _T因此，差分格式稳定的充分必要条件是Hk -C,定理9.3（稳定性必要条件）差分格式稳定的必要条件是存在与常数M,使谱半径，(H ) _1 M * * 定理9.4（稳定性充分条件）设H为正规矩阵，即HH =H式也是差分格式稳定的充分条件。下面讨论几种差分格式的稳定性。为便于讨论，引进N-1阶矩阵° 1 11 ° 1+ + +s=I1°1I1°uk满足方程(9.33 )(9.34 )，误差向量；k(9.35 )(9.36 )h, t无关的(9.37 )H ,则(9.37 )这个特殊矩阵的特征值为(9.38 )S =2cos匚 j =1,2, N -1jN例9-1古典显格式此时H =A=(1-2r)l+rS。利用(9.38)式和三角函数公式,可求得H的特征值为2 j 兀、.