偏微分方程數(shù)值解法

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

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

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

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

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

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

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

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

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

10、uk 1,j例1用五點(diǎn)菱形格式求解拉普拉斯(其中(x, y)0 x, yuk,j 1uk,j 14uk,j0Laplace )方程第一邊值問(wèn)題_1(0,0) (1,0)(2,0)(3,0)解網(wǎng)格中有四個(gè)內(nèi)點(diǎn)，均為正則內(nèi)點(diǎn)。由五點(diǎn)菱形格式，得方 程組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萬(wàn),40 lg3U2,10.4603488u1,20.3467842, u2,20.5080467當(dāng)h 時(shí)，對(duì)1卜2 (Uk 1,jUk 1,j Uk,j 1 Uk,j 1 4uk,j) fk,j禾用點(diǎn)(k,j), (k 1, j 1), (k 1, j 1)構(gòu)造的差分格式，稱為五點(diǎn)矩形格式，簡(jiǎn)記為. 2 口 Uk,jfk,jh其中口uk,j uk 1,j 1uk 1,j 1uk 1,j 1uk 1,j 14uk,j ,其截?cái)嗾`差為2五點(diǎn)菱形格式與矩形格式的截?cái)嗾`差均為O(h ),稱它們具有二階精度。如果用更多

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

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

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

15、(h)u(2)用中心差商近似，即x則得邊界條件的差分近似為2其截?cái)嗾`差為o(h )0誤差的階數(shù)提高了，但出現(xiàn)定解區(qū)域外的節(jié)點(diǎn)(1, j)和(n 1, j),這就需要將解拓展到定解區(qū)域外。可以通過(guò)用內(nèi)節(jié)點(diǎn)上的U值插值求出u 1,j和un 1,j ,也可以假定熱傳導(dǎo)方程在邊界上也成立，將差分方程擴(kuò)展到邊界節(jié)點(diǎn)上，由此消去 u i,j和 u n 1, j。(三)幾種常用的差分格式以熱傳導(dǎo)方程的初邊值問(wèn)題 為例給出幾種常用的差分格式。(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將其與初始條件及第一類邊界條件結(jié)合，我們得到求解此問(wèn)題的一種差分格式由于第0層(j0)上節(jié)點(diǎn)處的u值已知(uk,0k),由此即可算出u在第一層(j1)上節(jié)點(diǎn)處的近似值U k,1。重復(fù)使用此式，可以逐層計(jì)算出所有的Uk,j ,因此此差分格式稱為古典顯式格式。又因式中只出現(xiàn)相鄰兩個(gè)時(shí)間層的節(jié)點(diǎn)，故此式是二層顯式格式。(2)古典隱式格式將式Uk,j Uk,j 1Uk 1,j a2Uk, jUk 1,jh2整理并與初始條件及第一類邊界條件式聯(lián)立，得差分格式如下a其中r 卜2 。雖然第0層上的U值仍為已知，但不能由上式直接計(jì)算以上各層節(jié)點(diǎn)上的值 Uk,j ,必須通過(guò)解下列線性方

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

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

19、 j計(jì)算第j 1層的值uk,j 1時(shí)，需求解三對(duì)角方程組此方程組的系數(shù)矩陣嚴(yán)格對(duì)角占優(yōu)，故仍可用追趕法求解。例2用古典顯式格式求初邊值問(wèn)題 TOC o 1-5 h z 的數(shù)值解，取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é)果代入又可求得如此逐層計(jì)算，得全部節(jié)點(diǎn)上的數(shù)值解為 4雙曲型方程的差分解法對(duì)二階波動(dòng)方程_u_u如果令v1f , v2x ,則方程可化成一階線性雙曲型方程組記v (Vi, V2 )T ,則方程組可表成矩陣形式矩陣 A有兩個(gè)不同的特征值a,故存在非奇異矩陣 P,使得作變換w Pv (wi,w2),方程組可化為w wtx方程組由二個(gè)獨(dú)立的一階雙曲型方程聯(lián)立而成。因此本節(jié)主要討論一階雙曲型方程的差分解法。4.1