Numerical simulation of the soliton solutions for a complex modified Korteweg-de Vries equation by a finite difference method

Tao Xu,Guowei Zhang,Liqun Wang,Xiangmin Xu and Min Li

1 State Key Laboratory of Heavy Oil Processing,China University of Petroleum,Beijing 102249,China

2 College of Science,China University of Petroleum,Beijing 102249,China

3 Department of Mathematics and Physics,North China Electric Power University,Beijing 102206,China

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

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]

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 singlehump (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],meshfree 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 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]

in which the block matrices FN×L,GN×Land HN×Lare given by

with

Taking N=1 in solution (3),the one-soliton solution can be written as

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

Figure 1.Three types of one-soliton solutions of equation (2): (a) Aan SH soliton with β1=0,γ1=1 and b) a DH soliton with β1=1,γ1=0 and c) an MH soliton with β1=0.5,γ1=1 an

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 βkand γ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 lim∣η∣→∞p=0 can be approximated by

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

Meanwhile,based on the Taylor series expansion,we give the following three useful formulas:

where a,b,c are three positive integers.For simplicity,we introduce the finite difference operators

withbeing 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

Based on the above equations,we obtain the following finite difference scheme of equation (8) as follows:

where the CN scheme is used for the temporal derivative,the central finite differences is used for the spatial derivatives,andMeanwhile,the boundary condition (9) is discretized as

and the initial condition (10) is discretized as

Remark 3.1.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:

Lemma 3.2.The local truncation error of the difference equation (20) to equation (2) is of O(τ2+ h2).

Proof.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.

Lemma 3.3.If we define the discrete mass as

then the CN-type finite difference method (20)-(22)conserves the mass in the discrete level,i.e.

Proof.Multiplying equation (20) byand adding the resulted equation with its complex conjugate,we obtain

Taking the summation of equation (26) from j=1 to N - 1 and noticing (24),we have

This,together with the boundary conditions (21) and (23),implies that=0for 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,

where (jk,nk) (k=1,2) are the indices such that

Then,we perform the von Neumann analysis for the corresponding linear equation and obtain the following result.

Lemma 3.4.The CN-type finite difference scheme (20)-(22)to equation (2) is unconditionally linearly stable.

Proof.Replacing the nonlinear terms of equation(20)like the way in (28),we have the following linear equation:

Then,we assume=vneikjhand substitute it into equation (31),yielding

that is,

From the above equation,the amplification factor can be obtained as

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:

with

From tnto tn+1,the iteration stops when

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

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 i .

with

where the coefficients are given by

At each step,we need to compute the coefficient matricesfrom 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 solutionand the analytical solution p(ηj,ξn) as

to quantify the numerical solutions.

4.1.Numerical simulation of one-soliton evolution

Figure 3.Time evolution of error norms for three cases of onesoliton solutions in figures 2(a)-(c).

Figure 4.Time evolution of discrete masses for three cases of onesoliton solutions in figures 2(a)-(c).

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.Transverse plots of the two-soliton collision associated with figure 6(a) at (a) ξ=0,(b) ξ=17,(c) ξ=23,and (d) ξ=40.

Figure 7.Time evolution of error norms for three cases of twosoliton solutions in figures 5(a)-(c).

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-11for the iterative computation of equation (35).With different values of β1,γ1and λ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 M1is a conserved quantity and agrees well with the exact value as given in (7).

Figure 8.Time evolution of discrete masses for three cases of twosoliton solutions in figures 5(a)-(c).

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

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 and the other parameters are the same as those in figure 9.

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.

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.

Table 1.Spatial error analysis of three one-soliton solutions with time step τ=0.01 and different mesh sizes.

Table 2.Time error analysis of three one-soliton solutions with mesh step h=0.001 and different time steps.

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 CNtype 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β1= 0.5,γ1= 1,λ1=- 0 .3 +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.

For the study of linear stability,we perturb the initial solution p1in the form

Regarding the above three one-soliton cases,figure 10 makes a comparison of the exact analytical solution and numerical solutions at ξ=10 within (41).It can be seen that the soliton solutions enjoy a good linear stability against the initial perturbations.

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).