Numerical simulation of the soliton solutions for a complex modified Korteweg-de Vries equation by a
本站小编 Free考研考试/2022-01-02
Tao Xu,1,2, Guowei Zhang2, Liqun Wang2, Xiangmin Xu2, Min Li31State Key Laboratory of Heavy Oil Processing, China University of Petroleum, Beijing 102249, China 2College of Science, China University of Petroleum, Beijing 102249, China 3Department of Mathematics and Physics, North China Electric Power University, Beijing 102206, China
First author contact:*Author to whom should be corresponding address. Received:2020-08-22Revised:2020-10-13Accepted:2020-12-6Online:2021-01-28
Abstract In this paper, a Crank-Nicolson-type finite difference method is proposed for computing the soliton solutions of a complex modified Korteweg-de Vries (MKdV) equation (which is equivalent to the Sasa-Satsuma equation) with the vanishing boundary condition. It is proved that such a numerical scheme has the second-order accuracy both in space and time, and conserves the mass in the discrete level. Meanwhile, the resulting scheme is shown to be unconditionally stable via the von Nuemann analysis. In addition, an iterative method and the Thomas algorithm are used together to enhance the computational efficiency. In numerical experiments, this method is used to simulate the single-soliton propagation and two-soliton collisions in the complex MKdV equation. The numerical accuracy, mass conservation and linear stability are tested to assess the scheme’s performance. Keywords：complex modified Korteweg-de Vries equation;finite difference method;soliton solutions
PDF (848KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Tao Xu, Guowei Zhang, Liqun Wang, Xiangmin Xu, Min Li. Numerical simulation of the soliton solutions for a complex modified Korteweg-de Vries equation by a finite difference method. Communications in Theoretical Physics, 2021, 73(2): 025005- doi:10.1088/1572-9494/abd0e5
1. Introduction
Up to now, a number of integrable nonlinear evolution equations have been derived to describe various nonlinear wave phenomena, like the solitons, breathers and rogue waves [1-6]. As a universe integrable model, the nonlinear Schrödinger equation (NLSE) appears in optics [7], fluids [8], plasmas [9], quantum field theory [10], Bose-Einstein condensates [11], etc. Different from the Korteweg-de Vries equation, the NLSE field is complex and thus it admits the envelope solitons which are formed on the balance between the group velocity dispersion and cubic nonlinearity [7]. By including some higher-order dispersion and nonlinear terms, the NLSE has also been generalized to describe the propagation of ultrashort optical pulses [12, 13] and internal waves [14]. When the third-order dispersion is considered, two important integrable versions are the Hirota equation [15] and Sasa-Satsuma (SS) equation [16].
Through the Galilean-phase transformations, the Hirota equation and SS equation can be respectively written in the simple form [16-18]$\begin{eqnarray}{\partial }_{\xi }p+{\partial }_{\eta }^{3}p+6\,| p{| }^{2}{\partial }_{\eta }p=0,\end{eqnarray}$$\begin{eqnarray}{\partial }_{\xi }p+{\partial }_{\eta }^{3}p+6\,| p{| }^{2}{\partial }_{\eta }p+3\,p\,{\partial }_{\eta }(| p{| }^{2})=0,\end{eqnarray}$which are usually referred to be two integrable versions of the complex modified Korteweg-de Vries (MKdV) equation, where ξ and η are respectively the normalized time and space variables, p is the complex-valued function of ξ and η. Both equations (1) and(2) are solvable by the inverse scattering transform (IST) [16, 19-21], and their integrable properties have been detailed from different aspects, including the Lax pair [16, 22], Painlevé property [22], conservation laws [23, 24], Hamiltonian structures [24, 25], Darboux transformation (DT) [18, 26-28], bilinear representation [15, 29, 30], and squared eigenfunctions [31]. Meanwhile, a lot of efforts have been devoted to the exact analytical localized-wave solutions and their dynamical evolution behaviors (e.g. see [18-22, 26-30, 32-43]).
Similarly to the NLSE, equation (1) has the bright soliton with the vanishing background and exhibits the standard elastic soliton collisions [15, 17, 18, 26, 32]. Quite differently, because of the additional nonlinear term, equation (2) admits the single-hump (SH) and double-hump (DH) solitons as well as the multi-hump (MH) soliton having periodical oscillating behavior with the vanishing background [16, 20]. Moreover, it has been found that equation (2) can exhibit the unusual shape-changing soliton collisions, that is, the MH soliton changes into an SH or DH one upon a collision, and vice versa [27]. From the viewpoint of physical applications, the DH and MH solitons may be used to design new error preventable line-coding scheme which is robust against the impairments arising from intrachannel collisions [44], and the shape-changing soliton collisions may have applications in all-optical information processing, optical switching, and routing of optical signals [45, 46]. In addition, we point that the symmetric and asymmetric DH solitons also exist in other integrable systems [47-49].
However, we note that the exact analytical solutions of equation (2) may become weak in practical physical systems. This is because one must consider the possibility that the initial waves are not matched to the exact analytical solutions, or that there is a slight violation from the fixed ratio of the last three terms in equation (2) [7, 50]. Therefore, it is necessary and physically important to develop some numerical methods to simulate the soliton solutions of equation (2). In the past decades, a non-integrable complex MKdV equation was studied by the split-step Fourier spectral method [51], mesh-free collocation method [52], Petrov-Galerkin method [53], and multisymplectic method [54, 55]. But for the integrable cases, only equation (1) was numerically solved by the integrable IST scheme which does not suffer from any instabilities [56], whereas equation (2) has received little attention in the numerical aspect.
In this paper, we will propose a conservative finite difference method to equation (2) with the vanishing boundary condition, and use it to numerically simulate the propagation and collisions for three types of solitons. The rest of this paper is organized as follows. In section 2, we review the exact analytical soliton solutions as well as their propagation and collision properties. In section 3, we present a Crank-Nicolson (CN)-type finite difference scheme for equation (2), and make numerical analysis for the accuracy, mass conservation and linear stability. In section 4, we perform numerical experiments for the one- and two-soliton solutions, and check the numerical accuracy, mass conservation and linear stability against an initial perturbation. In section 5, we address the conclusions and discussions of this paper.
2. Exact analytical soliton solutions
By using the N-fold DT and choosing the seed solution p=0, we obtain the N-soliton solutions of equation (2) in the determinant form [27]$\begin{eqnarray}{p}_{N}=-2\,\displaystyle \frac{{{\boldsymbol{\chi }}}_{2N+1,2N-1,2N}}{{{\boldsymbol{\chi }}}_{2N,2N,2N}},\end{eqnarray}$with ${{\boldsymbol{\chi }}}_{L,{L}^{{\rm{{\prime} }}},2N}$ ($L\mbox{'}=4N-L$) defined as the following three-component Wronskian:$\begin{eqnarray}\begin{array}{c}{{\boldsymbol{\chi }}}_{L,L\text{\unicode{x02019}},2N}=\left|\begin{array}{ccc}{{\boldsymbol{F}}}_{N\times L} & -{{\boldsymbol{G}}}_{N\times L\text{\unicode{x02019}}} & -{{\boldsymbol{H}}}_{N\times N}\\ {{\boldsymbol{F}}}_{N\times L}^{\ast } & -{{\boldsymbol{H}}}_{N\times L\text{\unicode{x02019}}}^{\ast } & -{{\boldsymbol{G}}}_{N\times N}^{\ast }\\ {{\boldsymbol{G}}}_{N\times L}^{\ast } & {{\boldsymbol{F}}}_{N\times L\text{\unicode{x02019}}}^{\ast } & {\boldsymbol{O}}\\ {{\boldsymbol{H}}}_{N\times L} & {{\boldsymbol{F}}}_{N\times L\text{\unicode{x02019}}} & {\boldsymbol{O}}\\ {{\boldsymbol{H}}}_{N\times L}^{\ast } & {\boldsymbol{O}} & {{\boldsymbol{F}}}_{N\times N}^{\ast }\\ {{\boldsymbol{G}}}_{N\times L} & {\boldsymbol{O}} & {{\boldsymbol{F}}}_{N\times N}\end{array}\right|\,(L=2N,2N+1),\end{array}\,\end{eqnarray}$in which the block matrices FN×L, GN×L and HN×L are given by$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{F}}}_{N\times L} & = & {\boldsymbol{A}}{{\boldsymbol{\Theta }}}_{+}{{\boldsymbol{\Lambda }}}_{N\times L}^{+},\quad {{\boldsymbol{G}}}_{N\times L}={\boldsymbol{B}}{{\boldsymbol{\Theta }}}_{-}{{\boldsymbol{\Lambda }}}_{N\times L}^{-},\\ {{\boldsymbol{H}}}_{N\times L} & = & {\boldsymbol{C}}{{\boldsymbol{\Theta }}}_{-}{{\boldsymbol{\Lambda }}}_{N\times L}^{-},\end{array}\end{eqnarray}$with$\begin{eqnarray*}\begin{array}{l}{\boldsymbol{A}}={\boldsymbol{I}},\quad {\boldsymbol{B}}=\mathrm{diag}({\beta }_{1},\ldots ,{\beta }_{N}),\quad {\boldsymbol{C}}=\mathrm{diag}({\gamma }_{1},\ldots ,{\gamma }_{N}),\\ {{\boldsymbol{\Theta }}}_{+}=\mathrm{diag}\left({{\rm{e}}}^{{\theta }_{1}},\ldots ,{{\rm{e}}}^{{\theta }_{N}}\right),\quad {{\boldsymbol{\Theta }}}_{-}=\mathrm{diag}\left({{\rm{e}}}^{-{\theta }_{1}},\ldots ,{{\rm{e}}}^{-{\theta }_{N}}\right),\\ {{\boldsymbol{\Lambda }}}_{N\times L}^{+}={\left({\lambda }_{n}^{m-1}\right)}_{N\times L},\\ {{\boldsymbol{\Lambda }}}_{N\times L}^{-}={\left[{\left(-{\lambda }_{n}\right)}^{m-1}\right]}_{N\times L},\quad {\theta }_{k}={\lambda }_{k}\eta -4\,{\lambda }_{k}^{3}\xi ,\end{array}\end{eqnarray*}$where ${\{{\beta }_{k},{\gamma }_{k},{\lambda }_{k}\}}_{k=1}^{N}$ are 3 N parameters in ${\mathbb{C}}$.
Taking N=1 in solution(3), the one-soliton solution can be written as$\begin{eqnarray}\begin{array}{l}{p}_{1}\,=\,\displaystyle \frac{4{\lambda }_{1{\rm{R}}}\left[{\rm{i}}{\lambda }_{1{\rm{I}}}\left({\bar{\beta }}_{1}{\bar{\lambda }}_{1}{{\rm{e}}}^{2{\theta }_{1}}-{\gamma }_{1}{\lambda }_{1}{{\rm{e}}}^{2{\bar{\theta }}_{1}}\right)+{\delta }_{1}\left({\bar{\beta }}_{1}{\lambda }_{1}{{\rm{e}}}^{-2{\bar{\theta }}_{1}}-{\gamma }_{1}{\bar{\lambda }}_{1}{{\rm{e}}}^{-2{\theta }_{1}}\right)\right]}{{\lambda }_{1{\rm{I}}}^{2}{{\rm{e}}}^{4{\theta }_{1{\rm{R}}}}+| {\delta }_{1}{| }^{2}{{\rm{e}}}^{-4{\theta }_{1{\rm{R}}}}-2{\lambda }_{1{\rm{R}}}^{2}\left({\beta }_{1}{\gamma }_{1}{{\rm{e}}}^{-4{\rm{i}}{\theta }_{1{\rm{I}}}}+{\bar{\beta }}_{1}{\bar{\gamma }}_{1}{{\rm{e}}}^{4{\rm{i}}{\theta }_{1{\rm{I}}}}\right)+2| {\lambda }_{1}{| }^{2}\left(| {\beta }_{1}{| }^{2}+| {\gamma }_{1}{| }^{2}\right)},\end{array}\end{eqnarray}$where ${\theta }_{1}={\lambda }_{1}\eta -4\,{\lambda }_{1}^{3}\xi $, ${\delta }_{1}=| {\beta }_{1}{| }^{2}{\bar{\lambda }}_{1}-| {\gamma }_{1}{| }^{2}{\lambda }_{1}$, i is the imaginary unit, the subscripts R and I denote the real and imaginary parts of the involved parameters, respectively. Depending on the parameters β1, γ1 and λ1, solution(6) can display three different soliton profiles [16, 20, 27]: (i) For β1 or γ1=0 and $| {\lambda }_{1{\rm{R}}}| \leqslant \sqrt{3}| {\lambda }_{1{\rm{I}}}| $, ∣p1∣ has only one maximum and hence exhibits the usual SH soliton (see figure 1(a)); (ii) With β1 or γ1=0 and $| {\lambda }_{1{\rm{R}}}| \gt \sqrt{3}| {\lambda }_{1{\rm{I}}}| $, ∣p1∣ has two identical maxima and thus represents the DH soliton (see figure 1(b)); (iii) If β1 and γ1≠0, ∣p1∣ behaves as a localized multi-hump entity which oscillates periodically both in ξ and η (see figure 1(c)). For all the three cases, the soliton center trajectories lie in the line ${ \mathcal C }:\,2\,{\theta }_{1{\rm{R}}}=\mathrm{ln}\sqrt{\tfrac{\left|{\delta }_{1}\right|}{\left|{\lambda }_{1{\rm{I}}}\right|}}$, so the soliton propagates with the velocity $v=4\left({\lambda }_{1{\rm{R}}}^{2}-3\,{\lambda }_{1{\rm{I}}}^{2}\right)$. In addition, it can be obtained that the soliton mass is dependent on λ1R, that is,$\begin{eqnarray}M={\int }_{-\infty }^{+\infty }| p{| }^{2}{\rm{d}}\eta =4\,| {\lambda }_{1{\rm{R}}}| .\end{eqnarray}$
Figure 1.
New window|Download| PPT slide Figure 1.Three types of one-soliton solutions of equation (2): (a) Aan SH soliton with β1=0, γ1=1 and ${\lambda }_{1}=\tfrac{3}{5}+\tfrac{2}{5}{\rm{i}};$ (b) a DH soliton with β1=1, γ1=0 and ${\lambda }_{1}=\tfrac{3}{5}+\tfrac{1}{10}{\rm{i}};$ (c) an MH soliton with β1=0.5, γ1=1 and ${\lambda }_{1}=\tfrac{1}{2}+\tfrac{1}{5}{\rm{i}}$.
With N=2, solution(3) describes the elastic two-soliton collisions in the sense that the soliton velocities and masses are conserved before and after collision. According to the profiles of asymptotic solitons, the soliton collisions can be classified into three different cases [27]: (i) If (β1,β2), (γ1,γ2), (β1,γ2) or (β2,γ1)=0, the two interacting solitons are either the SH or DH ones. (ii) If none of βk and γk (k=1,2) is 0, the collisions occur between two MH solitons with internal oscillation. (iii) If only one of βk,γk (k=1,2) is 0, one soliton is always of the MH type, but the other one may experience the change of shape upon a collision, that is, the SH or DH soliton changes into an oscillating MH one, or the MH soliton into an SH or DH one. It should be noted that the shape-changing soliton collision is a peculiar property of equation (2), which has never been reported in other (1+1)-dimensional single-component integrable equations.
3. CN-type finite difference scheme and numerical analysis
By choosing an interval Ω=[L1,L2] properly and large enough, equation (2) with the vanishing boundary condition ${\mathrm{lim}}_{| \eta | \to \infty }p=0$ can be approximated by$\begin{eqnarray}\begin{array}{l}{\partial }_{\xi }p+{\partial }_{\eta }^{3}p+6\,| p{| }^{2}{\partial }_{\eta }p+3\,p\,{\partial }_{\eta }(| p{| }^{2})=0,\\ {L}_{1}\lt \eta \lt {L}_{2},\,\xi \gt 0,\end{array}\end{eqnarray}$$\begin{eqnarray}p({L}_{1},\xi )=p({L}_{2},\xi )=0,\qquad \xi \geqslant 0,\end{eqnarray}$$\begin{eqnarray}p(\eta ,\xi =0)={p}_{0}(\eta ),\qquad {L}_{1}\leqslant \eta \leqslant {L}_{2}.\end{eqnarray}$In order to obtain the numerical method for discretizing the problem(8)-(10), we choose the time step τ>0 and mesh size h=(L2−L1)/N (with N as an even positive integer), and denote the grid points and time steps as$\begin{eqnarray}{\eta }_{j}:= {L}_{1}+{jh},\,j=0,\cdots ,\,N;\quad {\xi }_{n}:= n\tau ,\,n=0,1,\cdots .\end{eqnarray}$Meanwhile, based on the Taylor series expansion, we give the following three useful formulas:$\begin{eqnarray}\begin{array}{l}p\left({\eta }_{j\pm 1},{\xi }_{n+\tfrac{1}{2}}\right)=p({\eta }_{j},{\xi }_{n+\tfrac{1}{2}})\\ \quad +\,\displaystyle \sum _{k=1}^{a}\displaystyle \frac{{\left(\pm h\right)}^{k}}{k!}{\left[\displaystyle \frac{{\partial }^{k}p}{\partial {\eta }^{k}}\right]}_{j}^{n+\tfrac{1}{2}}+O\left({h}^{a+1}\right),\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}p\left({\eta }_{j},{\xi }_{n+\tfrac{1}{2}\pm \tfrac{1}{2}}\right)=p({\eta }_{j},{\xi }_{n+\tfrac{1}{2}})\\ \quad +\,\displaystyle \sum _{k=1}^{b}\displaystyle \frac{{\left(\pm \tau /2\right)}^{k}}{k!}{\left[\displaystyle \frac{{\partial }^{k}p}{\partial {\xi }^{k}}\right]}_{j}^{n+\tfrac{1}{2}}+O\left({\tau }^{b+1}\right),\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}{\left|p\left({\eta }_{j\pm 1},{\xi }_{n+\tfrac{1}{2}}\right)\right|}^{2}={\left|p\left({\eta }_{j},{\xi }_{n+\tfrac{1}{2}}\right)\right|}^{2}\\ \quad +\,\displaystyle \sum _{k=1}^{c}\displaystyle \frac{{\left(\pm h\right)}^{k}}{k!}{\left[\displaystyle \frac{{\partial }^{k}| p{| }^{2}}{\partial {\eta }^{k}}\right]}_{j}^{n+\tfrac{1}{2}}+O\left({h}^{c+1}\right),\end{array}\end{eqnarray}$where a,b,c are three positive integers. For simplicity, we introduce the finite difference operators$\begin{eqnarray}\begin{array}{rcl}{\delta }_{\xi }^{+}{p}_{j}^{n} & = & \displaystyle \frac{{p}_{j}^{n+1}-{p}_{j}^{n}}{\tau },\quad {\delta }_{\eta }{p}_{j}^{n}=\displaystyle \frac{{p}_{j+1}^{n}-{p}_{j-1}^{n}}{2h},\\ {\delta }_{\eta }^{3}{p}_{j}^{n} & = & \displaystyle \frac{{p}_{j+2}^{n}-2{p}_{j+1}^{n}+2{p}_{j-1}^{n}-{p}_{j-2}^{n}}{2{h}^{3}},\end{array}\end{eqnarray}$with ${p}_{j}^{n}$ being the numerical approximation of p at the point (ηj,ξn) (j=0,⋯, N;n=0,1,⋯).
With proper integers for a, b and c, equations (12)-(14) give rise to$\begin{eqnarray}\displaystyle \frac{p\left({\eta }_{j},{\xi }_{n+1}\right)-p\left({\eta }_{j},{\xi }_{n}\right)}{\tau }={\left[\displaystyle \frac{\partial p}{\partial \xi }\right]}_{j}^{n+\tfrac{1}{2}}+O\left({\tau }^{2}\right),\end{eqnarray}$$\begin{eqnarray}\displaystyle \frac{p\left({\eta }_{j+1},{\xi }_{n+\tfrac{1}{2}}\right)-p\left({\eta }_{j-1},{\xi }_{n+\tfrac{1}{2}}\right)}{2h}={\left[\displaystyle \frac{\partial p}{\partial \eta }\right]}_{j}^{n+\tfrac{1}{2}}+O\left({h}^{2}\right),\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{{\left|p\left({\eta }_{j+1},{\xi }_{n+\tfrac{1}{2}}\right)\right|}^{2}-{\left|p\left({\eta }_{j-1},{\xi }_{n+\tfrac{1}{2}}\right)\right|}^{2}}{2h}\\ \quad =\,{\left[\displaystyle \frac{\partial | p{| }^{2}}{\partial \eta }\right]}_{j}^{n+\tfrac{1}{2}}+O\left({h}^{2}\right),\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}\displaystyle \frac{p\left({\eta }_{j+2},{\xi }_{n+\tfrac{1}{2}}\right)-2p\left({\eta }_{j+1},{\xi }_{n+\tfrac{1}{2}}\right)+2p\left({\eta }_{j-1},{\xi }_{n+\tfrac{1}{2}}\right)-p\left({\eta }_{j-2},{\xi }_{n+\tfrac{1}{2}}\right)}{2{h}^{3}}=\,{\left[\displaystyle \frac{{\partial }^{3}p}{\partial {\eta }^{3}}\right]}_{j}^{n+\tfrac{1}{2}}+O\left({h}^{2}\right),\end{array}\end{eqnarray}$where $[\,\cdot \,{]}_{j}^{n+\tfrac{1}{2}}$ denotes the exact value of the involved function at the point $({\eta }_{j},{\xi }_{n+\tfrac{1}{2}})$ (j=0,⋯, N;n=0,1,⋯).
Based on the above equations, we obtain the following finite difference scheme of equation (8) as follows:$\begin{eqnarray}\begin{array}{l}{\delta }_{\xi }^{+}{p}_{j}^{n}+{\delta }_{\eta }^{3}{p}_{j}^{n+\tfrac{1}{2}}+3{p}_{j}^{n+\tfrac{1}{2}}\left({\bar{p}}_{j+1}^{n+\tfrac{1}{2}}+{\bar{p}}_{j-1}^{n+\tfrac{1}{2}}\right){\delta }_{\eta }{p}_{j}^{n+\tfrac{1}{2}}\\ \quad +\,3{p}_{j}^{n+\tfrac{1}{2}}{\delta }_{\eta }\left({\left|{p}_{j}^{n+\tfrac{1}{2}}\right|}^{2}\right)=0,\end{array}\end{eqnarray}$where the CN scheme is used for the temporal derivative, the central finite differences is used for the spatial derivatives, and ${p}_{j}^{n+\tfrac{1}{2}}=\tfrac{1}{2}({p}_{j}^{n+1}+{p}_{j}^{n})$. Meanwhile, the boundary condition(9) is discretized as$\begin{eqnarray}{p}_{0}^{n}={p}_{N}^{n}=0,\quad n\geqslant 0,\end{eqnarray}$and the initial condition(10) is discretized as$\begin{eqnarray}{p}_{j}^{0}={p}_{0}({\eta }_{j}),\quad j=0,1,\cdots ,\,N.\end{eqnarray}$
In addition to the boundary condition(21), we should also impose the numerical derivatives pη=0 at two end points of the interval Ω. Hence, we introduce the artificial grids j=−1 and j=N+1, and take the following conditions:$\begin{eqnarray}{p}_{-1}^{n}={p}_{N+1}^{n}=0,\quad n\geqslant 0.\end{eqnarray}$
The local truncation error of the difference equation (20) to equation (2) is of O(τ2+h2).
The result can be directly obtained by using equations (16)-(19). □
For the CN-type finite difference scheme(20)-(22), we can get the mass conservation in the discrete level as follows.
If we define the discrete mass as$\begin{eqnarray}{M}_{1}^{n}:= h\displaystyle \sum _{j=1}^{N-1}| {p}_{j}^{n}{| }^{2},\end{eqnarray}$then the CN-type finite difference method(20)-(22) conserves the mass in the discrete level, i.e.$\begin{eqnarray}{M}_{1}^{n}\equiv {M}_{1}^{0},\quad n\geqslant 0.\end{eqnarray}$
Multiplying equation (20) by ${\bar{p}}_{j}^{n+\tfrac{1}{2}}=\tfrac{1}{2}\left({\bar{p}}_{j}^{n+1}+{\bar{p}}_{j}^{n}\right)$ and adding the resulted equation with its complex conjugate, we obtain$\begin{eqnarray}\begin{array}{rcl}{\delta }_{\xi }^{+}| {p}_{j}^{n}{| }^{2} & = & -\left({\bar{p}}_{j}^{n+\tfrac{1}{2}}{\delta }_{\eta }^{3}{p}_{j}^{n+\tfrac{1}{2}}+{p}_{j}^{n+\tfrac{1}{2}}{\delta }_{\eta }^{3}{\bar{p}}_{j}^{n+\tfrac{1}{2}}\right)\\ & & -6{\left|{p}_{j}^{n+\tfrac{1}{2}}\right|}^{2}{\delta }_{\eta }{\left|{p}_{j}^{n+\tfrac{1}{2}}\right|}^{2}.\end{array}\end{eqnarray}$Taking the summation of equation (26) from j=1 to N−1 and noticing(24), we have$\begin{eqnarray}\begin{array}{l}{M}_{1}^{n+1}-{M}_{1}^{n}=h\displaystyle \sum _{j=1}^{N-1}\left(| {p}_{j}^{n+1}{| }^{2}-| {p}_{j}^{n}{| }^{2}\right)\\ =\,-h\tau \displaystyle \sum _{j=1}^{N-1}\left({\bar{p}}_{j}^{n+\tfrac{1}{2}}{\delta }_{\eta }^{3}{p}_{j}^{n+\tfrac{1}{2}}+{p}_{j}^{n+\tfrac{1}{2}}{\delta }_{\eta }^{3}{\bar{p}}_{j}^{n+\tfrac{1}{2}}\right)\\ \ \ -\,6\tau \displaystyle \sum _{j=1}^{N-1}{\left|{p}_{j}^{n+\tfrac{1}{2}}\right|}^{2}\left({\left|{p}_{j+1}^{n+\tfrac{1}{2}}\right|}^{2}-{\left|{p}_{j-1}^{n+\tfrac{1}{2}}\right|}^{2}\right)\\ =\,-\displaystyle \frac{\tau }{2{h}^{2}}\left(2{\bar{p}}_{1}^{n+\tfrac{1}{2}}{p}_{0}^{n+\tfrac{1}{2}}+2{\bar{p}}_{0}^{n+\tfrac{1}{2}}{p}_{1}^{n+\tfrac{1}{2}}-{\bar{p}}_{-1}^{n+\tfrac{1}{2}}{p}_{1}^{n+\tfrac{1}{2}}-{\bar{p}}_{1}^{n+\tfrac{1}{2}}{p}_{-1}^{n+\tfrac{1}{2}}\right.\\ \ \ -\,{\bar{p}}_{2}^{n+\tfrac{1}{2}}{p}_{0}^{n+\tfrac{1}{2}}-{\bar{p}}_{0}^{n+\tfrac{1}{2}}{p}_{2}^{n+\tfrac{1}{2}}+{\bar{p}}_{N}^{n+\tfrac{1}{2}}{p}_{N-2}^{n+\tfrac{1}{2}}+{\bar{p}}_{N-2}^{n+\tfrac{1}{2}}{p}_{N}^{n+\tfrac{1}{2}}\\ \ \ \left.+\,{\bar{p}}_{N+1}^{n+\tfrac{1}{2}}{p}_{N-1}^{n+\tfrac{1}{2}}+{\bar{p}}_{N-1}^{n+\tfrac{1}{2}}{p}_{N+1}^{n+\tfrac{1}{2}}-2{\bar{p}}_{N}^{n+\tfrac{1}{2}}{p}_{N-1}^{n+\tfrac{1}{2}}-2{\bar{p}}_{N-1}^{n+\tfrac{1}{2}}{p}_{N}^{n+\tfrac{1}{2}}\right)\\ \ \ -\,6\tau \left({\left|{p}_{N-1}^{n+\tfrac{1}{2}}\right|}^{2}{\left|{p}_{N}^{n+\tfrac{1}{2}}\right|}^{2}-{\left|{p}_{1}^{n+\tfrac{1}{2}}\right|}^{2}{\left|{p}_{0}^{n+\tfrac{1}{2}}\right|}^{2}\right),\qquad n\geqslant 0.\end{array}\end{eqnarray}$This, together with the boundary conditions(21) and(23), implies that ${M}_{1}^{n+1}-{M}_{1}^{n}=0$ for all n≥0. □
Next, we examine the linear stability of the CN-type finite difference scheme(20)-(22) by the von Neumann method. Following the way in [57], we freeze some variables in the nonlinear terms of equation (20), that is,$\begin{eqnarray}\begin{array}{rcl}{{\mathfrak{M}}}_{1} & = & {p}_{{j}_{1}}^{{n}_{1}+\tfrac{1}{2}}\left({\bar{p}}_{{j}_{1}+1}^{{n}_{1}+\tfrac{1}{2}}+{\bar{p}}_{{j}_{1}-1}^{{n}_{1}+\tfrac{1}{2}}\right),\\ {{\mathfrak{M}}}_{2} & = & {p}_{{j}_{2}}^{{n}_{2}+\tfrac{1}{2}}\left(\left|{p}_{{j}_{2}+1}^{{n}_{2}+\tfrac{1}{2}}\right|+\left|{p}_{{j}_{2}-1}^{{n}_{2}+\tfrac{1}{2}}\right|\right),\end{array}\end{eqnarray}$where (jk,nk) (k=1,2) are the indices such that$\begin{eqnarray}\begin{array}{l}\left|{p}_{{j}_{1}}^{{n}_{1}+\tfrac{1}{2}}\left({\bar{p}}_{{j}_{1}+1}^{{n}_{1}+\tfrac{1}{2}}+{\bar{p}}_{{j}_{1}-1}^{{n}_{1}+\tfrac{1}{2}}\right)\right|\\ =\,\mathop{\max }\limits_{0\leqslant j\leqslant N;n\geqslant 0}\left|{p}_{j}^{n+\tfrac{1}{2}}\left({\bar{p}}_{j+1}^{n+\tfrac{1}{2}}+{\bar{p}}_{j-1}^{n+\tfrac{1}{2}}\right)\right|,\end{array}\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}\left|{p}_{{j}_{2}}^{{n}_{2}+\tfrac{1}{2}}\left(\left|{p}_{{j}_{2}+1}^{{n}_{2}+\tfrac{1}{2}}\right|+\left|{p}_{{j}_{2}-1}^{{n}_{2}+\tfrac{1}{2}}\right|\right)\right|\\ =\,\mathop{\max }\limits_{0\leqslant j\leqslant N;n\geqslant 0}\left|{p}_{j}^{n+\tfrac{1}{2}}\right|\left(\left|{p}_{j+1}^{n+\tfrac{1}{2}}\right|+\left|{p}_{j-1}^{n+\tfrac{1}{2}}\right|\right).\end{array}\end{eqnarray}$Then, we perform the von Neumann analysis for the corresponding linear equation and obtain the following result.
The CN-type finite difference scheme(20)-(22) to equation (2) is unconditionally linearly stable.
Replacing the nonlinear terms of equation (20) like the way in(28), we have the following linear equation:$\begin{eqnarray}\begin{array}{l}{p}_{j}^{n+1}-{p}_{j}^{n}+\displaystyle \frac{\tau }{4{h}^{3}}\left({p}_{j+2}^{n+1}-2{p}_{j+1}^{n+1}+2{p}_{j-1}^{n+1}\right.\\ \quad \left.-\,{p}_{j-2}^{n+1}+{p}_{j+2}^{n}-2{p}_{j+1}^{n}+2{p}_{j-1}^{n}-{p}_{j-2}^{n}\right)\\ \quad +\displaystyle \frac{3\tau }{4h}{{\mathfrak{M}}}_{1}\left({p}_{j+1}^{n+1}-{p}_{j+1}^{n}+{p}_{j-1}^{n+1}-{p}_{j-1}^{n}\right)\\ \quad +\,\displaystyle \frac{3\tau }{4h}{{\mathfrak{M}}}_{2}\left(\left|{p}_{j+1}^{n+1}+{p}_{j+1}^{n}\right|-\left|{p}_{j-1}^{n+1}+{p}_{j-1}^{n}\right|\right)=0.\end{array}\end{eqnarray}$Then, we assume ${p}_{j}^{n}={v}^{n}{{\rm{e}}}^{{\rm{i}}{kjh}}$ and substitute it into equation (31), yielding$\begin{eqnarray}\begin{array}{l}{v}^{n+1}-{v}^{n}+\displaystyle \frac{\tau }{4{h}^{3}}\left({v}^{n+1}{{\rm{e}}}^{2{\rm{i}}{kh}}-2{v}^{n+1}{{\rm{e}}}^{{\rm{i}}{kh}}+2{v}^{n+1}{{\rm{e}}}^{-{\rm{i}}{kh}}\right.\\ \left.-\,{v}^{n+1}{{\rm{e}}}^{-2{\rm{i}}{kh}}+{v}^{n}{{\rm{e}}}^{2{\rm{i}}{kh}}-2{v}^{n}{{\rm{e}}}^{{\rm{i}}{kh}}+2{v}^{n}{{\rm{e}}}^{-{\rm{i}}{kh}}-{v}^{n}{{\rm{e}}}^{-2{\rm{i}}{kh}}\right)\\ +\,\displaystyle \frac{3\tau {{\mathfrak{M}}}_{1}}{4h}\left({v}^{n+1}{{\rm{e}}}^{{\rm{i}}{kh}}+{v}^{n}{{\rm{e}}}^{{\rm{i}}{kh}}-{v}^{n+1}{{\rm{e}}}^{-{\rm{i}}{kh}}-{v}^{n}{{\rm{e}}}^{-{\rm{i}}{kh}}\right)=0,\end{array}\end{eqnarray}$that is,$\begin{eqnarray}\begin{array}{l}\left(1+\displaystyle \frac{\tau }{4{h}^{3}}{{\rm{e}}}^{2{\rm{i}}{kh}}-\displaystyle \frac{\tau }{2{h}^{3}}{{\rm{e}}}^{{\rm{i}}{kh}}+\displaystyle \frac{\tau }{2{h}^{3}}{{\rm{e}}}^{-{\rm{i}}{kh}}\right.\\ \quad \left.-\,\displaystyle \frac{\tau }{4{h}^{3}}{{\rm{e}}}^{-2{\rm{i}}{kh}}+\displaystyle \frac{3\tau {{\mathfrak{M}}}_{1}}{4h}{{\rm{e}}}^{{\rm{i}}{kh}}-\displaystyle \frac{3\tau {{\mathfrak{M}}}_{1}}{4h}{{\rm{e}}}^{-{\rm{i}}{kh}}\right){v}^{n+1}\\ =\,\left(1-\displaystyle \frac{\tau }{4{h}^{3}}{{\rm{e}}}^{2{\rm{i}}{kh}}+\displaystyle \frac{\tau }{2{h}^{3}}{{\rm{e}}}^{{\rm{i}}{kh}}-\displaystyle \frac{\tau }{2{h}^{3}}{{\rm{e}}}^{-{\rm{i}}{kh}}+\displaystyle \frac{\tau }{4{h}^{3}}{{\rm{e}}}^{-2{\rm{i}}{kh}}\right.\\ \quad \left.-\,\displaystyle \frac{3\tau {{\mathfrak{M}}}_{1}}{4h}{{\rm{e}}}^{{\rm{i}}{kh}}+\displaystyle \frac{3\tau {{\mathfrak{M}}}_{1}}{4h}{{\rm{e}}}^{-{\rm{i}}{kh}}\right){v}^{n}.\end{array}\end{eqnarray}$From the above equation, the amplification factor can be obtained as$\begin{eqnarray}v=\displaystyle \frac{1-{\rm{i}}\left[\tfrac{\tau }{2{h}^{3}}\sin (2{kh})-\tfrac{\tau }{{h}^{3}}\sin ({kh})+\tfrac{3\tau }{2h}{{\mathfrak{M}}}_{1}\sin ({kh})\right]}{1+{\rm{i}}\left[\tfrac{\tau }{2{h}^{3}}\sin (2{kh})-\tfrac{\tau }{{h}^{3}}\sin ({kh})+\tfrac{3\tau }{2h}{{\mathfrak{M}}}_{1}\sin ({kh})\right]},\end{eqnarray}$where ∣v∣≡1 implies that the finite difference scheme(20)-(22) is unconditionally linearly stable. □
Although the above finite difference scheme has the second-order accuracy both in space and time and enjoys the unconditional linear stability, one has to solve a set of fully nonlinear coupled algebraic equations at every time step, which is time consuming and tedious in programming. In order to improve the computational efficiency and make the algorithm easy to program, we adopt the iterative method [58, 59] and the Thomas algorithm together to solve equation (20). To be specific, the iterative scheme can be described as follows:$\begin{eqnarray}\begin{array}{l}\,-\displaystyle \frac{\tau }{4{h}^{3}}{p}_{j-2}^{n+1,l+1}+\left(\displaystyle \frac{\tau }{2{h}^{3}}-6{A}_{j}^{n,l}\right){p}_{j-1}^{n+1,l+1}+\,(1+3{B}_{j}^{n,l}){p}_{j}^{n+1,l+1}+\left(6{A}_{j}^{n,l}-\displaystyle \frac{\tau }{2{h}^{3}}\right){p}_{j+1}^{n+1,l+1}+\displaystyle \frac{\tau }{4{h}^{3}}{p}_{j+2}^{n+1,l+1}\\ \,=\,\displaystyle \frac{\tau }{4{h}^{3}}{p}_{j-2}^{n}+\left(6{A}_{j}^{n,l}-\displaystyle \frac{\tau }{2{h}^{3}}\right){p}_{j-1}^{n}+\,(1-3{B}_{j}^{n,l}){p}_{j}^{n}\\ \,+\,\left(\displaystyle \frac{\tau }{2{h}^{3}}-6{A}_{j}^{n,l}\right){p}_{j+1}^{n}-\displaystyle \frac{\tau }{4{h}^{3}}{p}_{j+2}^{n},\quad 1\leqslant j\leqslant N-1,\,n\geqslant 0,\end{array}\end{eqnarray}$with$\begin{eqnarray*}\begin{array}{rcl} & & {A}_{j}^{n,l}=\displaystyle \frac{\tau }{32h}\left({p}_{j}^{n+1,l}+{p}_{j}^{n}\right)\left({\bar{p}}_{j+1}^{n+1,l}+{\bar{p}}_{j+1}^{n}+{\bar{p}}_{j-1}^{n+1,l}+{\bar{p}}_{j-1}^{n}\right),\\ & & {B}_{j}^{n,l}=\displaystyle \frac{\tau }{16h}\left({\left|{p}_{j+1}^{n+1,l}+{p}_{j+1}^{n}\right|}^{2}-{\left|{p}_{j-1}^{n+1,l}+{p}_{j-1}^{n}\right|}^{2}\right),\\ & & {p}_{j}^{n+1,0}={p}_{j}^{n},\quad \mathop{\mathrm{lim}}\limits_{l\to +\infty }{p}_{j}^{n+1,l+1}={p}_{j}^{n+1}.\end{array}\end{eqnarray*}$From tn to tn+1, the iteration stops when$\begin{eqnarray}\mathop{\max }\limits_{1\leqslant j\leqslant N-1}\left|{p}_{j}^{n+1,l+1}-{p}_{j}^{n+1,l}\right|\lt \epsilon ,\end{eqnarray}$where l=0,1,2,⋯ denotes the iterative time, ε is a given error bond. For simplicity, the iterative equation (35) can be written in the matrix form$\begin{eqnarray}{{\boldsymbol{C}}}_{+}^{n,l}{{\boldsymbol{p}}}^{n+1,l}={{\boldsymbol{C}}}_{-}^{n,l}{{\boldsymbol{p}}}^{n,l},\end{eqnarray}$with$\begin{eqnarray}{{\boldsymbol{C}}}_{\pm }^{n,l}=\left(\begin{array}{cccccc}{a}_{1}^{\pm } & {b}_{1}^{\pm } & \pm c & & & 0\\ -{b}_{2}^{\pm } & \ddots & \ddots & \ddots & & \\ \mp c & \ddots & \ddots & \ddots & \ddots & \\ & \ddots & \ddots & \ddots & \ddots & \pm c\\ & & \ddots & \ddots & \ddots & {b}_{N-2}^{\pm }\\ 0 & & & \mp c & -{b}_{N-1}^{\pm } & {a}_{N-1}^{\pm }\end{array}\right),\quad {{\boldsymbol{p}}}^{n,l}=\left(\begin{array}{c}{p}_{1}^{n,l}\\ {p}_{2}^{n,l}\\ \vdots \\ \vdots \\ {p}_{N-2}^{n,l}\\ {p}_{N-1}^{n,l}\end{array}\right),\end{eqnarray}$where the coefficients are given by$\begin{eqnarray}\begin{array}{rcl}{a}_{j}^{\pm } & = & 1\pm 3{B}_{j}^{n,l},\quad {b}_{j}^{\pm }=\pm \left(6{A}_{j}^{n,l}-\displaystyle \frac{\tau }{2{h}^{3}}\right),\\ c & = & \displaystyle \frac{\tau }{4{h}^{3}},\quad 1\leqslant j\leqslant N-1.\end{array}\end{eqnarray}$At each step, we need to compute the coefficient matrices ${{\bf{C}}}_{\pm }^{n,l}$ from pn,l. Thus, the penta-diagonal system of equation (37) can be solved via the Thomas algorithm, and the computational cost per time step is O(N).
4. Numerical results
In this section, we will numerically simulate the soliton evolution and collisions in equation (2) with the vanishing background, and test the numerical accuracy and linear stability of the adopted finite difference method. Throughout this section, we use the l∞-norm of error between the numerical solution ${p}_{j}^{n}$ and the analytical solution p(ηj,ξn) as$\begin{eqnarray}{e}_{\infty }({t}_{n}):=\mathop{\max }\limits_{1\leqslant j\leqslant N-1}| p({\xi }_{n},{\eta }_{j})-{p}_{j}^{n}| ,\end{eqnarray}$to quantify the numerical solutions.
4.1. Numerical simulation of one-soliton evolution
First, we simulate the one-soliton evolution by choosing solution(6) at ξ=−15 as the initial data. In our simulations, we take the interval Ω=[−30,30] so that the boundaries do not affect the soliton propagation. Meanwhile, we set h=0.05,τ=0.001 and error bond ε=10−11 for the iterative computation of equation (35). With different values of β1, γ1 and λ1, figures 2(a)-(c) respectively illustrate the time evolution of the SH, DH and MH solitons, which shows the same dynamical behavior as the exact analytical solutions behave [27]. More precisely, we calculate the error norms between the numerical and analytical solutions for ξ∈[0,30], as shown in figure 3. It can be seen that the error norms for three cases in figures 2(a)-(c) are within the order 10−3. Besides, we give the time evolution of the discrete masses defined in(24) from ξ=0 to 30 (see figure 4), which shows that M1 is a conserved quantity and agrees well with the exact value as given in(7).
Figure 2.
New window|Download| PPT slide Figure 2.Time evolution for three different one-soliton solutions of equation (2): (a) an SH soliton with β1=0, γ1=2 and λ1=−0.25+0.2i; (b) a DH soliton with β1=2, γ1=0 and λ1=−0.3+0.06i; (c) an MH soliton with β1=0.5, γ1=2 and ${\lambda }_{1}=-0.2+\tfrac{\pi }{16}\,{\rm{i}}$.
Figure 3.
New window|Download| PPT slide Figure 3.Time evolution of error norms for three cases of one-soliton solutions in figures 2(a)-(c).
Figure 4.
New window|Download| PPT slide Figure 4.Time evolution of discrete masses for three cases of one-soliton solutions in figures 2(a)-(c).
We should mention that the time-splitting pseudo-spectral (TSPS) method is very effective in computing the envelope solitons of NLSE [60]. In comparison, the TSPS method performs better in the accuracy and costs less CPU time than the finite difference method. However, the former requires a judicious choice of the time step and mesh size. For example, with the same τ,h and initial values as adopted in figures 2(a)-(c), the TSPS method shows a good performance in simulating the SH and DH solitons, but there occurs a severe distortion for the MH soliton when the evolution time is less than 10.
4.2. Numerical simulation of two-soliton collisions
Next, we simulate the collisions between two solitons for equation (2) by the same finite difference method. To do so, we take solution(3) at ξ=−20 as the initial data, and solve the problem numerically on the interval Ω=[−40,40] with the mesh size h = 0.01, time step τ=0.001 and error bond ε=10−11. In figures 5(a) and (b), we demonstrate the shape-preserving collisions between an SH soliton and a DH soliton and between two MH solitons, respectively. Ignoring the internal oscillation of MH solitons, all the interacting solitons in such two cases retain their individual shapes and velocities, and experience only a small phase shift after the collision. Differently, figure 5(c) depicts a shape-changing two-soliton collision, that is, an SH soliton changes into an MH one when interacting with another MH soliton. However, this type of collision is also elastic in the sense that the masses of interacting solitons are conserved, which can be confirmed in figure 8. All the collisions are clean and display no dispersed radiation (for example, see figure 6), indicating the elastic nature.
Figure 5.
New window|Download| PPT slide Figure 5.Time evolution of three different two-soliton solutions of equation (2): (a) shape-preserving collision between an SH soliton and a DH soliton with β1=β2=1, γ1=γ2=0, λ1=−0.32−0.1i and λ2=0.32−0.3i; (b) shape-preserving collision between two MH solitons with β1=β2=1, γ1=γ2=1, λ1=−0.375−0.12i and λ2=0.375−0.3i; (c) shape-changing collision between an MH soliton and a DH soliton with β1=β2=1, γ1=0, γ2=1, λ1=0.4−0.08i and λ2=−0.3−0.3i.
Figure 6.
New window|Download| PPT slide Figure 6.Transverse plots of the two-soliton collision associated with figure 6(a) at (a) ξ=0, (b) ξ=17, (c) ξ=23, and (d) ξ=40.
Moreover, we give the time evolution of error norms and discrete masses from ξ=0 to 40 in figures 7 and 8. The error norms for three cases in figures 5(a)-(c) are within the order 10−3, but the last two cases exhibit a stronger oscillation with the time ξ because of the periodical behavior of MH solitons. From figure 8, we find that there is very little dissipation for the discrete masses and the values are very close to 4(∣λ1R∣+∣λ2R∣), which implies that the total masses of two solitons are conserved before and after the collision.
Figure 7.
New window|Download| PPT slide Figure 7.Time evolution of error norms for three cases of two-soliton solutions in figures 5(a)-(c).
Figure 8.
New window|Download| PPT slide Figure 8.Time evolution of discrete masses for three cases of two-soliton solutions in figures 5(a)-(c).
4.3. Numerical accuracy and linear stability
In this subsection, we use the one-soliton solution in(6) to numerically check the accuracy and linear stability of the CN-type finite difference method.
By setting a fixed time step but different mesh sizes, table 1 shows the error norms at the time ξ=2 for three different one-soliton solutions, where the parameters are selected as β1=1,γ1=0,λ1=−0.25+0.2i for the SH soliton, β1=1,γ1=0,λ1=−0.3+0.075i for the DH soliton, and ${\beta }_{1}=0.5,{\gamma }_{1}=1,{\lambda }_{1}=-0.3+\tfrac{\pi }{12}\,{\rm{i}}$ for the MH soliton, respectively. Likewise, table 2 presents the error norms at the time ξ=2 with a fixed mesh size but different time steps, where the parameters for the three one-soliton cases are the same as those in table 1. Evidently, the numerical experiments verify that the CN-type finite difference method is second-order accurate both in space and time, which is consistent with the statement in lemma 3.2.
Table 1. Table 1.Spatial error analysis of three one-soliton solutions with time step τ=0.01 and different mesh sizes.
For the study of linear stability, we perturb the initial solution p1 in the form$\begin{eqnarray}p={p}_{1}+\varepsilon {{\rm{e}}}^{-{\left(\eta -{\eta }_{0}\right)}^{2}},\end{eqnarray}$where ${\eta }_{0}=4\left({\lambda }_{1{\rm{R}}}^{2}-3\,{\lambda }_{1{\rm{I}}}^{2}\right){\xi }_{0}-\tfrac{\mathrm{ln}\sqrt{\tfrac{| {\lambda }_{1{\rm{I}}}| }{| {\delta }_{1}| }}}{2{\lambda }_{1{\rm{R}}}}$ represents the central position of soliton envelope at the initial time ξ=ξ0. The problem is solved numerically on the interval Ω=[−30,30] with the mesh size h = 0.005, time step τ=0.001 and error bond ε=10−11. By selecting different values for β1,γ1,λ1 and taking $\varepsilon =\tfrac{1}{50},\tfrac{1}{100},\tfrac{1}{200}$, we depict the time evolution of error norms between the numerical solutions of equation (2) with the initial data(41) and its exact analytical solution(6), as seen in figure 9.
Figure 9.
New window|Download| PPT slide Figure 9.Time evolution of the error norms from ξ=0 to 10 for three types of one-soliton solutions with different initial perturbations, where the parameters are chosen as follows: (a) β1=0, γ1=2 and λ1=−0.25+0.2i, (b) β1=2, γ1=0 and λ1=−0.3+0.06i, (c) β1=0.5, γ1=2 and ${\lambda }_{1}=-0.2+\tfrac{\pi }{16}\,{\rm{i}}$.
Regarding the above three one-soliton cases, figure 10 makes a comparison of the exact analytical solution and numerical solutions at ξ=10 with $\varepsilon =\tfrac{1}{100}$ in(41). It can be seen that the soliton solutions enjoy a good linear stability against the initial perturbations.
Figure 10.
New window|Download| PPT slide Figure 10.Comparison of the exact analytical solution and numerical solutions at ξ=10 for (a) SH soliton, (b) DH soliton, and (c) MH soliton, where $\varepsilon =\tfrac{1}{100}$ and the other parameters are the same as those in figure 9.
5. Conclusions
Recently, the soliton dynamics of the complex MKdV equation (2) (which is an integrable equation equivalent to the SS equation) has attracted much interest both in mathematics and physics. In this paper, we have proposed a CN-type finite difference method for computing the soliton solutions of equation (2) with the vanishing boundary condition. It has been shown that the resulting numerical scheme admits the second-order accuracy both in space and time, conserves the discrete mass, and enjoys the unconditional linear stability in the sense of the von Nuemann analysis. For solving a fully nonlinear coupled algebraic equations obtained from equation (20), we have adopted an iterative method and the Thomas algorithm to enhance the computational efficiency. In numerical experiments, we have simulated the single-soliton propagation and two-soliton collisions, which shows a good performance in the numerical accuracy, mass conservation and linear stability against with the initial perturbations. Therefore, our proposed numerical method is efficient in simulating the soliton solutions of equation (2) with the vanishing boundary condition. In the future, considering that equation (2) admits more abundant and complicated solitonic behaviors with the nonvanishing background [37, 40, 43], one may explore the finite difference scheme to equation (2) with the Neumann boundary condition [60]. However, it will be a challenging work to maintain the discrete mass conservation.
Acknowledgments
This work was partially supported by the Natural Science Foundation of Beijing Municipality (Grant No. 1212007), and by the Science Foundations of China University of Petroleum, Beijing (Grant Nos. 2 462 020YXZZ004 and 2 462 020XKJS02).
AblowitzM JClarksonP A1992Solitons, Nonlinear Evolution Equations and Inverse Scattering Cambridge Cambridge University Press [Cited within: 1]
DaiC QFanYZhangN2019 Re-observation on localized waves constructed by variable separation solutions of (1+1)-dimensional coupled integrable dispersionless equations via the projective Riccati equation method Appl. Math. Lett. 96 2026 DOI:10.1016/j.aml.2019.04.009
WangY YChenLDaiC QZhengJFanY2017 Exact vector multipole and vortex solitons in the media with spatially modulated cubic-quintic nonlinearity Nonlinear Dyn. 90 12691275 DOI:10.1007/s11071-017-3725-5
LiMYueX LXuT2020 Multi-pole solutions and their asymptotic analysis of the focusing Ablowitz-Ladik equation Phys. Scr. 95 055222 DOI:10.1088/1402-4896/ab4503
XuTChenYLiMMengD X2019 General stationary solutions of the nonlocal nonlinear Schrödinger equation and their relevance to the PT-symmetric system Chaos 29 123124 DOI:10.1063/1.5121776
LiMShuiJ JHuangY HWangLLiH J2018 Localized-wave interactions for the discrete nonlinear Schrödinger equation under the nonvanishing background Phys. Scr. 93 115203 DOI:10.1088/1402-4896/aae213 [Cited within: 1]
AgrawalG P2002Nonlinear Fiber Optics San Diego, CA Academic [Cited within: 3]
XuTTianBQiF-H2012 Bright N-soliton solution to the vector Hirota equation from nonlinear optics with symbolic computation Z. Naturforsch. A 67 3949 DOI:10.5560/zna.2011-0055 [Cited within: 4]
ZhangG QChenS YYanZ Y2020 Focusing and defocusing Hirota equations with non-zero boundary conditions: Inverse scattering transforms and soliton solutions Commun. Nonlinear Sci. Numer. Simul. 80 104927 DOI:10.1016/j.cnsns.2019.104927 [Cited within: 1]
MihalacheDTornerLMoldoveanuFPanoiuN CTrutaN1993 Inverse-scattering approach to femtosecond solitons in monomode optical fibers Phys. Rev. E 48 4699 DOI:10.1103/PhysRevE.48.4699 [Cited within: 2]
YangBChenY2019 High-order soliton matrices for Sasa-Satsuma equation via local Riemann-Hilbert problem Nonlinear Anal.-Real 45 918941 DOI:10.1016/j.nonrwa.2018.08.004 [Cited within: 1]
SwatersG EDosserH VSutherlandB R2012 Conservation laws, Hamiltonian structure, modulational instability properties and solitary wave solutions for a higher-order model describing nonlinear internal waves Stud. Appl. Math. 128 159182 DOI:10.1111/j.1467-9590.2011.00533.x [Cited within: 2]
SergyeyevADemskoiD2007 Sasa-Satsuma (complex modified Korteweg-de Vries II) and the complex sine-Gordon II equation revisited: Recursion operators, nonlocal symmetries, and more J. Math. Phys. 48 042702 DOI:10.1063/1.2710552 [Cited within: 1]
ZhangH QYuanS S2017 Dark soliton solutions of the defocusing Hirota equation by the binary Darboux transformation Nonlinear Dyn. 89 531538 DOI:10.1007/s11071-017-3469-2 [Cited within: 3]
TaoY SHeJ S2012 Multisolitons, breathers, and rogue waves for the Hirota equation generated by the Darboux transformation Phys. Rev. E 85 026601 DOI:10.1103/PhysRevE.85.026601
AnkiewiczASoto-CrespoJ MAkhmedievN2010 Rogue waves and rational solutions of the Hirota equation Phys. Rev. E 81 046602 DOI:10.1103/PhysRevE.81.046602
LiuCRenYYangZ YYangW L2017 Superregular breathers in a complex modified Korteweg-de Vries system Chaos 27 083120 DOI:10.1063/1.4999916
LiMZhangX FXuTLiL L2020 Asymptotic analysis and soliton interactions of the multi-pole solutions in the Hirota equation J. Phys. Soc. Japan. 89 054004 DOI:10.7566/JPSJ.89.054004
Soto-CrespoJ MDevineNHoffmannN PAkhmedievN2014 Rogue waves of the Sasa-Satsuma equation in a chaotic wave field Phys. Rev. E 90 032902 DOI:10.1103/PhysRevE.90.032902
MuGQinZ Y2016 Dynamic patterns of high-order rogue waves for Sasa-Satsuma equation Nonlinear Anal.-Real 31 179209 DOI:10.1016/j.nonrwa.2016.01.001
BandelowUAkhmedievN2012 Sasa-Satsuma equation: soliton on a background and its limiting cases Phys. Rev. E 86 026606 DOI:10.1103/PhysRevE.86.026606
ChenS H2013 Twisted rogue-wave pairs in the Sasa-Satsuma equation Phys. Rev. E 88 023202 DOI:10.1103/PhysRevE.88.023202
XuTLiMLiL2015 Anti-dark and Mexican-hat solitons in the Sasa-Satsuma equation on the continuous wave background Europhys. Lett. 109 30006 DOI:10.1209/0295-5075/109/30006 [Cited within: 2]
MarutaAInoueTNonakaYYoshikaY2002 Bi-soliton propagating in dispersion-managed system and its application to high-speed and long-haul optical transmission IEEE J. Sel. Top. Quantum Electron. 8 640650 DOI:10.1109/JSTQE.2002.1016368 [Cited within: 1]
ShirakiENishizawaNItohK2011 Ultrafast all-optical signal regenerator using pulse trapping in birefringent fibers J. Opt. Soc. Am. B 28 26432649 DOI:10.1364/JOSAB.28.002643 [Cited within: 1]
NishizawaNGotoT2003 Ultrafast all optical switching by use of pulse trapping across zero-dispersion wavelength Opt. Express 11 359365 DOI:10.1364/OE.11.000359 [Cited within: 1]
Cisneros-AkeL APradoH PLópez VillatoroD JCarretero-GonzálezR2018 Multi-hump bright solitons in a Schrödinger-mKdV system Phys. Lett. A 382 837845 DOI:10.1016/j.physleta.2018.01.031 [Cited within: 1]
PradoH PCisneros-AkeL A2020 The direct method for multisolitons and two-hump solitons in the Hirota-Satsuma system Phys. Lett. A 384 126471 DOI:10.1016/j.physleta.2020.126471
PradoH PCisneros-AkeL A2019 Multi-hump bright and dark solitons for the Schröinger-Korteweg-de Vries coupled system Chaos 29 053133 DOI:10.1063/1.5092985 [Cited within: 1]
KivsharYu SAgrawalG P2003Optical Solitons: From Fibers to Photonic Crystals San Diego Academic [Cited within: 1]
FengB FMitsuiT1998 A finite difference method for the Korteweg-de Vries and the Kadomtsev-Petviashvili equations J. Comput. Appl. Math. 90 95116 DOI:10.1016/S0377-0427(98)00006-5 [Cited within: 1]
WangT CGuoB LXuQ B2013 Fourth-order compact and energy conservative difference schemes for the nonlinear Schrödinger equation in two dimensions J. Comput. Phys. 243 382399 DOI:10.1016/j.jcp.2013.03.007 [Cited within: 1]
LiS CLiX GCaoJ JLiW B2017 High-order numerical method for the derivative nonlinear Schrödinger equation Int. J. Model. Simul. Sci. Comput. 8 1750017 DOI:10.1142/S1793962317500179 [Cited within: 1]
BaoW ZTangQ LXuZ G2013 Numerical methods and comparison for computing dark and bright solitons in the nonlinear Schrödinger equation J. Comput. Phys. 235 423445 DOI:10.1016/j.jcp.2012.10.054 [Cited within: 2]
,1,2, Guowei Zhang2, Liqun Wang2, Xiangmin Xu2, Min Li31State Key Laboratory of Heavy Oil Processing, 
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
New window|Download| PPT slide
