求解三維熱傳導(dǎo)方程的高效精確算子分裂方法

武莉莉

(1. 寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院, 寧夏 銀川 750021; 2. 寧夏師范學(xué)院 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院, 寧夏 固原 756000)

本文考慮如下三維熱傳導(dǎo)方程的初邊值問題：

為了討論方便,假設(shè)空間域Ω=[0,1]×[0,1]×[0,1].T表示溫度;ν為熱擴(kuò)散系數(shù);α(x,y,z)是光滑函數(shù),g0,g1,d0,d1,l0和l1是常數(shù).

熱傳導(dǎo)方程是一類非常重要的偏微分方程,能夠描述很多自然環(huán)境和工程設(shè)備中的物理現(xiàn)象,比如超導(dǎo)體中的熱傳導(dǎo),燃燒室內(nèi)的溫度分布等.只有一小部分問題能夠獲得精確解,絕大多數(shù)問題要得到精確解是非常困難的,甚至是不可能的,所以尋求其數(shù)值求解方法具有非常重要的理論價值和現(xiàn)實(shí)意義.

求解該方程的數(shù)值方法有很多,如有限差分法[1-3], 有限體積法[4], 有限元法[5-6]以及配點(diǎn)法和譜方法[7-8], 等等. 這些方法中有限差分方法由于格式構(gòu)造簡單容易理解引起了廣大學(xué)者的關(guān)注. 對于高維問題而言,一個重要的方面是如何減低計(jì)算成本和如何提高計(jì)算效率. 對于高維熱傳導(dǎo)方程而言,顯式差分格式由于受到苛刻的穩(wěn)定性條件的限制很難應(yīng)用,而全隱格式由于在每個時間層上都需要迭代求解,也存在計(jì)算量大、計(jì)算效率低下的顯著缺點(diǎn).為了克服這兩方面的缺點(diǎn),二十世紀(jì)五六十年代算子分裂方法應(yīng)用而生[9-10].目前有兩種算子分裂方法應(yīng)用最廣,其中一種是交替方向隱式(ADI )方法[8,11-16],另一種是局部一維化(LOD)方法[17-21]. 他們均可以將高維問題轉(zhuǎn)化為多個一維問題來求解,而一維問題往往可以采用快速的三對角矩陣算法 (TDMA). 在目前已有的關(guān)于算子分裂方法的研究中,很多空間高精度的算法已經(jīng)被提出來了, 但是時間方向的精度仍然是低精度的,這是由于傳統(tǒng)的Peaceman-Rachford方法[11-13], Crank-Nicolson方法[8,12,16-17,20], 或者其他低階精度的顯式或者隱式方法[19,21]被用于時間的離散和處理.

本文構(gòu)造時間和空間均具有高精度的數(shù)值算法,采用算子分裂的局部一維化方法將方程(1)表示為三個一維方程,再利用高精度有限差分格式和高精度Padé近似來分別離散方程(1)中的空間導(dǎo)數(shù)和時間導(dǎo)數(shù), 從而得到兩種高精度有限差分格式.從理論上嚴(yán)格證明格式的穩(wěn)定性,給出數(shù)值實(shí)驗(yàn)結(jié)果.

1 算子分裂方法

根據(jù)局部一維化方法的原理,方程(1)從形式上可以分裂為三個一維方程:

為了從第n時間層獲得第n+1時間層未知函數(shù)的值, 可以假設(shè)方程(6)可以實(shí)現(xiàn)由第n時間層獲得n+1/3時間層的值, 方程(7)可以實(shí)現(xiàn)由第n+1/3時間層獲得第n+2/3時間層的值, 方程(8)可以實(shí)現(xiàn)由第n+2/3時間層獲得n+1時間層的值. 下一步將針對上述方程構(gòu)造高精度格式.

2 高精度緊致差分格式

以τ表示時間步長,空間取等間距網(wǎng)格,步長用h表示。網(wǎng)格點(diǎn)為(xi,yj,zk,tn),xi=ih,yj=jh,zk=kh,tn=nτ,i,j,k=0,1,…,N,h=1/N,n≥0. 為了建立方程(6)的高精度緊致差分格式,利用如下緊致差分公式[22]:

(9)

(10)

式(10)具有空間四階精度.添加上初始和邊界條件可以得到如下的常微分方程的初值問題：

(11)

由此,很容易獲得方程(11)的解為

T(t)=-B-1F+etA-1B[T(0)+B-1F]

(12)

(13)

eP=(12-6P+P2)-1(12+6P+P2)+O(P4)

(14)

(15)

式中:τ為時間步長.將方程 (15) 代入方程(13)可得

(16)

(17)

把方程 (17) 記為FOD格式.應(yīng)用該格式能夠?qū)崿F(xiàn)由tn到tn+1的計(jì)算.由推導(dǎo)過程易知該格式時間和空間均具有四階精度.

下面構(gòu)造求解方程(1～4)的六階精度差分格式.利用文獻(xiàn)[1]中結(jié)論, 可以得到二階導(dǎo)數(shù)的六階精度離散格式：

(18)

(20)

為此,可以得到方程(6)在x方向的常微分方程系統(tǒng):

(21)

式中

T(t),F和B的定義如前面所述.方程(21)的解析解為

T(t)=-B-1F+etD-1B[T(0)+B-1F]

(22)

離散方程 (22) 可以得

(23)

利用六階精度的(3,3) Padé 格式[23]計(jì)算時間t的導(dǎo)數(shù),可得

eP=(120-60P+12P2-P3)-1×

(120+60P+12P2+P3)+O(P6)

(24)

(25)

將式(25)代入式(23), 整理后可得

(26)

利用類似的方程處理y和z方向, 最終可以得到求解方法(6～8)的新的有限差分格式, 表示如下：

(27)

將式(27)記為SOD格式.由推導(dǎo)過程可知,其截?cái)嗾`差為O(τ6+h6).

3 穩(wěn)定性分析

引理1假設(shè)λ是(N-1)×(N-1)階矩陣A-1B的特征值,x是其相應(yīng)的特征向量.則λ是實(shí)數(shù),并且滿足λ≤0.

證明假設(shè)λ是特征值,x是矩陣A-1B對應(yīng)的特征向量.他們滿足如下的關(guān)系:

A-1Bx=λx

或者

xTBx=λxTAx

由于

因此

xTBx<0

(28)

且有

因此

xTAx>0

上述結(jié)果即可表明λ是實(shí)數(shù),并且滿足：λ≤0.

引理2假設(shè)μ是矩陣D-1B的特征值,x是其相對應(yīng)的特征向量.則μ是實(shí)數(shù),并且滿足μ≤0.

證明對于矩陣D有

利用不等式

2xy≥ -x2-y2

和

2xy≤x2+y2

可以得到

結(jié)合式(28)的結(jié)果,可以得到μ≤0.

定理1FOD格式 (17) 和SOD格式(27)均是無條件穩(wěn)定的.

證明令λi(i=1,2,…,N-1)表示A-1B的特征值, 利用引理1,可以得到

由此可得

式中

令μi(i=1,2,…,N-1)為矩陣D-1B的特征值,利用引理2, 可得

因此

式中

ρ(H1)和ρ(H2)是矩陣H1和H2的譜半徑.至此定理得證.

推論1FOD格式(17)和SOD格式(27)均是無條件收斂的.

證明由FOD格式(17)和SOD格式(27)的建立過程很容易得出這兩種格式的截?cái)嗾`差分別為O(τ4+h4)和O(τ6+h6),即格式(17)和格式(27)與原微分方程(1)均是相容的.又由定理1可得格式(17)和格式(27)均是無條件穩(wěn)定的,因此根據(jù)Lax等價性定理,格式 (17) 和格式(27) 均是無條件收斂的.

4 數(shù)值算例

采用兩個算例來驗(yàn)證本文方法的性能,分別利用本文所提的FOD和SOD格式進(jìn)行計(jì)算,并且為了與低精度方法進(jìn)行比較,還采用了Douglas-Gunn格式[10]進(jìn)行了計(jì)算.全部的計(jì)算都是在PC機(jī)上,采用Matlab語言編程完成的.

格式的收斂階定義如下：

式中:L2-err表示L2范數(shù)誤差；h1和h2分別為兩個不同空間網(wǎng)格步長.L2-err定義如下：

算例1

0≤x,y,z≤1,t>0

該問題的分析解為

u(x,y,z,t)=e-3π2tsin(πx)sin(πy)sin(πz)

為了衡量時間精度,設(shè)定h=0.02,讓時間步長τ發(fā)生變化.表1分別給出了Douglas-Gunn格式[10]、本文的FOD格式 (17)和SOD格式 (27)的計(jì)算結(jié)果.可以看出Douglas-Gunn格式時間僅有二階精度,FOD格式具有四階精度,SOD格式具有六階精度.為了衡量空間精度,設(shè)定τ=0.001,讓空間步長h發(fā)生變化.表2分別給出了三種格式的計(jì)算結(jié)果.結(jié)果顯示本文FOD格式和SOD格式空間上也分別達(dá)到了四階和六階精度,而Douglas-Gunn格式精度僅有二階.為了驗(yàn)證格式的整體精度,設(shè)定h=τ=1/10,1/20,1/30,1/40, 表3給出了t=0.5時刻三種方法的數(shù)值計(jì)算結(jié)果.可以看出三種格式均達(dá)到了各自理論上的精度階.FOD格式和SOD格式整體具有四階和六階精度,而Douglas-Gunn格式[10]整體上僅具有二階精度.

表1 算例1當(dāng)h=0.02,t=0.5時,不同τ的L2范數(shù)誤差及收斂階Tab.1 L2-error and accuracy order with h=0.02,t=0.5 for Example 1

表2 算例1當(dāng)τ=0.001,t=0.5時,不同h的L2范數(shù)誤差及收斂階Tab.2 L2-error and accuracy order with τ=0.001,t=0.5 for Example 1

表3 算例1當(dāng)τ=h,t=0.5時的L2范數(shù)誤差及收斂階Tab.3 L2-norm error and accuracy order with τ=h,t=0.5 for Example 1

算例2:

0≤x,y,z≤1,t>0

該問題的分析解為

u(x,y,z,t)=e-4π2tsin(4πx)sin(4πy)sin(4πz)

針對算例2,為了評估本文格式的時間、空間精度以及整體精度.表4～表6給出了上述三種格式的計(jì)算結(jié)果.其中表4取定空間步長h=0.012 5,讓時間步長τ變化,考察了時間精度；表5取定時間步長τ=0.005,讓空間步長h變化,考察了空間精度；表6讓時間步長和空間步長同時變化,考察了格式的整體精度.通過三種格式的數(shù)值結(jié)果比較,仍然可以得到與理論分析相一致的結(jié)論.

表4 算例2當(dāng)h=0.012 5,t=0.5時,不同τ的L2范數(shù)誤差及收斂階Tab.4 L2-error and accuracy order with h=0.012 5 at t=0.5 for Example 2

表5 算例2當(dāng)τ=0.005, t=0.5時,不同h的L2范數(shù)誤差及收斂階Tab.5 L2-error and accuracy order with τ=0.005 at t=0.5 for Example 2

表6 算例2當(dāng)τ=h, t=0.5時的L2范數(shù)誤差及收斂階Tab.6 L2-norm error and accuracy order with τ=h at t=0.5 for Example 2

5 結(jié)論

本文提出了數(shù)值求解三維熱傳導(dǎo)方程的兩種局部一維化算子分裂方法,該方法把三維問題轉(zhuǎn)化為三個一維問題來求解.空間方向分別采用四階和六階緊致差分離散,而時間方向分別采用(2,2) Padé式和(3,3) Padé公式進(jìn)行離散,最終獲得了整體具有四階精度和六階精度的兩種差分格式,其截?cái)嗾`差分別為O(τ4+h4)和O(τ6+h6).通過矩陣?yán)碚摲治鰢?yán)格證明了兩種格式均是無條件穩(wěn)定的.通過數(shù)值實(shí)驗(yàn)驗(yàn)證了兩種格式均可達(dá)到各自的理論精度.本文方法較已有文獻(xiàn)中方法[8,11-13,16-17,19-21]的優(yōu)勢在于時間方向也達(dá)到了高精度,并且由于采用了局部一維化方法,因此同時保證了算法的高效性.

致謝：本文得到寧夏師范學(xué)院校級科研項(xiàng)目(NXSFZDA2001)的資助,在此表示感謝.