Lattice Boltzmann model for a generalized Gardner equation with time-dependent variable coefficients

Accepted Manuscript

Lattice Boltzmann model for a generalized Gardner equation with time-dependent variable coefficients Wen-Qiang Hu, Yi-Tian Gao, Zhong-Zhou Lan, Chuan-Qi Su, Yu-Jie Feng PII: DOI: Reference:

S0307-904X(17)30067-7 10.1016/j.apm.2017.01.061 APM 11559

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

23 June 2016 11 January 2017 17 January 2017

Please cite this article as: Wen-Qiang Hu, Yi-Tian Gao, Zhong-Zhou Lan, Chuan-Qi Su, Yu-Jie Feng, Lattice Boltzmann model for a generalized Gardner equation with time-dependent variable coefficients, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.01.061

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A lattice Boltzmann model for variable-coefficient Gardner equation is derived. • The effects of the free parameters induced in our model are discussed in detail. • Different types solutions for vc-Gardner equation have been simulated.

AC

CE

PT

ED

M

AN US

CR IP T

• Numerical solutions agree well with analytical solutions.

1

ACCEPTED MANUSCRIPT

Lattice Boltzmann model for a generalized Gardner equation with time-dependent variable coefficients

CR IP T

Wen-Qiang Hu, Yi-Tian Gao ∗, Zhong-Zhou Lan, Chuan-Qi Su, Yu-Jie Feng Ministry-of-Education Key Laboratory of Fluid Mechanics and National Laboratory for Computational Fluid Dynamics, Beijing University of

AN US

Aeronautics and Astronautics, Beijing 100191, China

Abstract

In this paper, lattice Boltzmann model for a generalized Gardner equation with timedependent variable coefficients, which can provide some more realistic models than their constant-coefficient counterparts, is derived through selecting equilibrium distribution func-

M

tion and adding the compensate function, appropriately. Effects and approximate value range of the free parameters, which are introduced to adjust the single relaxation time and

ED

equilibrium distribution function, are discussed in detail, as well as the impact of the lattice space step and velocity. Numerical simulations in different situations of this equation are conducted, including the propagation and interaction of the solitons, the evolution of the

PT

non-propagating soliton and the propagation of the double-pole solutions. It is found that the numerical results match well with the analytical solutions, which demonstrates that the

AC

CE

current lattice Boltzmann model is a satisfactory and efficient algorithm.

Keywords: Lattice Boltzmann model; Generalized Gardner equation with time-dependent variable coefficients; Numerical simulations; Soliton solutions

∗

Corresponding author, with e-mail address as [email protected]

2

ACCEPTED MANUSCRIPT

1. Introduction

AC

CE

PT

ED

M

AN US

CR IP T

As encountered in some branches of science and engineering, e.g., hydrodynamics, plasmas, elastic media and optical communication [1–4], nonlinear evolution equations (NLEEs) have played a crucial role in the study of varieties of nonlinear phenomena, especially to analyze the propagation of waves. Since it can provide much physical information and more insight into the physical aspects of the problems and thus lead to further applications, the investigation for the solutions of the NLEEs plays an important role in nonlinear science fields [5–7]. However, due to the complexity of the NLEEs, there is no unified method to find all solutions of the NLEEs, which makes the investigation of the solving methods for the NLEEs as an attractive research undertaking [8]. For decades, various powerful methods have been proposed to investigate the analytical and numerical solutions of the NLEEs, such as the Hirota bilinear method [9] and Darboux transformation [10] for analytical solutions, and finite difference method [11], finite element method [12] and spectral method [13] for numerical solutions. Among multifarious numerical algorithms, the lattice Boltzmann method (LBM), as an effective numerical approach, has achieved success in simulating some complex fluid flows [14–17]. Recently, LBM has been extended successfully to simulate some NLEEs, such as anisotropic dispersion equation [18], convection-diffusion equation [19], Burgers equation [20], Korteweg-de Vries (KdV) equation [21], modified Korteweg-de Vries (mKdV) equation [22], combined KdVmKdv equation [23], KdV-Burgers equation [24], the ion- and electron-acoustic solitary waves in beam-plasma system [25], and the generalized nonlinear wave equations [8, 26, 27]. Unlike traditional numerical methods which discretize the governing equations in time and space, LBM is based on kinetic theory which tracks the dynamics of microcosmic particle ensembles [16]. Through the particle distribution function and equilibrium distribution function, the macroscopic variables are educed and the macroscopic equations are restored exactly [27]. The above mentioned numerical studies [8, 18–25] for NLEEs via the LBM show that the LBM is an effective and accurate numerical algorithm for NLEEs. However, these studies are all considering about the constant-coefficient NLEEs, while the variable-coefficient NLEEs can be used to describe some complex and realistic phenomena with the consideration of inhomogeneous media and nonuniform boundaries [28, 29], which have not been investigated by LBM, to our knowledge. Hence, in this paper, we focus on a generalized Gardner equation with time-dependent variable coefficients, which has been proposed to model the propagation of weakly nonlinear long waves in a KdV-type medium that is characterized by varying dispersive and nonlinear coefficients [30, 31], or describe the atmospheric blocking phenomenon [32], ut − 6f0 (t)uux − 6f1 (t)u2 ux + f2 (t)uxxx − f3 (t)ux + f4 (t)(Cu + xux ) = 0 ,

(1)

where u is the wave-amplitude function of the scaled space coordinate x and time coordinate t, C is a nonzero constant, f0 (t), f1 (t) and f2 (t) correspond to the quadratic nonlinear, cubic 3

ACCEPTED MANUSCRIPT

nonlinear and dispersive coefficients, respectively, while −f3 (t) + f4 (t)x and f4 (t)C denote the dissipative and line-damping terms [33]. A group classification problems can be represented by Eq. (1) through different assignments of the continuous function. Some special cases of Eq. (1) have been acquired as follows, (I) When f0 (t) = − α6 , f2 (t) = β, f1 (t) = f3 (t) = f4 (t) = 0, Eq. (1) is transformed into Korteweg-de Vries (KdV) equation for the surface waves in the shallow water [34], (2)

CR IP T

ut + αuux + βuxxx = 0.

(II) When f0 (t) = − α3 , f1 (t) = β2 , f2 (t) = 1, f3 (t) = f4 (t) = 0, Eq. (1) is transformed into the constant-coefficient Gardner equation, or the combined KdV-modified KdV equation, for the wave propagation of the bound particle, sound wave and thermal pulse [35], ut + 2αuux − 3βu2 ux + uxxx = 0.

(3)

AN US

Another type of the constant-coefficient Gardner equation is of the following form with f0 (t) = − α6 , f1 (t) = β6 , f2 (t) = δ, f3 (t) = f4 (t) = 0, for the description of internal solitary waves in shallow seas [36], ut + αuux + βu2 ux + δuxxx = 0.

(4)

PT

ED

M

Ref. [23] has established a lattice BGK model for Eq. (4), whose tunable parameters in ChapmanEnskog expansion of the local equilibrium distribution function are determined by the coefficient of Eq. (4). (III) When f0 (t) = −K0 (t), f1 (t) = 0, f2 (t) = K0 (t), f3 (t) = −4K1 (t), f4 (t) = −h(t) and C = 2, Eq. (1) is transformed into the variable-coefficient and nonisospectral KdV equation, or the h-t-KdV equation, for the interfacial waves in a two-layer liquid and Alfv´en waves in a collisionless plasma [37], (5)

CE

ut + K0 (t)(uxxx + 6uux ) + 4K1 (t)ux − h(t)(2u + xux ) = 0.

AC

, f1 (t) = − f (t) , f2 (t) = g(t), f3 (t) = f4 (t) = 0 and f (t)g(t) 6= 0, (IV) When f0 (t) = − k(t) 6 6 Eq. (1) is transformed into a type of variable-coefficient Gardner equations, ut + k(t)uux + f (t)u2 ux + g(t)uxxx = 0,

(6)

which have been investigated for its exhaustive classification by using the groups of equivalence transformations [38]. Further developments of the exhaustive classification for Eq. (6) have been studied in Ref. [39]. The remaining part of this paper will be structured as follows. Lattice Boltzmann model for Eq. (1) will be derived in Section 2. The effects of the free parameters induced in our lattice Boltzmann model will be discussed in Section 3. Detailed numerical simulations on preceding 4

ACCEPTED MANUSCRIPT

cases will be performed in order to examine the accuracy and the stability of our model in Section 4. Finally, conclusions will be summarized in Section 5.

2. A lattice Boltzmann model for Eq. (1) Define S(x, t) = −f4 (t)(C − 1)u as a source item, Eq. (1) has been transformed into a concise modality as follows, (7)

CR IP T

ut − 6f0 (t)uux − 6f1 (t)u2 ux + f2 (t)uxxx − f3 (t)ux + f4 (t)(u + xux ) = S(x, t).

AN US

Without losing generality, we adopt the Bhatnagar-Gross-Krook (BGK) collision operator [16] with the D1Qb velocity model which uses 1-dimensional space and b velocity directions. In this − paper, we use the five velocity model (b = 5) which can be defined as → e = {e0 , e1 , e2 , e3 , e4 } = {0, c, −c, 2c, −2c}, where c is a scale factor constant, which represents the lattice velocity and defined as ∆x , (8) c= ∆t where ∆x and ∆t are represented the lattice space and time step, respectively. The evolution law of the particle distribution function for Eq. (1) can be written as  1 fα (x, t) − fα(eq) (x, t) + Υα (x, t) , (9) τ where fα (x, t) and fαeq (x, t) (α = 0, 1, . . . , 4) are the local particle distribution function and the local equilibrium distribution function respectively, τ is the single relaxation time, and Υα (x, t) is appended as a compensation in the process for recovering the macroscopic equations. We assume that the compensate function Υα (x, t) in the following form,   Υα (x, t) =  hα (x, t) + sα (x, t) , (10)

PT

ED

M

fα (x + eα ∆t, t + ∆t) − fα (x, t) = −

AC

CE

where hα (x, t) is the collisionless function, sα (x, t) representes the component of the source item S(x, t) impacting on the evolution equation at different velocity direction and  is the Knudsen number given by the ratio of the particle mean free path to the characteristic length of fluids. Here, we assume  is small enough and set  = ∆t. For the sake of simplicity and without loss of generality, we assume sα at different velocity direction has the same value. Hence it can be defined by

f4 (t)(C − 1)u S(x, t) =− , (11) b 5 where b is the number of different velocity direction. Applying Taylor expansion to Eq. (9), and retaining the terms up to O(∆t4 ), we have sα (x, t) =

∆t2 ∆t3 (∂t + eα ∂x )2 fα + (∂t + eα ∂x )3 fα + O(∆t4 ) 2 6  1 (eq) = − fα (x, t) − fα (x, t) + (hα + sα ). τ

∆t(∂t + eα ∂x )fα +

5

(12)

ACCEPTED MANUSCRIPT

Applying multi-scale Chapman-Enskog expansion [40] up to the two-order in time, we have ∂t = ∂t1 + 2 ∂t2 , fα =

∞ X

(13a)

n fα(n) = fα(0) + fα(1) + 2 fα(2) + 3 fα(3) + . . . ,

(13b)

n=0

sα = ∆t · ς1 = ς1 .

(13c)

O(0 ) : −

CR IP T

Substituting Eqs. (13a)-(13c) into Eq. (12), we can obtain series of differential equations for the first three orders of ,  1  (0) fα − fα(eq) = 0, i.e. fα(0) = fα(eq) , τ

(14a)

M

AN US

  1 O(1 ) : ∂x eα fα(0) = − fα(1) + hα , (14b) τ  1  1 O(2 ) : (∂t1 + eα ∂x )fα(0) + ∂x2 e2α fα(0) = − fα(2) + ς1 , (14c) 2 τ     1 2  2 (1)  1 3  3 (0)  2 (0) O(3 ) : ∂t2 fα(0) + ∂t1 fα(1) + ∂x eα fα(2) + ∂x,t e f + ∂x eα fα + ∂x eα fα α α 1 2 6 (14d) 1 (3) = − fα . τ Similar with general LBM, we define the macroscopic physical quantity u as distribution function X u= fα (x, t). (15)

ED

α

(eq)

PT

Due to the conservation law of local mass, fα and fα X X fα(eq) (x, t). fα (x, t) = u= α

(16)

α

CE

From Eq. (14a), we have X X fα(0) (x, t) = u, fα(n) (x, t) = 0, n > 0. α

satisfy

(17)

α

AC

To recover Eq. (7), some constraints on equilibrium distribution function and collisionless function are imposed as follows, which are enlightened by [27, 41], X eα fα(0) = 0 , (18a) α

X

e2α fα(0) = 0 ,

(18b)

α

X α

e3α fα(0) =

(τ 2

f2 (t)u , − τ + 61 )∆t2

(18c) 6

ACCEPTED MANUSCRIPT

X

hα = 0 ,

(18d)

α

X α

X

eα hα = −

3f0 (t)u2 + 2f1 (t)u3 + f3 (t)u − f4 (t)xu , τ ∆t

(18e)

e2α hα = 0.

(18f)

Through Eq. (14b) and noting Eqs. (18a)-(18f), we have hX i X X (1) 2 (0) e α fα = τ eα hα − ∂x e α fα α

α

α

3f0 (t)u2 + 2f1 (t)u3 + f3 (t)u − f4 (t)xu =− , ∆t

X

e2α fα(1) = τ

hX α

α

e2α hα − ∂x

X

e3α fα(0)

α

τ f2 (t)ux =− 2 . (τ − τ + 61 )∆t2

i

AN US

and

CR IP T

α

(19)

(20)

Similarly, through Eq. (14c) and noting Eqs. (18a)-(18f), we have

τ ( 1 − τ )f2 (t)uxx = − 22 . (τ − τ + 16 )∆t2

M

α

hX i 1 hX n hX i o i X e2α fα(1) ) + ∂x2 eα fα(2) = −τ ∂t1 eα fα(0) + ∂x e3α fα(0) − eα ς1 2 α α α α

ED

X

(21)

Summing Eq. (14c) over α and using Eqs. (15)-(19), we have f4 (t)(u + xux ) 6f0 (t)u + 6f1 (t)u2 + f3 (t) ux + = 5ς1 . ∆t ∆t

PT

∂t1 (u) −

(22)

CE

Similarly, summing Eq. (14d) over α and using Eqs. (15)-(21), we have ∂t2 (u) +

f2 (t)uxxx = 0. ∆t2

(23)

AC

The macroscopic equation Eq. (7) is exactly recovered by taking Eq. (22)×∆t+Eq. (23)×∆t2 : ut − 6f0 (t)uux − 6f1 (t)u2 ux + f2 (t)uxxx − f3 (t)ux + f4 (t)(u + xux ) = S(x, t).

(24)

Substituting the form of the source term S(x, t), the generalized Gardner equation with timedependent variable coefficients (1) can be obtained. In order to derive the equilibrium distribution function, we should introduce one new constraint other than Eqs. (16)-(18c). Considering the symmetry of the lattice velocity and regarding

7

ACCEPTED MANUSCRIPT

velocity-0 as the dominant term, we assume that the equilibrium distribution function satisfies the following constraint a small free parameter σ, (0) σf0

=

4 X

(0)

(0)

(0)

(0)

fα(0) = f1 + f2 + f3 + f4 ,

(25)

α=1

CE

PT

ED

M

AN US

CR IP T

where σ is introduced to adjust equilibrium distribution function. Coupling Eq. (25) and Eqs. (16)-(18c), we can derive the equilibrium distribution function as follows,  1  u, α=0 ,   1 + σ      2σ κ    − u, α=1,   3(1 + σ) 6       2σ κ (0) + u, α=2 , fα = (26) 3(1 + σ) 6        σ κ  − − u, α=3 ,   6(1 + σ) 12        σ κ  − + u, α=4 , 6(1 + σ) 12   where κ = f2 (t)/ (τ 2 − τ + 61 )c3 ∆t2 . Coupling Eqs. (18d)-(18f) and assuming the collisionless function in the directions of velocity1 and velocity-2 are the same, we can derived the collisionless function as follows,  1   H, α=0 ,   2     1   − H , α=1 ,     3 1 (27) hα = − H , α=2 ,  3    1   H, α=3 ,    3      − 1 H , α=4 , 6

AC

where H = [3f0 (t)u2 + 2f1 (t)u3 + f3 (t)u − f4 (t)xu] / (τ c∆t). Considering the time-dependent variable coefficients in Eq. (1), we think the single relaxation time τ should not be a constant with fixed value in the specific numerical simulations. Recalled Eq. (26), we can simply set the parameter κ as a given constant with fixed value. Therefore, due to the stability of the equation which requires τ > 0.5 [42], we can obtain the following relation, r 1 ∆t 1 τ (t) = + + f2 (t). (28) 2 12 κ∆x3 Considering all parameters are real, one should guarantee 1 ∆t 12∆tf2 (t) + f (t) > 0 ⇒ κ > − . 2 12 κ∆x3 ∆x3 8

(29)

ACCEPTED MANUSCRIPT

3. Discussions for the free parameters

CR IP T

We can analyze the order of the above inequality (29). In our simulations, the lattice space step ∆x = 10−2 ∼ 10−4 and the lattice time step ∆t = 10−2 ∼ 10−6 . We can easily obtain ∆t/∆x3 = 1 ∼ 1010 . If the function f2 (t) is negative, the parameter κ may be so large that exceed the limit of calculations. Hence, the time-dependent variable coefficient f2 (t) should be positive, and the parameter κ can be an arbitrary positive constant, which no longer consider the limit (29).

max |u(xj , t) − u∗ (xj , t)| ,

(30b)

PT

E∞ =

ED

M

AN US

As mentioned above, the constant σ is introduced to provide compensatory constraint of the equilibrium distribution functions to adjust them into appropriate relations. And the parameter κ is used to adjust the single relaxation time τ . Hence, the parameters σ and κ are all adjustment factors, which are important for the numerical simulations. Besides, the lattice space grid ∆x and lattice velocity c are also important for the numerical computation. Hereby, in this section, we will discuss the effects of κ, σ, ∆x and c, in our current lattice Boltzmann model, which can guide us to choose suitable values of the free parameters and obtain nice numerical simulation results. The validation of our work is confirmed by comparing the analytical solutions and the numerical results. The root-mean-square error E2 , the max error E∞ , and the global relative error GRE are used to measure the accuracy of our model [43]: v u N X 1u (30a) [u(xj , t) − u∗ (xj , t)]2 , E2 = t N j=1

(30c)

j=1,2,...,N

|u(xj , t) − u∗ (xj , t)| , PN ∗ (x , t)| |u j j=1

j=1

CE

GRE =

PN

AC

where u(xj , t) and u∗ (xj , t) are the numerical result and the analytical solution, respectively, and N is the number of grid points. A practical method to calculate the rate of convergence for a discretization method is to implement the following formula [43]: log Enew − log Eold , p≈ (31) log Hnew − log Hold where p is the convergent order, E represents the error, H represents the sequence, and the subscript new and old represent the new and old statuses, respectively. n Denote fi,j = fi (xj , tn ), unj = u(xj , tn ), xj = xmin + j∆x, tn = n∆t, where (i = 0, 1, 2, 3, 4; j = 1, 2, . . . , N ), n denotes the n-th layer time, j is spatial grid, and xmin is the left border of 9

ACCEPTED MANUSCRIPT

calculation section. Then, we can rewrite the evolution law (9) by classical finite difference notation: 1 n 1 n,(0) n = (1 − )fi,j + fi,j + ∆t(hi,j + si,j ), gi,j τ τ

(32a)

n+1 n fi,j = gi,j−E i

(32b)

CR IP T

where hi,j = hi (xj , tn ) in the expression (27), si,j = si (xj , tn ) in the expression (11), and E = {0, 1, −1, 2, −2}. We name Eq. (32a) as the collision process, and Eq. (32b) as the stream process.

Set Value

AN US

Calculate Coefficients

Initialize

IterNum=1 Calculate Coefficients

M

Collision

IterNum=IterNum+1

ED

Stream

CE

PT

Boundary Processing

If IterNum>=TMax

No

Yes Output Data

AC

Fig. 1. The flowchart of our numerical simulation procedures.

The non-equilibrium extrapolation scheme proposed in Ref. [44] is used for boundary condition. We can write specific formulas of boundary treatment by classical finite difference notation

10

ACCEPTED MANUSCRIPT

as follows, n,(0)

n,(0)

n+1 n+1 f3,2 = f3,2 + (f3,3 − f3,3 ), n,(0)

n,(0)

n+1 n+1 f4,N −1 = f4,N −1 + (f4,N −2 − f4,N −2 ),

do I = 0, 1, 2, 3, 4 n,(0)

(33)

n,(0)

n+1 n+1 fI,1 = fI,1 + (fI,2 − fI,2 ), n,(0)

n,(0)

n+1 n+1 fI,N = fI,N + (fI,N −1 − fI,N −1 ),

CR IP T

end do.

AN US

Furthermore, we give the flowchart of our numerical simulation procedures, as shown in Fig. 1. We think one can conduct the numerical simulation based on our LBM via this flowchart easily. We note that all parameters should be recalculated at the beginning of every time steps, due to the time-dependent variable coefficients. Without loss of generality, we consider one-soliton solutions for a KdV equation using our LBM model. Let f0 (t) = −1, f2 (t) = 1, f1 (t) = f3 (t) = f4 (t) = 0 and C = 1, Eq. (1) is transformed into the most common form of the KdV equation as follows, ut + 6uux + uxxx = 0.

(34)

u(x, t) =

et+x . 5(et + 0.1ex )2

M

An analytical one-soliton solution for Eq. (34) is given as,

ED

3.1. Effects of the parameter κ

AC

CE

PT

Firstly, we investigate the effects of the the parameter κ. Let σ = 0.015, ∆x = 0.01 and c = 100, we select the case of κ = 0.1, κ = 1.0 and κ = 1.5 for comparison. We give the E2 , E∞ and GRE of numerical result and one-soliton solutions for Eq. (34) at different times with different values of κ, which is presented in Tab. 1-3. We find that the E2 , E∞ and GRE of κ = 1.5 is less than the case of κ = 0.1 and κ = 1.0 at t = 1. Nevertheless, with the evolution of Eq. (34) continuing, the numerical results with κ = 1.5 diverge from the analytical solutions after a period of time (at t = 2), while the numerical results with κ = 0.1 and κ = 1.0 can couple with the analytical solutions for a longer time. We can also observe that the numerical results with κ = 0.1 less than the analytical solutions. Hence, we make a conjecture and consider the free parameter 1/κ as the damping. When the damping κ is selected as a large value, i.e. 1/κ as a small value, the numerical perturbations, e.g., truncation error, in the numerical simulation can not be restrained and get amplification as the evolution continuing, which lead to the divergency results. When the damping κ is selected as a small value, i.e. 1/κ as a large value, the inhibiting effects are strong to make the numerical results less than the analytical solutions.

11

ACCEPTED MANUSCRIPT

κ = 0.1

t=0.5

t=1.0

t=1.5

t=2.0

E2 E∞ GRE

3.1166e-005 2.1621e-002 4.1105e-002

1.1686e-004 3.7771e-002 7.9564e-002

2.4255e-004 5.2222e-002 1.1459e-001

3.9611e-004 6.5545e-002 1.4956e-001

Table 1. The error of numerical result and one-soliton solution for Eq. (34) at different times with κ = 0.1.

t=0.5

t=1.0

t=1.5

E2 E∞ GRE

4.9703e-007 2.4583e-003 5.1203e-003

1.8430e-006 4.3774e-003 1.0324e-002

3.7800e-006 6.1965e-003 1.4982e-002

t=2.0

CR IP T

κ = 1.0

6.1008e-006 7.7320e-003 1.9424e-002

Table 2. The error of numerical result and one-soliton solution for Eq. (34) at different times with κ = 1.0.

t=0.5

t=1.0

t=1.5

t=2.0

E2 E∞ GRE

2.2355e-007 1.6218e-003 3.4900e-003

8.3297e-007 2.9475e-003 6.9624e-003

1.9144e-004 4.4626e-002 1.1054e-001

? ? ?

AN US

κ = 1.5

M

Table 3. The error of numerical result and one-soliton solution for Eq. (34) at different times with κ = 1.5. (? represents the divergency.)

ED

3.2. Effects of the constant σ

AC

CE

PT

Secondly, we consider the constant σ induced in Eq. (26). Let κ = 1.0, ∆x = 0.01 and c = 100, we select the case of σ = −0.001, σ = 0.015 and σ = 0.5 for comparison. Detailed numerical simulations are shown in Fig. 2. We find that the numerical results with σ = 0.5 are less than the analytical solutions, while the numerical results with σ = −0.001 overtop the analytical solutions and have already oscillated in the middle part. Hence, we can consider σ as an equilibrium parameter, which can adjust the amplitude in the process of the evolution. In addition, the σ should be a positive value in our investigation.

12

ACCEPTED MANUSCRIPT

0.6

0.6 Analytical solution at t=2 Numerical result at t=2

0.5

Analytical solution at t=2 Numerical result at t=2

0.5

0.3

0.3

0.2

0.2

0.1

0.1

0

0

-0.1

-10

-5

0

5

10

15

-0.1

20

-10

-5

CR IP T

U

0.4

U

0.4

0

X

5

10

15

20

X

(a)

(b)

0.6 Analytical solution at t=2 Numerical result at t=2

0.5

U

0.3

0.2

0.1

0

-0.1

-10

-5

0

AN US

0.4

5

10

15

20

M

X

(c)

ED

Fig. 2. Numerical result and analytical solution at t = 1 with (a) σ = −0.001; (b) σ = 0.015; (c) σ = 0.5.

3.3. Effects of the lattice space step ∆x and the lattice velocity c

AC

CE

PT

Finally, we take the lattice space step ∆x and the lattice velocity c into account. We set κ = 1.0 and σ = 0.015. The errors and numerical order for different groups of ∆x and c are listed in Tab. 4-6.From the space-time evolution graphes of the numerical results, which are omitted here, we find that he numerical results of all these above cases in Tab. 4-6 match the analytical solutions well. Meanwhile, from Tab. 4-6, we can find that the convergence rate order of the approach of increasing the lattice velocity c is about 1.98 (based on the root-mean-square errors E2 ), while the approach of decreasing the lattice space step ∆x can hardly aggrandize the accuracy of the numerical results. Nevertheless, the stability of the numerical simulations may be reduced, which means that we should give new values of other free parameters to obtain nice numerical simulation results. ∆x = 0.02

∆x = 0.01

13

∆x = 0.005

ACCEPTED MANUSCRIPT

t=1

E2

Order

c=50 c=100 c=200

7.1745e-006 1.8228e-006 4.8907e-007

 1.9767 1.8980

E2

Order

7.2398e-006  1.8276e-006 1.9860 4.8401e-007 1.9168

E2

Order

7.2844e-006  1.8439e-006 1.9820 4.8659e-007 1.9220

Table 4. The E2 of numerical result and one-soliton solution for Eq. (34) and the convergence rates at t = 1.

c=50 c=100 c=200

∆x = 0.01 E∞ Order

 1.0631 0.9645

9.1508e-003  4.3587e-003 1.0700 2.2144e-003 0.9770

9.1235e-003 4.3664e-003 2.2376e-003

∆x = 0.005 E∞ Order

CR IP T

t=1

∆x = 0.02 E∞ Order

9.1739e-003  4.3748e-003 1.0683 2.2173e-003 0.9804

1.9628e-002 1.0276e-002 5.4035e-003

 0.9336 0.9273

∆x = 0.01 GRE Order

∆x = 0.005 GRE Order

1.9691e-002  1.0291e-002 0.9362 5.3900e-003 0.9330

1.9738e-002  1.0328e-002 0.9344 5.4037e-003 0.9345

M

∆x = 0.02 GRE Order

t=1 c=50 c=100 c=200

AN US

Table 5. The E∞ of numerical result and one-soliton solution for Eq. (34) and the convergence rates at t = 1.

3.4. Brief summary

ED

Table 6. The GRE of numerical result and one-soliton solution for Eq. (34) and the convergence rates at t = 1.

AC

CE

PT

From above discussion, we can summarize as follows: (1) The parameter κ is an adjustment factor of the single relaxation time τ , whose reciprocal effects the evolution of the NLEEs as the damping, which should be positive, i.e., κ > 0 ; (2) The constant σ is an adjustment factor of the equilibrium functions, which can adjust the amplitude in the process of the evolution and always be positive, i.e., σ > 0 ; (3) The lattice space step ∆x and the lattice velocity c can impact on the accuracy and the stability of the numerical simulation results. In addition, via plentiful computational examples for the preceding cases based on our current LBM, we can give the approximate value range of the parameters κ and σ, 0.0 < κ < 2.0,

0 < σ < 0.5,

which may be a guidance to choose the appropriate free parameters. Setting the initial value of the free parameters as arbitrary values in above scopes, adjusting the value of the free parameters, and comparing with the analytical solutions in the meantime 14

ACCEPTED MANUSCRIPT

continually, we can obtain the appropriate free parameters for Eq. (1) at this set of C, fi (t) (i = 0, 1, . . . , 4). After that, we can use our current LBM to figure out the evolutions and search some new solutions of Eq. (1) at this set with some other specific initial conditions, which have no analytical solution at present.

4. Numerical simulations

Example 4.1

The soliton solution of the constant-coefficient Gardner equation in Ref. [35]

is given by 6γ p √ , 2α + 4α2 − 18βγ cosh γ(x − γt)

AN US

u(x, t) =

CR IP T

In this section, following the guidance of the discussions in Section 3, we will select appropriate free parameters and present some numerical simulations for preceding cases, which implies that our current lattice Boltzmann model is a satisfactory and efficient algorithm.

(35)

where f0 (t) = − α3 , f1 (t) = β2 , f2 (t) = 1 and f3 (t) = f4 (t) = 0. In this simulation, we set α = 5, β = 0.5, γ = 2. Hence, Eq. (3) has the form as follows, (36)

M

ut + 10uux − 1.5u2 ux + uxxx = 0.

AC

CE

PT

ED

Some other parameters are κ = 1.1135, σ = 0.01571, ∆x = 0.02, c = 100 (i.e., the corresponding time step is ∆t = 0.0002) and the computation domain is fixed on I = [−10, 10]. We present the comparison between detailed numerical results and analytical solutions. Fig. 3 shows the two-dimensional visual comparisons at some different times. The space-time evolution graph of the numerical results is shown in Fig. 4. Tab. 7 gives the E2 , E∞ and GRE of the numerical results at some different times. All of them manifest that the numerical results agree with the analytical solutions well.

(a)

(b)

Fig. 3. Numerical result(a) and analytical solution(b) for the propagation of the soliton from t = 0 to t = 2.

15

ACCEPTED MANUSCRIPT

0.8 Numerical result at t=0.4

0.7

Analytical solution at t=0.4 Numerical result at t=1.2

0.6

Analytical solution at t=1.2 Numerical result at t=2.0

0.5

Analytical solution at t=2.0

U

0.4 0.3

0.1 0 −0.1

−10

−8

−6

−4

−2

0 X

2

4

CR IP T

0.2

6

8

10

Figs. 4. Numerical result and analytical solution at t = 0.4, 1.2, 2.0

E2 E∞ GRE

t=1.2

1.1848e-005 1.1215e-003 1.9623e-003

t=2.0

AN US

t=0.4

6.7327e-005 2.4707e-003 5.2572e-003

1.3272e-004 3.3827e-003 7.8095e-003

M

Table 7. Comparison of numerical result and soliton solution for the constant-coefficient Gardner equation at different times.

Example 4.2

ED

A nonpropagating one-soliton solution for nonisospectral and variable coefficient KdV equation with time varying nonvanishing boundary condition in Ref. [37] is given by

(37b)

f (t) = 8

(37c)

t

0

(37a)

0

0

  3 f (t) η K0 (t)η 2 + K0 (t)L(t) + K1 (t) dt, ϑ = − ηx, 2 2

CE

Z

PT

u(x, t) = L(t) + 2η 2 sech2 (ϑ), Z Z   t   t h(t)dt , 2h(t)dt , η = η(0)exp L(t) = L(0)exp

AC

where f0 (t) = −K0 (t), f1 (t) = 0, f2 (t) = K0 (t), f3 (t) = −4K1 (t), f4 (t) = −h(t) and C = 2. In this simulation, we set K0 (t) = 0.1, K1 (t) = 0 and h(t) = 3(1+t22t)(2+t2 ) . Hence, Eq. (5) has the form as follows, ut + 0.01(uxxx + 6uux ) −

2t (2u + xux ) = 0. 3(1 + t2 )(2 + t2 )

(38)

1 Using Eqs. (37a)-(37c) and setting L(0) = 30 , η(0) = 0.1, the solution of Eq. (38) can be calculated in the following form, " ) r  2 2/3 (  # 2+1 t +1 t t 5 3 u(x, t) = 0.02 2 sech2 −0.1 2 x + 0.0024t − 0.001697 arctan √ + . t +2 t +2 3 2

16

ACCEPTED MANUSCRIPT

M

AN US

CR IP T

Some other parameters are κ = 1.0998, σ = 0.05274, ∆x = 0.05, c = 100 (i.e., the corresponding time step is ∆t = 0.0005) and the computation domain is fixed on I = [−50, 50]. Detailed comparison between the numerical results and the analytical solutions can be found in Figs. 5, Figs. 6 and Tab. 8, which show good agreement. Furthermore, we recall Tab. 7, and compare it with Tab. 8. Clearly, we can find that the value in Tab. 8 are far less than them in Tab. 7. We can notice that the value of the variable coefficients for the nonlinear terms of Eq. (1) in Example 4.2 are much less than them in Example 4.1. We conjecture that the nonlinear terms generate the disturbance in the space-time evolution based on our LBM. Large values of the nonlinear terms cause large errors between the numerical results and the analytical solutions. This assumption is being researched in our current works, which will be presented in our future papers.

(a)

(b)

ED

Figs. 5. Numerical result(a) and analytical solution(b) for the propagation of the soliton from t = 0 to t = 5.

Numerical result at t=1 Analytical solution at t=1 Numerical result at t=2 Analytical solution at t=2 Numerical result at t=3 Analytical solution at t=3 Numerical result at t=4 Analytical solution at t=4 Numerical result at t=5 Analytical solution at t=5

PT

0.06 0.055

CE

0.05

AC

U

0.045 0.04

0.035 0.03 0.025 0.02

−50

−40

−30

−20

−10

0 X

10

20

30

40

Figs. 6. Numerical result and analytical solution at t = 1, 2, 3, 4, 5.

17

50

ACCEPTED MANUSCRIPT

E2 E∞ GRE

t=1.0

t=2.0

t=3.0

t=4.0

t=5.0

2.0973e-012 5.0900e-006 4.0772e-005

3.8743e-012 6.9100e-006 2.6913e-005

5.4865e-012 7.6600e-006 3.5734e-005

6.4729e-012 1.0370e-005 4.0676e-005

6.7293e-012 8.0800e-06 4.1870e-005

CR IP T

Table 8. Comparison of numerical result and nonpropagating one-soliton solution for the nonisospectral and variable coefficient KdV equation at different times.

Example 4.3

Let f0 (t) = 1, f1 (t) = −0.01, f2 (t) = 0.01, f3 (t) = t, f4 (t) = 0, C = 1. Using the solution (16) in Ref. [32], a one-soliton solution for Eq. (1) can be derived as follows, 2

ED

M

AN US

0.020202e0.01t+0.5t +x u(x, t) = − . 1.02041e0.02t + 2.0202e0.01t+0.5t2 +x + et2 +2x

(a)

(b)

Figs. 7. Numerical result(a) and analytical solution(b) for the propagation of the soliton from t = 0 to t = 3. −3

x 10

CE

0

PT

1

−1

U

AC

−2

−3

Analytical solution at t=1 Numerical result at t=1 Analytical solution at t=2 Numerical result at t=2 Analytical solution at t=3 Numerical result at t=3

−4

−5

−6 −20

−15

−10

−5

0 X

5

10

15

Figs. 8. Numerical result and analytical solution at t = 1.0, 2.0, 3.0.

18

20

ACCEPTED MANUSCRIPT

In this simulation, parameters are set as follows: κ = 0.1500, σ = 0.01523, ∆x = 0.02, c = 100 (i.e., the corresponding time step is ∆t = 0.0002), and the computation domain is fixed on I = [−20, 20]. Figs. 7, Figs. 8 and Tab. 9 present the two-dimensional visual comparisons at some different times, the space-time evolution graph of the numerical results, and the E2 , E∞ and GRE of the numerical results at some different times, respectively. All of them manifest that the numerical results agree with the analytical solutions well. t=2.0

6.1920e-013 3.5510e-006 5.5832e-004

9.6446e-012 1.4366e-005 2.2225e-003

t=3.0

CR IP T

E2 E∞ GRE

t=1.0

4.9164e-011 3.3299e-005 5.0255e-003

Table 9. Comparison of numerical result and one-soliton solution for the variables-coefficient Gardner equation at different times.

Example 4.4

M

AN US

Setting f0 (t) = 1, f1 (t) = −2, f2 (t) = 2, f3 (t) = 2, f4 (t) = 0, C = 1 and using the double-pole solution (19) for Eq. (1) in Ref. [32], we can derive the analytical solution as follows, by the assignment of each coefficient, i h 3t 3t x 3t x e 4 + 2 e 2 +x (t + 2x − 6) + 16e 4 + 2 − 2(t + 2x + 6) u(x, t) = 3t +x . 9t 3t x 3x e 2 (t2 + 4tx + 4x2 + 4) + 4e 4 + 2 (t + 2x − 2) + 2e3t+2x − 2e 4 + 2 (t + 2x + 2) + 8

AC

CE

PT

ED

Some other parameters are κ = 1.5370, σ = 0.1130, ∆x = 0.02, c = 100 (i.e., the corresponding time step is ∆t = 0.0002) and the computation domain is fixed on I = [−20, 20]. Detailed numerical results can be found in Figs. 9, Figs. 10 and Tab. 10. All of them manifest that the numerical results agree with the analytical solutions well.

(a)

(b)

Figs. 9. Numerical result(a) and analytical solution(b) for the propagation of the soliton from t = 0 to t = 1.

19

ACCEPTED MANUSCRIPT

1.5

1.5

1

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1

-1.5

-20

-15

-10

-5

0

5

10

15

-1.5

20

-20

-15

-10

CR IP T

Analytical solution at t=0.6 Numerical result at t=0.6

U

U

Analytical solution at t=0.2 Numerical result at t=0.2

-5

X

0

5

10

15

20

X

(a)

(b)

1.5 Analytical solution at t=1.0 Numerical result at t=1.0

1

0

-0.5

-1

-1.5

-20

-15

-10

-5

AN US

U

0.5

0

5

10

15

20

M

X

(c)

t=0.2

t=0.6

t=1.0

2.4733e-004 9.6198e-003 3.7464e-003

6.6804e-004 1.2548e-002 8.2515e-003

3.6825e-003 4.2899e-002 1.6670e-002

CE

PT

E2 E∞ GRE

ED

Figs. 10. Numerical result and analytical solution at (a) t = 0.2, (b) t = 0.6, (c) t = 1.0.

Table 10. Comparison of numerical result and double-pole solution for vc-Gardner equatio at different times.

AC

From the front four examples, it can be seen clearly that the numerical results have a commendable accuracy to the analytical solutions when we take the appropriate κ and σ, which implies that the present lattice Boltzmann model is a satisfactory and efficient algorithm. Therefore, it is expected that the lattice Boltzmann model can be used to figure out the evolutions and search some new solutions of these NLEEs which can be represented by Eq. (1) with some specific initial conditions.

5. Conclusions In this paper, by using the Chapman-Enskog expansion, a lattice Boltzmann model for the 20

ACCEPTED MANUSCRIPT

AN US

CR IP T

generalized Gardner equation with time-dependent variable coefficients, i.e., Eq. (1), have been derived through selecting equilibrium distribution function and the compensate function, appropriately. The D1Q5 velocity model has been used in numerical simulations with different forms of Eq. (1). Two free parameters κ and σ are introduced to adjust the single relaxation time τ and the equilibrium distribution functions, respectively. Effects and approximate value range of these two free parameters have been discussed in detail, as well as the selection approach. Effects of the lattice space step ∆x and the lattice velocity c are also discussed. Soliton solution for the constant-coefficient Gardner equation (36), the nonpropagating onesoliton solution for h-t-KdV equation (38), the one-soliton solution and the double-pole solution for variable-coefficient Gardner equation (1) have been simulated by taking appropriate free parameters, i.e., κ, σ, ∆x and c. It has been found that the numerical results agree well with the analytical solutions, which indicates that the current lattice Boltzmann model is a satisfactory and efficient algorithm. From the derivation of lattice Boltzmann model in Section 2 of this paper, as the derivation in Refs. [8, 18–26], we can find that the one dimensional NLEEs, which can be modeled by the lattice Boltzmann methods in these references, as well as Eq. (1), have the following specific forms, N

∂u X ∂Ψ(u) αk + = F (u), ∂t k=1 ∂xk

(39)

ED

M

where u is the wave-amplitude function of the scaled coordinate x and t, N is the highest order of the derivatives, αk are constants or known functions, Ψ(u) and F (u) are all polynomial functions of u. Function F (u) is called the source item, which can be zero. In another word, a NLEE with the form (39) may derive its lattice Boltzmann model. Whereas, some important NLEEs do not have the form (39), such as the Camassa-Holm equation in shallow water [45],

PT

ut + 2ux − uxxt + 3uux = 2ux uxx + uuxxx ,

(40)

AC

CE

which can not derive its own lattice Boltzmann model, in our current derivation techniques. Forms (39) limit the application of the lattice Boltzmann method for numerically simulating NLEEs. We will investigate this restriction in our future works. Furthermore, the physical meanings of the two free parameters κ and σ are not clear. We only consider the lattice Boltzmann model for a generalized Gardner equation with time-dependent variable coefficients, not investigating the space-dependent variable coefficients Gardner equation in this paper. Extension of our lattice Boltzmann model for some higher derivative NLEEs are not given in this paper. All of these inadequacies will be replenished in our further study.

Acknowledgements We express our sincere thanks to the Editors, Referees and members of our discussion group for their valuable suggestions. This work has been supported by the National Natural Science 21

ACCEPTED MANUSCRIPT

Foundation of China under Grant No.11272023, and by the Fundamental Research Funds for the Central Universities under Grant No.50100002016105010.

References [1] R. Grimshaw, D. Pelinovsky, E. Penlinovsky and T. Talipova, Wave group dynamics in weakly nonlinear long-wave models, Phys. D 159, 35 (2001).

CR IP T

[2] A. V. Slyunyaev, Dynamics of localized waves with large amplitude in a weakly dispersive medium with a quadratic and positive cubic nonlinearity, J. Exp. Theor. Phys. 92, 529 (2001).

AN US

[3] R. Grimshaw, E. Penlinovsky, T. Talipova and A. Kurkin, Simulation of the transformation of internal solitary waves on oceanic shelves, J. Phys. Oceanogr. 34, 2774 (2004). [4] W. Q. Hu, Y. T. Gao, C. Zhao, Y. J. Feng and C. Q. Su, Oscillations in the interactions among multiple solitons in an optical fibre, Z. Naturforsch. A 71, 1079 (2016). [5] B. B. Kadomtsev and V. I. Petviashvili, On the stability of solitary waves in weakly dispersing media, Sov. Phys. Dokl. 15, 539 (1970).

M

[6] H. Zhang, A note on exact complex travelling wave solutions for (2 + 1)-dimensional B-type Kadomtsev-Petviashvili equation, Appl. Math. Comput. 216, 2771 (2010).

PT

ED

[7] W. Q. Hu, Y. T. Gao, S. L. Jia, Q. M. Huang and Z. Z. Lan, Periodic wave, breather wave and travelling wave solutions of a (2 + 1)-dimensional B-type Kadomtsev-Petviashvili equation in fluids or plasmas, Eur. Phys. J. Plus 131, 390 (2016).

CE

[8] H. L. Lai and C. F. Ma, Lattice Boltzmann model for generalized nonlinear wave equations, Phys. Rev. E 84, 046708 (2011).

AC

[9] R. Hirota, Exact solution of the Kortewegde Vries equation for multiple collisions of solitons, Phys. Rev. Lett. 27, 1192 (1971). [10] V. G. Bagrov and B. F. Samsonov, Darboux transformation, factorization, and supersymmetry in one-dimensional quantum mechanics, Theor. Math. Phys. 104, 1051 (1995). [11] A. R. Mitchell and D. F. Griffiths, The finite difference method in partial differential equations, John Wiley (1980). [12] C. Johnson, Numerical solution of partial differential equations by the finite element method, Courier Corporation (2012).

22

ACCEPTED MANUSCRIPT

[13] D. G. Dommermuth, D. K. P. Yue, A high-order spectral method for the study of nonlinear gravity waves, J. Fluid Mech. 1987, 184: 267-288. [14] S. Y. Chen and G. D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30, 329 (1998). [15] Y. H. Qian, S. Succi and S. A. Orszag, Recent advances in lattice Boltzmann computing, Annu. Rev. Comput. Phys. 3, 195 (1995).

CR IP T

[16] R. Benzi, S. Succi and M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222, 145 (1992). [17] H. D. Chen, S. Kandasamy, S. Orszag, R. Shock, S. Succi and V. Yakhot, Extended Boltzmann kinetic equation for turbulent flows, Science 301, 633 (2003).

AN US

[18] I. Ginzburg, Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation, Adv. Water Resour. 28, 1171 (2005). [19] B. C. Shi and Z. L. Guo, Lattice Boltzmann model for nonlinear convection-diffusion equations, Phys. Rev. E 79, 016701 (2009).

M

[20] J. Y. Zhang and G. W. Yan, Lattice Boltzmann method for one and two-dimensional Burgers equation, Phys. A 387, 4771 (2008). [21] G. W. Yan, A lattice Boltzmann equation for waves, J. Comput. Phys. 161, 61 (2000).

ED

[22] H. B. Li, P. H. Huang, M. R. Liu and L. J. Kong, Simulation of the MKdV equation with lattice Boltzmann method, Acta. Phys. Sini. 50, 837 (2001).

PT

[23] C. F. Ma, A lattice BGK model for simulating solitary waves of the combined KdV-MKdV equation, Int. J. Mod. Phys. B 25, 589 (2011).

CE

[24] C. Y. Zhang, H. L. Tan, M. R. Liu and L. J. Kong, A lattice Boltzmann model and simulation of KdV-Burgers equation, Commun. Theor. Phys. 42, 281 (2004).

AC

[25] H. M. Wang, A lattice Boltzmann model for the ion- and electron-acoustic solitary waves in beam-plasma system, Appl. Math. Compu. 279, 62 (2016). [26] B. C. Shi, N. Z. He and Z. L. Guo, Lattice Boltzmann Model for High-Order Nonlinear Partial Differential Equations, arXiv preprint arXiv: 0907.1720 (2009). [27] Z. H. Chai, B. C. Shi and L. Zheng, A unified lattice Boltzmann model for some nonlinear partial differential equations, Chaos, Solitons and Fractals 36, 874 (2008). [28] Y. Nakamura and I. Tsukabayashi, Observation of modified Kortewegde Vries solitons in a multicomponent plasma with negative ions, Phys. Rev. Lett. 52, 2356 (1984). 23

ACCEPTED MANUSCRIPT

[29] R. A. Szoeke, An effect of the thermobaric nonlinearity of the equation of state: A mechanism for sustaining solitary Rossby waves, J. Phys. Oceanogr. 34, 2042 (2004). [30] H. Triki, T. R. Taha and A. M. Wazwaz, Solitary wave solutions for a generalized KdVCmKdV equation with variable coefficients, Math. Comput. Simul. 17, 80 (1867).

CR IP T

[31] G. Q. Meng, Y. T. Gao, X. Yu, Y. J. Shen and Y. Qin, Painlev´e analysis, Lax pair, B¨acklund transformation and multi-soliton solutions for a generalized variable-coefficient KdV-mKdV equation in fluids and plasmas, Phys. Scr. 85, 055010 (2012). [32] P. Wang, B. Tian, W. J. Liu, Y. Jiang and Y. S. Xue, Interactions of breathers and solitons of a generalized variable-coefficient Korteweg-de Vries-modified Korteweg-de Vries equation with symbolic computation, Eur. Phys. J. D 66, 233 (2012).

AN US

[33] X. Yu, Y. T. Gao and Z. Y. Sun, Solitonic propagation and interaction for a generalized variable-coefficient forced Korteweg-de Vries equation in fluids, Phys. Rev. E 83, 056601 (2011). [34] T. Brugarino, Painlev´e property, auto-B¨acklund transformation, Lax pairs, and reduction to the standard form for the Korteweg-De Vries equation with nonuniformities, J. Math. Phys. 30, 1013 (1989).

M

[35] A. M. Wazwaz, New solitons and kink solutions for the Gardner equation, Commun. Non. Sci. Numer. Simul. 12, 1395 (2007).

ED

[36] K. W. Chow, R. H. J. Grimshaw and E. Ding, Interactions of breathers and solitons in the extended Korteweg-de Vries equation, Wave Motion 43, 158 (2005).

PT

[37] W. L. Chan and K. S. Li, Nonpropagating solitons of the variable coefficient and nonisospectral KortewegCde Vries equation, J. Math. Phys. 30, 2521 (1989).

AC

CE

[38] O. Vaneeva, O. Kuriksha, C. Sophocleous, Enhanced group classification of Gardner equations with time-dependent coefficients, Commun. Nonlinear Sci. Numer. Simul. 22, 1243 (2015). [39] R. de la Rosa, M. L. Gandarias, M. S. Bruz´on, Equivalence transformations and conservation laws for a generalized variable-coefficient Gardner equation, Commun. Nonlinear Sci. Numer. Simul. 40, 71 (2016). [40] S. Chapman and T. G. Cowling, The mathematical theory of non-uniform gases 3rd ed, Cambridge University Press (1970). [41] F. Liu, W. P. Shi and F. F. Wu, A lattice Boltzmann model for the generalized Boussinesq equation, Appl. Math. Comp. 274, 331 (2016). 24

ACCEPTED MANUSCRIPT

[42] J. D. Sterling, S. Y. Chen, Stability analysis of lattice Boltzmann methods, J. Comp. Phys. 123, 196 (1996). [43] R. L. Burden and J. D. Faires, Numerical Analysis 7th ed, Brooks/Cole (2001). [44] Z. L. Guo, C. G. Zheng and B. C. Shi, Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method, Chin.Phys. 11, 366 (2002).

AC

CE

PT

ED

M

AN US

CR IP T

[45] R. Camassa and D. D. Holm, An Integrable Shallow Water Equation with Peaked Solitons, Phys. Rev. Lett. 71, 1661 (1993).

25