偏微分方程差分方法

1、第9章偏微分方程的差分方法含有偏導(dǎo)數(shù)的微分方程稱為偏微分方程。由于變量的增多和區(qū)域的復(fù)雜性，求偏微分方程的精確解一般是不可能的，經(jīng)常采用數(shù)值方法求方程的近似解。偏微分方程的數(shù)值方法種類較多，最常用的方法是差分方法。差分方法具有格式簡單，程序易于實現(xiàn)，計算量小等優(yōu)點(diǎn)，特別適合于規(guī)則區(qū)域上偏微分方程的近似求解。本章將以一些典型的偏微分方程為例，介紹差分方法的基本原理和具體實現(xiàn)方法。9.1 橢圓型方程邊值問題的差分方法9.1.1 差分方程的建立最典型的橢圓型方程是Poisson(泊松)方程：2u：2u一一Au三一(r+:)=f(x,y),(x,y)eG(9):xcyG是x,y平面上的有界區(qū)域，其邊界

2、r為分段光滑的閉曲線。當(dāng)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)都是已知的函數(shù)，k(x,y)>0O滿足方程(9.1)和上述三種邊值條件之一的光滑函數(shù)u(x,y)稱為橢圓型方程邊值問題的解。用差分方法求解偏微分方程，就是要求出精確解u(x,y)在區(qū)域G的一些離散節(jié)點(diǎn)(xy。上的近似值口

3、一"(為，y)。差分方法的基本思想是，對求解區(qū)域G做網(wǎng)格剖分，將偏微分方程在網(wǎng)格節(jié)點(diǎn)上離散化，導(dǎo)出精確解在網(wǎng)格節(jié)點(diǎn)上近似值所滿足的差分方程，最終通過求解差分方程，通常為一個線性方程組，得到精確解在離散節(jié)點(diǎn)上的近似值。設(shè)G=0<x<a,0<y<b為矩形區(qū)域，在x,y平面上用兩組平行直線x=ih1,i=0,1,N1,h1=a/N1y=jh2,j=0,1,N2,h2=b/N2將G剖分為網(wǎng)格區(qū)域，見圖9-1。h1,h2分別稱為x方向和y方向的剖分步長，網(wǎng)格交點(diǎn)(xi,yO稱為剖分節(jié)點(diǎn)(區(qū)域內(nèi)節(jié)點(diǎn)集合記為Gh=(X,y);(xi,y。CG),網(wǎng)格線與邊界r的交點(diǎn)稱為邊界

4、點(diǎn)，邊界點(diǎn)集合記為現(xiàn)在將微分方程（yi）處，方程（9.1）為2u/、-（Xi,yi）：x需進(jìn)一步離散（點(diǎn)函數(shù)值u（Xi,差商表達(dá)式:：2u-2X,：2u-2y9.1）在每一個內(nèi)節(jié)點(diǎn)（Xi,yi）上進(jìn)行離散。在節(jié)點(diǎn)（Xi,二2u2(為7)=f(Xi,y)(x/)Ghy(9.5)9.5）中的二階偏導(dǎo)數(shù)。為簡化記號，簡記節(jié)點(diǎn)（Xi,yi）=（i,j）,節(jié)yi）=u（i,j）o利用一元函數(shù)的Taylor展開公式，推得二階偏導(dǎo)數(shù)的12(i,j)=77u(i1,j)-2u(i,j)u(i-1,j)0(h2)A(i,j)=u(i,j1)-2u(i.j)u(i,j-1)0(hf)h2代入（9.5）式中，得到方

5、程（9.1）在節(jié)點(diǎn)（i,j）處的離散形式1-7u(i1,j)-2u(i1j)u(i-1,j)-h10(h12h2),(i,j)Gh1-2u(i,j1)-2u(i,j)u(i,j-1)其中=f（為，）。舍去高階小項0（h2+h；）,就導(dǎo)出了u（i,j）的近似值ui,j所滿足的差分方程11.-7Tui+,j-2ui,j+ui4.,j-77ui,j4r-2加+ui,j=fi,j,(i,j)Gh(9.6)h1h2在節(jié)點(diǎn)（i,j）處方程（9.6）逼近偏微分方程（9.1）的誤差為O（h2+h；）,它關(guān)于剖分步長是二階的。這個誤差稱為差分方程逼近偏微分方程的截斷誤差，它的大小將影響近似解的精度。在差分方程（

6、9.6）中，每一個節(jié)點(diǎn)（i,j）處的方程僅涉及五個節(jié)點(diǎn)未知量u，j,ui+1jui-1,j,ui,j+1,ui,j-1,因此通常稱（9.6）式為五點(diǎn)差分格式，當(dāng)h1=h2=h時，它簡化為h2ui1,jUi,jui,j1Ui,j-4Ui,j=fi,j,(i,j)Gh差分方程（9.6）中，方程個數(shù)等于內(nèi)節(jié)點(diǎn)總數(shù)，但未知量除內(nèi)節(jié)點(diǎn)值Ui,j,（i,j）£Gh外，還包括邊界點(diǎn)值。例如，點(diǎn)（1,j）處方程就含有邊界點(diǎn)未知量U0,j。因此，還要利用給定的邊值條件補(bǔ)充上邊界點(diǎn)未知量的方程。對于第一邊值條件式（9.2）,可直接取Ui,j=a（Xi,yi）,（i,j）Th（9.7）(03*圖9-2對于

7、第三（k=0時為第二）邊值條件式（9.4）,以左邊界點(diǎn)（1,j）為例，見圖9-2,利用一階差商公式.:uu（0,j）-u（1,j）（0,j）二1O（1）一nh1則得到邊界點(diǎn)（0,j）處的差分方程u0,j-也h1-k0,ju0,j=r0,j(9.8)聯(lián)立差分方程（9.6）與（9.7）或（9.8）就形成了求解Poisson方程邊值問題的差分方程組，它實質(zhì)上是一個關(guān)于未知量Ui,j的線性代數(shù)方程組，可采用第2,3章介紹的方法進(jìn)行求解。這個方程組的解就稱為偏微分方程的差分近似解，簡稱差分解。考慮更一般形式的二階橢圓型方程UU、c：uU-(A-)+丁(Bt)+Ct+D丁+Eu=f(x,y),(x,y產(chǎn)G

8、(9.9)二x二x二y二y二x二y1.其中A(x,y)>Amn>0,B(x,y)>B/>0,E(x,y)>0。引進(jìn)半節(jié)點(diǎn)x1=xi±-h1i二22,1y1=yi士一h2利用一階中心差商公式，在下點(diǎn)（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衛(wèi)u(i-7-1)2),x2h1:u:u.對（B）,類似處理，就可推得求解方程（9.9）的差分方程.y二y二y-ai1,j

9、ui,jai,jui,j-ai,jdui,j1'ai,jdui,j1aai,jui,j1y9.10)=f(i,j),(i,j)Gh其中nh1_=幾(Ai1Ci,j)2，j2,_2,工_、=hi(AiG,j)i工"2ai,j1ai,j1士九B.ii,j-2”九Bii，j-22Di，j)-")(9.ii)B.i)Ei,ji、_2_2_利=%(AA,)h2(B.ii'一"i一，ji,j-222顯然，當(dāng)系數(shù)函數(shù)A（x,y）=B（x,y）=i,C（x,y）=D（x,y）=E（x,y）=0時，橢圓型方程（9.9）就成為Poisson方程（9.1）,而差分方程（

10、9.10）就成為差分方程（9.6）。容易看出，差分方程（9.10）的截斷誤差為O（h；+h；）階。9.1.2一般區(qū)域的邊界條件處理前面已假設(shè)G為矩形區(qū)域，現(xiàn)在考慮G為一般區(qū)域情形，這里主要涉及邊界條件的處理?？紤]Poisson方程第一邊值問題(9.12)M=f(x,y),(x,y)wGU=s(x,y),(x,y)w其中G可為平面上一般區(qū)域，例如為曲邊區(qū)域。仍然用兩組平行直線：x=x0+ihi,y=yO+jh2,i,j=0,±1,對區(qū)域G進(jìn)行矩形網(wǎng)格剖分，見圖9-3。如果一個內(nèi)節(jié)點(diǎn)（i,j）的四個相鄰節(jié)點(diǎn)（i+1,j）,（i-1,j）,（i,j+1）和（i,j-1）屬于G=G，j，則稱

11、其為正則內(nèi)點(diǎn)，見圖9-3中打“?！碧栒?；如果一個節(jié)點(diǎn)（i,j）屬于G且不為正則內(nèi)點(diǎn)，則稱其為非正則內(nèi)點(diǎn)，見圖9-3中打".”號者。記正則內(nèi)點(diǎn)集合為G非正則內(nèi)點(diǎn)集合為"。顯然，當(dāng)G為矩形區(qū)域時，Gh=Gh,甲=曙成立。（如圖9-4中E點(diǎn)）上的u的值作為Ubh1二(A)6uC,h1、C,"Xb-Xa在正則內(nèi)點(diǎn)（i,j）處，完全同矩形區(qū)域情形，可建立五點(diǎn)差分格式1.八r1.-S-Ui,jU-J-h/jLMC,OaGh©13)在方程（9.13）中，當(dāng)（i,j）點(diǎn)臨近邊界時，將出現(xiàn)非正則內(nèi)點(diǎn)上的未知量，因此必須補(bǔ)充非正則內(nèi)點(diǎn)處的方程。若非正則內(nèi)點(diǎn)恰好是邊界點(diǎn)，如圖

12、9-4中D點(diǎn)，則利用邊界條件可取Ud=oc（D）對于不是邊界點(diǎn)的非正則內(nèi)點(diǎn)，如圖9-4中點(diǎn)，一般可采用如下兩種處理方法。a.直接轉(zhuǎn)移法.取與點(diǎn)B距離最近的邊界點(diǎn)u（B）的近似值UB,即Ub=u（E）="（E）直接轉(zhuǎn)移法的優(yōu)點(diǎn)是簡單易行，但精度較低，只為一階近似。b.線性插值法.取B點(diǎn)的兩個相鄰點(diǎn)（如圖9-4中邊界點(diǎn)A和正則內(nèi)點(diǎn)C作為插值節(jié)點(diǎn)對u（B）進(jìn)行線性插值u(B)Xc-Xbu(A)Xb-Xau(C)O(h12)Xc-XaXc-Xa則得到點(diǎn)B處的方程線性插值法精度較高，為二階近似。對每一個非正則內(nèi)點(diǎn)進(jìn)行上述處理，將所得到的方程與（9.13）式聯(lián)立，就組成了方程個數(shù)與未知量個數(shù)相一

13、致的線性代數(shù)方程組。求解此方程組就可得到一般區(qū)域上邊值問題（9.12）的差分近似解。對于一般區(qū)域上二階橢圓型方程（9.9）的第一邊值問題，可完全類似處理。第二、三邊值條件的處理較為復(fù)雜，這里不再討論。9.2拋物型方程的差分方法本節(jié)介紹拋物型方程的差分方法，重點(diǎn)討論差分格式的構(gòu)造和穩(wěn)定性分析。9.2.1 一維問題作為模型，考慮一維熱傳導(dǎo)的初邊值問題2*=a-2+f(x,t),0cxMl,0<t<T(9.14).:tjxu(x,0)=*(x),0cx<l(9.15)u（0,t）=g（t）,u（l,t）=g2（t）,0<t<T（9.16）其中a是正常數(shù)，f（x,t）,中

14、（x）,g（t）和g2（t）都是已知的連續(xù)的函數(shù)?，F(xiàn)在討論求解問題（9.14）-（9.18）的差分方法。首先對求解區(qū)域G=0WxWl,0<t<T進(jìn)行網(wǎng)格剖分。取空間步長h=l/N,時間步長t=T/M,其中N,M是正整數(shù)，作兩族平行直線x=xj=jh,j=0,1,Ntf=k,k=0,1,M將區(qū)域G剖分成矩形網(wǎng)格，見圖9-5,網(wǎng)格交點(diǎn)（為，t。稱為節(jié)點(diǎn)。0123_N-1圖95用差分方法求解初邊值問題（9.14）-（9.16）就是要求出精確解u（x,t）在每k個下點(diǎn)（xj,tk）處的近似值uju（xj,tk）。為簡化記號，簡記節(jié)點(diǎn)（xj,tk）=u（j,k）。利用一元函數(shù)的Taylor展

15、開公式，可推出下列差商表達(dá)式(9.17)(9.18)(9.19)(9.20)衛(wèi)u(j,k1)-u(j,k)0()Zt衛(wèi)此水)川IDo()由衛(wèi)(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.古典顯格式在區(qū)域G的內(nèi)節(jié)點(diǎn)(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島1 h2j其中f；=f(Xi,tk)o舍去高階小項O(w+h2),就得到節(jié)點(diǎn)近似值(差分解)uk

16、所滿足的差分方程kkk"j(9.21)uj1-25Ujk2 a：2fjh顯然，在節(jié)點(diǎn)(j,k)處，差分方程(9.21)逼近偏微分方程(9.14)的誤差為O(T+h2),這個誤差稱為截斷誤差，它反映了差分方程逼近偏微分方程的精度?，F(xiàn)將(9.21)式改寫為便于計算的形式，并利用初邊值條件(9.15)與(9.16)補(bǔ)充上初始值和邊界點(diǎn)方程，則得到k+k4/c、kk-kuj=ruj由十(12r)uj+ruj+Tfjj=1,2,，N1,k=0,1，M-1/、,'',(9.22uj=5(Xj),j=1,2,，N-1kkjo=g1(tk),un=g2(tk),k=0,1,m其中r=

17、a=稱為網(wǎng)比。h2與時間相關(guān)問題差分方程的求解通常是按時間方向逐層進(jìn)行的。對于差分方程kk1(9.22),當(dāng)?shù)趉層節(jié)點(diǎn)值山已知時，可直接計算出第k+1層節(jié)點(diǎn)值uj。這樣，從第0層已知值u0=邛(為)開始，就可逐層求出各時間層的節(jié)點(diǎn)值。差分方程(9.22)的求解計算是顯式的，無須求解方程組，故稱為古典顯格式。此外，在式(9.22)中，每個內(nèi)節(jié)點(diǎn)處方程僅涉及k和k+1兩層節(jié)點(diǎn)值，稱這樣的差分格式為雙層格式。差分方程(9.22)可表示為矩陣形式kk+八ku=AuFkk=0,1,M-1(9.23)其中1 -2rr-r1-2rr+A=,'+，二rr1-2rkkk.kTU二(U1,UN)Fk=(f

18、1k-rg1(tk),f2k，f：2f；rg2(tk)T：=(：(X1)J-，；(XnJt2 .古典隱格式在區(qū)域G的內(nèi)節(jié)點(diǎn)(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),則得到如下差分方程kuj一ujkuj1二a.fjk(9.24)它的截斷誤差為O(T+h2),逼近精度與古典顯格式相同。改寫(9.24)式為便于計算的形式，并補(bǔ)充上初始值與邊界點(diǎn)方程，則得到kkkk上士卜一ruj由十(1十2r)ujruj=uj+f(9.25)j=1,2,，

19、N-1,k=0,1，Mu0=&Xj),j=1,2,，N-1k.kU0=g1(tk),un=g2(tk),k=0,1,M與古典顯格式不同，在差分方程（9.25）的求解中，當(dāng)?shù)趉-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=、丁+''-r-r1+2r向量uk,FkW同（9.23）式中定義。從（9.26）式看到，古典隱格式在

20、每一層計算時，都需求解一個三對角形線性方程組，可采用追趕法求解。3.Crank-Nicolson格式（六點(diǎn)對稱格式）利用一元函數(shù)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+）點(diǎn)離散偏微分方程（9.14）,然后利用（9.20）式進(jìn)2一步離散二階偏導(dǎo)數(shù)，則可導(dǎo)出差分方程k"1k"1k1kkk1（9.27）11Uj1-2UjUj4uj1-2UjUjka二a22fj2h2h2j其截斷誤差為O。2+h2）,在時間方向的逼近階較顯

21、格式和隱格式高出一階。這個差分格式稱為Crank-Nicolson格式，有時也稱為六點(diǎn)對稱格式，它顯然是雙層隱式格式。改寫（9.27）式，并補(bǔ)充初始值和邊界點(diǎn)方程得到k書，、k書k書-ruj書+2(1+r)Uj-rujj=ruj書+2(1-r)uj+ruj+24j2j=1,2,，N1,k=0,1，M-1u0=%xj),j=1,2,，N-1k，、k，、,一，一jo=g1(tk),un=g2(tk),k=0,1，m它的矩陣形式為(9.28)1_kHk_k(IB)uk1=(IA)ukF2,k=1,2,Mu0=,k-0,1,M-1(9.29)在每層計算時，仍需求解一個三對角形方程組。4.Richard

22、son格式利用公式（9.19）和（9.20）,可導(dǎo)出另一個截斷誤差為O。2+h2）階的差分方程k1k1Uj-ujkkkUj1-2%Uj二afjh稱之為Richardson格式?？筛膶憺閡：1二u2r(u：1-2u：U：D2f；(9.32)這是一個三層顯式差分格式。在逐層計算時，需用到u：工和u：兩層值才能得到k+1層彳tuk+o這樣，從第0層已知值U0=（Xj）開始，還須補(bǔ)充上第一層值u1,才能逐層計算下去。可采用前述的雙層格式計算U1。除上述四種差分格式外，還可構(gòu)造出許多逼近偏微分方程（9.14）的差分格式，但并不是每個差分格式都是可用的。一個有實用價值的差分格式應(yīng)具有如下性質(zhì)：（1）收斂性

23、。對任意固定的節(jié)點(diǎn)（Xj,tk）,當(dāng)剖分步長E,hT0時，差分解U：應(yīng)收斂到精確解U（Xj,tk）。（2）穩(wěn)定性。當(dāng)某一時間層計算產(chǎn)生誤差時，在以后各層的計算中，這些誤差的傳播積累是可控制而不是無限增長的。理論上可以證明，在一定條件下，穩(wěn)定的差分格式必然是收斂的。因此，這里主要研究差分格式的穩(wěn)定性。作為例子，先考查Richardson格式的穩(wěn)定性。設(shè)口：是當(dāng)計算過程中帶有誤差時，按Richardson格式（9.30）得到的實際計算值。u：是理論值，誤差e：=U；-U：。假定右端項f：的計算是精確的，網(wǎng)比r=1,則e：滿足j2jk1k1kkkej=ej+（0書一20+e）（9.31）設(shè)前k-1層

24、計算時精確的，誤差只是在第k層j0點(diǎn)發(fā)生，即-kik-kej=0,ej°=與白=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-68s89

25、s-68s31s-8££k+5-10S49s-144£273s-388s273s-144£49s-10Sk+671£-260s641£-1096s1311£-1096s641£-260s71£1從表中看出，誤差是逐層無限增長的。表中的計算雖然是就網(wǎng)比r=進(jìn)行的,2實際上對任何r>0都會產(chǎn)生類似現(xiàn)象，所以Richardson格式是不穩(wěn)定的。利用誤差傳播圖表方法考查差分格式的穩(wěn)定性雖然直觀明了，但只能就具體取定的r值進(jìn)行，并且也不適用于隱式差分格式。9.2.2差分格式的穩(wěn)定性前節(jié)構(gòu)造的幾種雙層差分格式都

26、可以表示為如下的矩陣方程形式kk4k(9.32)u=HuFu°二其中H稱為傳播矩陣。對于顯格式H=A,隱格式H=B-1,六點(diǎn)對稱格式H=（I+B）（I+A）。一般的三層格式也可以轉(zhuǎn)化為雙層格式。為了討論方便，設(shè)在初始層產(chǎn)生誤差總°,且假定右端項Fk的計算是精確的。用uk表示當(dāng)初始層存在誤差/時，由差分格式（9.32）得到的計算解，則k口kku=Hu+F-0A,0u=牛+名記誤差向量J=Uk_u;則Sk滿足方程7k=H”，k=1,2,1°為初始誤差定義9.1稱差分格式（9.32）是穩(wěn)定的，如果對任意初始誤差s0在某種范數(shù)下滿足6klMc卜0|,k>0,0<

27、;T<T0其中C為與h,°無關(guān)的常數(shù)。這個定義表明，當(dāng)差分格式穩(wěn)定時，它的誤差傳播是可控制的。從（9.34）式遞推得到kk0T；=H；,0-k-因此，差分格式穩(wěn)定的充分必要條件是Hk<C,0<k<-定理9.3（穩(wěn)定性必要條件）差分格式穩(wěn)定的必要條件是存在與常數(shù)M,使譜半徑：（H）-1M、.-.*定理9.4（穩(wěn)定性充分條件）設(shè)H為正規(guī)矩陣，即HH=Huk滿足方程(9.33)(9.34)，誤差向量飛(9.35)(9.36)h,t無關(guān)的(9.37)H,則(9.37)式也是差分格式穩(wěn)定的充分條件。下面討論幾種差分格式的穩(wěn)定性。為便于討論，引進(jìn)-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)式和三角函數(shù)公式，可求得H的特征值為2 ji=1-4rsin(),j=1,2,N-1j2N1為使穩(wěn)定性條件式(9.39)成立，必須且只須rM2。由于H=八為實對稱矩陣，所1以古典顯格式穩(wěn)定的充分