有限差分方法

有限差分方法

正文

一種求偏微分（或常微分）方程和方程組定解問題的數值解的方法，簡稱差分方法。
微分方程的定解問題就是在滿足某些定解條件下求微分方程的解。在空間區域的邊界上要滿足的定解條件稱為邊值條件。如果問題與時間有關，在初始時刻所要滿足的定解條件，稱為初值條件。不含時間而只帶邊值條件的定解問題，稱為邊值問題。與時間有關而只帶初值條件的定解問題，稱為初值問題。同時帶有兩種定解條件的問題，稱為初值邊值混合問題。
定解問題往往不具有解析解，或者其解析解不易計算。所以要採用可行的數值解法。有限差分方法就是一種數值解法，它的基本思想是先把問題的定義域進行格線剖分，然後在格線點上，按適當的數值微分公式把定解問題中的微商換成差商，從而把原問題離散化為差分格式，進而求出數值解。此外，還要研究差分格式的解的存在性和唯一性、解的求法、解法的數值穩定性、差分格式的解與原定解問題的真解的誤差估計、差分格式的解當格線大小趨於零時是否趨於真解(即收斂性)，等等。
有限差分方法具有簡單、靈活以及通用性強等特點，容易在計算機上實現。
偏微分方程初值問題的差分法　許多物理現象隨著時間而發生變化、如熱傳導過程、氣體擴散過程和波的傳播過程都與時間有關。描述這些過程的偏微分方程具有這樣的性質：若初始時刻t=t0的解已給定,則t>t0時刻的解完全取決於初始條件和某些邊界條件。利用差分法解這類問題，就是從初始值出發，通過差分格式沿時間增加的方向，逐步求出微分方程的近似解。
雙曲型方程的差分方法　最簡單的雙曲型方程的初值問題是：

式中嫓(x)為已知初值函式。這初值問題的解是：

u(x,t)＝嫓(x－at)。　(2)

由(2)可見，(1a)(1b)的解(2)當a＞0時代表一個以有限的速度a沿特徵線x-at=常數向右傳播的波,而解u(x,t)在P(慜,慙)點的值完全由嫓(x)在x軸上的點A(慜－а慙,0)的值決定。A點就是雙曲型方程(1a)在P點的依賴域(圖1)。現以初值問題(1)為例介紹初值問題差分方法的基本思想。

①剖分格線　用格線覆蓋(1a),(1b)的定解區域,如圖2所示,在x，t平面的上半部作兩族平行於坐標軸的直線：

x=xj=jΔx，　j＝0,±1,±2,…

　 　

t=tn=nΔt，　n＝0,1,2,…

並稱之為格線線。Δx,Δt分別稱為空間步長和時間步長。格線線的交點(jΔx，nΔt)稱為格點。
②建立差分格式　以下除特別聲明外,總設a＞0,由泰勒公式，有：

，　 (3a)

和

-1<θ2<0(3b)

解出，代入(1a)得：

(4)

式中

(5)

E是微分方程(1a)用它的解在相鄰三個格點(見圖2)上的值的差分來表示的形式。 略去(4)中關於Δx,Δt的高階項E,得到一個較簡單的差分方程,但微分方程的解u(jΔx,nΔt)不再是這方程的解；設這個方程的解是u，u滿足的方程是：

；(6)

式(6)還可寫成：　

；　(6′)

初值條件(1b)此時就是：
u=嫓(xj),　j=0,±1,±2,…。 (7)
差分方程(6)和相應的初值條件(7)合稱差分格式，利用這些格式可逐步算出t=Δt,2Δt,…各時間層的u,u,…，等等。這個把微分方程化為近似的差分方程的過程常稱為離散化。
③差分格式的截斷誤差和相容性　(5)中的E是把微分方程充分光滑的解代入差分方程(6)的結果,它說明微分方程(1a)和差分方程(6)的區別，稱為差分格式(6)的截斷誤差。式(6)的截斷誤差對Δt和Δx都是一階的，寫成O(Δx＋Δt)，因此稱差分格式(6)為一階相容格式。一般說，如果Δx，Δt趨於零，截斷誤差也趨於零,則差分方程與微分方程是相容的。不相容的格式的解不能作為原微分方程的近似解，因而是無用的。方程(1a)的離散化過程也不是唯一的。例如取數值微分公式：

代替微分方程(1a)中的

，

可得另一個差分方程：

，　(9)

它的截斷誤差是O(Δx2＋Δt)階的,也是相容的差分格式，再若用數值微分公式
代替(1a)中的,又得到截斷誤差為O(Δx＋Δt)的相容差分格式：　

。　(11)

但是，並不是每個相容格式都有用。

④差分格式的收斂性　設P(慜,慙)是求解區域中的一點，取步長Δx,Δt使慜＝jΔx，慙＝nΔt,用差分格式算出u,如果當Δx,Δt→0時u－u(慜，慙)→0,便可用步長充分小時的u作為微分方程的解u(jΔx，nΔt)的近似，這種差分格式便是收斂的。
雙曲型微分方程的解，對求解區域內一點(慜，慙)而言，在初值區域內有一個依賴域,差分方程也是如此,對於差分方程(6),點(jΔx,nΔt)的依賴域是初值線上區間【(j－n)Δx,jΔx】。如令Δt/Δx=r=常數,慜=jΔx,慙=nΔt，則差分方程(6)在點(慜,慙)的依賴域為【慜－a慙/r，慜】,並且步長比r固定時，依賴域與Δx，Δt無關。
差分方程(9)在(慜,慙)的依賴域是【慜-a慙/r,慜+a慙/r】，而差分方程(11)的依賴域則是【慜，慜＋│a│慙/r】,R.庫朗等人曾經證明，差分格式收斂的一個必要條件是差分方程的依賴域應包含微分方程的依賴域，這個條件叫作“庫朗條件”。從圖3中可以看到，對於差分方程(6)，這個條件是慜-a慙/r≤慜-a慙≤慜，即對於格式(9)，庫朗條件是

,

兩者不同。對於格式(11),庫朗條件是
　慜≤慜-a慙≤慜+│a│慙/r；
在a＞0時,顯然不能成立,所以格式(11)當a＞0時不收斂， 因而也是無用的。格式(6)a＞0在而庫朗條件滿足時，的確是收斂的。因為的離散化誤差

適合

由此可知：

又因差分格式與微分方程的初值相同，於是可知：

這說明條件

滿足時，格式(6)收斂。
如果a＜0，格式(6)不收斂。但當時，　 格式(11)收斂。這兩個格式稱為“迎風格式”,因為a＞0時，用向後差商代替,往上風取近似值；當a＜0時則用向前差商代替,也是往上風取近似值。可見作(1)的差分格式時，要考慮波的傳播方向。

⑤差分格式的穩定性　用一個差分格式計算時，初值的誤差必然要影響到以後各層。通常希望這誤差的影響不會越來越大，以致完全歪曲了差分方法的真解，這便是穩定性問題。討論時，常把問題化簡，設初值有誤差ε，而以後的計算並不產生誤差,由於誤差ε，使變成了+ε， 但+ε仍滿足所適合的差分格式。定義一種衡量 t=tn層格點上ε的大小的所謂範數‖ε‖，若有常數K＞0 使當 Δt、Δx→0 而0≤t＝nΔt≤T 時，恆有‖ε‖≤K‖ε‖,則稱此差分格式是穩定的。 以格式(6)為例，ε適合差分方程：

可證當時，取，則有

這說明，用格式(6)計算時，若步長比合於庫朗條件,則初值誤差的影響不增長，即使Δt縮小，算到t＝T時,也不再增大，因而格式是穩定的。
對於線性偏微分方程組的穩定性理論,J.von諾伊曼曾用傅立葉分析作了系統研究，把差分方程的解表成諧波的疊加，考察其中一個諧波

　(12)

的增長情況，式中k為實數；G＝G(k,Δt)稱為增長因子。若對於一切諧波，(12)的振幅一致有界，即對一切合於0≤nΔt≤T的n和充分小的Δt都有|Gn|≤K,K為常數,則此差分格式是穩定的。具體地說，對格式(6),把(12)代入(6)，得：

而

故當時,│G│≤1，解的振幅不增加,所以格式(6)是穩定的。
相容性和庫朗條件都不能保證穩定性，例如對格式(9)，把(12)代入，得：

而

故當sinkΔx厵0時，恆有|G|＞1,解的振幅逐層增加,所以雖然格式(9)是相容的格式，並且適合庫朗條件,但它仍是不穩定的，因而也是無用的。
P.D.拉克斯1956年曾證明：對於線性偏微分方程組的適定的初值問題，一個與之相容的線性差分格式是收斂的格式的充分必要條件是這格式的穩定性。
非線性問題沒有相應的等價定理。
拋物型方程的差分方法　拋物型方程的定解問題是初值問題或初值邊值問題。為了說明拋物型方程差分方法的某些特點，考慮熱傳導方程的初值、邊值問題：　

式中a＞0為常數；┃(x)為給定的連續函式。這裡也是用直線
　x=xj=jΔx，　j＝0,1,…,Μ,
　t=tn=nΔt，　n＝0,1,2,…,N(＝T/Δt)，
剖分求解區域為矩形格線(見圖4),式中Δx＝1/Μ,Μ為正整數。利用數值微分公式：

中的(14a)及(14b)，從微分方程(13)可得差分格式：

這裡計算u囃時只用到前一層的u囐，u怹及u囐，是一個顯示格式。要是用(14a)和(14c)，則得到一個隱式差分格式：

此時計算u囃必須解一個線性代數方程組。 用諾伊曼方法，可以證明：當常數時,對於任何r值，格式(16)是無條件穩定的;格式(15)則是條件穩定的,穩定性條件是。顯然格式(15)的時間步長要受到較大的限制。隱式格式歸結為解線性代數方程組。最簡便的方法是追趕法,其要點如下:對於每個時間步n＋1，把(16)改寫成：　

用消去法可導出如下兩套遞推公式：

和　

用公式(18a)從j＝0齣發，逐步先求出αj和βj，然後用(18b)從j＝Μ-1逐步求出u囄,u囅,…,u囀,u囆來。當Aj＞0，Bj＞0，Cj≥Aj+Bj時，這兩個遞推過程都是穩定的。

偏微分方程邊值問題的差分法　物理上的定常問題，如彈性力學中的平衡問題， 亞聲速流、 不可壓粘性流、電磁場及引力場等可歸結為橢圓型方程。其定解問題為各種邊值問題， 即要求解在某個區域D內滿足微分方程，在邊界上滿足給定的邊界條件。橢圓型方程的差分解法可歸結為選取合理的差分格線，建立差分格式，求解代數方程組以及考察差分格式的收斂性等問題。
泊松方程是橢圓型方程的典型例子，它的第一邊值問題為：

式中D為x,y平面上某個封閉區域；дD為它的邊界（圖5);┃(x,y),g(x,y)為連續函式；u(x,y)為未知解。以下簡單介紹其差分解法的基本思想。
通常可將定解區域剖分成矩形格線或三角形格線。三角形格線對不規則區域較為方便。 為簡便起見，設D為單位正方形，x 和y方向均取為等距步長h，並用直線xi=ih，yj＝jh（i，j＝0，1，2，…，N）將此正方形D={0≤x≤y≤1}剖分成正方形格線。

在格點(i,j)上,微商uxx、uyy分別用x、y方向的二階中心差商來代替，得到差分格式：

的截斷誤差是二階的。將u嗎按順序u1,1,u1,2,…，u,u2,1,u2,2,…，排列，並用矢量形式表示為：
u=(u1,1,u1,2,…,u,u2,1,u2,2,…,u)TT代表轉置，則(20)可寫成方程組　

，　(21)

　其中係數矩陣A為三對角塊狀方陣

,

I為(N－1)×(N－1)的單位矩陣,S亦為(N－1)×(N－1)矩陣。
偏微分方程邊值問題的差分方程組的特點是係數矩陣中非零元素很少，即是稀疏矩陣。近年來由於稀疏矩陣技術的發展，解差分方程組時，直接法受到了較多的重視。疊代法是用逐次逼近的方式得到差分方程組的解，它的存儲量小，程式簡單，因此常用於橢圓型差分方程組的求解。疊代方法很多，最基本的有三種：
①同時位移法（也稱雅可比法）：　

，(23)

n代表疊代的次數。
②逐個位移法（也稱賽德耳法）：　

，(24)

當已算出時，取，否則。 ③鬆弛法：

，(25)

當已算出時，，否則， ，式中ω為鬆弛因子；0＜ω＜2,ω＜1為低鬆弛法,ω＞1為超鬆弛法。三個方法中超鬆弛法收斂最快，是常用的方法之一。
差分方法的發展和套用　前面闡述了兩個自變數，線性方程的差分法。實際問題常會遇到多個自變數，非線性的方程或方程組；它們還可能是混合型的偏微分方程（如機翼的跨聲速繞流），其解包含著各種間斷（如激波間斷、按觸間斷等）。非線性問題的差分法求解是十分困難的。隨著電子計算機的發展，在解決各種非線性問題中，差分法得到了很快的發展，並且出現了許多新的思想和方法，如守恆差分格式，時間相關法、分步法等。
守恆差分格式　數學物理偏微分方程通常代表某種物理、力學中的守恆律（如質量守恆、動量守恆、能量守恆、粒子守恆等）。原始問題的差分格式，若能保持同樣的守恆性質，就稱為守恆差分格式。守恆性反映出物理問題的整體性質，用它來檢驗差分格式的好壞是合理的。對於間斷的問題，守恆格式特別重要。從積分守恆關係式出發，利用積分插值方法容易得到守恆格式。這時對於複雜的求解區域、各種類型邊界條件、間斷係數等複雜情況都可以處理。
時間相關法　把定常的微分問題用一個相應的非定常問題來代替，然後用差分法解後者的初值問題，要求當t→∞時,它的穩定解為原來問題的解，這類方法叫作時間相關法。實踐上，當計算時間足夠大時，就能得到滿足給定精度的近似解。例如拉普拉斯方程第一邊值問題：

可以用熱傳導方程的初邊值問題：

(27)

來代替。若用顯式格式計算(27)，可避免解大型代數方程組。特別是當微分方程的類型在定解區域內發生變化時，可只用一種類型來算，而使問題大大化簡。這種方法在定常問題中廣泛使用。缺點是達到定常解的計算時間較長，有待改進。
分步法　把複雜的問題的每一時間步分解成幾個中間步，例如把多維問題按坐標分解為幾個一維問題，然後用差分法解這些比較簡單的各中間步，最後得到原始問題的近似解，這類方法叫作分步法。交替方向法、預估－修正法、時間分裂法、因式分解法等都屬此類。以二維拋物型方程定解問題：

為例，用顯式格式求解，時間步長受穩定性條件：

的限制，用隱式格式,則歸結為大型線性代數方程組,解起來比較麻煩。1955年皮斯曼－拉什福德提出交替方向隱式格式：

　(i＝1,2,…,N－1;j＝1,2,…,Μ－1;n＝0,1,2,…)δ為中心差分算符，第一步x方向取隱式，y方向取顯式，第二步則相反。兩步合成無條件穩定的格式。由於每一步可用追趕法求解，大大簡化了解法。交替方向法出現後，進一步發展了各種形式的分步格式，並可推廣到任何維數的方程或方程組的情形，困難在於邊界條件的處理。
有限差分方法已成為解各類數學物理問題的主要數值方法，也是計算力學中的主要數值方法之一。有些解偏微分問題的方法（如特徵線法、直線法）實質上也是差分方法的一種形式。在固體力學中，有限元方法出現以前，主要採取差分方法；在流體力學中，差分方法仍然是主要的數值方法。當然，對於某些具有複雜的幾何形狀及複雜的流動現象的實際問題，差分方法還有待進一步發展。
參考書目
　馮康等編：《數值計算方法》，國防工業出版社，北京，1978。
　胡祖熾編:《計算方法》,高等教育出版社,北京,1959。
　清華大學、北京大學《計算方法》編寫組編：《計算方法》，科學出版社，北京，1980。
　朱幼蘭等著：《初邊值問題差分法及繞流》，科學出版社，北京，1980。
　R.D.里奇特邁爾著，何旭初等譯：《初值問題差分方法》， 科學出版社， 北京， 1966。 (R.D.Richtmyer,Difference　Methods　for　Initial-Value Problems, Interscience Pub., New York,1957.)
　R.D.Richtmyer, K.W.Morton, Difference　Methods for Initial-Value Problems, 2nd ed.,Interscience Pub.,New York,1967.

配圖

相關連線