“good”Boussinesq 方程的平均向量場方法

黃榮芳，孫建強，蔣朝龍

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

考慮如下的“good”Boussinesq 方程

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

其中vx=ut.

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

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

1 平均向量場方法

對給定的Hamilton 系統

可知Hamilton 系統具有能量守恒的特性.

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

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

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

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

由式(8)，可知

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

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

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

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

即

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

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

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

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

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

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

定義2

為插值空間，其中gj(x)可以被顯示表示為

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

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

下面用uj來表示導數的值

同理可得

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

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

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

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

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

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

3 數值實驗

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

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

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

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

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

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

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

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

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

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

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

4 小 結

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

［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］. 應用數學和力學，2002 ，23(27):743 －748.

［7］胡偉鵬，鄧子辰.廣義Boussinesq 方程的多辛算法［J］. 應用數學和力學，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.