黃榮芳，孫建強，蔣朝龍

(海南大學(xué) 信息科學(xué)技術(shù)學(xué)院，海南 ?？?70228)

考慮如下的“good”Boussinesq 方程

方程(1)描述了非線性淺水波在2 個方向的傳播，具有孤立子相互作用的互動機制［1－2］.孤立波只存在有限范圍的速度，孤子可以保持其形狀和速度碰撞后小振幅孤子.然而對于大型振幅孤子，孤子可能發(fā)展成所謂的反孤立子.在周期或者零邊界條件下，方程(1)有如下的守恒特性

其中vx=ut.

對“good”Boussinesq 方程已經(jīng)有大量的理論和計算方法研究:El－Zoheriry［3］構(gòu)造了“good”Boussinesq方程的有限差分格式，并做了穩(wěn)定性分析，保結(jié)構(gòu)算法在求解“good”Boussinesq 方程具有顯著的優(yōu)勢;Aydina 和Kara?zen［4］構(gòu)造了“good”Boussinesq 方程的辛和多辛LO－BATTO 格式;曾文平［5］和Huang［6］等構(gòu)造了“good”Boussinesq 方程的多辛preissman 格式;胡偉鵬［7］等研究了廣義Boussinesq 方程的多辛格式;蔡家祥［8］等構(gòu)造了“good”Boussinesq 方程的局部保結(jié)構(gòu)算法，并取得了很好的數(shù)值結(jié)果.然而，在數(shù)值求解時，能精確保持“good”Boussinesq 能量守恒特性的數(shù)值算法很少. 因此在數(shù)值求解時精確保持“good”Boussinesq 方程的能量對正確地模擬方程具有重要的意義.

最近在保結(jié)構(gòu)算法領(lǐng)域內(nèi)，Quispel［9］和Mclachlan［10］提出了精確地保持Hamilton 系統(tǒng)能量的平均向量場方法.平均向量場方法已被應(yīng)用于KdV 方程，麥克斯韋方程等的求解［11］.筆者利用平均向量場方法求解“good”Boussinesq 方程，并數(shù)值模擬孤立波在不同振幅下的演化行為和能量守恒特性.

1 平均向量場方法

對給定的Hamilton 系統(tǒng)

可知Hamilton 系統(tǒng)具有能量守恒的特性.

對式(4)在時間方向進行離散

其中 H(zn+1，zn)是Hamilton 系統(tǒng)能量函數(shù)的離散梯度為反對稱矩陣.

定理1 離散梯度格式(8)保持Hamilton 系統(tǒng)能量守恒.

證明 由離散梯度的定義可知

由式(8)，可知

對于給定的Hamilton 系統(tǒng)，Quispel 和Mclachlan 給出了精確保Hamilton 系統(tǒng)能量守恒的二階平均向量場方法

其中h 是時間步長［3］.

定理2 平均向量場能夠精確保持Hamilton 系統(tǒng)的離散能量守恒.

證明 方程(11)可以改寫成

即

由微積分基本定理，可以得到(H(zn+1)－H(zn))/h=0.

因此，平均向量場方法可以在每個時間層上保持Hamilton 系統(tǒng)能量守恒.

2 “good”Boussinesq 方程的平均向量場格式

方程(1)可以寫成無限維Hamilton 系統(tǒng)形式

其中z=(u，v)T，Hamilton 函數(shù)為

假設(shè)空間積分區(qū)間Ω=［a，b］，空間長度L=a－b.將區(qū)間Ω=［a，b］分為N 等分，其中N 為正偶數(shù)，h=L/N 為空間步長.xj=a+hj，j=0，…，N－1 為空間配置點，uj和vj是對u(x，t)和v(x，t)在配置點xj的近似.

定義2

為插值空間，其中g(shù)j(x)可以被顯示表示為

定義如下的插值算子IN，對于任意

由于式(18)的正交性，可知

下面用uj來表示導(dǎo)數(shù)的值

同理可得

用D1近似?x.在空間上，對式(15)進行譜離散，可以得到“good”Boussinesq 方程的Fourier 擬譜半離散形式

由式(21)和(22)可知式(25)和(26)等價于

其中A 是一階譜微分矩陣D1和二階譜微分矩陣D2的乘積，式(27)和(28)可以寫成如下的Hamilton 形式

其中Z=(u，v)T，相應(yīng)的離散Hamilton 能量函數(shù)為

在時間方向上，將式(11)應(yīng)用到式(29)得

對式(31)和(32)積分，可得到?！癵ood”Boussinesq 方程離散能量的二階平均向量場格式

3 數(shù)值實驗

為了驗證上述理論分析的有效性，利用平均向量場格式數(shù)值模擬“good”Boussinesq 方程在不同振幅下的孤立波的演化行為.定義相對能量誤差如下

其中H(Z0)是初始離散能量，H(Zn)是在t=nΔt 時的離散能量.

首先，考慮利用平均向量場格式模擬孤立子演化過程.初始條件為

其中x0和σ 是實參數(shù)，A 是振幅，并且A=3σ2/2.

?。?5≤x≤75，x0=0，τ =0.02，N =240 和周期邊界條件，對孤立波在振幅A =1.48 時，進行數(shù)值模擬，如圖1 和圖2 所示.

圖1 孤立波在A=1.48 時的演化行為

圖2 在t [0,60]時的相對能量誤差圖

圖1 為“good”Boussinesq 方程在振幅A =1.48 時，孤子波在t ［0，60］時的演化行為.圖1 表明孤立波演化過程中，分裂成2 列波，2 列波分別向相反的方向運動，這與文獻［3］的結(jié)果是一致的.圖2 表明孤立波演化過程中的相對能量誤差圖.從圖2 可知，孤立波演化過程中相對能量誤差可以到達10－13，系統(tǒng)能量守恒.由此可知，平均向量場格式不僅能很好地模擬孤立波的演化行為，并且能精確保持方程的能量.

圖3 孤立波在A=1.55 時的演化行為

圖4 在t [0,8]時的相對能量誤差圖

圖3 為“good”Boussinesq 方程在振幅A=1.55 時，孤立波在t ［0，8］時的演化行為.從圖3 中可以知道當在振幅A=1.55 時，孤立波演化一段時間之后出現(xiàn)了爆破(blows－up)現(xiàn)象.這與文獻［3］的結(jié)果也是一致的. 圖4 表明孤立波演化過程中方程相對能量誤差圖. 從圖4 中可以看出孤立波演化過程中，即便出現(xiàn)了孤立波的爆破(blows－up)現(xiàn)象，但是方程的能量相對誤差可以達到10－13，系統(tǒng)能量守恒.

4 小 結(jié)

本文構(gòu)造了“good”Boussinesq 方程的平均向量場格式，利用平均向量場格式研究了孤立波在不同振幅條件下的演化行為. 數(shù)值結(jié)果表明，平均向量場格式不僅能很好地模擬孤立子波的演化行為，并且能精確保持方程能量. 平均向量場算法為數(shù)值模擬具有能量守恒的微分方程提供了新的選擇.

［1］Manoranjan V S，Mitchell A R ，Morris J L L. Numerical solution of“good”Boussinesq equation［J］. SIMA J. Sci. Stat.Comput，1984，5(4):946 －957.

［2］Frutos DE J，Ortega T，Sanz－Serna J M. Pseudo－spectral method for the“good”Boussinesq equation［J］. Math. Comput，1991，57(195):109 －122.

［3］El－Zoheiry H. Numerical investigation for the solitary waves interaction of the“good”Boussinesq equation［J］. Appl. Num.Math，2003，45(2/3):161 －173.

［4］Aydin A ，Karas?zen B.Symplectic and multi－symplectic Lobatto methods for the ”good”Boussinesq equation［J］. J. Math.Phy，2008，49(8):083509.

［5］Huang L Y，Zeng W P，Qin M Z. A new multi－symplectic scheme for nonlinear“good”Boussinesq equation［J］. J. Comput.Math.，2003，21(6):703 －714.

［6］曾文平，黃浪揚，秦孟兆. “good”Boussinesq 方程的多辛算法［J］. 應(yīng)用數(shù)學(xué)和力學(xué)，2002 ，23(27):743 －748.

［7］胡偉鵬，鄧子辰.廣義Boussinesq 方程的多辛算法［J］. 應(yīng)用數(shù)學(xué)和力學(xué)，2008，29(7):839 －845.

［8］Cai J X，Wang Y S. Local structure－preserving algorithms for the“good”Boussinesq equation［J］. J. Comput. Phy，2013(239):72 －89.

［9］Quispel G R W ，Mclachlan D I. A new class of energy－preserving numerical integration methods［J］.J Phys. A:Math. Theor，2008，41(4):1 －7.

［10］Mclachlan D I，Quispel G R W，Robidoux N. Geometric integration using discrete gradients［J］.Phil. Trans . Roy. Soc.A. 1999，357(1754);1 021 －1 045.

［11］Celledoni E，Grimm R I，Mclachlan D I，et al. Preserving energy resp. dissipation in numerical PDEs using the average vector field method［J］. J. Comput. Phys.，2012，231(20):6 770 －6 789.