1、第9章 偏微分方程的差分方法含有偏導數的微分方程稱為偏微分方程。由于變量的增多和區域的復雜性，求偏微分方程的精確解一般是不可能的，經常采用數值方法求方程的近似解。偏微分方程的數值方法種類較多，最常用的方法是差分方法。差分方法具有格式簡單，程序易于實現，計算量小等優點，特別適合于規則區域上偏微分方程的近似求解。本章將以一些典型的偏微分方程為例，介紹差分方法的基本原理和具體實現方法。9.1橢圓型方程邊值問題的差分方法9.1.1 差分方程的建立最典型的橢圓型方程是Poisson（泊松）方程 （9.1）G是x,y平面上的有界區域，其邊界為分段光滑的閉曲線。當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)都是已知的函數，k(x,y)0。滿足方程（9.1）和上述三種邊值條件之一的光滑函數u(x,y)稱為橢圓型方程邊值問題的解。用差分方法求解偏微分方程，就是要求出精確解u(x,y)在區域G的一些離散節點（xi，yi）上的近似值ui,j(xi，yi）。差分方法的基本思想是，對求解區域G做網格剖分，將偏微分方程在網格節點上離散化，導出精確解在網格節點上近似值所滿足的差分方程，最終通過求解

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

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

5、i-1,j，ui,j+1，ui,j-1，因此通常稱（9.6）式為五點差分格式，當h1= h2=h時，它簡化為差分方程（9.6）中，方程個數等于內節點總數，但未知量除內節點值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）聯立差分方程（9.6）與（9.7）或

6、（9.8）就形成了求解Poisson方程邊值問題的差分方程組，它實質上是一個關于未知量ui,j的線性代數方程組，可采用第2，3章介紹的方法進行求解。這個方程組的解就稱為偏微分方程的差分近似解，簡稱差分解。考慮更一般形式的二階橢圓型方程 （9.9）其中A(x,y)Amin>0, B(x,y) Bmin >0, E(x,y) 0。引進半節點利用一階中心差商公式，在節點（i,j）處可有對類似處理，就可推得求解方程（9.9）的差分方程 （9.10）其中 （9.11）顯然，當系數函數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 一般區域的邊界條件處理前面已假設G為矩形區域，現在考慮G為一般區域情形，這里主要涉及邊界條件的處理。考慮Poisson方程第一邊值問題 （9.12）其中G可為平面上一般區域，例如為曲邊區域。仍然用兩組平行直線：x=x0+ih1,y=y0+jh2,i,j=0,±1,對區域G進行矩形網格剖分，見圖9-3。 如果一個內節點（i,j）的四個相鄰節點（i+1,j），（i-1,j），（i,j+1）和（i,j-1）屬于，則稱其為正則內點，見圖9-3中打“。”號者；如

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

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

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

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

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

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

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

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

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

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