許 潔， 謝樹森

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

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

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

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

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

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

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

1 問題及記號

上述SRLW方程可寫成如下等價形式：

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

πH:a=x1/2

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

在粗網(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)等價寫成如下形式

(2)

(3)

(4)

基本步驟可以概括為：

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

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

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

利用截斷技術(shù)和逆估計可以證明如下誤差估計：

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

且滿足條件

Δt=o(H1/4),

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

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

注3由于差分格式是非線性的，利用文獻[19]中的截斷技巧和逆估計可以證明上述結(jié)論成立。

4 數(shù)值算例

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

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

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

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

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

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

選取初始條件為：

且可得方程的精確解為

(5)

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)語

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