8.1有關物理問題與數(shù)學物理方程8.2有限差分原理8.3矩形域中泊松方程的有限差分法8.4差分方程的迭代解法

8.5非矩形邊界區(qū)域泊松方程的有限差分法8.6一維擴散方程的有限差分法

8.7二維擴散方程的有限差分法8.8一維波動方程的有限差分法習題八下篇計算物理學

8.1有關物理問題與數(shù)學物理方程8.1.1方程的導出

一個連續(xù)體，如氣體、液體或固體，以及一個場，如電磁場、溫度場等,其狀態(tài)可用一時空函數(shù)u(x,y,z,t)來描寫。在這個連續(xù)體或場中發(fā)生的物理過程遵循特定的自然規(guī)律，其數(shù)學描述就是u隨時空變化的方程，它們常以偏微分方程的形式出現(xiàn)。

1．描述穩(wěn)定過程的泊松方程(橢圓型方程)

例如在靜電場中，其狀態(tài)可用標勢函數(shù)u(x,y,z)來描寫，根據(jù)高斯定理有

·(εE)=ρ

(8.1)

其中,ε為介質(zhì)的介電常數(shù)，ρ為自由電荷密度，

。而電場強度E又可以表示為E=-

u，將其代入上式可得-·(ε

u)=ρ。對于均勻介質(zhì)，ε為常數(shù)，于是有

(8.2)

該方程為靜電場中的泊松方程，其中表示

又稱為拉普拉斯方程。而在靜磁場中，同靜電場類似，有

A=-μj

(8.3)

該方程為靜磁場中的泊松方程，其中A是磁矢勢，μ是介質(zhì)的磁導率，j為傳導電流密度。

2．描述輸運過程的擴散方程(拋物型方程)

以物質(zhì)輸運為例,設ρ(x,y,z,t)為某一流體的密度函數(shù)，由于它的不均勻而發(fā)生擴散。按照擴散定律，物質(zhì)的擴散流密度j正比于密度的梯度，兩者反向，即有j=-D

ρ，其中D為擴散系數(shù)。由質(zhì)量守恒可導出連續(xù)性方程

(8.4)若D是均勻的，則有

(8.5)

該方程為擴散方程。同樣，在連續(xù)體中，由于溫度T的不均勻而發(fā)生熱傳導，對于這種方程可用類似方程描寫為

(8.6)

其中K為溫度傳導率(假設均勻)，式(8.6)描述的方程即為熱傳導方程。

3．描述振動傳播過程的波動方程(雙曲型方程)

對于交變電磁場，場強E、B

同場勢A、u的關系為

(8.7)

將上式代入Maxwell方程，可以得到場勢滿足的波動方程

(8.8)

其中c=1/為波速。對于在連續(xù)介質(zhì)中機械振動的傳播過程，有如下類似的波動方程

(8.9)

其中,u是介質(zhì)的位移，r是楊氏模量，a是波速，f是外力。8.1.2方程的分類

以上我們按照數(shù)理方程所描述的物理過程將它們分成三類，而在數(shù)學理論中，則按照這些方程的結(jié)構進行分類。二階線性偏微分方程的一般形式可以寫成(設自變量只有兩個，即x和y)

(8.10)

1．橢圓型方程(B2-4AC<0)

如二維泊松方程屬于這一類型。事實上，在這種情況下,方程(8.2)和(8.3)可寫成

(8.11)

對比式(8.10)可知，B=0，A=C=1。

2.拋物型方程(B2-4AC=0)

如一維擴散方程或熱傳導方程屬于這一類型，方程(8.5)和(8.6)可以寫成

(8.12)

對比式(8.10)可知，B=C=0,A=1。

3．雙曲型方程(B2-4AC>0)

如一維波動方程屬于這一類型，方程(8.8)和(8.9)可以寫成

(8.13)

對比式(8.10)可知，B=0，A=1，C=-1。8.1.3邊界條件和初始條件

從物理上講，以上給出的數(shù)理方程僅適用于描寫在一連續(xù)體或場的內(nèi)部發(fā)生的物理過程。但僅靠這些方程還不足以完全確定物理過程的具體特征，因為物理過程的具體特征還與連續(xù)體或場的初始狀態(tài)和邊界受到的外界影響有關。

從數(shù)學上講，一個偏微分方程會有無限多個解，偏微分方程加上邊界條件和初始條件，才能構成一個定解問題。

1．第一類邊界條件

若u代表方程中的未知函數(shù)，用Γ表示方程適用區(qū)域D的邊界。第一類邊界條件為

u|Γ=u0(rb,t)

(8.14)

其中,u0(rb,t)是定義在Γ上的已知函數(shù)，rb是相應邊界點的位矢。在這種邊界條件下邊界上連續(xù)體或者場的狀態(tài)是已知的。

2．第二類邊界條件

表達式為

(8.15)

其中,n表示Γ的外法線，q0(rb,t)是定義在Γ上的已知函數(shù)。若u是電磁場的勢，則代表場強；若u是密度(或溫度、或位移)，則代表流量(或熱流，或應力)。當q0=0時稱為第二類齊次邊界條件。

3.第三類邊界條件

表達式為

(8.16)其中a0、b0

和c0是定義在Γ上的已知函數(shù)。這一條件的物理情況較為復雜。對于熱傳導問題，它對應著連續(xù)體通過表面與外界發(fā)生輻射或?qū)α鞯确绞浇粨Q熱量，表面熱流正比于表面溫度u與外界溫度u0之差，即=k(u-u0)。

對于給定的物理問題，其邊界條件可能是很復雜的。例如在區(qū)域D的一部分邊界Γ1上有第一類條件，一部分邊界Γ2上有第二類條件，而其余部分則滿足第三類條件，或者是上述三式不能表示的另外的條件，需要具體分析邊界的物理情況，然后予以確切的數(shù)學描述。

4．初始條件

與時間坐標t相聯(lián)系，給出初始瞬間待求函數(shù)u在場域各處的值，即

u|t=0=f1(r)

(8.17)

以及初始瞬間場域各處u對時間的變化率，即

(8.18)通常僅含初始條件的定解問題稱為初值問題(柯西問題)；沒有初始條件而只有邊界條件的定解問題稱為邊值問題；既有初始條件又有邊界條件的定解問題，則稱為混合問題(也稱初邊值問題)。對于由拉普拉斯方程構成的第一、第二和第三類邊值問題，常分別稱為狄利克雷、諾伊曼和洛平問題。數(shù)學物理方程有許多是線性方程，對于這一類特定的物理問題往往有特定的解析方法，如分離變量法、積分變換法、格林函數(shù)法等等。解有時能用各種初等函數(shù)和超越的特殊函數(shù)來表達,但這些只限于比較典型的情況。更多的實際物理問題只能用非線性方程或方程組來描述，其求解方法更為復雜，只有少數(shù)問題有解析解,這時需要借助計算機來求偏微分方程的數(shù)值解。變分解法、有限元法、有限差分法、邊界元等方法是數(shù)值解法中應用廣、精度較高的幾種常見解法，其共同特點是將求解偏微分方程的問題轉(zhuǎn)化為求解線性代數(shù)方程組的問題。8.2.1差商公式

將微分方程中的微商用差商代替(離散化)，是差分法求解偏微分方程的基礎。設u是坐標x的函數(shù)，取Δx=h等分坐標軸，結(jié)點坐標為xi，相應的有ui(i=1,2,…)，則有泰勒展開

(8.19)8.2有限差分原理

(8.20)

(8.21)

(8.22)

1.誤差為O(h)的差商公式

由式(8.19)可得一階向前差商公式

(8.23)

由式(8.20)可得一階向后差商公式

(8.24)

式(8.21)-2×式(8.19),可得二階向前差商公式

(8.25)

式(8.22)-2×式(8.20),可得二階向后差商公式

(8.26)

2．誤差為O(h2)的差商公式

將式(8.25)代入式(8.19)，可得一階向前差商公式

(8.27)

將式(8.26)代入式(8.20),可得一階向后差商公式

(8.28)［式(8.19)-式(8.20)］/(2h),可得一階中心差商公式

(8.29)

［式(8.19)＋式(8.20)］/h2，可得二階中心差商公式

(8.30)8.2.2差分格式的收斂性和穩(wěn)定性

用某種近似差商代替相應的微商，叫做一種差分格式。所謂差分格式的收斂性，是指當步長h→0時，差分方程的解是否收斂于微分方程的解。顯然只有那些具有收斂性的差分格式才是所需要的。

所謂差分格式的穩(wěn)定性，是指誤差Δu在運算過程中不會失控，即累計誤差是否會無限增加。任何初值擾動對差分數(shù)值解的影響隨時間推移不再增加(強穩(wěn)定)或在一段時間內(nèi)有界(弱穩(wěn)定)。

有關差分法的數(shù)學理論對上述收斂性和穩(wěn)定性都作了詳細討論。以下僅作結(jié)論性說明，不予以數(shù)學證明。8.3.1五點差分格式

在二維情況下，泊松方程形式為

(8.31)

若f(x,y)=0，則式(8.31)為拉普拉斯方程。8.3矩形域中泊松方程的有限差分法取Δx=Δy=h的正方形網(wǎng)格覆蓋x-y平面(如圖8.1所示)，結(jié)點坐標為(xi,yj)(i=1,2,…,N;j=1,2,…,M)。結(jié)點處的函數(shù)為u(xi,yj)=uij。在(i,j)點，利用中心差商公式(8.30)，則式(8.31)變?yōu)?/p>

(8.32)式(8.32)中涉及到五個點，即中央點(i,j)及周圍四點，如圖8.1所示。這種差分格式的誤差為O(h2)。差分法的數(shù)學理論已經(jīng)證明，由上述五點格式得到的差分方程，在給定的邊界條件下具有唯一的解，且當h→0時趨于微分方程的解。圖8.18.3.2矩形域的拉普拉斯方程

(8.33)

【例8.1】

用有限差分法求解拉普拉斯方程，邊界條件如圖8.2所示。圖8.2若取h=5，如圖8.2所示有三個內(nèi)點，相應的u值記為u1、u2、u3。根據(jù)式(8.32)，可列出關于三個內(nèi)點的差分方程組

它的矩陣形式為

解得。若取h=2.5，精度會提高，當然計算量會增大。有限差分法計算結(jié)果與本題精確解的比較如表8.1所示。表8.1有限差分法計算結(jié)果與精確解的比較根據(jù)(8.32)式可以得到泊松方程的差分方程組

(8.34)8.4差分方程的迭代解法

該方程組中方程個數(shù)等于網(wǎng)格結(jié)點數(shù)，并且每個方程的左邊最多有五項。故該方程組的系數(shù)矩陣具有大型、稀疏和帶狀的特點。為節(jié)省內(nèi)存和機時，一般不采用消元法，而是采用迭代法，包括如下三種方法。

(1)同步法。

將式(8.34)化為

(8.35)上式中右邊的u用第k步的值，則左邊得到第k+1步的值。這種迭代方法稱為同步法，它需要兩套內(nèi)存，一套存第k步的值，一套存第k+1步的值，收斂較慢。

(2)異步法。

由于u的計算是按照i、j

由小到大進行的，迭代到某一步在計算新的uij時，新的ui，j-1和ui-1，j已被算出，所以式(8.35)可以改寫為

(8.36)

這種迭代方法稱為異步法，它只需一套內(nèi)存，收斂較快。

(3)逐次超松弛迭代法(SOR法)。

在異步法的基礎上，加權平均，得到SOR法。引入超松弛因子ω，將式(8.36)改寫為

(8.37)當ω=1時，式(8.37)就是式(8.36)。當在ω∈(1,2)內(nèi)取一適當值時，可獲得較快的收斂速度。如在例8.1中，當h=2.5并要求|u(k+1)-u(k)|<10-3時，取不同的ω值，所需的迭代次數(shù)k也不同，如表8.2所示。表8.2從表8.2可以看出，ω值的選取對迭代次數(shù)k影響較大。顯然，并非ω越大或越小時，迭代次數(shù)k越小。那么如何選取最佳的ω值使得迭代次數(shù)k越小呢？通常采用如下方法:

(1)對于矩形區(qū)域，常數(shù)邊界條件下超松弛因子ω常取為

(經(jīng)驗公式)其中，n+1、m+1分別為x軸和y軸上的等分結(jié)點數(shù)。而對于正方形區(qū)域，ω取為

其中n+1為每邊上的等分結(jié)點數(shù)。

(2)在(1,2)內(nèi)，依次取幾個ω值，用少數(shù)幾步迭代，找出最佳值。

SOR法迭代過程中除了超松弛因子ω的選取外，關于初值的選取也很重要。事實上，關于上述各種迭代法的收斂性問題，可以證明，當h→0或k→∞時，差分方程的解將趨近于微分方程的解。收斂性與初值無關，但收斂速度與初值有關。關于初值的選取問題，若f(x,y)=0，可取u|Γ平均值為初始值，或任取初值。若f(x,y)≠0，可取u=0為初值。另外，關于收斂的判別方法，可事先指定出誤差的范圍ε，當時，迭代終止。

采用有限差分方法求解偏微分方程時，還需注意對邊界條件的處理。對于規(guī)則矩形邊界而言，研究區(qū)域的網(wǎng)格結(jié)點落在邊界上，第一類邊界條件顯然無需再做處理，而對于第二、三類邊界條件，可采用以下差分格式。根據(jù)以下第三類邊界條件

(8.38)為找到誤差為O(h2)的邊界條件的差分式，在圖8.3所示的xd和yd處用式(8.27)中的一階向前差商代替該處的，而在xu和yu處用式(8.28)中的一階向后差商代替該處的，解出邊界值，其結(jié)果為

(8.39)圖8.3一般情況下，a和b是x、y的函數(shù)，上式應作相應變化。顯然，當a=0時，上式即為第二類邊界條件的差分格式。

通常邊界是不規(guī)則的，邊界條件處理方法如下：對于第一類邊界條件u|Γ=u0，Γ為邊界。若結(jié)點在邊界上，直接代入。若結(jié)點不在邊界上，采用不對稱網(wǎng)格方法，如圖8.4所示，邊界與網(wǎng)格線交點為D和B，對于泊松方程,則有

(8.40)

(8.41)

其中，。此時泊松方程在邊界上的差分格式為

(8.42)

對于第二類邊界條件，若結(jié)點在邊界上(如圖8.5所示)，則有

(8.43)其中，。其差分格式為

(8.44)圖8.4圖8.5若結(jié)點不在邊界上，過P點向邊界作垂線，交邊界于P′。用P′的法向單位矢n做為P的法向矢，實質(zhì)上是把P近似為P′，差分格式處理同式(8.44)。

對于第三類邊界條件

的處理方法是前兩類邊界條件處理方法的組合。

【例8.2】

試采用有限差分法編程求解拉普拉斯方程在網(wǎng)格結(jié)點處的值。

圖8.6如圖8.6所示，取h=1，ω=1.25，迭代終止誤差為ε=10-4。

計算程序:

programmain

integeri,j

realx(5),y(4),u(5,4),t,p,q

real,parameter::w=1.25,eps=1e-4,h=1,pi=3.1415926!w，

eps分別表示ω和ε

open(1,file=′example.dat′)

doi=1,5!輸入邊界條件

x(i)=(i-1)*h

u(i,1)=sin(pi*x(i)/4)

u(i,4)=0.0

enddo

doj=1,4

y(j)=(j-1)*h

u(1,j)=y(j)*(y(j)-3)

u(5,j)=0.0

enddo

doi=2,4!SOR法賦初值

doj=2,3

u(i,j)=0.0

enddo

enddo

10p=0.0

doi=2,4

doj=2,3

t=u(i,j)

u(i,j)=w*(u(i,j-1)+u(i-1,j)+u(i,j+1)+u(i+1,j))/4+(1-w)*u(i,j)

q=abs(u(i,j)-t)

if(q.gt.p)p=q !p為每次迭代的最大誤差

enddo

enddo

if(p.gt.eps)goto10doi=1,5 !輸出結(jié)果

doj=1,4

write(*,*)i,j,u(i,j)

write(1,*)i,j,u(i,j)

enddo

enddo

end計算結(jié)果:

iju(i,j)

110.000000e+00

12-2.000000

13-2.000000

140.000000e+00

217.071068e-01

22-4.403373e-01

23-6.375487e-01

240.000000e+00

311.000000

321.690345e-01

33-1.098504e-01

340.000000e+00

417.071069e-01

422.263135e-01

432.911735e-02

440.000000e+00

510.000000e+00

520.000000e+00

530.000000e+00

540.000000e+008.5.1圓形域中泊松方程的有限差分解

1．泊松方程的差分格式

這里仍以二維泊松方程為例，該方程表示為

2u=f(r,φ)

(0<r<r0,0≤φ≤2π)

(8.45)8.5非矩形邊界區(qū)域泊松方程的有限差分法

該方程在極坐標系下可以寫為

(8.46)圖8.7用一組等角線和等距離同心圓分割區(qū)域，將半徑N-1等分，圓周M-1等分(如圖8.7所示)，即有

(8.47)因此有結(jié)點坐標

(8.48)用誤差為O(Δr2)和O(Δφ2)的中心差商代替(8.46)式中的微商，可以得到泊松方程的差分格式為

(8.49)

其中，α0=［Δφ(i-1)］-2，α1=1-(2i-2)-1,α2=1+(2i-2)-1

。

2．周期性條件

u(r,φ)=u(r,φ+2π)

(8.50)

即滿足ui1=uiM，ui2=ui,

M+1。

3.圓心處差分格式

方程(8.49)不適合于圓心(i=1)。如圖8.8所示，利用直角坐標系下的五點差分格式，有

(8.51)圖8.8在一般情況下，可將上式寫為

(8.52)

其中是半徑為h的圓周上各結(jié)點的平均值，即

(8.53)

4.邊界條件的差分格式

(8.54)

用誤差為O(h2)的一階向后差商代替，則有

(8.55)

【例8.3】

同軸線內(nèi)導線半徑為a，外導線半徑為b,兩導線間填充均勻介質(zhì)，相對介電常數(shù)為εr，兩導線間電位差為V0，求同軸線介質(zhì)間的電勢分布。圖8.9對于內(nèi)外導線間的電場分析，可理想化為二維場問題，其間電勢函數(shù)u(r,φ)應滿足拉普拉斯方程，極坐標下拉普拉斯方程形式為

根據(jù)有限差分原理，同樣作五點中心差商，其離散結(jié)點如圖8.9所示。對上式進行化簡，可得

引入超松弛因子，進行超松馳迭代

這里松弛因子ω的計算式，采用經(jīng)驗公式

，其中p、q分別為半徑r方向上和角φ方向上的結(jié)點總數(shù)。

在利用上述迭代公式求得同軸線介質(zhì)間的電勢u分布后，可以利用

求出電場。而根據(jù)電勢u分布沿角度方向的對稱性，電勢u僅沿法向有變化，

這里利用誤差為O(Δr2)的一階向前差商公式

就可以求出內(nèi)表面的電場。利用

和

可以求出同軸線的特征阻抗，而同軸線的特征阻抗精確解為

本例計算中,取a=4.0cm，b=6.909cm，Δr=0.1cm，Δφ=，V0=18V。圖8.10給出了采用有限差分方法的計算結(jié)果，即同軸線沿半徑方向上的電勢分布曲線，其中電勢分布的解析解為u(r,φ)=。我們還分別計算了相對介電參數(shù)εr的值為1、2、3、4、5時的特征阻抗，表8.3中，Z0、Z1分別表示特征阻抗的精確解和數(shù)值解。從圖8.10和表8.3中可以看出數(shù)值解和精確解誤差較小，從而說明了有限差分方法的可行性。圖8.10同軸線電勢沿徑向分布結(jié)果表8.3特征阻抗數(shù)值解與精確解對比及誤差8.5.2軸對稱場區(qū)域泊松方程的有限差分解

采用柱坐標系(r,,z)。由于是軸對稱場，u(r,z)與無關。在柱坐標系中泊松方程為

(8.56)圖8.11取Δx=Δz=h，用正方形網(wǎng)格覆蓋r-z平面(如圖8.11所示)，有

(8.57)用誤差為O(h2)的中心差商代替微商，且u(ri,zj)=uij，則有

(8.58)其中，α1=［1-1/(2i-2)］，α2=［1+1/(2i-2)］(i=2,3,…,N-1;j=2,3,…,M-1)。上式不適用于軸上各點(i=1)。應用羅畢達法則有

(8.59)將上式代入式(8.56)，所以在軸上有

(8.60)軸上的差分公式可以寫為

(8.61)

對于其他形狀的邊界，可用網(wǎng)格邊界近似代替實際邊界，或者用插值近似方法，但由于程序編制過于復雜，通常不采用差分法，而是利用有限元方法或邊界元方法。8.6.1隱式六點差分格式(C－N格式)

以下介紹一維擴散方程或熱傳導方程的有限差分解法?？紤]一維擴散方程的定解問題

(8.62)8.6一維擴散方程的有限差分法圖8.12

取Δx=h，Δt=τ進行離散化，如圖8.12所示，結(jié)點坐標為

(8.63)結(jié)點處的函數(shù)為u(xi,tk)=。在點，用中心差商，用(i,k)和(i,k+1)兩點的中心差商的平均值代替，則式(8.62)中的偏微分方程變?yōu)?/p>

(8.64)

引入P=τD/h2，，將上式中的含u

(k+1)項移至等號左邊，將含u(k)項移至等號右邊，式(8.64)變?yōu)?/p>

(8.65)上式表明由k時的u值可求得k+1時的u值，但要解聯(lián)立方程組，所以這種差分格式是隱式的。整個方程涉及到六個結(jié)點處的u值，所以稱為隱式六點差分格式，又稱為CrankNicolson格式，簡稱C－N格式，誤差為O(τ2)+O(h2)，是無條件穩(wěn)定的。8.6.2邊界條件的差分格式

由式(8.62)知，一維擴散方程的邊界條件為

(8.66)

圖8.13在x軸上設置兩個虛格點i=0和i=N+1(見圖8.13)。

用中心差商代替式(8.66)中的，則得

(8.67)由式(8.67)解出u0和un+1，分別代入i=1和i=N的式(8.65)，得到

(8.68(a))

(8.68(b))8.6.3差分方程組及其求解

把式(8.65)和式(8.68(a))、(8.68(b))結(jié)合起來，構成差分方程組，其形式為

AU=R(8.69)其中,U=(u1,u2,…,un)是未知量組成的矢量。系數(shù)矩陣A是三對角的，而R是由前一時刻的u值組成的矢量R=(R1,R2,…,Rn)。該方程組可利用5.3節(jié)中的追趕法進行求解。由式(8.65)和式(8.68(a))、(8.68(b))可知

(8.70)

(8.71)對于三對角矩陣A，用ai2表示對角元，用ai1表示左旁元，用ai3表示右旁元，則差分方程組的矩陣形式為

(8.72)利用高斯消元法，消去左旁元可得

(8.7３)其中

(8.74)由(8.73)式的最后一個方程可以解出un=Rn′/an2′，再用回代法可以解出

(8.75)8.6.4計算程序

1．主程序流圖

dimensionx(N),u(N),R(N),A(N,3)

(1)輸入D,a0,a1,a2,b1,b2,c1,c2,tmax,h,τ。

輸入結(jié)點坐標xi=(i-1)h(i=1,2,…,N)

輸入初值ui(x,0)=u0(xi)

(2)調(diào)用計算原始和消元后的系數(shù)矩陣A。

callcoef(N,a1,b1,a2,b2,h,P1,A)

(3)

do50k=1,kmax

t=k*tao!tao對應于τ

R(1)=(b1*P2-h*a1)*u(1)+b1*u(2)+2*h*c1

R(N)=(b2*P2-h*a2)*u(N)+b2*u(N-1)+2*h*c2

do30i=2,N-1

30

R(i)=u(i-1)+2*P2*u(i)+u(i+1)

(4)消元法解方程AU=R。

callsolve(N,A,R)

do40i=1,N

x(i)=(i-1)*h

u(i)=R(i)

40

write(*,*)t,x(i),u(i)

50

continue

end

2．子程序

(1)系數(shù)矩陣(原始的和消元后的)子程序。

subroutinecoef(N,a1,b1,a2,b2,h,P1,A)

dimensionA(N,3)

A(1,2)=P1*b1+h*a1

A(1,3)=-b1

A(N,2)=P1*b2+h*a2

A(N,1)=-b2

do10i=2,N-1

A(i,1)=-1,A(i,2)=2*P1

10

A(i,3)=-1

do20i=2,N

20

A(i,2)=A(i,2)-A(i,1)*A(i-1,3)/A(i-1,2)

return

end

(2)消元法解三對角系數(shù)矩陣方程的子程序。

subroutinesolve(N,A,R)

dimensionA(N,3),R(N)

do10i=2,N

10

R(i)=R(i)-A(i,1)*R(i-1)/A(i-1,2)

R(N)=R(N)/A(N,2)

do20i1=1,N-1

I=N-i1

20R(i)=(R(i)-A(i,3)*R(i+1))/A(i,2)

return

end8.7.1交替方向隱式差分格式(ADI格式)

考慮以下二維擴散方程的定解問題

(8.76)

邊界條件見后。8.7二維擴散方程的有限差分法求解(8.76)方程的要點是把二維問題化為一維問題，并使差分方程組的系數(shù)矩陣變?yōu)槿龑切问健?/p>

取Δx=Δy=h的正方形網(wǎng)格覆蓋x-y平面，并取Δt=τ。結(jié)點坐標為

(8.77)結(jié)點處的函數(shù)為。

在(i,j,k+1/2)點，用中心差商，用k+1時的中心差商，而用k時的中心差商代替，則式(8.76)中的偏微分方程變?yōu)?/p>

(8.78(a))這種時間上的不一致引起的偏差，可用下一個時刻的偏差來補償。在(i,j,k+3/2)點，用中心差商，用k+1時的中心差商，而用k+2時的中心差商代替，則式(8.76)變?yōu)?/p>

(8.78(b))與8.6節(jié)類似，引入P=τD/h2，P1=

。經(jīng)整理，式(8.78(a))和式(8.78(b))分別變?yōu)?/p>

(8.79(b))(8.79(a))

比較可知，式(8.79(a))和(8.79(b))的形式與式(8.65)相同，所以求解方法也相同。應注意，k為奇數(shù)時用式(8.79(b))沿y方向算，k為偶數(shù)時用式(8.79(a))沿x方向計算，并且只輸出后一結(jié)果，以便減少因為時間上的不一致引起的偏差。這種格式稱為交替方向隱式格式，簡稱ADI格式。顯然，其誤差為O(τ2)+O(h2)，并且可以證明是無條件穩(wěn)定的。8.7.2邊界條件的差分格式

式(8.76)中的偏微分方程滿足如下邊界條件

(8.80(a))

(8.80(b))

如圖8.14所示。圖8.14邊界條件可以用兩種方法給出誤差為O(h2)的邊界條件差分式。一種方法設置虛結(jié)點，用中心差商代替，可得虛結(jié)點處的u。沿y方向在邊界①、②上有

(8.81(a))

(8.81(b))沿x方向在邊界③、④上有

(8.81(c))

(8.81(d))

另一種方法是分別用向前和向后代替①、②和③、④邊界處的，可得邊界點處u(見式(8.39))

(8.82(a))

(8.82(b))(8.82(c))(8.82(d))

將式(8.81(c))、(8.81(d))與(8.79(b))結(jié)合起來，消去ui0和ui,m+1，可得到沿y方向計算的方程組。對于給定的i，其矩陣式為

Auy=R

其中，uy(j)=uij(j=1,2,…,m),R(j)=ui-1,j+2P2uij+ui+1,j。而R(1)←R(1)b3＋2hc3(xi,t),R(M)←R(M)b4+2hc4(xi,t)。系數(shù)矩陣A具有三角形式：A(1,2)=b3P1+ha3，A(1,3)=-b3，A(M,2)=b4P1+ha4,A(M,1)=-b4,A(j,1)=-1,A(j,2)=P1，A(j,3)=-1(j=2,…,M-1)。

在沿y方向計算時，是對每個給定的i(i=2,…,N-1)逐列進行的。至于i=1和N這兩列邊界的u值，利用式(8.82(a))、(8.82(b))進行計算。對于沿x方向的計算，有類似的情況。8.7.3計算程序流圖

計算程序流程為

dimensionx(N),y(M),u(N,M),R(N)(設N>M),A(N,3),B(M,3)

輸入D,a0,b0,tmax,h,τ,(a,b)1,2,3,4，P=2τD/h2,P1=1/P+1,P2=1/P-1，N=1+a0/h，M=1+b0/h，kmax=tmax/τ。原始的和消元后的系數(shù)矩陣A和B(若與t有關，應置于k的循環(huán)之外)

callcoef(N,a1,b1,a2,b2,h,P1,A)

callcoef(M,a3,b3,a4,b4,h,P2,B)

do8k=1,kmax,2

沿y方向算。t=k*τ

do②i=2,N-1

do①j=1,M①Rj=ui-1,j+2*P2*uij+ui+1,j

R1=R1*b3+2*h*c3(xi,t)

Rm=Rm*b4+2*h*c4(xi,t)

callsolve(M,B,R)

do②j=1,M

②Vij=Rj

邊界值①和②上的u值為

do③j=1,M

V1j=［2*h*c1(yj,t)+4*b1*V2j-b1*V3j］/(2*h*a1+3*b1)③Vnj=［2*h*c2(yj,t)+4*b2*Vn-1,j-b2*Vn-2,j］/(2*h*a2+3*b2)

沿x方向算。t=t+τ

do⑤j=2,M-1

do④i=1,N

④Ri=Vij-1+2*P2*Vij+Vij+1

R1=R1*b1+2*h*c1(yj,t)

Rn=Rn*b2+2*h*c2(yj,t)

callsolve(N,A,R)

do⑤i=1,N⑤uij=Ri

邊界③和④上的u值為

do⑥i=1,N

ui1=［2*h*c3(xi,t)+4*b3*ui2-b3*ui3］/(2ha3+3b3)

⑥uim=［2*h*c4(xi,t)+4*b4*ui,m-1-b4*ui,m-2］/(2*h*a4+3*b4)8輸出t,u

需要說明的是，子程序coef和solve與8.6節(jié)的相同。8.7.4二維顯式格式

對于式(8.76)，在(i,j,k)點，用向前差商，

用中心差商代替，則可以得到

(8.83)引入P=τD/h2，可把式(8.83)化為

(8.84)

利用式(8.84)可由k時的u直接求k+1時的u，不必解聯(lián)立方程，故稱顯式格式。誤差為O(τ)+O(h2)。顯式格式雖然簡單，但它的穩(wěn)定性是有條件的，要求

(8.85)此外，在設計程序時需用兩套u的數(shù)組，一套存k時刻的，另一套存k+1時刻的。

對于邊界條件式(8.80)，將用誤差O(h)的向前(①、③邊界)或向后(②、④邊界)差商代替，就得到邊界條件的差分式。計算程序流程為

dimensionx(N),y(M),u(N,M),V(N)

輸入D,a0,b0,tmax,h,P,(a,b,c)1,2,3,4，τ=h2P/D,N=1+a0/h,M=1+b0/h,kmax=tmax/τ。

結(jié)點坐標與初值i=1,…,N,xi=(i-1)h；j=1,…,M,yj=(j-1)h，uij=u0(xi,yj)時間演化

do40k=1,kmax

t=k*τ

do10i=2,N-1

do10j=2,M-1

10Vij=(1-P)*uij+P*(ui,j-1+ui,j+1+ui-1,j+ui+1,j)

do20i=2,N-1

do20j=2,M-1

20uij=Vij

邊界值

do25j=1,M

u1j=(b1*u2j+h*c1)/(h*a1+b1)①邊界

25unj=(b2*un-1,j+h*c2)/(h*a2+b2)②邊界

do30i=1,N

ui1=(b3*ui2+h*c3)/(h*a3+b3)③邊界

30uim=(b4*ui,m-1+h*c4)/(h*a4+b4)④邊界

40輸出t,u8.8.1顯式差分格式

考慮如下一維波動方程的定解問題

(8.86)8.8一維波動方程的有限差分法取Δx=h，Δt=τ進行離散化，節(jié)點坐標為

(8.87)節(jié)點處的函數(shù)為u(xi,tk)=。在(i,k)點，用中心差商代替,則式(8.86)中的偏微分方程變?yōu)?/p>

(8.88)引入P=(τc/h)2，上式變?yōu)?/p>

(8.89)

由k-1和k時刻的u可直接求k+1時刻的u，不必解聯(lián)立方程組，故這種差分方程格式是顯式的。誤差為O(τ2)+O(h2)，可以證明，當P<1時，這種格式是收斂和穩(wěn)定的。8.8.2初值、邊界條件的差分格式

初值條件用前向差商代替，有