偏微分方程數值解法

1、十 章 偏 微 分方 程 數 值 解 法偏微分方程問題，其求解十分困難。除少數特殊情況外，絕大多數情況均難以求出精確解。因此，近似解法就顯得更為重要。本章僅介紹求解各類典型偏微分方程定解問題的差分方法。 1差分方法的基本概念1.1幾類偏微分方程的定解問題橢圓型方程：其最典型、最簡單的形式是泊松( Poisson )方程特別地，當 f(x,y) 0時，即為拉普拉斯(Laplace )方程，又稱為調和方程Poisson方程的第一邊值問題為其中 為以 為邊界的有界區(qū)域，為分段光滑曲線，稱為定解區(qū)域，f (x, y), (x, y)分別為 ,上的已知連續(xù)函數。第二類和第三類邊界條件可統(tǒng)一表示為其中n為

2、邊界的外法線方向。當0時為第二類邊界條件，0時為第三類邊界條件。拋物型方程：其最簡單的形式為一維熱傳導方程方程可以有兩種不同類型的定解問題：初值問題初邊值問題其中 (x), g1(t), g2(t)為已知函數，且滿足連接條件邊界條件 u(0,t) g1(t),u(l,t) g2(t)稱為第一類邊界條 件。第二類和第三類邊界條件為其中 1(t) 0, 2(t) 0。當 1(t)2(t)0 時，為第二類邊界條件，否則稱為第三類邊界條件。雙曲型方程：最簡單形式為一階雙曲型方程物理中常見的一維振動與波動問題可用二階波動方程描述，它是雙曲型方程的典型形式。方程的初值問題為邊界條件一般也有三類，最簡單的初

3、邊值問題為1.2差分方法的基本概念差分方法又稱為有限差分方法或網格法，是求偏微分方程定 解問題的數值解中應用最廣泛的方法之一。它的基本思想是：先對求解區(qū)域作網格剖分，將自變量的連 續(xù)變化區(qū)域用有限離散點(網格點)集代替；將問題中出現(xiàn)的連 續(xù)變量的函數用定義在網格點上離散變量的函數代替；通過用網 格點上函數的差商代替導數，將含連續(xù)變量的偏微分方程定解問 題化成只含有限個未知數的代數方程組(稱為差分格式)。如果 差分格式有解，且當網格無限變小時其解收斂于原微分方程定解 問題的解，則差分格式的解就作為原問題的近似解(數值解)。 因此，用差分方法求偏微分方程定解問題一般需要解決以下問題:(1)選取網格

4、；(2)對微分方程及定解條件選擇差分近似，列出差分格式；(3)求解差分格式；(4)討論差分格式解對于微分方程解的收斂性及誤差估計。下面，用一個簡單的例子來說明用差分方法求解偏微分方程 問題的一般過程及差分方法的基本概念。設有一階雙曲型方程初值問題。(1)選取網格：-2h-h0h2h3h首先對定解區(qū)域D (x,t) x ,t 0作網格剖分，最簡單常用一種網格是用兩族分別平行于 x軸與t 軸的等距直線x Xk kh,t tj j (k 0, 1, 2L ,j 0,1,2,L u等D 分成許 多小矩形區(qū)域。這些直線稱為網格線，其交點稱為網格點，也稱為節(jié)點,h和 分別稱作 X 方向和t方向的步長。這種

5、網格稱為矩形網格(2)對微分方程及定解條件選擇差分近似，列出差分格式： 如果用向前差商表示一階偏導數，即其中01, 21方程處可表小為其中 u(Xk,0)(Xk)(k 0,1,2,L)。由于當h, 足夠小時，在式中略去R(Xk,tj),就得到一個與方程相近似的差分方程此處，Uk, j可看作是問題的解在節(jié)點(x k , t j )處的近似值。同初值條件結合，就得到求問題的數值解的差分格式。式O( qh p),則稱差分方稱為差分方程的截斷誤差。如果一個差分方程的截斷誤差為 R程對t是q階精度，對x是p階精度的。顯然，截斷誤差的階數 越大，差分方程對微分方程的逼近越好。若網格步長趨于0時，差分方程的

6、截斷誤差也趨于0,則稱 差分方程與相應的微分方程是相容u 這是用差分方法求解偏微 分方程問題的必要條件。如果當網格步長趨于0時，差分格式的解收斂到相應微分方 程定解問題的解，則稱這種差分格式是收斂的。 2橢圓型方程第一邊值問題的差分解法本節(jié)以Poisson方程為基本模型討論第一邊值問題的差分方法。差分格式的建立考慮Poisson方程的第一邊值問題取h, 分別為x方向和y方向的步長，如圖所示，以兩族平行線x xk kh, y yj j (k,j 0, 1, 2,L )將定R (xk, yj)解區(qū)域剖分成矩形網 格。節(jié)點的全體記為xk kh, yjj , k, j為整數 0定解區(qū)域內部的節(jié)點稱為內

7、點，記內點集 R 為 h 。邊界 與網格 線的交點稱為邊界點，邊界點全體記為h o與節(jié)點(x, yj)沿x方向或y方向只差一個步長的點(xk i,yjM(xk,yj .稱為節(jié)點(Xk,yj)的相鄰節(jié)點。如果一個內點的四個相鄰節(jié)點均屬于，稱為正一 (i)一 一則內點，正內點的全體記為，至少有一個相鄰干點不屬于-一“一 “一，(2)、目工的內點稱為非正則內點，非正則內點的全體記為。問題是要求出第一邊值問題在全體內點上的數值解。為簡便，記(k,j) (Xk,yj), u(k, j) u(4,yj),fk,j f(xk，yj)。對正則(1),人，,內點(k,j),由二階中心差冏公式2u2Poisson

8、 方程 x2了 f(x，V)在點(k, j)處可表小為 其中為其截斷誤差表示式，略去R(k, j),即得與方程相近似的差分方程式中方程的個數等于正則內點的個數，而未知數Uk,j則除了包含正則內點處解 U 的近似值外，還包含一些非正則內點處U的近似值，因而方程個數少于未知數個數。在非正則內點處Poisson方程的差分近似不能按上式給出，需要利用邊界條件得到。邊界條件的處理可以有各種方案，下面介紹較簡單的兩種。 (1)直接轉移用最接近非正則內點的邊界點上的 u值作為該點上 u 值的 近似，這就是邊界條件的直接轉移。例如，點 P(k, j)為非正則內點，其最接近的邊界點為 Q點，則有上式可以看作是用

9、零次插值得到非正則內點處 u的近似值，容易求出，其截斷誤差為 O(h )。 將上式代入，方程個數即與未知數 個數相等。(2)線性插值這種方案是通過用同一條網格線上與點P相鄰的邊界點與內點作線性插值得到非正則內點 P(k, j)處U值的近似。由點R與T的線性插值確定U(P)的近似值U k, j ,得2其中d RP ,其截斷誤差為O(h )。將其與方程相近似的差分程聯(lián)立，得到方程個數與未知數個數相等的方程組，求解此方程組可得Poisson方程第一邊值問題的數值解。上面所給出的差分格式稱為五點菱形格式，實際計算時經常取h ,此時五點菱形格式可化為簡記為 uk,jfk,j其中 uk,j uk 1,j

10、uk 1,j例1用五點菱形格式求解拉普拉斯(其中(x, y)0 x, yuk,j 1uk,j 14uk,j0Laplace )方程第一邊值問題_1(0,0) (1,0)(2,0)(3,0)解網格中有四個內點，均為正則內點。由五點菱形格式，得方 程組h2(U2,1U0,1U1 ,2U1, 04u1,1 )5 (U31U1,1 U2,2 U2,04U2,1 )0h2(U3,2U0,2 U1，3 U1,14U1,2)0 代入邊界條件U1,2 U2,3 U2,14u2,2 )0U1,0U0,1U1,3U3,1其解為,16 ig I, u 9,10 ig 一, u 9,25|gm,u,37ig 一, u

11、 92,00,22,33,2U1,10.2756919 ,ig25,13喝,341g萬,40 lg3U2,10.4603488u1,20.3467842, u2,20.5080467當h 時，對1卜2 (Uk 1,jUk 1,j Uk,j 1 Uk,j 1 4uk,j) fk,j禾用點(k,j), (k 1, j 1), (k 1, j 1)構造的差分格式，稱為五點矩形格式，簡記為. 2 口 Uk,jfk,jh其中口uk,j uk 1,j 1uk 1,j 1uk 1,j 1uk 1,j 14uk,j ,其截斷誤差為2五點菱形格式與矩形格式的截斷誤差均為O(h ),稱它們具有二階精度。如果用更多

12、的點構造差分格式，其截斷誤差的階數可 以提高，如利用菱形格式及矩形格式所涉及的所有節(jié)點構造出的 九點格式就是具有四階精度的差分格式。 3拋物型方程的差分解法以一維熱傳導方程為基本模型討論適用于拋物型方程定解問題的幾種差分格式。3.1差分格式的建立首先對xt平面進行網格剖分。分別取h,為x方向與t方向的步長，用兩族平行直線x xkkh(k 0, 1, 2,L ),ttjj (j 0,1,2 )，將 xt 平面剖分成矩形網格，節(jié)點為(xk,tj)(k 0, 1, 2,L ,j 0,1,2L )。為簡便，記(k, j) (xk,tj) , u(k, j) u(xk,tj),k (xk), g1jgt

13、j),g2jg2 (tj ),1 j1 (tj ),2j 2(tj ) o(一)微分方程的差分近似在網格內點(k, j)處，對 一J分別采用向前、向后及中心差商公式一維熱傳導方程可分別表示為由此得到一維熱傳導方程的不同差分近似2、上述差分方程所用到的節(jié)點各不相同。其截斷誤差分別為O( h ),222、0( h )和0( h )。因此，它們都與一維熱傳導方程相容。 如果將式1()中的u k, j用2 (Uk 1Uk，J 1 )代替，則可得到又一種差分近似差分方程用到四個節(jié)點。由Taylor公式容易得出,2. 2 x故其的截斷誤差為0( h )h 2。因而不是對任意的h,0,此差分方程都能逼近熱傳

14、導方程(a 0)僅當 o(h)時，才成立。綜上可知，用不同的差商公式可以得到微分方程的不同的差 分近似。構造差分格式的關鍵在于使其具有相容性、收斂性和穩(wěn) 定性。前面三個方程都具有相容性，而此方程則要在一定條件下 才具有相容性。(二)初、邊值條件的處理為用差分方法求解定解問題初值問題初邊值問題還需對定解條件進行離散化。對初始條件及第一類邊界條件，可直接得到 TOC o 1-5 h z lT苴中n, m 一八十h對第二、三類邊界條件需用差分近似。下面介紹兩種較簡單的處理方法。u(1)在左邊界(x 0 )處用向前差商近似偏導數在右邊界(xl )處用向后差商近似則得邊界條件的差分近似為其截斷誤差為 0

15、(h)u(2)用中心差商近似，即x則得邊界條件的差分近似為2其截斷誤差為o(h )0誤差的階數提高了，但出現(xiàn)定解區(qū)域外的節(jié)點(1, j)和(n 1, j),這就需要將解拓展到定解區(qū)域外。可以通過用內節(jié)點上的U值插值求出u 1,j和un 1,j ,也可以假定熱傳導方程在邊界上也成立，將差分方程擴展到邊界節(jié)點上，由此消去 u i,j和 u n 1, j。(三)幾種常用的差分格式以熱傳導方程的初邊值問題 為例給出幾種常用的差分格式。(1)古典顯式格式a令r . 2 ,貝Uuk,j 1 uk,j uk 1, j 2uk, j uk 1, ja2 0h2可改寫成uk,j 1 ru k 1,j(1 2r)

16、uk,j ru k 1,j將其與初始條件及第一類邊界條件結合，我們得到求解此問題的一種差分格式由于第0層(j0)上節(jié)點處的u值已知(uk,0k),由此即可算出u在第一層(j1)上節(jié)點處的近似值U k,1。重復使用此式，可以逐層計算出所有的Uk,j ,因此此差分格式稱為古典顯式格式。又因式中只出現(xiàn)相鄰兩個時間層的節(jié)點，故此式是二層顯式格式。(2)古典隱式格式將式Uk,j Uk,j 1Uk 1,j a2Uk, jUk 1,jh2整理并與初始條件及第一類邊界條件式聯(lián)立，得差分格式如下a其中r 卜2 。雖然第0層上的U值仍為已知，但不能由上式直接計算以上各層節(jié)點上的值 Uk,j ,必須通過解下列線性方

17、程組才能由Uk,j計算Uk,j 1 ,故此差分格式稱為古典隱式格式。此方程組是三對角方程組，且系數矩陣嚴格對角占優(yōu)，故解存在唯一。(3) Richardson 格式Richardson格式是將式整理后與初始條件及第一類邊界條件式聯(lián)立。其計算公式為這種差分格式中所涉及的節(jié)點出現(xiàn)在j 1, j, j 1三層上，故為三層顯式格式。Richardson格式是一種完全不穩(wěn)定的差分格式，因此它在實際 計算中是不能采用的。(4)杜福特-弗蘭克爾(DoFort-Frankel )格式DoFort-Frankel 格式也是三層顯式格式，它是由式與初始條件及第一類邊界條件式結合得到的。具體形式如下：用這種格式求解

18、時，除了第0層上的值U k,0由初值條件得到，必須先用二層格式求出第1層上的值Uk,1 ,然后再按上式逐層計算 Uk,j(j 2,3, ,m)。(5)六點隱式格式對二階中心差商公式2在點(k,j 1)與點(k, j) 處的二階中心差商的平如果用 Y 2 xu均值近似2u2x(i 1)(k，J 2)處的值，即(k i )同時在點(k, j 2)處的值也用中心差商近似，即這樣又得到熱傳導方程的一種差分近似22、其截斷誤差為O(h ),將上式與初始條件及第一類邊界條件式聯(lián)立并整理，得差分格式此格式涉及到六個節(jié)點，它又是隱式格式，故稱為六點隱式格式。與古典隱式格式類似，用六點格式由第 j層的值u k,

19、 j計算第j 1層的值uk,j 1時，需求解三對角方程組此方程組的系數矩陣嚴格對角占優(yōu)，故仍可用追趕法求解。例2用古典顯式格式求初邊值問題 TOC o 1-5 h z 的數值解，取h 1,0.5 0a 1,r ar0.5(x)x2解這里 h h2，(x)x ,g) 0,g2(t) 9。由格式可得到將初值u k ,0代入上式，即可算出將邊界條件u0,10,u3,19及上述結果代入又可求得如此逐層計算，得全部節(jié)點上的數值解為 4雙曲型方程的差分解法對二階波動方程_u_u如果令v1f , v2x ,則方程可化成一階線性雙曲型方程組記v (Vi, V2 )T ,則方程組可表成矩陣形式矩陣 A有兩個不同的特征值a,故存在非奇異矩陣 P,使得作變換w Pv (wi,w2),方程組可化為w wtx方程組由二個獨立的一階雙曲型方程聯(lián)立而成。因此本節(jié)主要討論一階雙曲型方程的差分解法。4.1