求解一維對(duì)流方程的高精度緊致差分格式

侯波,葛永斌

(寧夏大學(xué)數(shù)學(xué)統(tǒng)計(jì)學(xué)院,寧夏 銀川750021)

1.引言

對(duì)流方程在生物數(shù)學(xué)、能源開(kāi)發(fā)、空氣動(dòng)力學(xué)等許多領(lǐng)域都具有十分廣泛的應(yīng)用,因此求解該類(lèi)方程具有非常重要的理論價(jià)值和實(shí)際意義.然而,由于實(shí)際問(wèn)題通常十分復(fù)雜,往往難以求得精確解,因此研究其精確穩(wěn)定的數(shù)值解法是十分必要的.

針對(duì)對(duì)流方程國(guó)內(nèi)外很多學(xué)者提出了很多的數(shù)值方法.如張?zhí)斓潞蛯O傳灼[1]針對(duì)一維對(duì)流方程采用待定系數(shù)法,得到了兩層四點(diǎn)格式和四階六點(diǎn)格式,并且是無(wú)條件穩(wěn)定的,該方法適用于在點(diǎn)數(shù)確定的前提下,得到精度高的差分格式; 于志玲和朱少紅[2]針對(duì)一維問(wèn)題建立了中間層為兩個(gè)節(jié)點(diǎn)的三層顯格式,其截?cái)嗾`差為O(τ2+h2); 曾文平[3]針對(duì)一維對(duì)流方程推導(dǎo)出了一種兩層半顯式格式，其截?cái)嗾`差為O(τ2+h2),該格式是無(wú)條件穩(wěn)定的.姚朝輝等人[5]將二階的迎風(fēng)格式和中心差分格式進(jìn)行加權(quán)得到了WSUC格式,該格式是無(wú)條件穩(wěn)定的;但該格式時(shí)間方向和空間方向僅有二階精度.湯寒松等人[6]通過(guò)立方插值擬質(zhì)點(diǎn)方法(CIP 方法),給出了一些保單調(diào)的CIP格式; Erdogan[9]針對(duì)一維的對(duì)流方程推導(dǎo)出了一種指數(shù)擬合的差分格式,其截?cái)嗾`差為O(τ2+h2); Bourchtein[10]構(gòu)造了對(duì)流方程的三層五點(diǎn)中心型蛙跳格式,該格式的截?cái)嗾`差為O(τ4+h4); 即該格式時(shí)間和空間均具有四階精度,但是該格式是三層的,空間方向需要五個(gè)點(diǎn),并且是條件穩(wěn)定的; Kim[11]構(gòu)造了多層無(wú)耗散的迎風(fēng)蛙跳格式,即時(shí)間和空間分別具有二階、四階、六階精度,但格式為三層甚至是四層的,并且六階格式空間方向最多需要五個(gè)點(diǎn),給靠近邊界的內(nèi)點(diǎn)的計(jì)算帶來(lái)困難.

綜上所述,文獻(xiàn)中已經(jīng)有的數(shù)值計(jì)算方法大多為低階精度的,而高精度方法涉及多個(gè)時(shí)間層,需要一個(gè)或多個(gè)時(shí)間啟動(dòng)步,或者空間方向的網(wǎng)格節(jié)點(diǎn)多于三個(gè),這都給計(jì)算造成困難或不便.為此本文將構(gòu)造一種緊致格式,這里緊致格式的定義為對(duì)時(shí)間導(dǎo)數(shù)項(xiàng)的離散采用不超過(guò)三個(gè)時(shí)間層,而對(duì)空間導(dǎo)數(shù)項(xiàng)的離散采用不超過(guò)三個(gè)網(wǎng)格點(diǎn),時(shí)間和空間即可以達(dá)到高階精度（三階及三階以上）的格式.本文擬構(gòu)造的格式時(shí)間方向僅用到兩個(gè)時(shí)間層上的函數(shù)值,在每個(gè)時(shí)間層上僅涉及到三個(gè)空間網(wǎng)格點(diǎn),格式時(shí)間和空間具有整體的四階精度.該格式的優(yōu)點(diǎn)是無(wú)須啟動(dòng)步的計(jì)算,并且在對(duì)靠近邊界點(diǎn)的計(jì)算時(shí),不會(huì)用到計(jì)算域以外的網(wǎng)格節(jié)點(diǎn).此外該格式為無(wú)條件穩(wěn)定的,可以采用比較大的時(shí)間步長(zhǎng)進(jìn)行計(jì)算.最后通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證本文格式的精確性和穩(wěn)定性.

2.差分格式的建立

考慮如下一維對(duì)流方程：

給定初始條件為:

給定周期性邊界條件為:

其中,u(x,t)為未知函數(shù),f為非齊次項(xiàng),a為對(duì)流項(xiàng)系數(shù),φ(x)為已知函數(shù).

將求解區(qū)域[b,c]等距剖分為N個(gè)子區(qū)間：b=x0,x1,··· ,xN?1,xN=c,并且定義時(shí)間也采用等距剖分,步長(zhǎng)用τ表示.在本文中,我們利用分別表示u在(xi,tn),點(diǎn)處的函數(shù)值.

其中

同理可得:

將(2.7)和(2.8)代入(2.4)整理可得:

為了使該格式在時(shí)間方向和空間方向上均達(dá)到四階精度,須對(duì)(2.9)式中的項(xiàng)進(jìn)行二階的離散,同時(shí)為了保證本文格式的緊致性,即空間方向不超過(guò)三個(gè)網(wǎng)格點(diǎn),我們對(duì)(2.1)式進(jìn)行如下變形:

從而可得:

將(2.13)代入(2.11)得:

即,

即,

由推導(dǎo)過(guò)程可知,該格式的截?cái)嗾`差為O(τ4+τ2h2+h4),即格式(2.15) 在時(shí)間和空間上均可達(dá)到四階精度.我們注意到,格式為兩層格式,并且格式每層僅用到三個(gè)網(wǎng)格點(diǎn),形成的代數(shù)方程組系數(shù)矩陣為循環(huán)三對(duì)角矩陣,可采用追趕法進(jìn)行求解[8],同時(shí)由于要求未知時(shí)間層上(第n+1層)中間點(diǎn)的系數(shù)不能等于0,即因此

3.穩(wěn)定性分析

下面采用von Neumann方法分析本文所推導(dǎo)的差分格式(2.15)的穩(wěn)定性.

對(duì)于(2.15)式,舍掉非齊次項(xiàng)f,即假設(shè)f項(xiàng)精確成立,令uni=ηneIσxi,其中,η為振幅,σ為波數(shù),為虛數(shù)單位,有

兩邊同時(shí)約掉eIσxi,并整理可得：

利用Euler公式,即eIσh=cosσh+Isinσh,e?Iσh=cosσh ?Isinσh,可得:

對(duì)上式進(jìn)行化簡(jiǎn)整理有

從而可得格式(2.15)的誤差放大因子為:

由von Numann穩(wěn)定性定理可知當(dāng)|G| ≤1 時(shí),格式是穩(wěn)定的,由(3.5) 可得|G|=1,因此,格式(2.15) 是無(wú)條件穩(wěn)定的.

4.數(shù)值實(shí)驗(yàn)

為了驗(yàn)證本文格式(2.15)的精確性和穩(wěn)定性,現(xiàn)考慮以下三個(gè)具有精確解的初邊值問(wèn)題.分別采用Crank-Nicolson(C-N)格式,文[7] 中格式和本文格式(2.15) 進(jìn)行計(jì)算; 其中,最大絕對(duì)誤差及收斂階的定義為:

L∞(h1)和L∞(h2) 為空間網(wǎng)格步長(zhǎng)分別為h1和h2時(shí)的最大絕對(duì)誤差.

問(wèn)題1[7]:

該問(wèn)題的精確解為:u(x,t)=sin[π(x ?t)].

表1 問(wèn)題1 當(dāng)λ=τ/h=0.5,t=1 時(shí)刻的最大絕對(duì)誤差及收斂階

表2 問(wèn)題1 當(dāng)τ=λh,t=2 時(shí)刻的最大絕對(duì)誤差

圖1 問(wèn)題1 當(dāng)N=32,τ=0.03125,t=0.2 時(shí)刻的數(shù)值解與精確解

表1給出了針對(duì)問(wèn)題1三種格式在不同空間步長(zhǎng)h下,當(dāng)λ=τ/h=0.5,t=1時(shí)的最大絕對(duì)誤差和收斂階.我們發(fā)現(xiàn)C-N格式在時(shí)間和空間上都為二階精度,由于文[7]格式時(shí)間具有二階精度,空間具有四階精度,因此當(dāng)取τ=O(h) 時(shí),格式空間僅有二階精度,而本文格式時(shí)間和空間均為四階精度.圖1給出N=32,τ=0.03125,t=0.2數(shù)值解與精確解對(duì)比圖,可以看出數(shù)值解與精確解吻合的很好.

表2給出了當(dāng)h=1/16和h=1/32時(shí),τ=λh,t=2時(shí)刻對(duì)問(wèn)題1采用三種格式計(jì)算的最大絕對(duì)誤差.可以看出網(wǎng)格比λ最大取到12.8,計(jì)算仍然是穩(wěn)定的,因此本文格式是無(wú)條件穩(wěn)定的.并且本文格式在所有參數(shù)下,其計(jì)算結(jié)果比C-N 格式和文[7]格式計(jì)算結(jié)果更加精確.

問(wèn)題2[7]:

該問(wèn)題的精確解為:u(x,t)=ecos[π(x?t)].

表3 問(wèn)題2 當(dāng)λ=τ/h=0.5,t=1 時(shí)刻的最大絕對(duì)誤差及收斂階

表4 問(wèn)題2 當(dāng)τ=λh,t=2 時(shí)刻的最大絕對(duì)誤差

表3和表4給出了針對(duì)問(wèn)題2利用本文格式和C-N格式以及文[7]格式的計(jì)算結(jié)果.表3考察了格式的精度,表4驗(yàn)證了格式的穩(wěn)定性.可以看出本文格式在時(shí)間和空間上均可達(dá)到四階精度,并且是無(wú)條件穩(wěn)定的.

問(wèn)題3

該問(wèn)題的精確解為:u(x,t)=cos[π(x ?t)].

表5 問(wèn)題3 當(dāng)λ=τ/h=0.5,a=0.5,t=1 時(shí)刻的最大絕對(duì)誤差及收斂階

問(wèn)題3為非齊次問(wèn)題,由于文[7]的方程模型為齊次方程,不能計(jì)算非齊次問(wèn)題,因此該問(wèn)題我們采用本文格式和C-N進(jìn)行計(jì)算和比較,表5給出了兩種格式在不同空間步長(zhǎng)h下,當(dāng)t=1時(shí)的最大絕對(duì)誤差和收斂階.可以看出當(dāng)λ=τ/h=0.5,a=0.5 時(shí),C-N格式在時(shí)間和空間上都為二階精度,而本文格式時(shí)間和空間均為四階精度.

5.結(jié)論

本文針對(duì)一維對(duì)流方程提出了一種兩層隱式高精度緊致差分格式,時(shí)間和空間均采用泰勒級(jí)數(shù)展開(kāi)法以及截?cái)嗾`差余項(xiàng)修正法進(jìn)行處理,格式截?cái)嗾`差為O(τ4+τ2h2+h4),即該格式在時(shí)間和空間上均可達(dá)到四階精度.并通過(guò)von Neumann方法分析得到該格式為無(wú)條件穩(wěn)定的.最后通過(guò)三個(gè)數(shù)值算例驗(yàn)證了格式的精確性和穩(wěn)定性.通過(guò)上述研究,我們可以得出如下結(jié)論:

1.文獻(xiàn)(如[10-11])中的高精度格式往往是時(shí)間多層格式,需要另外構(gòu)造啟動(dòng)步的計(jì)算格式,如果采用低精度格式啟動(dòng),必然會(huì)影響以后時(shí)間層的計(jì)算精度.而本文格式僅為兩層格式,無(wú)須啟動(dòng)步的計(jì)算,時(shí)間即可達(dá)到四階精度.

2.文獻(xiàn)(如[1,10-11])中的高精度格式空間方向上往往超過(guò)三個(gè)網(wǎng)格節(jié)點(diǎn),導(dǎo)致靠近邊界的內(nèi)點(diǎn)計(jì)算困難,需要采用特殊處理,而本文格式僅用到三個(gè)網(wǎng)格節(jié)點(diǎn),可以有效避免這一問(wèn)題.

3.盡管本文格式要求2,這是本文格式的一個(gè)缺陷,但是由于本文格式是無(wú)條件穩(wěn)定的,從理論上講可以采用任意網(wǎng)格比,因此可以很容易避開(kāi)2的條件限制,使得這一缺陷并不太影響格式的使用.

5.本文方法可直接推廣到二維和三維問(wèn)題中去,我們將另文報(bào)道.