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

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

2、ace(拉普拉斯）方程。橢圓型方程的定解條件主要有如下三種邊界條件 第一邊值條件 （9.2） 第二邊值條件 （9.3） 第三邊值條件 （9.4）這里，n表示上單位外法向，(x,y),(x,y),(x,y)和k(x,y)都是已知的函數(shù)，k(x,y)0。滿足方程（9.1）和上述三種邊值條件之一的光滑函數(shù)u(x,y)稱為橢圓型方程邊值問題的解。用差分方法求解偏微分方程，就是要求出精確解u(x,y)在區(qū)域G的一些離散節(jié)點（xi，yi）上的近似值ui,j(xi，yi）。差分方法的基本思想是，對求解區(qū)域G做網(wǎng)格剖分，將偏微分方程在網(wǎng)格節(jié)點上離散化，導出精確解在網(wǎng)格節(jié)點上近似值所滿足的差分方程，最終通過求解

3、差分方程，通常為一個線性方程組，得到精確解在離散節(jié)點上的近似值。設G=0<x<a, 0<y<b為矩形區(qū)域，在x,y平面上用兩組平行直線x=ih1, i=0,1,N1, h1=a/N1 y=jh2, j=0,1,N2, h2=b/N2將G剖分為網(wǎng)格區(qū)域，見圖9-1。h1，h2分別稱為x方向和y方向的剖分步長，網(wǎng)格交點（xi，yi）稱為剖分節(jié)點（區(qū)域內(nèi)節(jié)點集合記為Gh=（xi，yi）; （xi，yi）G）,網(wǎng)格線與邊界的交點稱為邊界點，邊界點集合記為h。 現(xiàn)在將微分方程（9.1）在每一個內(nèi)節(jié)點（xi，yi）上進行離散。在節(jié)點（xi，yi）處，方程（9.1）為 （9.5）需進

4、一步離散（9.5）中的二階偏導數(shù)。為簡化記號，簡記節(jié)點（xi，yi）=(i,j)，節(jié)點函數(shù)值u（xi，yi）=u(i,j)。利用一元函數(shù)的Taylor展開公式，推得二階偏導數(shù)的差商表達式代入（9.5）式中，得到方程（9.1）在節(jié)點(i,j)處的離散形式 其中。舍去高階小項，就導出了u(i,j)的近似值ui,j所滿足的差分方程 （9.6）在節(jié)點(i,j)處方程（9.6）逼近偏微分方程（9.1）的誤差為，它關(guān)于剖分步長是二階的。這個誤差稱為差分方程逼近偏微分方程的截斷誤差，它的大小將影響近似解的精度。在差分方程（9.6）中，每一個節(jié)點(i,j)處的方程僅涉及五個節(jié)點未知量ui,j，ui+1,j，u

5、i-1,j，ui,j+1，ui,j-1，因此通常稱（9.6）式為五點差分格式，當h1= h2=h時，它簡化為差分方程（9.6）中，方程個數(shù)等于內(nèi)節(jié)點總數(shù)，但未知量除內(nèi)節(jié)點值ui,j ,(i,j)Gh外，還包括邊界點值。例如，點(1,j)處方程就含有邊界點未知量u0,j。因此，還要利用給定的邊值條件補充上邊界點未知量的方程。對于第一邊值條件式（9.2），可直接取ui,j=（xi，yi）， (i,j)h （9.7）對于第三（k=0時為第二）邊值條件式（9.4），以左邊界點(1,j)為例，見圖9-2,利用一階差商公式 則得到邊界點(0,j)處的差分方程 （9.8）聯(lián)立差分方程（9.6）與（9.7）或

6、（9.8）就形成了求解Poisson方程邊值問題的差分方程組，它實質(zhì)上是一個關(guān)于未知量ui,j的線性代數(shù)方程組，可采用第2，3章介紹的方法進行求解。這個方程組的解就稱為偏微分方程的差分近似解，簡稱差分解?？紤]更一般形式的二階橢圓型方程 （9.9）其中A(x,y)Amin>0, B(x,y) Bmin >0, E(x,y) 0。引進半節(jié)點利用一階中心差商公式，在節(jié)點（i,j）處可有對類似處理，就可推得求解方程（9.9）的差分方程 （9.10）其中 （9.11）顯然，當系數(shù)函數(shù)A(x,y)=B(x,y)=1, C(x,y)=D(x,y)=E(x,y)=0時，橢圓型方程（9.9）就成為P

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

8、果一個節(jié)點（i,j）屬于且不為正則內(nèi)點，則稱其為非正則內(nèi)點，見圖9-3中打“.”號者。記正則內(nèi)點集合為，非正則內(nèi)點集合為。顯然，當G為矩形區(qū)域時，成立。 在正則內(nèi)點（i,j）處，完全同矩形區(qū)域情形，可建立五點差分格式 （9.13）在方程（9.13）中，當（i,j）點臨近邊界時，將出現(xiàn)非正則內(nèi)點上的未知量，因此必須補充非正則內(nèi)點處的方程。若非正則內(nèi)點恰好是邊界點，如圖9-4中D點，則利用邊界條件可取uD=(D)對于不是邊界點的非正則內(nèi)點，如圖9-4中B點，一般可采用如下兩種處理方法。a.直接轉(zhuǎn)移法.取與點B距離最近的邊界點（如圖9-4中E點）上的u的值作為u(B)的近似值uB，即uB=u(E)=

9、(E)直接轉(zhuǎn)移法的優(yōu)點是簡單易行，但精度較低，只為一階近似。b.線性插值法.取B點的兩個相鄰點（如圖9-4中邊界點A和正則內(nèi)點C作為插值節(jié)點對u(B)進行線性插值則得到點B處的方程 線性插值法精度較高，為二階近似。對每一個非正則內(nèi)點進行上述處理，將所得到的方程與（9.13）式聯(lián)立，就組成了方程個數(shù)與未知量個數(shù)相一致的線性代數(shù)方程組。求解此方程組就可得到一般區(qū)域上邊值問題（9.12）的差分近似解。對于一般區(qū)域上二階橢圓型方程（9.9）的第一邊值問題，可完全類似處理。第二、三邊值條件的處理較為復雜，這里不再討論。9.2 拋物型方程的差分方法本節(jié)介紹拋物型方程的差分方法，重點討論差分格式的構(gòu)造和穩(wěn)定

10、性分析。9.2.1 一維問題 作為模型，考慮一維熱傳導的初邊值問題 （9.14） （9.15） （9.16）其中a是正常數(shù)，都是已知的連續(xù)的函數(shù)?，F(xiàn)在討論求解問題（9.14）-(9.18)的差分方法。首先對求解區(qū)域G=0xl, 0tT進行網(wǎng)格剖分。取空間步長h=l/N,時間步長=T/M,其中N,M是正整數(shù)，作兩族平行直線 將區(qū)域G剖分成矩形網(wǎng)格，見圖9-5，網(wǎng)格交點（xj，tk）稱為節(jié)點。用差分方法求解初邊值問題（9.14）-（9.16）就是要求出精確解u(x,t)在每個節(jié)點（xj，tk）處的近似值。為簡化記號，簡記節(jié)點（xj，tk）=u(j,k)。 利用一元函數(shù)的Taylor展開公式，可推出

11、下列差商表達式 （9.17） （9.18） （9.19） （9.20）1.古典顯格式 在區(qū)域G的內(nèi)節(jié)點(j,k)處，利用公式（9.17）和（9.20），可將偏微分方程（9.14）離散為其中。舍去高階小項，就得到節(jié)點近似值（差分解）所滿足的差分方程 （9.21）顯然，在節(jié)點(j,k)處，差分方程（9.21）逼近偏微分方程（9.14）的誤差為，這個誤差稱為截斷誤差，它反映了差分方程逼近偏微分方程的精度?，F(xiàn)將（9.21）式改寫為便于計算的形式，并利用初邊值條件（9.15）與（9.16）補充上初始值和邊界點方程，則得到 （9.22）其中稱為網(wǎng)比。與時間相關(guān)問題差分方程的求解通常是按時間方向逐層進行的。

12、對于差分方程（9.22），當?shù)趉層節(jié)點值已知時，可直接計算出第k+1層節(jié)點值。這樣，從第0層已知值開始，就可逐層求出各時間層的節(jié)點值。差分方程（9.22）的求解計算是顯式的，無須求解方程組，故稱為古典顯格式。此外，在式（9.22）中，每個內(nèi)節(jié)點處方程僅涉及k和k+1兩層節(jié)點值，稱這樣的差分格式為雙層格式。差分方程（9.22）可表示為矩陣形式 （9.23）其中 2. 古典隱格式 在區(qū)域G的內(nèi)節(jié)點(j,k)處，利用公式（9.18）和（9.20），可將偏微分方程（9.14）離散為舍去高階小項，則得到如下差分方程 （9.24）它的截斷誤差為，逼近精度與古典顯格式相同。改寫（9.24）式為便于計算的形式

13、，并補充上初始值與邊界點方程，則得到 （9.25）與古典顯格式不同，在差分方程（9.25）的求解中，當?shù)趉-1層值已知時，必須通過求解一個線性方程組才能求出第k層值，所以稱（9.25）式為古典隱格式，它也是雙層格式。差分方程（9.25）的矩陣形式為 （9.26）其中 向量同（9.23）式中定義。從（9.26）式看到，古典隱格式在每一層計算時，都需求解一個三對角形線性方程組，可采用追趕法求解。3.Crank-Nicolson格式（六點對稱格式）利用一元函數(shù)Taylor展開公式可得到如下等式 使用這兩個公式，在點離散偏微分方程（9.14），然后利用（9.20）式進一步離散二階偏導數(shù)，則可導出差分方

14、程 （9.27）其截斷誤差為，在時間方向的逼近階較顯格式和隱格式高出一階。這個差分格式稱為Crank-Nicolson格式，有時也稱為六點對稱格式，它顯然是雙層隱式格式。改寫（9.27）式，并補充初始值和邊界點方程得到 （9.28）它的矩陣形式為 （9.29）在每層計算時，仍需求解一個三對角形方程組。4. Richardson格式利用公式（9.19）和（9.20），可導出另一個截斷誤差為階的差分方程稱之為Richardson格式?？筛膶憺?(9.32)這是一個三層顯式差分格式。在逐層計算時，需用到和兩層值才能得到k+1層值。這樣，從第0層已知值開始，還須補充上第一層值，才能逐層計算下去。可采用

15、前述的雙層格式計算。除上述四種差分格式外，還可構(gòu)造出許多逼近偏微分方程（9.14）的差分格式，但并不是每個差分格式都是可用的。一個有實用價值的差分格式應具有如下性質(zhì)：（1）收斂性。對任意固定的節(jié)點（xj，tk）,當剖分步長時，差分解應收斂到精確解u（xj，tk）。（2）穩(wěn)定性。當某一時間層計算產(chǎn)生誤差時，在以后各層的計算中，這些誤差的傳播積累是可控制而不是無限增長的。理論上可以證明，在一定條件下，穩(wěn)定的差分格式必然是收斂的。因此，這里主要研究差分格式的穩(wěn)定性。作為例子，先考查Richardson格式的穩(wěn)定性。設是當計算過程中帶有誤差時，按Richardson格式（9.30）得到的實際計算值。是

16、理論值，誤差。假定右端項的計算是精確的，網(wǎng)比，則滿足 （9.31）設前k-1層計算時精確的，誤差只是在第k層點發(fā)生，即。則利用（9.31）式可得到誤差的傳播情況，見表9-1。表9-1 r=1/2時Richardson格式的誤差傳播 j kj0-4j0-3j0-2j0-1j0j0+1j0+2j0+3j0+4k00000000k+1000-2000k+200-47-400k+30-617-2417-60k+4-831-6889-6831-8k+5-1049-144273-388273-14449-10k+671-260641-10961311-1096641-26071 從表中看出，誤差是逐層無限

17、增長的。表中的計算雖然是就網(wǎng)比進行的，實際上對任何r>0都會產(chǎn)生類似現(xiàn)象，所以Richardson格式是不穩(wěn)定的。利用誤差傳播圖表方法考查差分格式的穩(wěn)定性雖然直觀明了，但只能就具體取定的r值進行，并且也不適用于隱式差分格式。9.2.2 差分格式的穩(wěn)定性前節(jié)構(gòu)造的幾種雙層差分格式都可以表示為如下的矩陣方程形式 （9.32）其中H稱為傳播矩陣。對于顯格式H=A, 隱格式H=B-1,六點對稱格式H=(I+B) -1 (I+A)。一般的三層格式也可以轉(zhuǎn)化為雙層格式。為了討論方便，設在初始層產(chǎn)生誤差，且假定右端項的計算是精確的。用表示當初始層存在誤差時，由差分格式（9.32）得到的計算解，則滿足方

18、程 （9.33）記誤差向量，則滿足方程 （9.34）定義9.1 稱差分格式（9.32）是穩(wěn)定的，如果對任意初始誤差，誤差向量在某種范數(shù)下滿足 （9.35）其中C為與h,無關(guān)的常數(shù)。這個定義表明，當差分格式穩(wěn)定時，它的誤差傳播是可控制的。從（9.34）式遞推得到因此，差分格式穩(wěn)定的充分必要條件是 （9.36）定理9.3 （穩(wěn)定性必要條件）差分格式穩(wěn)定的必要條件是存在與h,無關(guān)的常數(shù)M,使譜半徑 （9.37）定理9.4 （穩(wěn)定性充分條件）設H為正規(guī)矩陣，即,則（9.37）式也是差分格式穩(wěn)定的充分條件。下面討論幾種差分格式的穩(wěn)定性。為便于討論，引進N-1階矩陣 這個特殊矩陣的特征值為 （9.38）例9-1古典顯格式 此時H=A=(1-2r)I+rS。 利用（9.38）式和三角函數(shù)公式，可求得H的特征值為 為使穩(wěn)定性條件式（9.39）成立，必須且只須。由于H=A為實對稱矩陣，所以古典顯格式穩(wěn)定的充分必要條件是網(wǎng)比例9-2 古典隱格式 此時H=B-1,B=(1+2r)I-rS。利用（9.38）式可求得H的特征值為 顯然，對任意r>0，條件（9.37）成立。注意，H=B-1仍為實對稱矩陣，所以古典隱格式對任何網(wǎng)比r>0都是穩(wěn)定的，稱為絕對穩(wěn)定。例9-3 六點對稱格式 此時H=(I+B)-1(I+A)，利用矩陣A和B的特征值可得到矩陣H的特征值為 則對任意r>