二維拋物方程的有限差分法

本文格式為Word版，下載可任意編輯——二維拋物方程的有限差分法華北電力大學本科畢業(yè)設計（論文）

有限元方法(FiniteElementMethods)的基礎是變分原理和分片多項式插值。該方法的構(gòu)造過程包括以下三個步驟。首先，利用變分原理得到偏微分方程的弱形式(利用泛函分析的知識將求解空間擴大)。其次，將計算區(qū)域劃分為有限個互不重疊的單元(三角形、四邊形、周邊體、六面體等)。再次，在每個單元內(nèi)選擇適合的節(jié)點作為求解函數(shù)的插值點，將偏微分方程中的變量改寫成由各變量或其導數(shù)的節(jié)點值與所選用的分片插值基函數(shù)組成的線性表達式，得到微分方程的離散形式[1]。

有限體積法(FiniteVolumeMethods)又稱為控制體積法。其基本思路是：將計算區(qū)域劃分為一系列互不重疊的控制體，并使每個網(wǎng)格點周邊有一個控制體；將待求解的微分方程對每一個控制體積積分，便得出一組離散方程。該方法的未知量為網(wǎng)格點上的函數(shù)值。為了求出控制體積的積分，須假定函數(shù)值在網(wǎng)格點控制體邊界上的變化規(guī)律。從積分區(qū)域的選取方法看來，有限體積法屬于有限元方法中檢驗函數(shù)取分片常數(shù)插值的子區(qū)域法；從未知量的近似方法看來，有限體積法屬于采用局部近似的離散方法[1]。

有限單元法的基本求解思想是把計算域劃分為有限個互不重疊的單元，在每個單元內(nèi)，選擇一些適合的節(jié)點作為求解函數(shù)的插值點，借助于變分原理或加權余量法，將微分方程離散求解[1]。無網(wǎng)格方法的近似函數(shù)建立在一些離散的節(jié)點上，不需要借助網(wǎng)格，在涉及到網(wǎng)格畸變，網(wǎng)格移動等問題中顯示出了明顯的優(yōu)勢[1]。

1.2.2有限差分方法的發(fā)展

隨著人們對顯示格式的進一步研研究，提出了好多新的高精度，絕對穩(wěn)定的差分格式。文獻2探討了變系數(shù)熱傳導方程的兩層絕對穩(wěn)定差分格式。研究了二維變系數(shù)非齊次熱傳導方程的兩層絕對穩(wěn)定的差分格式問題。首先運用Pade迫近導出了差分格式,給出了差分格式的截斷誤差;探討了差分格式的絕對穩(wěn)定性和收斂性,且收斂階為O(k2?h4);最終給出了數(shù)值例子,數(shù)值結(jié)果和理論結(jié)果是吻合的。

文獻3第一部分用待定系數(shù)法對P維拋物型方程(P?1,2,3,4)構(gòu)造出了高精度(截斷誤差達O(k2?h4))能顯式計算，穩(wěn)定性較好(r?1/2)的三層(特別狀況下可以是兩層)顯式差分格式，在空間變量更多的狀況下，指出了構(gòu)造高精度差分格式的一般方法。

文獻3其次部分用算子方法對二維和三維拋物型方程構(gòu)造出了高精度(截斷誤差達

O(k2?h4))的絕對穩(wěn)定的交替方向隱式差分格式，在空間變量更多的狀況下，也指出了構(gòu)造高精度交替方向隱式差分格式的一般方法。并附有數(shù)值例子，它說明理論分析的正確性和所建立的差分格式的有效性。

文獻4對二階拋物型方程構(gòu)造了一族含參數(shù)高精度三層差分格式。當參數(shù)滿足一定的條件時,差分格式絕對穩(wěn)定，其局部截斷誤差階數(shù)最高可達O(k2?h4)。適當?shù)卣{(diào)理參數(shù)，可以得到一個七點顯式差分格式和一個兩層六點隱格式。數(shù)值例子說明,對穩(wěn)定性所作的分析是正確的。

文獻5多維拋物型方程和雙曲方程的差分解法對一般m維熱傳導方程，波動方程的初、

2

華北電力大學本科畢業(yè)設計（論文）

邊值問題，提出若干個交替方向的差分格式。

1.3差分格式建立的基礎[6]

1.3.1區(qū)域剖分[14]

構(gòu)造一維線性拋物方程

（x,t）??u??u?u?(（ax,t）)?（bx,t）?（cx,t）u(1-2)?t?t?x?x的有限差分迫近，首先將求解區(qū)域?用兩組平行于x軸和t軸的直線構(gòu)成網(wǎng)格覆蓋，網(wǎng)格邊長在x方向為?x?h，t方向為?t?k，網(wǎng)格節(jié)點為(xm,tn)，xm?mh,tn?nk，m,n為整數(shù)[6]。網(wǎng)格如圖1-1。

圖1-1一維拋物方程的網(wǎng)格剖分

二維拋物方程的區(qū)域剖分將剖分區(qū)域擴展到三維空間。首先將求解區(qū)域V用三組平行于x軸，y軸和t軸的直線構(gòu)成網(wǎng)格覆蓋，網(wǎng)格邊長在x方向為?x?h，y方向為

?y?h，t方向為?t?k，網(wǎng)格節(jié)點為(xl,ym,tn)，xl?lh,ym?mh,tn?nk，l,m,n為整數(shù)[3]。

網(wǎng)格如圖1-2[6]。

圖1-2二維拋物方程的網(wǎng)格剖分

1.3.2差商代替微商

差分方法的基本思想就是以差商代替微商?？紤]如下兩個Taylor公式：1113n?n1)?u(t)??u(t)?h??u(2)t?h???u()?t?h?n()u()t?hu(t?h(Oh)(1-3)

2!3!n!3

華北電力大學本科畢業(yè)設計（論文）

u(t?h)?u(t)?u?(t)h?111u??(t)h2?u???(t)h3???u(n)(t)hn?O(hn?1)(1-4)2!3!n!u?(ti)?u(ti?1)?u(ti)?O(h)h從式(1-3)得到：

從式(1-4)得到：

u?(ti)?

u(ti?1)?u(ti)?O(h)h從式(1-3)-式(1-4)得到：

u?(ti)?u(ti?1)?u(ti?1)?O(h2)2h從式(1-3)+式(1-4)得到：

u(ti?1)?2u(ti)?u(ti?1)2??u(t)??O(h)i2h下表1-1面介紹常用的幾種差商[6]。

表1-1：常用的幾種差商形式

有限差分近似公式ur?ur?1（向前差分）hur?1?ur（向后差分）hur?1/2?ur?1/2（中心差分）2h

誤差階O(h)O(h)O(h2)由式(1-1)可得：

n?1n(1-5)um?exp(kL)um1.3.3差商代替微商格式的誤差分析

為了分析分析數(shù)值方法的確切度，考察收斂性。往往在ui?u(ti)成立的假定下，估計誤差ei?1?u(ti?1)?ui?1這種誤差稱為“局部截斷誤差〞局部截斷誤差是以點ti的確切解u(ti)為出發(fā)值，用數(shù)值方法推進到下一個點ti?1而產(chǎn)生的誤差[8]。若如下：

h2h2u(ti?1)?u(ti)?hu?(ti)?u??(?)?u(ti)?hf(ti,u(ti))?u??(?)(1-6)

22h2ei?1?u(ti?1)?ui?1?u??(?)?O(h2)(1-7)

24

華北電力大學本科畢業(yè)設計（論文）

局部截斷誤差是關于hn的微小函數(shù)，則收斂[13]。

為了分析分析數(shù)值方法的確切度，考察收斂性。整體截斷誤差是以點t0的初始值的偏差：?i?1?u(ti?1)?ui?1稱為整體截斷誤差[8]。若如下：

為出發(fā)值，用數(shù)值方法推進i+1步到點ti?1，所得的近似值ui?1與確切值u(ti?1)

ui?1?ui?hf(ti,ui)

?i?1??i?h[f(ti,u(ti)?f(ti,ui)]?Dh2

h2ei?1?u(ti?1)?ui?1?u??(?)?O(h2)

2?i?1??i?h[f(ti,u(ti)?f(ti,ui)]?Dh2

?i?1??i?hL?i?Dh2

?i?1?(1?hL)i?1?0?Dh2[1?(1?hL)?(1?hL)2???(1?hL)i]

?(1?hL)i?1Dh2?0?[(1?hL)i?1?1][14]

hL整體截斷誤差不隨?0無限擴大，則收斂。

1.4本文主要研究內(nèi)容

本文主要研究二維拋物方程的有限差分方法，研究工作分為以下四個方面：1)差分格式的構(gòu)造方法

有限差分法法解二維拋物方程包括區(qū)域剖分和差商代替導數(shù)兩個過程。差分格式就是在網(wǎng)格結(jié)點上求出微分方程的近似解的一種方法，因此又稱為網(wǎng)格法。

2)差分格式的分析方法

數(shù)值方法是近似計算方法，需從收斂性，相容性，穩(wěn)定性方面考察。收斂性是考察步長足夠小時，數(shù)值解是否無限迫近解析解。穩(wěn)定性主要最初產(chǎn)生的小誤差在以后的計算中雖然會傳遞下去，但不會無限制的擴大。

3)顯式差分格式

顯式差分格式是差分方法中可逐層逐點分別求解的格式，是一種有限差分近似方法。顯式差分格式的優(yōu)點在于計算簡單，但它并不總是穩(wěn)定的。

4)隱式差分格式

與顯式差分格式不同，隱式差分格式中包括了（n+1）時間層上二個或二個以上結(jié)點

n?1n?1n?1處的未知值（例如Um，使用隱式差分格式和使用顯式差分格式求解完全不同。,U,U?1mm?1）

相對而言，使用隱式差分格式求解，每時間層包含有較多的計算工作量。從后面對差分格式的穩(wěn)定性分析可知，隱式格式的優(yōu)點在于，其穩(wěn)定性要求對步長比的限制大為放寬，而

5

華北電力大學本科畢業(yè)設計（論文）

這正是我們所期望的。

6

華北電力大學本科畢業(yè)設計（論文）

2顯式差分格式

現(xiàn)在，對二維拋物方程式(1-1)，從方程式(1-4)出發(fā)，構(gòu)造有限差分的顯式格式。由于一維熱傳導方程

2?u2?u?a(2-1)?t?x2

二維熱傳導方程

2?u?2u2?u?a(2?2)(2-2)?t?x?y在熱傳導，磁擴散等大量領域有重要的應用。實際導熱問題必然涉及邊值條件，在有限差分法中它們也必需差分化。因此，我們需要研究的不僅是差分方程本身，而且是包括全部內(nèi)部區(qū)域和所有邊界上的差分方程所組成的代數(shù)方程組。又稱為差分格式[3]。

2.1常系數(shù)熱傳導方程的古典顯式格式

首先考慮熱傳導方程的邊值問題

2??u?2u2?u??t?a（?x2??y2)???u(0,y,t)?u(L,y,t)?u(x,0,t)?u(x,L,t)?0?u(x,y,0)?f(x,y)???離散化

u0j，m?u(0,m,jk)?uLj，m?u(L,m,jk)?0

jjul，?u(l,0,jk)?u0l，M?u(M,m,jk)?0

ul0,m?u(lh,mh,0)?f(lh,mh)?fl,m[15]

2.1.1古典顯式格式格式的推導

現(xiàn)在對熱傳導方程式(2-1)推導其最簡單的顯隱式差分迫近——古典顯隱式格式。由

?122uln,m?exp(kDx?kDy)uln,m

故

7

華北電力大學本科畢業(yè)設計（論文）

11?12244uln,m?uln,m(1?kDx?kDy?k2Dx?k2Dy??)

2211222式中右邊假使僅保存二階導數(shù)項，且以2?x2替代Dx，2?y替代Dy，則得差分格式

hhkk?1Uln,m?Uln,m(1?2?x2?2?y2)

hh或者

?1Uln,m?(1-4r)Uln,m?r(Uln?1,m?Uln?1,m?Uln,m?1?Uln,m?1)(2-3)

r?h/k2

這是一個顯式格式(四點格式)，如圖(2-1)所示。

圖2-1：古典顯式格式

格式(2-3)應用在一維熱傳導問題中，得到古典顯式格式為：

n?1nnnUm?rUm?1?(1?2r)Um?rUm?1(2-3)

每一層各個節(jié)點上的值是通過一個方程組求解得到的。顯一隱格式區(qū)域分解方法就是以顯格式計算出相鄰子區(qū)域相交內(nèi)邊界的近似值的一種方法。顯一隱格式區(qū)域分解方法綜合了二者的優(yōu)點，借助前一層數(shù)值解的信息，用顯格式給出在這—層的子問題的未知內(nèi)邊界條件，把一個整體區(qū)域上的問題化為若干個子區(qū)域上的子問題，在每個子區(qū)域上用隱式方法求解，從而實現(xiàn)了并行。這可以從下面的計算過程看出來[8]。

2.1.3古典顯式格式的算法步驟

古典顯式格式，以m=1為例，列成矩陣有如下形式：

100000?u1,1?(1?4r)u1,1?r(u0,1?u2,1?u1,0?u1,2)?100000?u2,1?(1?4r)u2,1?r(u1,1?u3,1?u2,0?u2,2)?100000?u3,1?(1?4r)u3,1?r(u2,1?u4,1?u3,0?u3,2)???00000?u1L?1,1?(1?4r)uL?1,1?r(uL?2,1?uL,1?uL?1,0?uL?1,2)?系數(shù)矩陣為

8

華北電力大學本科畢業(yè)設計（論文）

r?1?4r???r1?4rr???????r1?4rr???r1?4r???為了得到二維問題的誤差估計，我們做一些改動。因此，這個算法從un到un?1的計算步驟為：

第一步，根據(jù)給出的初邊值條件，得出t=0時u0；

?1其次步，根據(jù)古典顯式格式的公式，Uln,m?(1-4r)Uln,m?r(Uln?1,m?Uln?1,m?Uln,m?1?Uln,m?1)，

由第un得出un?1[6]。

下文中將通過一個具體的數(shù)值例子說明白計算的方法，表達了這種格式的實用性和優(yōu)越性。

9

華北電力大學本科畢業(yè)設計（論文）

3隱式差分格式

與顯式差分格式不同，隱式差分格式中包括了(n+1)時間層上二個或二個以上結(jié)點處的

n?1n?1n?1未知值(例如Um,U,U?1mm?1)，使用隱式差分格式和使用顯式差分格式求解完全不同。相

對而言，使用隱式差分格式求解，每時間層包含有較多的計算工作量。從后面對差分格式的穩(wěn)定性分析可知，隱式格式的優(yōu)點在于，其穩(wěn)定性要求對步長比的限制大為放寬，而這正是我們所期望的[6]。

對二維拋物方程式(1-1)，從方程式(1-4)出發(fā)，構(gòu)造有限差分的顯式格式。由于熱傳導方程式(2-1)構(gòu)造差分格式。我們需要研究的不僅是差分方程本身，而且是包括全部內(nèi)部區(qū)域和所有邊界上的差分方程所組成的代數(shù)方程組。又稱為差分格式。

3.1古典隱式格式

現(xiàn)在對熱傳導方程推導其最簡單的隱式差分迫近——古典隱式格式。由

?122uln,m?exp(kDx?kDy)uln,m

故

22n?1nexp(?kDx?kDy)um?umun?1m124124?u(1?kD?kD?kDx?kDy??)22nm2x2y式中右邊假使僅保存二階導數(shù)項，且以則得差分格式

121212222???替代，且以替代，替代，DDDxxyxxy222hhhk2k2?x?2?y)h2hn?1nUm?Um(1?或者

?11n?1n?1n?1n(1?4r)Uln,m?r(Uln??1,m?Ul?1,m?Ul,m?1?Ul,m?1)?Ul,m(3-1)

r?h/k2

將格式(2-3)應用于一維熱傳導方程，古典顯式格式為：

n?1nnnUm?rUm?1?(1?2r)Um?rUm?1(3-2)

顯式與隱式的比較如下[6]：

(1)同一階數(shù)下，隱式的局部截斷誤差的系數(shù)的絕對值比顯式的要?。?2)顯式的計算工作量比隱式的?。?/p>

10

華北電力大學本科畢業(yè)設計（論文）

(3)隱式的穩(wěn)定范圍比顯式的大。

局部截斷誤差的階最高是3，式(3-2)是允許函數(shù)任意變化狀況下截斷誤差最小的二階方法。要再提高階就必需增加計算函數(shù)值的次數(shù)。上述式(3-2)又稱為古典隱式格式。

故格式用圖3-1表示，其截斷誤差階為O(k2?h2),與古典差分格式一致。這是一個四點差分格式，如圖3-1所示。

圖3-1：古典隱式格式

為了求得第(n+1)時間層上的Un?1的值，必需通過解線性代數(shù)方程組。

這是一個隱式差分格式，必需聯(lián)合其初邊值條件求解。格式(3-2)尋常稱為古典隱式格

22式。我們也可以通過直接用差分算子代替Dx,Dx,Dy,Dy的方即：

n?1n?un?1ul,m?ul,m()l,m??tkn?1n?1n?1?2un?1ul,m?1?2ul,m?ul,m?1(2)l,m??yh2n?1n?1n?1?2un?1ul?1,m?2ul,m?ul?1,m(2)l,m??xh2代入微分方程，得到格式(3-1)。

古典隱式格式的方程組和矩陣形式如下：

當知道第n層上的uij時，要確定第n+1層上各點值uij?1必需通過求解一個線性代數(shù)方程組。以m=1為例：

j?1j?1j?1j?1j?1j?(1?4r)u1,1?r(u0,1?u2,1?u1,0?u1,2)?u1,1?j?1j?1j?1j?1j?1j?(1?4r)u2,1?r(u1,1?u3,1?u2,0?u2,2)?u2?j?1j?1j?1j?1j?1j?(1?4r)u3,1?r(u2,1?u4,1?u3,0?u3,2)?u3,1???j?1j?1j?1j?1j?1j?(1?4r)uN?1,1?r(uN?2,1?uN,1?uN?1,0?uN?1,2)?uN?1,1?其系數(shù)矩陣如下：

11

華北電力大學本科畢業(yè)設計（論文）

?r?1?4r???r1?4r????????r?????r1?4r?r??r1?4r??3.2Crank-Nicolson隱式格式

Crank-Nicolson隱式差分格式是解熱傳導方程(2-1)的常用的差分格式，近年來，關于拋物方程的區(qū)域分解算法作為并行計算的一種有效工具，吸引了好多學者的注意。這類算法的主要困難在于；如何定義內(nèi)邊界點的值和在子區(qū)域上選取合理的計算解去近似。為了推導它，由式（1-4），有

11?1exp(?kL)uln,m?exp(kL)uln,m

22由

22L?Dx?Dy得

121211221122121211221122n?1[1?kDx?kDy?(kDx)?(kDy)??]um?[1?kDx?kDy?(kDx)?(kDy)222222222222??]uln,m(3-3)兩邊僅保存前二項，用

121222??y替代Dy代替，，則得差分格式Dxx22hh111?1(1?r?x2?r?y2)Uln,m?(1?r?y2)Uln.m

222這是一個隱式差分格式，稱為Crank-Nicolson差分格式，截斷誤差階為?(k2?h2)，也可寫為

11?1?1n?1n?1n?1nnnnn(1?2r)Uln,m?r(Uln,m?U?U?U)?(1?2r)U?r(U?U?U?U?1l,m?1l?1,ml?1,mml,m?1l,m?1l?1,ml?1,m)

44(3-4)

由于格式（3-4）中包括六個結(jié)點，故也可稱為六點格式（如圖3-2所示）。

12

華北電力大學本科畢業(yè)設計（論文）

圖3-2：Crank-Nicolson隱式格式

也可將

n?1nu?u?un?1l,ml,m()l,m2??tkn?1n?1n?1nnnu?2u?uu?2u?u?2un?11l?1,ml,ml?1,ml?1,ml,ml?1,m(2)l,m2?(?)22?x2hhn?1n?1n?1nnnu?2u?uu?2u?u?2un?11l,m?1l,ml,m?1l,m?1l,ml,m?1(2)m2?(?)22?y2hh代入微分方程(2-1)，得到Crank-Nicolson格式。解微分方程(2-1)，根據(jù)Crank-Nicolson格式得到的方程組：11?1?1n?1n?1n?1n(1?2r)Uln,m?r(Uln,m?U?U?U)?(1?2r)U?r(Uln,m?1?Uln,m?1?Uln?1,m?Uln?1,m)?1l,m?1l?1,ml?1,ml,m44其矩陣表達式為：

j?1j????uu1?r?r/21?rr/211??1，??1，???uj?1????uj???r/21?rr/2?11??r/21?r?r/2??2，??2，????????????????j?1????j?????r/21?r?r/2??uN?1，r/21?rr/2??uN?1，?11?j?1??????j????r1?2r?r/21?2r?????1?1??uN，?uN，3.3Douglas差分格式[6]

基于宛如Crank-Nicolson格式一樣的六個網(wǎng)格結(jié)點可獲得另一精度較高的差分格式，如

22在前式(3-3)中僅保存直到Dx，Dy的項，即有：

112n?1112n(1?kDx2?kDy)ul,m?(1?kDx2?kDy)ul,m

2222可令

Dx2uln,m?2nDyul,m1212n?(1??x)ul,m2xh12

1212n?2?y(1??y)ul,mh12則可得：

1212?1?(1??x)xh212

112Dy?2?y2(1??y2)?1h12Dx2?代入上式，則有如下差分格式

13

華北電力大學本科畢業(yè)設計（論文）

11111111?1[1?(r?)?x2?(r?)?y2]Uln,m?[1?(r?)?x2?(r?)?y2]Uln,m(3-5)

26262626它稱為Douglas差分格式，具有截斷誤差階O(k2?h4)。精度較高的差分格式(Douglas格式),并通過一個具體的數(shù)值例子說明白計算的方法,表達了這種格式的實用性和優(yōu)越性。

3.4加權六點隱式格式[9]

前面，我們已經(jīng)推導了熱傳導方程(2-1)的古典顯式顯示格式，古典隱式格式及Crank-Nicolson格式等。實際上，它們都可以作為本節(jié)推導的加權六點隱式格式的特別情形。由

?122uln,m?exp(kDx?kDy)uln,m

得

22?122exp(??kDx??kDy)uln,m?exp[(1??)kDx?(1??)kDx]uln,m,0???1

即

111124?124[1??kDx2??kDy??2k2Dx4??2k2Dy??]uln,m?[1?(1??)kDx2?(1??)kDy?(1??)2k2Dx?22224(1??)2k2Dy??]uln,m

兩邊去掉高于二階導數(shù)的項，且用

121222??代替，替代，則得差分格式DDxyxy22hh2?12(1??r?x2??r?y)Uln,m?[1?(1??)r?x2?(1??)r?y]Uln,m

或者

?1nnn?1n?1n?1n?1(1?2r?)Uln,m?r(1??)(Uln,m?1?Uln,m?1?Um,l?1?Um,l?1)?r?(Ul,m?1?Ul,m?1?Um,l?1?Um,l?1)?[1?

2r(1??)]Uln,m,0???1(3-5)這是一個六點差分格式（如圖3-3所示），稱為加權六點差分格式。

?1??

圖3-3：加權六點差分格式

14

華北電力大學本科畢業(yè)設計（論文）

3.5交替方向隱式格式

交替方向隱式格式是關于x方向的顯格式，關于y方向的隱格式。在n+1/2層上用古典顯式格式計算出過度值ujn?12，再在n+1層上用古典隱格式校正預計值。得到的穩(wěn)定性

和收斂性是相近的。但是，只探討了關于常系數(shù)熱傳導方程的在z和可方向都是一致的時間步長的情形，所以我們在這里試圖做一些擴展。我們首先考慮在內(nèi)邊界點用一個方向為顯格式，另一個方向為隱格式的情形，再考察在z和y方向都用顯格式的情形，而且我們探討在z和可方向用不同的時間步長的情形[6]。

3.5.1Peaceman-Rachford格式

r?h/k2考慮二維拋物方程(2-1)的差分解法，以上對顯式格式的穩(wěn)定性分析發(fā)現(xiàn)，

的要求比一維情形更苛刻。分析說明，維數(shù)越高時，要求時間步長越小，計算工作量越大[3]。

(P-R)ADI格式是由Peaceman-Rachford在1995年發(fā)表的，從第n層到n+1層，P-R格式分兩步進行，每一步只需解具有三對角系數(shù)的線性方程組[6]。P-R格式為

(1)?1/2Ul*,n?Uln,mmk/2?1?1/2Uln,m?Ul*,nm?12*n?12n(?U??Ul,m)(3-6)xl,my2h12*n?1/22n?1(?xUl,m??yUl,m)(3-7)2h(2)k/2?nnn，?yul,m?ul?1/2,m?ul?1/2,m?x，?y是中心差算子，?xuln,m?uln,m?1/2?unl,m?1/2P-R格式對任意r>0無條件穩(wěn)定。但預計P-R格式對三個變量空間問題，無條件穩(wěn)定性不再成立[3]。

3.5.2Rachford-Mitchell格式

Douglas和Rachford在1956年提出另一個隱式差分格式，即Douglas-Rachford格式

[6]

。D-R格式是第一個能被推廣到三維情形的交替方向隱式格式，二維D-R格式是

?1nUl*,nm?Ul,mk?1?1Uln,m?Ul*,nm??1*n?1*n?1Ul*?n1,m?2Ul,m?Ul?1,mh21n?1n?1Uln??1,m?2Ul,m?Ul?1,m?Uln,m?1?2Uln,m?Uln,m?1h2Uln,m?1?2Uln,m?Uln,m?1h2(3-7)

k?h2?(3-8)

3.5.3Mitchell-Fairweather格式

Mitchell和Fairweather在1964年推導了一個高精度ADI差分格式，稱為M-F格

15

華北電力大學本科畢業(yè)設計（論文）

式[6]。M-F格式截斷誤差達O(k2?h4)較之P-R格式和D-R格式更好，且無條件收斂[3]。

3.5.4交替方向隱式格式的算法步驟

ADI格式的算法步驟與古典顯式格式的算法步驟相像，只是ADI格式是隱式格式，需用追趕法求解。追趕法是適用于三對角矩陣的線性方程組求解的方法，并不適用于其他類型矩陣。

第一步，根據(jù)給出的初邊值條件，得出t=0時u0;

其次步，根據(jù)各個交替方向隱式格式中的第一個公式，由第un求得出un?1/2如式(3-6)，

1采用追趕法，確定追趕因子a=1+r/2,b=r/2,c=1+r/2,d=(1?2r)Uln,m?r(Uln,m?1?Uln,m?1?Uln?1,m,

4nnn0n+1/2n+1/2n+1/2=d-*;?Ul?1,m)wm=b-a*u,gm=(d-u*a)/w,然后根據(jù)追趕因子確定u。uM=gM,umgu?1mm第三步，再由追趕法由第un?1/2求得出un?1[6]。

由該算法，計算un的近似解的過程。下文中將通過一個具體的數(shù)值例子說明白計算的方法,表達了這種格式的實用性和優(yōu)越性。

16

華北電力大學本科畢業(yè)設計（論文）

華北電力大學本科畢業(yè)設計（論文）

4實例分析與結(jié)果分析

4.1算例

4.1.1已知有確切解的熱傳導問題

例1是二維Dirichlet邊值問題的熱傳導方程

??T?2T?2T?2?2(0?x?1;0?y?1;0?t?T)??t?x?y??(0?x?1;0?y?1)?T(x,y,0)?0?T(0,y,t)?T(1,y,t)?T(x,1,t)?T(x,0,t)?1(0?t?T)???方程的確切解為

exp(??2(k2?l2)t)T(x,y,t)?1?2??sin(k?x)sin(l?y)

?k?1,3,?l?1,3,?kl16??采用古典顯示格式，以求解得t?0.1，h?1/20，r?h/k2?1/8時，T關于x，y的函數(shù)如圖4-1所式。

2圖4-1：古顯格式解例1，T關于x，y的圖像（t?0.1，h?1/20，r?h/k?1/8）

分別以古典顯示格式，P-R(ADI)格式，以一致的步長(h?1/20，r?h/k2?1/8)求解

17

華北電力大學本科畢業(yè)設計（論文）

同一處(t?0.1)的T與解析解間的差值，列在誤差表中如表4-1，表4-1所示。

表4-1古顯格式解例1的誤差（t?0.1，h?1/20，r?h/k2?1/8）

誤差y=0y=0.1y=0.2y=0.3y=0.4x=000000x=0.10-0.00013-0.00025-0.00035-0.00041x=0.20-0.00025-0.00048-0.00066-0.00077x=0.30-0.00035-0.00066-0.0009-0.00106x=0.40-0.00041-0.00077-0.00106-0.00125x=0.50-0.00043-0.00081-0.00112-0.00131x=0.60-0.00041-0.00077-0.00106-0.00125x=0.70-0.00035-0.00066-0.0009-0.00106x=0.80-0.00025-0.00048-0.00066-0.00077x=0.90-0.00013

-0.00025

-0.00035

-0.00041

x=100000續(xù)表4-1古顯格式解例1的誤差（t?0.1，h?1/20，r?h/k2?1/8）

誤差y=0.6y=0.7y=0.8y=0.9y=1x=000000x=0.1-0.00041-0.00035-0.00025-0.000130x=0.2-0.00077-0.00066-0.00048-0.000250x=0.3-0.00106-0.0009-0.00066-0.000350x=0.4-0.00125-0.00106-0.00077-0.000410x=0.5-0.00131-0.00112-0.00081-0.000430x=0.6-0.00125-0.00106-0.00077-0.000410x=0.7-0.00106-0.0009-0.00066-0.000350x=0.8-0.00077-0.00066-0.00048-0.000250x=0.9-0.00041-0.00035-0.00025-0.000130x=100000

表4-2P-R(ADI)格式解例1的誤差（t?0.1，h?1/20，r?h/k2?1/8）誤差y=0y=0.1y=0.2y=0.3y=0.4x=000000x=0.108.86E-071.19E-061.47E-066.44E-07x=0.201.19E-062.05E-062.19E-069.66E-07x=0.301.47E-062.19E-061.48E-069.21E-07x=0.406.44E-079.66E-079.21E-07-8.6E-07x=0.501.18E-069.93E-073.37E-09-1.4E-06x=0.606.44E-079.66E-079.21E-07-8.6E-07x=0.701.47E-062.19E-061.48E-069.21E-07x=0.801.19E-062.05E-062.19E-069.66E-07x=0.908.86E-071.19E-061.47E-066.44E-07x=100

000

18

y=0.50-0.00043-0.00081-0.00112-0.00131-0.00138-0.00131-0.00112-0.00081-0.00043

0

y=0.501.18E-069.93E-073.37E-09-1.4E-06-1.4E-06-1.4E-063.37E-099.93E-071.18E-060

華北電力大學本科畢業(yè)設計（論文）

2續(xù)表4-2P-R(ADI)格式解例1的誤差（t?0.1，h?1/20，r?h/k?1/8）

誤差x=0x=0.1x=0.2x=0.3x=0.4x=0.5x=0.6x=0.7x=0.8x=0.9x=1

y=0.606.44E-079.66E-079.21E-07-8.6E-07-1.4E-06-8.6E-079.21E-079.66E-076.44E-070y=0.701.47E-062.19E-061.48E-069.21E-073.37E-099.21E-071.48E-062.19E-061.47E-060

y=0.801.19E-062.05E-062.19E-069.66E-079.93E-079.66E-072.19E-062.05E-061.19E-060

y=0.908.86E-071.19E-061.47E-066.44E-071.18E-066.44E-071.47E-061.19E-068.86E-070y=100000000000

4.1.2未知確切解的熱傳導問題

例2是二維Dirichlet邊值問題的熱傳導方程。

??u?2u?2u??2?2?t?x?y???u(x,y,0)?sin(x)?u(0,y,t)?u(?,y,t)?u(x,?,t)?u(x,?,t)?0???0?x,y??;0?t?1采用古典顯式格式，以求解得t?1，h??/20，r?h/k2?1/5時，u關于x，y的函數(shù)如圖4-2所式。

2圖4-2：古顯格式解例2，u關于x,y的圖像（t?1，h??/20，r?h/k?1/5）

19

華北電力大學本科畢業(yè)設計（論文）

采用P-R格式，t由0變化到1時，不同的t,u關于x，y的圖像（t?1，h??/20，

r?h/k2?1/8）對比圖如圖4-3。

2圖4-3：P-R格式解例2，t變化u關于x，y的圖像（h??/20，r?h/k?1/5）

4.2結(jié)果分析

由例1，可看出，古顯格式和P-R格式精度都較高，且P-R格式較古顯格式更接近解析解。由例2，可由結(jié)果圖4-2看出u的大致圖形，由圖4-3可看出u隨時間的推移，溫度漸漸靠近邊界的溫度。可見，數(shù)值方法也能一定程度上反應解的狀況。

20

華北電力大學本科畢業(yè)設計（論文）

5穩(wěn)定性探究與分析

5.1穩(wěn)定性問題的提出

考慮數(shù)值算例例1，例2，用古典顯示格式和P-R格式解時，不同的r可能導致影響格式的穩(wěn)定性。以例2分析，采用古典顯式格式，以求解得t?0.1，h?1/20，r?h/k2變化時，T關于x，y的函數(shù)如圖5-1所示。

圖5-1：古顯格式解例2，u關于x，y的圖像（t?0.1，h??/20）

下面我們先研究確鑿解和微分方程的解之差，時隨著時間層數(shù)n的增大而無限增大還是有所控制。假使這種區(qū)別是無限增加的，則差分格式不穩(wěn)定，不穩(wěn)定的格式是不可用的。

5.2幾種分析穩(wěn)定性的方法

隨著人們對差分格式的深入研究，發(fā)展了以下幾種穩(wěn)定性研究方法。

00

即把Um改成了Um+?，而在?-圖研究法是在假定在固定的某節(jié)點上引入一個誤差?，

00

這一層上其他節(jié)點的值保持不變，且假定計算時沒有引入誤差，我們探究此時誤差是否會無限增大的方法。穩(wěn)定性分析的矩陣方法，采用矩陣的2范數(shù)，來衡量穩(wěn)定性。

Fourior級數(shù)法引進Fourior級數(shù)，通過考察增長因子來考察穩(wěn)定性。Fourior級數(shù)法又稱為Von-Neumann方法。在判斷有限差分近似的穩(wěn)定性方法中，以Fourior方法使用較為

21

華北電力大學本科畢業(yè)設計（論文）

廣泛，它僅適用于線性常系數(shù)的有限差分近似。以一維熱傳導方程的顯式格式為例，過程如下：

第一步：首先，要研究的差分方程可寫為：

m?N1?aumn?1j?m??bmunj?m

m?N0其中N0=｛-1，0，1｝，N1=｛0｝

一維熱傳導方程的古典顯式格式則表示為，

?1nnnun?ru?(1?2r)u?rujj?1jj?1

其次步：其次，對uij進行變量分開：

nunj??Vpexp[ip2?pxj]ln第三步：將unj替代為Vexp[i?xj]代入所考察的有限差分方程。

Vn?1exp[i?xj]?rVnexp[i?xj?1]?(1?2r)Vnexp[i?xj]?rVnexp[i?xj?1]

Vn?1?rVnexp[i?h]?(1?2r)Vn?rVnexp[?i?h]

Vn?12?h?1?4rsin?r(cos?h?isin?h)?(1?2r)?r(cos?h?isin?h)?1?2r(1?cos?h)

2Vn當

Vn?1?1Vn即對所有的?1?1?4ssin2?h2?1，時，格式穩(wěn)定。由于0?sin2?h2?1，故r?1時一維熱傳2導方程的古典顯式格式穩(wěn)定。

運用Fourior級數(shù)法可推導出，一維熱傳導方程的古典隱式格式、Crank-Nicolson隱式格式無條件穩(wěn)定。即對任意的r,該格式都穩(wěn)定。加權六點隱式格式的穩(wěn)定性條件與?有關，如下：

假使r?假使??假使??1，則穩(wěn)定性條件為??0；211，則穩(wěn)定性條件為r?；22(1?2?)1，則無條件穩(wěn)定。21，古典隱式格式、Crank-Nicolson隱式、P-R格式、D-R格式、422

二維熱傳導方程的幾種差分格式的穩(wěn)定性也可采用Fourior級數(shù)法求得，得出古典顯式格式的穩(wěn)定性條件為r?華北電力大學本科畢業(yè)設計（論文）

M-F格式無條件穩(wěn)定。

5.3r變化對穩(wěn)定性的探究

5.3.1古典顯式格式的穩(wěn)定性

根據(jù)5.1中對古典顯式格式穩(wěn)定性的分析可知，古典顯式格式的穩(wěn)定性條件是

r?h/k2?1/4。采用古典顯式格式解例2，以求解得t?0.1，h?1/20，r?h/k2變化時，T關于x，y的函數(shù)如上圖5-1所示。采用古典顯式格式解例1，以求解得t?0.1，h?1/20，

r?h/k2變化時，數(shù)值解與解析解間誤差如下表5-1所示。

表5-1古顯格式解例1的誤差（t?0.1，h?1/20）

r=1/8r=2/8r=3/8000-0.00013-0.00026-2.4E+26-0.00048-0.00095-4.3E+26-0.0009-0.0018-5.5E+26-0.00125-0.00248-6.1E+26-0.00138-0.00273-6.3E+26-0.00125-0.002481.11E+27-0.0009-0.0018-1E+27-0.00048-0.00095-7.8E+26-0.00013-0.00026-4.3E+26000

r=5/809E+286.1E+299E+294.4E+291.3E+294.4E+291.7E+297.4E+299.3E+290

r=6/803.42E+324.72E+323.86E+322.39E+321.74E+322.39E+323.25E+322.58E+321.53E+320

r=7/804.17E+298.93E+299E+294.48E+294.84E+294.84E+298.16E+292.06E+295.41E+290

誤差x=0,y=0

x=0.1,y=0.1x=0.2,y=0.2x=0.3,y=0.3x=0.4,y=0.4x=0.5,y=0.5x=0.6,y=0.6x=0.7,y=0.7x=0.8,y=0.8x=0.9,y=0.9x=1.0,y=1.0r=4/80-5.9E+29-6.1E+29-5E+29-5.8E+29-7.7E+29-5E+29-5.4E+29-1E+29-7.8E+290r=104.48E+293.31E+285.95E+289.57E+298.93E+297.3E+291.54E+293.31E+288.16E+29

0

誤差x=0,y=0x=0.1,y=0.1x=0.2,y=0.2x=0.3,y=0.3x=0.4,y=0.4x=0.5,y=0.5x=0.6,y=0.6x=0.7,y=0.7x=0.8,y=0.8x=0.9