許 潔， 謝樹森

(中國海洋大學(xué)數(shù)學(xué)科學(xué)學(xué)院， 山東 青島 266100)

1984年，Seyler通過對(duì)弱非線性作用下等離子聲波和空間電荷波的傳播的分析，提出了正則長波方程的對(duì)稱描述，即對(duì)稱正則長波(Symmetric regularized long wave,SRLW)方程[1]：

一般情況下，通常無法得到此類非線性模型的解析解。因此，數(shù)值方法成為研究SRLW方程的一個(gè)重要手段。

關(guān)于SRLW方程的數(shù)值方法已經(jīng)有了一些相應(yīng)的研究。1987年，郭柏靈[2]用譜方法研究了一類SRLW方程的周期邊值問題。1989年，鄭家棟等[3]對(duì)SRLW方程提出了帶約束算子的Fourier擬譜方法，并證明了該方法的穩(wěn)定性和最優(yōu)誤差估計(jì)。1995年，任宗修[4]考慮了SRLW方程的Chebyshev擬譜方法，構(gòu)造了半離散和全離散Chebyshev擬譜格式，并得出了相應(yīng)的誤差估計(jì)。2003年，尚亞東和郭柏靈[5]分析了多維廣義SRLW方程的Chebyshev擬譜格式。2006年，孔令華等[6]用時(shí)間上的Euler中點(diǎn)格式和空間上的Fourier擬譜方法對(duì)SRLW方程構(gòu)造了一個(gè)多辛Fourier擬譜格式，并證明了該格式的離散守恒定律。

塊中心有限差分法(Block-centered finite differe-nce,BCFD)又稱為單元中心有限差分法，通過采用適當(dāng)?shù)臄?shù)值求積公式，可以認(rèn)為是最低階的Raviart-Thomas混合元法[7]。Weiser和Wheeler[8]研究了一維和二維矩形區(qū)域中帶Neumann邊界條件橢圓問題的BCFD方法，證明了在足夠光滑的條件下，對(duì)于所有非均勻網(wǎng)格，所求橢圓問題的近似解及其一階導(dǎo)數(shù)中張量積BCFD的離散L2范數(shù)誤差都是二階的。在文獻(xiàn)[9]中，Arbogast、Wheeler和Yotov提出了以張量系數(shù)為塊中心差分的橢圓問題的混合有限元方法。此外，文獻(xiàn)[10-14]分別提出了求解多尺度流動(dòng)問題、半導(dǎo)體器件問題、含時(shí)間變系數(shù)拋物方程、對(duì)流-擴(kuò)散問題和海水入侵問題的BCFD方法。BCFD方法在一般非均勻空間網(wǎng)格上具有超收斂性，對(duì)未知量及其導(dǎo)數(shù)都能保持二階空間精度。因此，它被廣泛應(yīng)用于邊界層和大梯度變形問題。

1996年，許進(jìn)超[15]提出二重網(wǎng)格(Two-grid,TG)法思想，這種方法可以看作是一種兩層預(yù)測校正方法。二重網(wǎng)格法不是在細(xì)網(wǎng)格上求解大規(guī)模非線性問題，而是在粗網(wǎng)格上求解小規(guī)模非線性系統(tǒng)，在細(xì)網(wǎng)格上求解大規(guī)模線性系統(tǒng)，這樣在保證良好穩(wěn)定性的同時(shí)，又減少計(jì)算量。因此，如文獻(xiàn)[16-18]所示，TG方法由于其高效性和有效性，特別是在非線性模型的大規(guī)模建模和仿真方面，受到了廣泛的關(guān)注。

本文給出了Crank-Nicolson二重網(wǎng)格塊中心有限差分(CN-TG-BCFD)的全離散格式。對(duì)模型進(jìn)行了數(shù)值實(shí)驗(yàn)，驗(yàn)證了該方法在均勻和非均勻網(wǎng)格上的二階收斂性。將CN-TG-BCFD格式得到的結(jié)果與Crank-Nicolson完全非線性塊中心差分(CN-BCFD)方法比較，兩者在同一離散范數(shù)下具有相同量級(jí)的誤差，但CN-TG-BCFD方法效率更高，且在計(jì)算大規(guī)模問題中更明顯。

1 問題及記號(hào)

上述SRLW方程可寫成如下等價(jià)形式：

(1)

其中(x,t)∈I×J:(a,b)×(0,T]。 邊界和初值條件分別為

ux(a,t)=ux(b,t)=ρ(a,t)=ρ(b,t)=0,t∈J,

其中u0(x)和ρ0(x)為兩個(gè)已知的光滑函數(shù)，其中ρ和u分別是無量綱電子電荷密度和流體速度[13]。該系統(tǒng)描述了弱非線性(1+1)維離子聲波和空間電荷波。

πH:a=x1/2

假設(shè)網(wǎng)格劃分是正則的，即存在一個(gè)正常量σ，使得

在粗網(wǎng)格上，定義如下離散內(nèi)積和范數(shù)

ΠHv(x)定義為v(x) 的分段線性插值：

x∈[xi,xi+1],1≤i≤N-1。

ΠHv(x1/2) 定義為：

定義 LHv(x)為v(x) 的分段線性插值：

x∈[xi-1/2,xi+1/2],1≤i≤N。

2 全離散二重網(wǎng)格塊中心差分格式

令p=-ux，則式(1)等價(jià)寫成如下形式

(2)

(3)

(4)

基本步驟可以概括為：

(1)利用初始條件，在粗網(wǎng)格上解非線性方程(3)，得到粗網(wǎng)格上的近似解。

(2)在細(xì)網(wǎng)格上解線性系統(tǒng)(4)，求出細(xì)網(wǎng)格點(diǎn)(yj,tm)處的數(shù)值解。

3 收斂性估計(jì)及穩(wěn)定性結(jié)果

利用截?cái)嗉夹g(shù)和逆估計(jì)可以證明如下誤差估計(jì)：

u,p,ρ∈C3([0,T];C3[a,b]),

且滿足條件

Δt=o(H1/4),

則存在一個(gè)與H和Δt無關(guān)的正常數(shù)C，使得

C(Δt2+h2+H7/2)。

注3由于差分格式是非線性的，利用文獻(xiàn)[19]中的截?cái)嗉记珊湍婀烙?jì)可以證明上述結(jié)論成立。

4 數(shù)值算例

在本節(jié)中，使用CN-TG-BCFD算法(3)和(4)分別在均勻和非均勻空間網(wǎng)格上進(jìn)行數(shù)值實(shí)驗(yàn)。此外，為了檢驗(yàn)二重網(wǎng)格方法的有效性，還將其與完全非線性CN-BCFD算法(3)在細(xì)網(wǎng)格上的結(jié)果進(jìn)行比較。

本文考慮的非均勻網(wǎng)格是均勻網(wǎng)格添加一個(gè)隨機(jī)擾動(dòng)后生成的。首先構(gòu)造一個(gè)均勻網(wǎng)格如下：

然后使用Matlab內(nèi)聯(lián)代碼對(duì)網(wǎng)格步長進(jìn)行隨機(jī)小擾動(dòng)如下所示：

Hi=H[1+0.01×ε×rand (i)],i=1,…,N-1。

定義非均勻網(wǎng)格點(diǎn)：

且HN=xN+1/2-xN-1/2。

選取初始條件為：

且可得方程的精確解為

(5)

模擬該孤立波在區(qū)間t∈[0,12]，x∈[-20,50]上的運(yùn)動(dòng)。

數(shù)值結(jié)果與精確解進(jìn)行比較，結(jié)果如表1～4和圖1～3所示。

表1和2分別為均勻網(wǎng)格下CN-TG-BCFD方法(3)、方法(4)和完全非線性CN-BCFD方法(3)的誤差、階數(shù)和CPU時(shí)間，誤差和CPU時(shí)間的曲線圖如圖2所示，其中h=H/5，Δt=h。 圖1顯示了Δt=h=0.1 時(shí)應(yīng)用二重網(wǎng)格法得到的數(shù)值解的運(yùn)動(dòng)，這與精確解(5)幾乎相同。觀察結(jié)果如下：(1)這兩種方法在離散范數(shù)中產(chǎn)生相同量級(jí)的誤差，并且在空間和時(shí)間上都具有二收斂性；(2)提出的二重網(wǎng)格方法更有效。例如，當(dāng)h=0.025時(shí)，非線性格式(3)花費(fèi)268 547 s(約為3 d)，而提出的二重網(wǎng)格方法(3)、方法(4)僅花費(fèi)不到17 min，對(duì)于較小的網(wǎng)格步長，差距更加明顯。因此，二重網(wǎng)格法是求解非線性SRLW方程的較好方法。

圖1 單孤立波數(shù)值解的運(yùn)動(dòng)(Δt=h=0.1)

(這里三角形正弦值為2。Here the tangent of the triangle is 2.)

表1 均勻網(wǎng)格下二重網(wǎng)格方法(3)和(4)的誤差和CPU時(shí)間

接下來，測試了大小為H′的非均勻網(wǎng)格上的收斂速度和CPU時(shí)間。設(shè)定隨機(jī)擾動(dòng)ε=0.05，h=H′/5和Δt=h/5，誤差和CPU時(shí)間顯示在表3和4中，并繪制在圖3中。在非均勻網(wǎng)格上仍然可以觀察到空間和時(shí)間上的二階精度。

(這里三角形正弦值為2。 Here the tangent of the triangle is 2.)

表2 均勻網(wǎng)格下非線性格式(3)的誤差和CPU時(shí)間

表3 非均勻網(wǎng)格下二重網(wǎng)格方法(3)和(4)的誤差和CPU時(shí)間

表4 非均勻網(wǎng)格下非線性格式(3)的誤差和CPU時(shí)間

5 結(jié)語

本文主要討論了對(duì)非線性SRLW方程運(yùn)用Crank-Nicolson二重網(wǎng)格塊中心有限差分方法進(jìn)行數(shù)值求解的問題。在空間上，采用了二重網(wǎng)格塊中心差分方法進(jìn)行離散，而對(duì)時(shí)間上，則采用了Crank-Nicolson方法進(jìn)行離散。同時(shí)，在均勻和非均勻網(wǎng)格上對(duì)孤立波問題進(jìn)行數(shù)值模擬，并將結(jié)果與完全非線性塊中心差分方法的數(shù)值結(jié)果進(jìn)行比較。數(shù)值實(shí)驗(yàn)表明，采用CN-TG-BC方法求解非線性孤立波方程是可行并有效的，并且這種方法與完全非線性差分方法相比更有優(yōu)勢(shì)。