Burgers-Fisher方程改進的交替分段Crank-Nicolson并行差分方法

潘悅悅,吳立飛,楊曉忠

(1.華北電力大學(xué) 控制與計算機工程學(xué)院,北京 102206;2.華北電力大學(xué) 數(shù)理學(xué)院,北京 102206)

§1 引言

本文研究數(shù)學(xué)物理中一類重要的偏微分方程,Burgers-Fisher(B-F)方程[1-2]:

B-F方程是刻畫擴散傳播,對流傳導(dǎo)作用的典型模型,可以大量用于氣體動力學(xué)和熱傳導(dǎo)的研究[3].其中,α ＞0,μ ＞0,β ＜0為常數(shù).當(dāng)β=0時,方程(1)為經(jīng)典的Burgers方程;當(dāng)α=0時,方程(1)為一維Fisher方程.

在Burgers-Fisher方程(1)的實際計算中,通常取有限的計算區(qū)間(a,b).本文考慮Burgers-Fisher方程(1)的如下初邊值問題,

初始條件:

邊界條件:

分別取確定的f(x),g1(t)和g2(t),

可得到解析解為

對于Burgers方程的數(shù)值研究,已有大量研究成果[4-8].但關(guān)于Burgers-Fisher方程的并行數(shù)值方法研究尚不多見.Kaya和El-Sayed(2004年)[9]研究了廣義Burgers-Fisher方程的分解格式,在不使用任何離散技術(shù)的情況下,得到數(shù)值解,并與精確解比較,證明該方法具有較高精度.Zhang等(2014年)[10]利用孤立波解提出了Burgers-Fisher方程的一個新精確差分格式和非標準有限差分格式,并證明了其精度.Namjoo等(2018年)[11]考慮了廣義Burgers-Fisher方程的非標準有限差分格式,討論了其正定性,一致性和有界性.然而,文獻中的大部分格式只考慮了計算精度,忽略了格式的計算效率.

隨著大型并行計算機的迅速發(fā)展,研究和構(gòu)造適用于大型并行計算機的計算方法成為主流.大規(guī)模并行計算的瓶頸之一是全局通訊耗時多,同步等待時間長.為了保證計算過程無需進行全局通訊,人們開始研究能適應(yīng)大型并行計算機的并行數(shù)值解法.Zhang和Li(1994年)[12]將具有理想并行性的顯格式和穩(wěn)定性好的隱格式,Crank-Nicolson格式結(jié)合,建立了兼具穩(wěn)定性和并行性的交替分段Crank-Nicolson格式.周毓麟(1997年)[13]建立了具有并行本性的差分方法的基本理論,針對擬線性拋物型方程組,構(gòu)造了非均勻網(wǎng)格上的具有并行本性的差分格式,證明了差分解的先驗估計,從而得到離散解的存在性,唯一性,收斂性和穩(wěn)定性.受到構(gòu)造并行本性差分方法的啟示,目前,許多發(fā)展方程都可以用此并行差分方法進行數(shù)值求解,并且可以直接和有效地在并行機上應(yīng)用[14-15].Sheng等(2007年)[16]提出兩個并行本性差分格式,通過使用前兩個時間層在邊界點的值,得到內(nèi)邊界條件,用全隱式計算子域中的值,然后更新邊界值,得到全局并行算法.曲富麗和王文洽(2007年)[17]設(shè)計了KdV方程的交替分段顯-隱差分格式,理論分析格式的線性絕對穩(wěn)定性,并通過數(shù)值試驗驗證了格式具有較高的計算精度.Guo和Liu(2013年)[18]構(gòu)造了四階拋物方程的一般交替差分格式,數(shù)值試驗以交替分組顯式,交替分段顯-隱式和交替分段Crank-Nicolson格式為例,驗證了格式的收斂性,穩(wěn)定性和計算精度.Xue和Feng(2018年)[19]提出了對流占優(yōu)擴散方程的一種新并行差分方法,在相同時間層使用兩種區(qū)域分解法,并取平均值得到新的數(shù)值解,數(shù)值算例證明算法具有較高的計算精度且適用于并行計算.Pan等(2020年)[20]構(gòu)造了Burgers-Fisher方程的交替分段純顯-隱和交替分段純隱-顯格式,數(shù)值試驗證明PASE-I(PASI-E)差分格式具有空間和時間二階精度,獲得了穩(wěn)定性和并行性兼顧的結(jié)果.并行本性差分格式可以直接應(yīng)用到分布式內(nèi)存的并行計算機系統(tǒng),使得處理機之間的通信量極小,容易平衡計算和通信,從而得到好的精度和并行計算的可擴展性.本文將以這類傳統(tǒng)并行差分格式為基礎(chǔ),構(gòu)造改進的交替分段Crank-Nicolson差分方法求解Burgers-Fisher方程.

本文結(jié)構(gòu)安排如下:§2將Crank-Nicolson格式做分段化處理,分段點處交替使用顯格式和隱格式,得到改進的交替分段Crank-Nicolson(IASC-N)格式.§3給出Burgers-Fisher方程IASC-N并行差分格式的理論分析.§4的數(shù)值試驗結(jié)果驗證理論分析的正確性和IASC-N格式計算效率的優(yōu)越性.

§2 IASC-N并行差分格式

對求解區(qū)域[0,L]×[0,T]進行網(wǎng)格剖分,取空間步長和時間步長,其中M,N是正整數(shù).記xi=ih(i=0,1,2……·M),tn=nτ(n=0,1,2……·),用表示解析解U(xi,tn)的數(shù)值解.為了構(gòu)造方程(1)的IASC-N格式,先分別給出以下三種離散格式.

(1) 古典顯式格式

(2) 古典隱式格式

(3) Crank-Nicolson格式

現(xiàn)在討論Burgers-Fisher方程的交替分段并行差分格式,本文的IASC-N格式將結(jié)合以上幾種差分格式使用.

設(shè)M -1=Jl,(J ≥3,l ≥3為正整數(shù),且J為奇數(shù)),將同一時間層上計算的點分為J段,記為S1,S2,……,SJ.以M=26,J=5,l=5為例,定義每段的內(nèi)邊界點及內(nèi)點,構(gòu)造IASC-N格式.

S1段有一個內(nèi)邊界點=6及4個內(nèi)點i1=2,3,4,5;S2段有兩個內(nèi)邊界點=7,=11及3個內(nèi)點i2=8,9,10;S3段有兩個內(nèi)邊界點=12,=16及3個內(nèi)點i3=13,14,15;S4段有兩個內(nèi)邊界點=17,=21及3個內(nèi)點i4=18,19,20;S5段有一個內(nèi)邊界點=22及4個內(nèi)點i5=23,24,25,26.

設(shè)n為偶數(shù),在第n層上的數(shù)值解為已知,現(xiàn)計算.n+1時間層上,在點xi(i=處使用古典顯式格式(4),在點處使用古典隱式格式(5),其余各點處使用經(jīng)典C-N格式(6).n+2時間層上,將n+1時間層顯式格式(4)的位置變?yōu)殡[式格式(5),隱式格式(5)的位置變?yōu)轱@式格式(4),其余各點仍使用經(jīng)典C-N格式(6).交替使用n+1層和n+2層上的格式,可得到IASC-N格式.具體參見圖1,□采用(4)式,■采用(5)式,*采用(6)式,其矩陣形式為

圖1 IASC-N格式構(gòu)造示意圖

式中

Un,F1,F2都是(M -1)維向量,A=-r1-r2,B=-r1+r2,I是(M -1)階單位矩陣.A1和A2均是(M -1)階對角矩陣,且滿足A1+A2=I.A1=diag(θ1,θ2,……,θM-2,θM-1),其中

§3 IASC-N并行差分方法的數(shù)值分析

3.1 IASC-N差分格式解的存在唯一性

為討論IASC-N格式解的存在唯一性,需要引入如下兩個引理.

引理1(Kellogg引理[21]) 設(shè)θ ＞0,矩陣(G+GT)是非負定的,則(θI+G)-1存在,并且有

引理2由IASC-N格式(7)定義的矩陣A1G和A2G是非負定矩陣.

證明引理2,只需證明A1G+(A1G)T和A2G+(A2G)T是非負定矩陣.作如下計算

考慮低速流情況,不妨設(shè)0＜u ＜1.由于β ＜0,r3＜0,則A1G+(A1G)T是對角占優(yōu)的三對角矩陣,其主對角元素為正實數(shù).顯然A1G+(A1G)T是非負定矩陣,同理A2G+(A2G)T也是非負定矩陣,故A1G和A2G是非負定矩陣.引理2證畢.

由Burgers-Fisher方程的初邊值條件可知,U0的值已知.假設(shè)第2n層已經(jīng)確定,求解第2n+1層差分方法的矩陣方程為

由引理1和引理2,等式右側(cè)已知,且(I+A1G)-1存在,所以方程(8)有唯一解.

同樣地,用IASC-N方法求解第2n+2層,矩陣方程為

同上分析可得,方程(9)有唯一解.綜上分析,有如下定理.

定理1Burgers-Fisher方程IASC-N格式(7)的解是存在且唯一的.

3.2 IASC-N差分格式的線性絕對穩(wěn)定性

引理3[21]設(shè)θ ＞0,矩陣(G+GT)是非負定的,則‖(θI -G)(θI+G)-1‖2≤1.

以下穩(wěn)定性分析中,假定(7)式中的系數(shù)uδ=a為常數(shù)且0＜a ＜1.消去Un+1,公式(7)可以改寫成Un=GUn-2.其中G為增長矩陣且

對任何偶數(shù)n,有

由引理1-3,對任意的r1,r2和r3,有

其中r3＜0.因此

其中C=1+4r1+2r2-r3,即‖Gn‖2≤C.定理得證.

定理2Burgers-Fisher方程的IASC-N格式(7)是線性絕對穩(wěn)定的.

3.3 IASC-N差分格式的計算精度

古典顯式格式:

古典隱式格式:

將以上兩格式中的各點分別在點un+1i處作Taylor級數(shù)展開,記截斷誤差分別為T1(τ,h),T2(τ,h),則有

由Burgers-Fisher方程IASC-N差分格式的分段化設(shè)計可以看出,在每段的“內(nèi)點”處使用經(jīng)典的C-N格式,其計算精度為二階.在每段的“內(nèi)邊界點”處分別使用古典顯式格式和古典隱式格式,在T1(τ,h)和T2(τ,h)的表達式中,分別包含了絕對值相同但符號相反的項,古典顯式和古典隱式不僅在同一時間層交替進行,在不同時間層也交替進行,部分誤差項將相互抵消,因此,IASC-N格式的計算精度為O(h2+τ2).

定理3Burgers-Fisher方程的IASC-N格式(7)計算精度為O(h2+τ2).

3.4 IASC-N差分格式的收斂性

考慮如下非線性方程:

假設(shè)u(k)(x)是方程(10)的第k個解,對于非線性方程(10),應(yīng)用擬線性過程,引入一個由以下遞推關(guān)系確定的線性方程序列.

其中k=0,1,2,……·是迭代索引.在函數(shù)空間,應(yīng)用Newton-Raphson-Kantorovich近似方法[22],選擇一個滿足初始條件f(x)的合理初始值u(0)(x).令u(k+1)=得到:

此外,假設(shè)函數(shù)a(x),b(x)和c(x)是空間方向上足夠光滑的函數(shù),且有

其中ν是終止計算的小規(guī)定值.是非線性邊值問題(10)的數(shù)值解.

在第n+1時間層上,令是由擬線性技術(shù)得到的序列.為了證明IASC-N格式的收斂性,考慮以下方程:

通過使用格林函數(shù),將以上方程轉(zhuǎn)化為積分函數(shù)

其中格林函數(shù)定義如下

取空間域上的最大范數(shù),簡化后得到

定理得證.

定理4Burgers-Fisher方程的IASC-N格式(7)是收斂的.

§4 數(shù)值試驗

數(shù)值試驗基于Intel Core i5-4200 CPU@2.50GHz,四核處理器,在Matlab R2014a環(huán)境下運行.

例1取α=0.1,μ=1,β=-0.25,δ=1.考慮如下Burgers-Fisher方程[23]:

滿足以下初邊值條件:

得到解析解為:

取時間層N=600,空間層M=101,分段數(shù)J=5,內(nèi)點數(shù),解析解曲面,C-N格式解曲面,IASC-N格式解曲面如下.圖2至圖4可以看出它們的形狀整體保持一致.在t=0.5處,將本文IASC-N格式的數(shù)值解與解析解,C-N格式解進行比較.計算結(jié)果如表1,可以看出IASC-N格式的數(shù)值解很好地逼近了解析解,且u的值在0到1之間,符合理論分析.把看作解析解,看作數(shù)值解,定義節(jié)點誤差(Node Error,NE),,圖5和圖6分別是C-N格式解和IASC-N格式解相對于解析解的節(jié)點誤差分布圖.從圖中可以看到,IASC-N格式的計算誤差與C-N格式的計算誤差相似,且最大誤差不超過6e-04.

圖5 C-N格式解的節(jié)點誤差分布圖

圖6 IASC-N格式解的節(jié)點誤差分布圖

表1 數(shù)值解與解析解的比較

圖2 解析解曲面

圖3 C-N格式解曲面

圖4 IASC-N格式解曲面

表2和表3分別是取M=601,當(dāng)N=100,200,300,400,500時,C-N格式解,IASC-N格式解相對于解析解的節(jié)點誤差.從表中數(shù)據(jù)可以看出,隨著時間層數(shù)的增加,當(dāng)0＜t ≤0.4 時,IASC-N格式解與解析解的節(jié)點誤差逐漸減小.當(dāng)0.4＜t ＜1時,IASC-N格式解與解析解的節(jié)點誤差雖然有少量增加,但都小于C-N格式解與解析解的節(jié)點誤差,且最大誤差均不超過6e-04.說明IASC-N并行差分方法是求解Burgers-Fisher方程的一種高精度差分方法.

表2 C-N格式解與解析解的節(jié)點誤差

表3 IASC-N格式解與解析解的節(jié)點誤差

時間收斂階的計算結(jié)果如表4,IASC-N格式在時間方向上的收斂階接近2階.從表5可以看到,Orderx的值接近于2.因此得到IASC-N格式的收斂階是O(h2+τ2).數(shù)值結(jié)果符合理論分析.

表4 兩種格式的數(shù)值誤差和時間收斂階

表5 兩種格式的數(shù)值誤差和空間收斂階

最后,取N=600,當(dāng)空間網(wǎng)格不斷加密時,比較和分析兩種格式的計算時間.圖7是IASCN格式和C-N格式計算時間的比較.從圖中可以看到,當(dāng)M ＜700時,兩格式的計算時間很接近,數(shù)據(jù)通信對循環(huán)產(chǎn)生的影響降低了IASC-N格式的計算效率,并行計算的優(yōu)越性不明顯.但隨著網(wǎng)格點數(shù)的增加,C-N格式的計算時間迅速增加,而IASC-N格式的計算時間增加緩慢,程序循環(huán)體執(zhí)行的影響遠大于數(shù)據(jù)通信的影響,并行差分格式高效率的優(yōu)勢越來越明顯.

圖7 IASC-N格式與C-N格式計算時間比較

比較C-N格式和IASC-N格式的計算速度,定義加速比Sp=T1/Tp和效率Ep=Sp/p(p=4)[26-27].固定時間層N=600,當(dāng)空間格點數(shù)分別為301,701,1101,1501,1901,2301時,兩種格式計算時間如表6.

表6 兩種格式計算時間及加速比和效率

由表6可以看出,當(dāng)空間格點數(shù)M ＞301時,IASC-N格式的計算時間一直比C-N格式計算時間少,且加速比Sp總是大于1.2.當(dāng)空間格點數(shù)M ＞701時,加速比和效率呈現(xiàn)遞增趨勢.隨著空間格點數(shù)的增加,IASC-N格式并行計算的特性更加突出.當(dāng)空間格點數(shù)M ＞1501時,和串行差分格式相比,IASC-N格式的計算時間大約能節(jié)約40%.

例2取L=1,T=1,.考慮Burgers-Fisher方程

解析解為

圖8是時間層N=600,空間層M=601,t=0.5時,三種格式解的比較,顯然IASC-N格式解曲線十分逼近解析解曲線,因此IASC-N并行差分方法是求解Burgers-Fisher方程的一種高精度差分方法.

圖8 三種格式解的比較

為了驗證IASC-N格式的穩(wěn)定性和計算精度,將解析解看作控制解,C-N格式解和IASCN格式解看作擾動解,定義時間層相對誤差和(the Sumof Relative Error for every Time level,SRET)和空間層能量誤差(Difference Total Energy,DTE):

取N=601,M=601,圖9可以看出,當(dāng)時間層數(shù)較小時,格式解的誤差稍大,但隨著時間層數(shù)的增加,兩種格式解的SRET迅速減小并趨向于0,且IASC-N格式解的SRET始終小于C-N格式解的SRET,說明IASC-N格式是計算穩(wěn)定的.圖10中兩種格式的DTE曲線形狀相似,且均不超過3e-06,說明IASC-N并行差分方法能夠有效求解Burgers-Fisher方程.

圖9 C-N格式解和IASC-N格式解的SRET

圖10 C-N格式解和IASC-N格式解的DTE

§5 結(jié)論

本文構(gòu)造的Burgers-Fisher方程IASC-N并行差分方法是一類典型的離散區(qū)域分裂法,既有網(wǎng)格區(qū)域分裂也有差分算子分裂.IASC-N并行差分方法是在傳統(tǒng)的有限差分理論基礎(chǔ)上發(fā)展起來的,便于直接地在多種類型的并行計算機上使用.

同時,本文用非線性方程線性化的方法證明IASC-N差分格式的線性絕對穩(wěn)定性和收斂性.特別當(dāng)空間分點數(shù)足夠大時,IASC-N差分格式在計算和通信方面都具有明顯的局部化特征,并行特性明顯.IASC-N并行差分方法對于求解Burgers-Fisher方程是高效可行的.