数学物理学报  2016, Vol. 36 Issue (6): 1092-1102   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
王俊杰
王连堂
一类DGH方程的多辛Fourier拟谱方法
王俊杰1, 王连堂2     
1. 普洱学院数学系 云南普洱 665000 ;
2. 西北大学数学系 西安 710127
摘要:DGH方程作为一类重要的非线性方程有着许多广泛的应用前景.通过正则变化,构造了DGH方程的多辛哈密尔顿系统.利用Fourier拟谱方法对此哈密尔顿系统进行数值离散,并构造了一种半隐式的多辛格式.数值算例结果表明该多辛离散格式具有较好的长时间数值稳定性.
关键词Hamilton系统     Fourier拟谱方法     多辛理论     DGH方程    
Multi-Symplectic Fourier Pseudospectral Method for a DGH Equation
Wang Junjie1, Wang Liantang2     
1. Mathematics Department of Pu Er University, Yunnan Pu'er 665000 ;
2. Mathematics Department of Northwest University, Xi'an 710127
Abstract: The DGH equation, an important nonlinear wave equation, has broad application prospect. With the canonical momenta, the multi-symplectic formulations for the DGH equation are presented. The multi-symplectic Fourier pseudospectral method is reviewed, and a semi-implicit scheme with certain discrete conservation laws is constructed to solve the DGH equation. The numerical experiments for DGH equation are given, showing that the multi-symplectic Fourier pseudospectral method is an efficient algorithm with excellent long-time numerical behaviors.
Key words: Hamilton system     Fourier pseudospectral method     Multi-symplectic theory     DGH equation    
1 引言

2001年, Dullin, Gottwald和Holm从Euler方程出发, 考虑浅水层自由表面的无旋不可压缩无粘水波运动时, 得到了一类带线性和非线性色散项的新型浅水波方程, 即DGH方程[1]

$ \begin{equation} u_t+2\omega u_x+3uu_x+ \gamma u_{xxx}-\alpha^2u_{xxt}-\alpha^2(uu_{xxx}+2u_xu_{xx})=0, \end{equation} $ (1.1)

其中$\omega, \gamma, \alpha$是常数. DGH方程被提出以后引起了广泛的关注.方程(1.1)具有一类奇异的孤立波解即尖峰孤立波$u(x, t)=c{\rm e}^{-|x-ct|}$, 尖峰孤立波在许多模型中相继被发现. DGH方程包含了两类相对独立的可积孤立波方程.当$\alpha\rightarrow 0$时, DGH方程变为KdV方程

$ \begin{equation} u_t+2\omega u_x+3uu_x+ \gamma u_{xxx}=0, \end{equation} $ (1.2)

KdV方程(1.2)具有光滑孤立波解

$ u(x, t)=u_0{\rm sech}^2(\sqrt{\frac{u_0}{\gamma}}(x-ct)), c=2\omega+u_0. $

$\gamma=0$时, DGH方程转化为另一类重要的可积方程, 即C-H方程

$ \begin{equation} u_t+2\omega u_x+3uu_x-\alpha^2u_{xxt}-\alpha^2(uu_{xxx}+2u_xu_{xx})=0. \end{equation} $ (1.3)

近几十年, 许多学者对方程(1.1)-(1.3)已经做了许多研究, 其中, Camassa和Holm在研究浅水波运动规律时, 推导出一类非线性水波方程(1.3).对任意的$ \omega $, C-H方程具有一个Lax对, 具有双哈密顿结构和无穷个守恒量[2-3].对$\omega=0 $, C-H方程有尖峰孤立波解, 更进一步研究表明, C-H方程具有多重尖峰孤波解.在文献[4-5]中, 作者对方程(1.3)进行了数值模拟, 且得到方程(1.3)的守恒性质.在文献[6]中, 作者对方程(1.3)进行了对称性和可积性研究.在文献[7]中, 作者对方程(1.3)的孤立子解和双孤立子解进行了研究.在文献[8]中, 作者对方程(1.1)的哈密尔顿结构和解的整体存在性及Blow-up现象进行了研究. Tian, Gui和Liu研究了DGH方程的Cauchy问题的局部适定性理论, 整体适定性理论, 解的极限行为, 孤立波的轨道稳定性, 散射理论, 新型尖峰孤立波解, 同时给出了DGH方程的散射数据.郭柏林和刘正荣通过使用平面的自治系统的定性分析研究方法, 研究了DGH方程的尖峰孤立波解.

由于问题(1.1)的非线性, 在实际应用中要求得方程(1.1)的初边值问题的精确解几乎不可能, 大部分情况下只能用数值方法来模拟方程(1.1), 但是对方程(1.1)的数值模拟的研究却很少.在文献[9]中, 给出了方程(1.1)数值方法的一些理论研究, 但是没有给出具体的数值方法.鉴于此, 本文的目的就是用数值方法来模拟方程(1.1)的初边值问题.如果采用一般的数值方法(如差分法、有限元法、谱方法), 这些方法都不是保结构的算法, 都会人为代入一定的耗散性, 都会破坏问题本来的一些物理性质, 在进行长期数值模拟时都会失真.由冯康先生开创的Hamilton系统的辛几何算法[10]在理论和应用上已日趋成熟并得到了蓬勃发展.现在, 一定程度上保持原问题的能力是评价一个数值算法优劣的标准.但用此方法对偏微分方程进行辛离散时具有局限性, 具体表现在此方法是个全局性的概念, 依赖于边界条件.为克服此局限性, Bridges和Reich引入了一个保持守恒型偏微分方程局部多辛结构的多辛积分的概念[11-13].本文的目的就是利用多辛算法来模拟方程(1.1), 这种算法可以进行长时间的数值模拟.由于这种算法不会破坏本身的物理性质, 所以本文的算法可以模拟孤立波解.

本文利用多辛Fourier拟谱方法对满足周期边界条件的DGH方程(1.1)进行了数值模拟.在第二节里通过引入正则动量, 验证了DGH方程(1.1)具有多辛哈密尔顿结构, 并证实此格式具有多辛守恒律、局部能量守恒和动量守恒.第三节给出了DGH方程(1.1)的离散多辛Fourier拟谱方法, 并证实此格式在离散格式下仍保持多辛守恒律.在第四节给出了数值模拟, 验证格式有长时间的稳定性.

2 DGH方程的多辛形式及守恒律

大量的偏微分方程可以表示为下列多辛哈密尔顿偏微分方程的形式[14-18]

$ \begin{equation} Mz_t+Kz_x=\nabla_zS(z), \end{equation} $ (2.1)

其中$M, K\in {\Bbb R}^{n\times n}(n\geq 3)$是反对称矩阵.$S:{\Bbb R}^n\rightarrow {\Bbb R}$是光滑函数, 称为哈密尔顿函数, $\nabla_zS(Z)$为函数$S(z)$的梯度. Bridges和Reich证明多辛哈密尔顿系统具有三个守恒律, 即多辛守恒律、局部能量守恒律、局部动量守恒律.

定理2.1  多辛哈密尔顿系统(2.1)满足多辛守恒律

$ \begin{equation} \frac{\partial}{\partial t}w+\frac{\partial}{\partial x}k=0, \end{equation} $ (2.2)

其中

$ w=\frac{1}{2}{\rm d}z\wedge M{\rm d}z, k=\frac{1}{2}{\rm d}z\wedge K{\rm d}z, $

分别对应于$t$$x$方向上的辛结构.

定理2.2  多辛哈密尔顿系统(2.1)具有局部能量守恒律

$ \begin{equation}\frac{\partial }{\partial t} E+\frac{\partial }{\partial x} F=0, \end{equation} $ (2.3)

其中

$ E=S(z)-\frac{1}{2}z^TKz_x, F=\frac{1}{2}z^TKz_t, $

局部动量守恒律

$ \begin{equation}\frac{\partial }{\partial t} I+\frac{\partial }{\partial x} G=0, \end{equation} $ (2.4)

其中

$ I=\frac{1}{2}z^TMz_x, G=S(z)-\frac{1}{2}z^TMz_t. $

$E$为能量密度; $F$为能量流; $I$为动量密度; $G$为动量流.

若引入矩阵分解

$ M=M_++M_-, K=K_++K_-, M_+^T=-M_-, K_+^T=-K_-, $

可以得到多辛守恒律、局部能量守恒律、局部动量守恒律的另一种表示, 其中

$ w={\rm d}z\wedge M_+{\rm d}z, k={\rm d}z\wedge K_+{\rm d}z, $
$ E=S(z)+z_x^TK_+z, F=-z^T_tK_+z, $
$ I=-z^T_xM_+z, G=S(z)+z^T_tM_+z. $

与辛结构相比, 多辛结构具有的守恒律, 在任意时空都成立, 不依赖于边界条件.大量的守恒方程都具有多辛结构, 例如KdV方程、KP方程、薛定谔方程等.下面我们建立DGH方程(1.1)的多辛结构.

引入正则动量

$ u_t=-2w_x, u=\phi_x, u_x=\psi, \frac{\alpha^2}{2}u_t-\frac{\gamma}{2}u_x=-\alpha^2u\psi+\Phi, $

系统(1.1)可以表示为下面一阶方程组

$ \begin{equation} \left\{ \begin{array}{lcl} \frac{1}{2}\phi_t-\frac{\alpha^2}{2}\psi_t+\frac{\gamma}{2}\psi_x-\Phi_x=-w-2\omega u-\frac{3}{2}u^2-\frac{\alpha^2}{2}\psi^2, \\ -\frac{1}{2}u_t+w_x=0, \\ -\phi_x=-u, \\ u_x=\psi, \\ \frac{\alpha^2}{2}u_t-\frac{\gamma}{2}u_x=-\alpha^2u\psi+\Phi. \end{array} \right.\end{equation} $ (2.5)

定义状态变量

$ z=(u, \phi, w, \Phi, \psi), $

系统(2.5)可以表示为多辛哈密尔顿偏微分方程的形式(2.1), 其中

$ M=\left[\begin{array}{ccccc} 0& ~~\frac{1}{2}~~&0&~~0~~&-\frac{\alpha^2}{2}\\ -\frac{1}{2}&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ \frac{\alpha^2}{2}&0&0&0&0 \end{array}\right], K=\left[\begin{array}{ccccc} 0&~~0~~&0&~~-1~~& \frac{\gamma}{2}\\ 0&0&1&0&0\\ 0&-1&0&0&0\\ 1&0&0&0&0\\ -\frac{\gamma}{2}&0&0&0&0 \end{array}\right], $

Hamilton函数

$ S(z)=-uw-\omega u^2-\frac{u^3}{2}-\frac{\alpha^2}{2}\psi^2u+\Phi\psi. $

系统(2.5)满足多辛守恒律(2.2), 其中

$ w=\frac{1}{2} ({\rm d}z\wedge M{\rm d}z)=\frac{1}{2}{\rm d}u\wedge {\rm d}\phi+\frac{\alpha^2}{2}{\rm d}\psi \wedge {\rm d}u, $
$ k=\frac{1}{2}({\rm d}z \wedge K {\rm d}z)= {\rm d}\Phi\wedge {\rm d}u +\frac{\gamma}{2} {\rm d}u\wedge {\rm d}\psi+{\rm d}\phi\wedge {\rm d}w. $

系统(2.5)具有能量守恒律(2.3), 其中

$ E=-\frac{1}{2}uw-\omega u^2-\frac{u^3}{2}-\frac{\alpha^2}{2}\psi^2u +\frac{1}{2}\Phi\psi-\frac{1}{2}(-\frac{\gamma}{2}\psi^2+\phi w_x-u\Phi_x+ \frac{\gamma}{2}u\psi_x), $
$ F=\frac{1}{2}(\Phi u_t-\frac{\gamma}{2}\psi u_t-w\phi_t+ \phi w_t-u\Phi_t+\frac{\gamma}{2}u\psi_t). $

系统(2.5)具有动量守恒律(2.4), 其中

$ I=-\frac{1}{4}\phi\psi+\frac{\alpha^2}{4}\psi^2+\frac{1}{4}u^2 -\frac{\alpha^2}{4}u\psi_x, $
$ G=-uw-\omega u^2-\frac{u^3}{2}-\frac{\alpha^2}{2}\psi^2u+\Phi\psi+\frac{1}{4}\phi u_t-\frac{\alpha^2}{4}\psi u_t-\frac{1}{4}u\phi_t+\frac{\alpha^2}{4}u\psi_t. $
3 DGH方程的多辛Fourier拟谱方法及离散守恒律

多辛是哈密尔顿偏微分方程的一个几何特征, 用数值方法模拟多辛哈密尔顿偏微分方程时, 自然希望能保持连续系统的这一特征.基于这个想法, 称一种能保持多辛守恒律的离散数值方法为多辛算法.

利用某种数值方法对多辛哈密尔顿系统(2.1)进行离散, 可以得到下面离散格式

$ \begin{equation} M\partial^{j, n}_tz_j^n+K\partial^{j, n}_xz_j^n=\nabla_zS(z_j^n), \end{equation} $ (3.1)

其中$z_j^n$为离散点$x=x_j, t=t_n$$z$的近似值, $ \partial_t^{j, n}, \partial_x^{j, n}$分别表示微分算子$\partial_t, \partial_x$在这种数值方法下离散点$x=x_j, t=t_n$处的离散值.

定义3.1  若一个数值方法满足离散多辛守恒律

$ \begin{equation} \partial^{j, n}_tw_j^n+\partial^{j, n}_xk_j^n=0, \end{equation} $ (3.2)

其中

$ w_j^n=\frac{1}{2}{\rm d}z_j^n\wedge M{\rm d}z_j^n, k_j^n=\frac{1}{2}{\rm d}z_j^n\wedge K{\rm d}z_j^n, $

称这样的算法为多辛算法.

Bridges和Reich结合传统谱方法在Fourier空间上提出了多辛Fourier谱离散格式, 基于这种思想, 陈景波和秦孟兆构造了发展方程的多辛Fourier拟谱格式.下面我们建立方程(1.1)的多辛Fourier拟谱方法.设$u(x, t)$是空间$[0, L]$上的周期函数, 即$u(0, t)=u(L, t)$.用$I_Nu(x, t)$表示为$u(x, t)$ 的插值函数, 配置点为$x_j=h\cdot j, j=0, 1, 2, \cdots N-1$, 其中$h =\frac{L}{N}$, N为偶数, 具体表达式为

$ \begin{equation} I_Nu(x, t)=\sum_{j=0}^{N-1}u_jg_j(x), \end{equation} $ (3.3)

其中$u_j=u(x_j, t)$, $g_j(x_k)=\delta_j^k$.$g_j(x)$$\frac{N}{2}$次正交三角多项式

$ \begin{equation} g_j(x)=\frac{1}{N}\sum_{l=-\frac{N}{2}}^\frac{N}{2}\frac{1}{C_l}{\rm e}^{{\rm i}l\mu(x-x_j)}, \end{equation} $ (3.4)

其中$C_l=1(|l|\neq \frac{N}{2}), C_{\frac{N}{2}}=C_{-\frac{N}{2}}=2, \mu=\frac{2\pi}{L}$.通过直接计算, 得到

$ \begin{equation} I_Nu(x, t)=\sum_{l=-\frac{N}{2}}^\frac{N}{2}\frac{1}{C_l}{\rm e}^{{\rm i}l\mu x}\frac{1}{N}\sum_{j=0}^{N-1}u_j{\rm e}^{-{\rm i}l\mu x_j}, \end{equation} $ (3.5)

其中

$ \begin{equation} u_j=I_Nu(x_j, t)=\sum_{l=-\frac{N}{2}}^\frac{N}{2}\frac{1}{C_l}{\rm e}^{{\rm i}l\mu x_j}\frac{1}{N}\sum_{k=0}^{N-1}u_k{\rm e}^{-{\rm i}l\mu x_k}. \end{equation} $ (3.6)

为了建立DGH方程(1.1)的Fourier拟谱多辛算法, 引入双线性定义

$ (u, v)_N=h\sum_{j=0}^{N-1}u(x_j, t)v(x_j, t). $

应用插值公式(3.3)得

$ (I_Nu(x_i, t), v(x_j, t))_N=(u(x_i, t), v(x_j, t))_N. $

Fourier拟谱方法的关键是计算$\frac{{\rm d}^k}{{\rm d}x^k}I_N(x, t)$$x=x_i$处的导数值, 对(3.5)式直接微分得到

$ \frac{{\rm d}^k}{{\rm d}x^k}I_N(x_i, t)=\frac{1}{N}\sum_{j=0}^{N-1}u_j \frac{{\rm d}^kg_j(x_i)}{{\rm d}x^k}=(D_ku)_i, $

其中$D_k$为谱微分矩阵

$ (D_1)_{i, j}=\left\{ \begin{array}{ll} \frac{1}{2}\mu(-1)^{i+j}\cot(\mu\frac{x_i-x_j}{2}), ~~ & i\neq j, \\ 0, & i= j, \end{array} \right.\\ (D_2)_{i, j}=\left\{ \begin{array}{ll} \frac{1}{2}\mu^2(-1)^{i+j+1}\frac{1}{\sin^2(\mu\frac{x_i-x_j}{2})}, ~~ & i\neq j, \\ -\mu^2\frac{N^2+2}{12}, & i= j. \end{array} \right. $

定义3.2  为了进行数值分析, 引入向量的Hadamard内积

$ u\circ v=(u_0v_0, u_1v_1, u_2v_2, \cdots u_{N-1}v_{N-1})^T, $

其中$u=(u_0, u_1, u_2, \cdots u_{N-1})^T$, $v=(v_0, v_1, v_2, \cdots v_{N-1})^T$.

定理3.1 为了进行数值分析, 引入向量的Hadamard内积微分形式为

$ D_N(u\circ v)=D_Nu\circ v+D_Nv\circ u, $

其中$u=(u_0, u_1, u_2, \cdots u_{N-1})^T$$v=(v_0, v_1, v_2, \cdots v_{N-1})^T$, $u_i=u(x_i, t), v_i=v(x_i, t), $$i=0, 1, 2 \cdots N-1.$

利用谱微分矩阵, 对$x$方向进行离散, 得到系统(2.5)的半离散系统

$ \begin{equation} \left\{ \begin{array}{ll} \frac{1}{2}\frac{{\rm d}\phi_i}{{\rm d}t}-\frac{\alpha^2}{2}\frac{\psi_i}{{\rm d}t}+\frac{\gamma}{2}(D_1\psi)_i-(D_1\Phi)_i=-w_i-2\omega u_i-\frac{3}{2}u_i^2-\frac{\alpha^2}{2}\psi_i^2, \\ -\frac{1}{2}\frac{{\rm d}u_i}{{\rm d}t}+(D_1w)_i=0, \\ -(D_1\phi)_i=-u_i, \\ (D_1u)_i=\psi_i, \\ \frac{\alpha^2}{2}\frac{{\rm d}u_i}{{\rm d}t}-\frac{\gamma}{2}(D_1u)_i=-\alpha^2u_i\psi_i+\Phi_i. \end{array} \right.\end{equation} $ (3.7)

定理3.2  系统(2.5)的多辛Fourier拟谱半离散格式(3.7)具有$N$个半离散多辛守恒律

$ \begin{equation} \frac{\rm d}{{\rm d}t}w_i+\sum_{j=0}^{N-1}(D_1)_{i, j}k_{i, j}=0, ~~i=1, 2, \cdots N-1 \end{equation} $ (3.8)

和整体辛守恒律

$ \frac{\rm d}{{\rm d}t}\sum_{i=0}^{N-1}w_i=0, $

其中

$ \left\{ \begin{array}{lcl} w_i=\frac{1}{2}({\rm d}z_i\wedge M {\rm d}z_i)=\frac{1}{2}{\rm d}u_i\wedge {\rm d}\phi_i+\frac{\alpha^2}{2} {\rm d}\psi_i\wedge {\rm d}u_i, \\ k_{i, j}=\frac{1}{2}({\rm d}z_i\wedge K {\rm d}z_j+{\rm d}z_j\wedge K {\rm d}z_i)=2{\rm d}\Phi_i\wedge {\rm d}u_j +\gamma {\rm d}u_i\wedge {\rm d}\psi_j+2{\rm d}\phi_i\wedge {\rm d}w_j. \end{array} \right. $

用隐式中点辛格式对半离散格式(3.7)时间方向进行离散, 得到系统(2.5)的全离散系统

$ \begin{equation} \left\{ \begin{array}{lcl} \frac{1}{2}D_t\phi_i^j-\frac{\alpha^2}{2}D_t\psi_i^j+\frac{\gamma}{2}(D_1A_t\psi^j)_i-(D_1A_t\Phi^j)_i\\ =-A_tw_i^{j}+2\omega A_tu_i^{j}-\frac{3}{2}(A_tu_i^{j})^2-\frac{\alpha^2}{2}(A_t\psi_i^{j})^2, \\ -\frac{1}{2}D_tu_i^j+(D_1A_tw^j)_i=0, \\ -(D_1A_t\phi^j)_i=A_tu_i^j, \\ (D_1A_tu^j)_i=-A_t\psi_i^j, \\ \frac{\alpha^2}{2}D_tu_i^j-\frac{\gamma}{2}(D_1A_tu^j)_i=-\alpha^2(A_t(u_i^j\psi_i^j))+A_t\Phi_i^j, \end{array} \right.\end{equation} $ (3.9)

其中

$ \begin{equation} D_tu^j=\frac{u^{j+1}-u^j}{\Delta t}, A_tu^j=\frac{u^{j+1}+u^j}{2}. \end{equation} $ (3.10)

定理3.3  系统(2.5)的多辛Fourier拟谱全离散格式(3.9)具有$N$个全离散多辛守恒律

$ \begin{equation} D_tw_i^j+\sum_{k=0}^{N-1}(D_1)_{i, k}A_tk_{i, k}^j=0, ~~i=1, 2, \cdots N-1 \end{equation} $ (3.11)

和整体辛守恒

$ \sum_{i=0}^{N-1}w_i^{j+1}=\sum_{i=0}^{N-1}w_i^{j}, $

其中

$ \left\{ \begin{array}{lcl} w_i^j=\frac{1}{2}({\rm d}z_i^j\wedge M {\rm d}z_i^j)=\frac{1}{2}{\rm d}u_i^j\wedge {\rm d}\phi_i^j+\frac{\alpha^2}{2} {\rm d}\psi_i^j\wedge {\rm d}u_i^j, \\ k_{i, k}^j=\frac{1}{2}({\rm d}z_i^j\wedge K {\rm d}z_k^j+{\rm d}z_k^j\wedge K {\rm d}z_i^j)=2{\rm d}\Phi_i^k\wedge {\rm d}u_j^k +\gamma {\rm d}u_i^k\wedge {\rm d}\psi_j^k+2{\rm d}\phi_i\wedge {\rm d}w_j^k. \end{array} \right. $

多辛Fourier拟谱方法不能满足离散局部能量守恒律, 为此我们考虑守恒律误差, 记

$ \begin{equation}(R_E)_i^j=D_tE_i^j+(D_1A_tF^j)_i, \end{equation} $ (3.12)

其中

$ \left\{ \begin{array}{lcl} E_i^j=S(Z_i^j)-\frac{1}{2}(Z_i^j)^TK(D_1Z^j)_i, \\ A_tF_i^j=\frac{1}{2}(A_tZ_i^j)^TK(D_tZ_i^j), \end{array} \right. $

$(R_E)_i^j$称为局部能量守恒律误差.

在周期边界条件下, 我们可以定义离散能量为

$ \begin{equation}\varepsilon^j=\Delta x\sum_{i=0}^{n-1}E_i^j, \end{equation} $ (3.13)

定理3.4  利用文献[19]类似的方法, 我们可以得到局部能量守恒律误差(3.12)和能量误差(3.13)满足

$ (R_E)_i^{j+\frac{1}{2}}=-\frac{1}{8\Delta t}(u_i^{j+1}-u_i^{j})^3-\frac{\alpha^2}{8\Delta t}(\psi_i^{j+1}-\psi_i^j)^2(u_i^{j+1}-u_i^j), $
$ \varepsilon^j=\varepsilon^0+\Delta x\sum_{j=0}^{n-1}\sum_{i=0}^{n-1}\bigg[-\frac{1}{8}(u_i^{j+1}-u_i^{j})^3-\frac{\alpha^2}{8}(\psi_i^{j+1}-\psi_i^j)^2(u_i^{j+1}-u_i^j)\bigg]. $

利用文献[19]类似的方法, 我们可得到如果$\Delta t$足够小, 则存在两个不依赖于$\Delta t$, $\Delta x$的两个常数$C_1, C_2, $且有

$ |(R_E)_i^{j+\frac{1}{2}}|\leq C_1\Delta t^2, |\varepsilon^{j+1}-\varepsilon^j|\leq C_2\Delta t^3. $
4 数值例子

为了说明多辛算法的诸多优点, 本文用多辛Fourier拟谱方法离散DGH方程(1.1)对应的多辛Hamilton方程(2.5), 消去辅助变量$\phi, w, \Phi, \psi$, 得到DGH方程(1.1)的多辛Fourier拟谱三层格式为

$ \begin{eqnarray*} &&A_tD_tU^j+2\omega A_tD_1U^j-\alpha^2 D_1^2D_tA_tU^j+\gamma A_tD_1^3U^j-\alpha^2(D_1^2A_t[U^{\frac{1}{2}}\cdot D_1U^{\frac{1}{2}}]\\ &&- \frac{1}{2}D_1A_t(D_1U^{\frac{1}{2}})^2)+\frac{3}{2}D_1A_t(U^{\frac{1}{2}})^2=0, \end{eqnarray*} $

在利用上面三层格式进行数值模拟时, 可用下面二层格式计算上面三层格式第二层上的初值

$ \begin{eqnarray*} &&D_tU^j+2\omega D_1U^j-\alpha^2 D_1^2D_tU^j+\gamma D_1^3U^j-\alpha^2(D_1^2[U^{\frac{1}{2}}\cdot D_1U^{\frac{1}{2}}]\\ &&- \frac{1}{2}D_1(D_1U^{\frac{1}{2}})^2)+\frac{3}{2}D_1(U^{\frac{1}{2}})^2=0, \end{eqnarray*} $

其中, $U=[u_0, u_1, u_2, \cdots u_{N-1}]^T, U^{\frac{1}{2}}=\frac{U^j+U^{j+1}}{2}.$

4.1 孤立波解

取参数$\alpha=1, \omega=1, \gamma=1$, 考虑下面DGH方程(1.1)的孤立波的初边值问题(周期边界条件)

$ \begin{equation} \left\{ \begin{array}{lcl} u_t+2 u_x+3uu_x+ u_{xxx}-u_{xxt}-(uu_{xxx}+2u_xu_{xx})=0, \\ u(x, 0)=-\frac{5}{2}+\cos(x), \\ u(0, t)=u(2\pi, t)=-\frac{5}{2}+\cos(t). \end{array} \right.\end{equation} $ (4.1)

由文献[20], 我们可以得到问题(4.1)有孤立波解

$ u(x, t)=-\frac{5}{2}+\cos(x+t). $

考虑时间步长$\Delta t=0.001, $空间步长$\Delta x=0.01$.利用多辛Fourier拟谱方法对DGH方程(1.1)的初边值问题(周期边界条件) (4.1)在区间$x\in[0, 2\pi], t\in [0, 100]$内模拟时, 得到DGH方程(1.1)的初边值问题(4.1)的孤立波的数值解.图 1给出了DGH方程的初边值问题(4.1)随时间的演化图, 表 1给出了精确解、多辛Fourier拟谱格式、有限差分法的比较.

图 1 DGH方程解的孤立波图

表 1 精确解、多辛Fourie拟谱格式、有限差分法的比较
4.2 尖峰孤立波解

取参数$\alpha=1, \omega=1, \gamma=1$, 考虑下面DGH方程(1.1)的尖峰孤立波的初边值问题(周期边界条件)

$ \begin{equation} \left\{ \begin{array}{lcl} u_t+2 u_x+3uu_x+ u_{xxx}-u_{xxt}-(uu_{xxx}+2u_xu_{xx})=0, \\ u(x, 0)=-\frac{3}{2}+\exp(-|x|), \\ u(-10, t)=u(10, t)=-\frac{3}{2}+\exp(-|10+t|). \end{array} \right.\end{equation} $ (4.2)

考虑时间步长$\Delta t=0.001, $空间步长$\Delta x=0.01$.利用多辛Fourier拟谱方法对DGH方程(1.1)的初边值问题(周期边界条件) (4.2)在区间$x\in[-10, 10], t\in[0, 100]$内模拟时, 得到DGH方程(4.2)的孤立波的数值解.图 2给出了DGH方程的初边值问题(4.2)随时间的演化图, 图 3给出了离散能量守恒律误差.

图 2 DGH方程解的尖峰孤立波图

图 3 局部能量守恒律误差
参考文献
[1] Camassa R, Holm D. An integrable shallow water equation with peaked solitions. Physical Review Letters, 1993, 11: 1661–1664.
[2] Fuchssteiner B, Fokas A. Symplectic strutures their Backlund transformation and herditary symmetries. Physica D, 1981, 4: 47–66. DOI:10.1016/0167-2789(81)90004-X
[3] Camassa R, Holm D, Hyman J. A new integrable shallow water equation. Advances in Applied Mechanics, 1994, 31: 1–33. DOI:10.1016/S0065-2156(08)70254-0
[4] Fisher M, Sehiff J. The Camassa Holme equations:conserved quantities and the initial value Problem. Physics Letters A, 1999, 259(3): 371–376.
[5] Clarkson P, Mansfield E, priestley T. Symetries of a class of nonlinear third-order partial differential equations. Mathematical and Computer Modelling, 1997, 25(8): 195–212.
[6] CooPer F, ShePard H. Solitons in the Camassa-Holm shallow watere equation. Physics Letters A, 1994, 194(4): 246–250. DOI:10.1016/0375-9601(94)91246-7
[7] Constantin A, Eseher J. Global existence and blow-up for a shallow watere equation. Annali Della Scuola Normale Superiore Di Pisa Classe Di Scienze, 1998, 26: 303–328.
[8] Tian L X, Shi Q. Boundary control of viscous Dullin-Gottwald-Holm equation. International Journal of Nonlinear Science, 2007, 4(1): 67–75.
[9] Feng K, Qin M Z. The symplectic methods for the computation of Hamiltonian equations//Lecture Notes in Math, 1297. Berlin:Springer, 1987:1-37
[10] Liu T T, Qin M Z. Multi-symplectic geometry and multi-symplectic preissmann scheme for the KP equation. Journal of Mathematical Physics, 2002, 43(8): 4060–4077. DOI:10.1063/1.1487444
[11] Tian Y M, Qin M Z, Zhang Y M, Ma T. The multi-symplectic numerical method for Gross-Pitaevskii equation. Computer Physics Communications, 2008, 176(6): 449–458.
[12] Wang Y S, Wang B, Qin M Z. Concatenating construction of multi-symplectic scheme for 2+1 dimensional sine-Gordon equation. Science in China (Series A), 2004, 47(1): 18–30. DOI:10.1360/02yc0163
[13] Escher J, Lechtenfeld O, Yin Z. Well-posedness and blow-up phenomena for the 2-component Camassa-Holm equation. Discrete and Continuous Dynamical Systems-Series A, 2007, 19: 493–513. DOI:10.3934/dcdsa
[14] Kong L H, Liu R X, Zheng X H. A survey on symplectic and multi-symplectic algorithms. Applied Mathematics and Computation, 2007, 186: 670–684. DOI:10.1016/j.amc.2006.08.012
[15] Leimkuhler B, Reich S. Simulating Hamiltonian Dynamics. Cambridge: Cambridge University Press, 2004.
[16] Islas A, Schober C. Backward error analysis for multisymplectic discretization of Hamiltonian PDEs. Mathematics and Computers in Simulation, 2005, 69: 290–303. DOI:10.1016/j.matcom.2005.01.006
[17] Moore B, Reich S. Backward error analysis for multi-symplectic integration methods. Numerische Mathematik, 2003, 95: 625–652. DOI:10.1007/s00211-003-0458-9
[18] Wang J. Multisymplectic forier pseudospectral method for the nonlinear schrodinger equation with wave operator. Journal of Computational Mathematics, 2007, 25(1): 31–48.
[19] 殷久利, 田立新. 一类非线性色散方程中的新型奇异孤立波. 物理学报 , 2009, 58 (6): 3632–3636.
Yin J L, Tian L X. New exotic solitary waves in one type of nonlinear dispersive equations. Acta Physica Sinica, 2009, 58(6): 3632–3636.