Sub-ODE method and solitary wave solutions for higher order nonlinear Schrödinger equation

Physics Letters A 363 (2007) 96–101 www.elsevier.com/locate/pla

Sub-ODE method and solitary wave solutions for higher order nonlinear Schrödinger equation Mingliang Wang a,b,∗ , Xiangzheng Li a , Jinliang Zhang a a College of Science, Henan University of Science & Technology, Luoyang 471003, PR China b Department of Mathematics, Lanzhou University, Lanzhou 730000, PR China

Received 11 April 2006; received in revised form 15 October 2006; accepted 27 October 2006 Available online 7 November 2006 Communicated by A.R. Bishop

Abstract With the aid of an ordinary differential equation (ODE) involving an arbitrary positive power of dependent variable and its various positive solutions, three types of solitary wave solution of the higher order nonlinear Schrödinger equation (NLSE) with non-Kerr terms have been found out, which are the bell type solitary waves, the kink type solitary waves and the algebraic solitary waves, provided that the coefficients of the higher order NLSE satisfy certain constraint conditions. © 2006 Elsevier B.V. All rights reserved. Keywords: Higher order nonlinear Schrödinger equation; Sub-ODE; Bell type solitary waves; Kink type solitary waves; Algebraic solitary waves

1. Introduction In this Letter, we consider the higher order nonlinear Schrödinger equation (NLSE) in the form   iuz + utt + 2|u|2 u + iαuttt + iβ |u|2 u t   + iγ |u|4 u t + δ|u|4 u = 0 (1.1) where α, β, γ and δ are constants, which was firstly proposed and investigated by Radhakrishnan, Kundu, Lakshmanan (RKL) [1], where u = u(z, t) is a complex function that represents a normalized complex amplitude of the pulse envelope in nonlinear optical fiber, the variable z is interpreted as the normalized propagation distance, t —retarded time. If the last four terms are omitted, Eq. (1.1) becomes the standard NLSE, which describes the propagation of picosecond pulses in optical fiber [2–4], and which is integrable. However, RKL model Eq. (1.1) describes the propagation of femtosecond optical pulses, the last four terms representing non-Kerr nonlinearity effects become important and non* Corresponding author.

E-mail address: [email protected] (M. Wang). 0375-9601/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2006.10.077

negligible, and RKL model Eq. (1.1) as well as the other higher order NLSE models proposed in the literature [5–11] are not integrable, because these models involve not only the higher order nonlinear terms but also the higher dispersion terms. The analytical and numerical solutions of these models have been actively investigated by many authors [5–12]. In recent years, various powerful methods have been presented to derive exact solutions of nonlinear partial differential equations (PDE) in mathematical physics, such as the inverse scattering transform [13], the Hirota’s bilinear operators [14], the homogeneous balance method [15–18], the Jacobi elliptic function expansion [19,20], the tanh-function expansion and its various extensions [21–26], the auxiliary (or subsidiary) ordinary differential equation method [27], the F-expansion method [28–32] and so on. The last three methods mentioned above belong to a class of method called sub-ODE method [33], the sub-ODEs which were often used are Riccati equation, Jacobi elliptic equation and projective Riccati equations, etc. In this Letter, a sub-ODE involving an arbitrary positive power of dependent variable is introduced as follows F  2 = AF 2 + BF 2+p + CF 2+2p ,

p > 0,

(1.2)

M. Wang et al. / Physics Letters A 363 (2007) 96–101

where A, B and C are constants, sub-ODE (1.2) admits the following positive solutions:  1/p 1 (a) F (ξ ) = (1.3) , σ < 1, √ cosh p Aξ − σ provided that   A > 0, B = 2σ A, C = σ 2 − 1 A; (1.4) (b)

1/p   A 1 1 p√ Aξ , F (ξ ) = ± tanh C 2 2 2

provided that A > 0,

 (c)

√ B = −2 AC;

C > 0,

F (ξ ) =

1 ξ2 + σ

1/p ,

σ > 0,

(1.7)

provided that 4 4 , C = − 2 σ. (1.8) p2 p The results from (1.2) to (1.8) can be examined directly. It should be noted that ODE (1.2) could admit other solutions, for example, the sinusoidal type solutions, the negative solutions for odd integer p, and so on. But for the sake of simplicity we neglect the cases here. The objective of the Letter is that with the aid of subODE (1.2) and its various solutions we look for solitary wave solutions for RKL model Eq. (1.1). In general, the presence of solitary wave solution depends on the coefficients α, β, γ and δ appearing in Eq. (1.1), so we shall also give the constraint condition of the coefficients, under which the solitary wave solutions of Eq. (1.1) can be found out. The rest of the Letter are organized as follows: in Section 2, by using traveling wave reduction Eq. (1.1) is reduced into a nonlinear ODE, which can be thought of as a special case of sub-ODE (1.2), provided that the coefficients of Eq. (1.1) satisfy certain constraint conditions; from Section 3 to Section 5 the bell type solitary wave solutions, the kink type solitary wave solutions and the algebraic solitary wave solutions as well as their corresponding constraint conditions are given respectively; in Section 6, a brief conclusion and a generalization are made for the results obtained in this Letter. A = 0,

B=

2. Traveling wave reduction of Eq. (1.1) Since u = u(z, t) is a complex function we suppose that u(z, t) = φ(ξ )eiη ,

ξ = kt + ωz,

η = λt + μz,

  Im: αk 3 φ  + ω + 2λk − 3αλ2 k φ  + 3βkφ 2 φ  + 5γ kφ 4 φ  = 0.

(2.1)

where k, ω, λ and μ are constants to be determined later, φ = φ(ξ ) > 0. Substituting (2.1) into Eq. (1.1), removing the common factor eiη we have the overdetermined ODEs for φ(ξ ) Re: (3αλ − 1)k 2 φ    = αλ3 − λ2 − μ φ + (2 − λβ)φ 3 + (δ − λγ )φ 5 , (2.2)

(2.3)

There are two cases to discuss. 2.1. When 3αλ − 1 = 0, Eq. (2.2) can be rewritten as φ  =

(1.5)

(1.6)

97

αλ3 − λ2 − μ 2 − λβ φ+ φ3 2 (3αλ − 1)k (3αλ − 1)k 2 δ − λγ + φ5. (3αλ − 1)k 2

(2.4)

Substituting (2.4) into (2.3) and setting the coefficients of φ  , φ 2 φ  and φ 4 φ  to zero respectively, yields α(αλ3 − λ2 − μ)k + ω + 2λk − 3αλ2 k = 0, 3αλ − 1 3α(2 − λβ)k + 3βk = 0, 3αλ − 1 5α(δ − λγ )k + 5γ k = 0. 3αλ − 1 From the equations above yields 1 1 − , 2α β     1 1 1 2 1 ω = 3α −2 − − 2α β 2α β

λ=

−

(2.5)

1 1 α[α( 2α − β1 )3 − ( 2α − β1 )2 − μ]

1 2

−

k,

3α β

(2.6)

and the constraint condition 2γ = βδ.

(2.7)

Substituting (2.5) into Eq. (2.4), multiplying the resulted equation by φ  , integrating it with respect to ξ once and taking the constant of integration to zero, yields φ  2 = A1 φ 2 + B1 φ 4 + C1 φ 6 ,

(2.8)

where A1 =

1 1 α( 2α − β1 )3 − ( 2α − β1 )2 − μ

( 12 −

3α 2 β )k

,

β γ , C1 = − , 2 2αk 3αk 2 μ and k are still to be determined later.

B1 = −

(2.9)

2.2. When 3αλ − 1 = 0, i.e., 1 , 3α then Eq. (2.2) becomes   3 αλ − λ2 − μ φ + (2 − λβ)φ 3 + (δ − λγ )φ 5 = 0.

λ=

(2.10)

Setting the coefficients of φ, φ 3 and φ 5 to zero respectively, from the resulted equations and using (2.10) yields μ=−

2 , 27α 2

(2.11)

98

M. Wang et al. / Physics Letters A 363 (2007) 96–101

and the constraint condition β = 6α

and 2γ = βδ,

which implies γ = 3αδ.

(2.12)

Substituting (2.10) and (2.12) into Eq. (2.3) and integrating it with respect to ξ once and taking the constant of integration to zero, after multiplying the resulted equation by φ  and integrating it again, yields φ  2 = A2 φ 2 + B2 φ 4 + C2 φ 6 ,

(2.13)

where the constant of integration is zero,   ω 1 3 + 2 2 , A2 = − B2 = − 2 , 3 αk 3α k k

k—arbitrary nonzero constant. C2 = −

δ , k2 (2.14)

ω and k are still to be determined later. Note that both ODE (2.8) and ODE (2.13) can be thought of as a special case of sub-ODE (1.2) when p = 2, so it is possible that we use the results described in Section 1 to solve (2.8) and (2.13), meanwhile to determine the rest constants k, ω and μ, and finally to obtain various exact solutions of RKL model Eq. (1.1). 3. Bell type solitary wave solutions 3.1. According to the solution (1.3) of Eq. (1.2) with coefficients (1.4) (p = 2), Eq. (2.8) admits a solution in the form 1/2  1 , σ < 1, φ(ξ ) = (3.1) √ cosh 2 A1 ξ − σ provided that 1 1 α( 2α − β1 )3 − ( 2α − β1 )2 − μ

( 12 −

3α 2 β )k

= A1 > 0,

 1 γ  1 β − 2 = σ 2 − 1 A1 . = 2σ A1 , 2k 2 α 3k α From (3.2) we derive three kind results which are

−

where σ1 is expressed by (3.3),       1 1 2 1 β 1 −2 − − + ξ = kt + 3α kz, 2α β α β 4σ1   1 1 − t η= 2α β       1 β 1 1 1 3 1 2 3 + − +α − − z, + − 4σ1 8σ1 α 2α β 2α β

(3.2)

3.1.2. When 2γ = βδ, αβ > 0 and αγ < 0, which implies δ < 0, we obtain   2 1 1 δ = σ2 (σ2 < −1), σ = δ− 1+ (3.6) 3 3 1 β A1 = − (3.7) > 0, 4σ2 k 2 α     1 β 1 1 3 1 2 1 3 + − , +α − − μ=− (3.8) 4σ2 8σ2 α 2α β 2α β k—arbitrary nonzero constant. Substituting (3.1) with (3.6)–(3.8) as well as (2.5)–(2.6) into (2.1), we have the bell type solitary wave solution of RKL model Eq. (1.1)  1/2 1 eiη , u2 (z, t) = 1 β cosh 2 − 4σ k 2 α ξ − σ2 2

where σ2 is expressed by (3.6),       1 1 2 1 β 1 −2 − − + kz, ξ = kt + 3α 2α β 2α β 4σ2   1 1 − t η= 2α β       1 β 1 1 3 1 2 1 3 + − +α − − z, + − 4σ2 8σ2 α 2α β 2α β k—arbitrary nonzero constant.

3.1.1. When 2γ = βδ, αβ > 0 and αγ > 0, which implies δ > 0, we obtain   2 1 1 δ = σ1 (−1 < σ1 < 0), σ = δ− 1+ (3.3) 3 3 1 β > 0, A1 = − (3.4) 4σ1 k 2 α     1 β 1 1 1 3 1 2 3 + − , +α − − μ=− (3.5) 4σ1 8σ1 α 2α β 2α β k—arbitrary nonzero constant. Substituting (3.1) with (3.3)–(3.5) as well as (2.5)–(2.6) into (2.1), we have the bell type solitary wave solution of RKL model Eq. (1.1)  1/2 1 eiη , u1 (z, t) = 1 β cosh 2 − 4σ k 2 α ξ − σ1 1

3.1.3. When 2γ = βδ, αβ < 0 and αγ > 0, which implies δ < 0, we obtain   2 1 1 σ = δ+ 1+ (3.9) δ = σ3 (0 < σ3 < 1), 3 3 1 β A1 = − (3.10) > 0, 4σ3 k 2 α     1 β 1 1 1 3 1 2 3 + − , (3.11) +α − − μ=− 4σ3 8σ3 α 2α β 2α β k—arbitrary nonzero constant. Substituting (3.1) with (3.9)–(3.11) as well as (2.5)–(2.6) into (2.1), we have the bell type solitary wave solution of RKL model Eq. (1.1)  1/2 1 eiη , u3 (z, t) = cosh 2 − 4σ1k 2 βα ξ − σ3 3

M. Wang et al. / Physics Letters A 363 (2007) 96–101

where σ3 is expressed by (3.9),       1 1 2 1 β 1 −2 − − + kz, ξ = kt + 3α 2α β 2α β 4σ3   1 1 − t η= 2α β       1 β 1 1 1 3 1 2 3 + − +α − − z, + − 4σ3 8σ3 α 2α β 2α β k—arbitrary nonzero constant. 3.2. Similarly, doing the analogous discussion for Eq. (2.13) with (2.14) as that we have done for Eq. (2.8) with (2.9) in Section 3.1, we can obtain the bell type solitary wave solutions of RKL model Eq. (1.1) as follows. 3.2.1. When 2γ = βδ, β = 6α and δ > 0, we have  1/2 1 eiη , u4 (z, t) = 3 cosh 2 − 2σ k 2 ξ − σ1 1

where

  2 1 1 σ1 = δ − 1 + (−1 < σ1 < 0), δ 3 3   1 1 2 3α − z, kz, η= t− ξ = kt + 2σ1 3α 3α 27α 2 k—arbitrary nonzero constant. 3.2.2. When 2γ = βδ, β = 6α and δ < 0, we obtain  1/2 1 eiη , u5 (z, t) = 3 cosh 2 − 2σ k 2 ξ − σ3 3

where

  2 1 1 σ3 = δ − 1 + (σ3 < −1), δ 3 3   1 1 2 3α − z, kz, η= t− ξ = kt + 2σ3 3α 3α 27α 2 k—arbitrary nonzero constant. 4. Kink type solitary wave solutions According to the solution (1.5) of ODE (1.2) with coefficients (1.6) (p = 2), we conclude that Eq. (2.8) with coefficients (2.9) and Eq. (2.13) with coefficients (2.14) admit a solution in the form  1/2   A0 1 1 A0 ξ , ± tanh φ(ξ ) = (4.1) C1 k 2 2 2 k2 provided that 1 1 − β1 )3 − ( 2α − β1 )2 − μ α( 2α

( 12

−

3α 2 β )k

=

A0 k2

(A0 > 0),

99

   A0 1 β 1 γ , − 2 = −2 − 2 2k α k2 3k α

C1 = −

1 γ >0 3k 2 α

(4.2) and   1 A0 ω + = 2 (A0 > 0), − 3 2 2 αk 3α k k 3 2 δ − 2 = − 2 −A0 δ, (4.3) C1 = − 2 > 0 k k k respectively. When 2γ = βδ, αβ > 0 and αγ < 0, which implies δ < 0, from (4.2) we obtain 3 β > 0, A0 = − 8δ α       1 1 3 β 1 3α 1 3 1 2 μ=α − + − − − . 2α β 2α β 8δ α 2 β

(4.4)

When 2γ = βδ, β = 6α and δ < 0, from (4.3) we obtain   9 1 9α A0 = − > 0, (4.5) ω= − k, 4δ 4δ 3α k—arbitrary non-zero constant. Substituting (4.1) with (4.4) as well as (2.5)–(2.6) into (2.1) we have a kink type solitary wave solutions of RKL model Eq. (1.1)   1/2  3 1 1 3 β ξ eiη , ± tanh − u6 (z, t) = − 2δ 2 2 8δ αk 2 where      

1 1 1 2 1 3β ξ = k t + 3α −2 − − + z , 2α β 2α β 8δ   1 1 − t η= 2α β        1 3 β 1 3α 1 3 1 2 1 − + − − − z. + α 2α β 2α β 8δ α 2 β Substituting (4.1) with (4.5) as well as (2.10)–(2.12) into (2.1) we have another kink type solitary wave solution of RKL model Eq. (1.1)   1/2  3 1 1 1 3 − 2ξ eiη , ± tanh u7 (z, t) = − 2δ 2 2 2 δk where     1 9α − z , ξ =k t + 4δ 3α

η=

1 2 z. t− 3α 27α 2

5. Algebraic solitary wave solutions According to the solution (1.7) of ODE (1.2) with coefficients (1.8) (p = 2), we conclude that Eq. (2.8) with (2.9) admits a solution in the form  1/2 1 φ(ξ ) = (5.1) , σ > 0, ξ2 + σ

100

M. Wang et al. / Physics Letters A 363 (2007) 96–101

provided that 1 1 α( 2α − β1 )3 − ( 2α − β1 )2 − μ

( 12 −

3α 2 β )k

= 0,

1 γ 1 β − 2 = −σ. − 2 = 1, (5.2) 2k α 3k α When 2γ = βδ, αβ < 0 and αγ > 0, which implies δ < 0, from (5.2) we obtain  δ β σ = − > 0, k=± − , 3 2α     1 1 3 1 2 1 − . − − μ=α (5.3) 2α β 2α β Substituting (5.1) with (5.3) as well as (2.5)–(2.6) into (2.1) we have algebraic solitary wave solutions of RKL model Eq. (1.1)  1/2 1 eiη , u8 (z, t) = ξ 2 − 3δ where



    

β 1 1 2 1 1 ξ =± − −2 t + 3α − − z , 2α 2α β 2α β         1 1 1 3 1 2 1 1 − − t+ α − − η= z. 2α β 2α β 2α β 6. Conclusion and generalization

Eq. (6.1) admits the bell type solitary wave solutions  1/p 1 u1 (z, t) = eiη , β 1 cosh p − (p+1)σ ξ − σ 1 k2 α 1

where 1p+2 δ− σ1 = 4p+1 either



ξ = kt +

 1+



1p+2 δ 4p+1

 1 β − kz, (p + 2)σ1 3α

2 (−1 < σ1 < 0),

η=

1 2 z t− 3α 27α 2

or

      1 1 1 2 1 β ξ = kt + 3α −2 − − + kz, 2α β 2α β (p + 2)σ1        1 1 1 1 3 1 2 1 − − t+ α − − η= 2α β 2α β 2α β   β 1 3α 1 + − z. (p + 2)σ1 α 2 β 6.2.2. When 2γ = βδ, αβ < 0 and αγ > 0, which implies δ < 0. (6.3) Eq. (6.1) admits the bell type solitary wave solutions  1/p 1 eiη , u2 (z, t) = β 1 cosh p − (p+1)σ k 2 α ξ − σ2 2

We have investigated RKL model Eq. (1.1), from which we briefly come to the conclusion as follow. 6.1. When 2γ = βδ, αβ > 0 and αγ > 0, which implies δ > 0, there exists the bell type solitary wave solutions for Eq. (1.1). When 2γ = βδ, αβ > 0 and αγ < 0, which implies δ < 0, there exist either bell type solitary wave solutions or kink type solitary wave solutions for Eq. (1.1). When 2γ = βδ, αβ < 0 and αγ > 0, which implies δ < 0, there exist either bell type solitary wave solutions or algebraic solitary wave solutions for Eq. (1.1). 6.2. In addition, as a generalization of Eq. (1.1), we consider the generalized higher order NLSE in the form   iuz + utt + 2|u|p u + iαuttt + iβ |u|p u t   + iγ |u|2p u t + δ|u|2p u = 0, p > 0. (6.1) Doing the analogous discussion for Eq. (6.1) as that we have done for Eq. (1.1), we also get the following results. 6.2.1. When 2γ = βδ, αβ > 0 and αγ > 0, which implies δ > 0. (6.2)

where

   1p+2 1p+2 2 (0 < σ2 < 1), δ+ 1+ δ σ2 = 4p+1 4p+1       1 1 2 1 β 1 −2 − − + kz, ξ = kt + 3α 2α β 2α β (p + 2)σ2        1 1 1 1 2 1 1 −2 − t+ α − − η= 2α β 2α β 2α β   1 β 1 3α + − z. (p + 2)σ2 α 2 β Under the constraint condition (6.3), Eq. (6.1) also admits the algebraic solitary wave solutions 1/p  1 eiη , u3 (z, t) = p+2 ξ 2 − 4(p+1) δ where



p2 β 2(p + 2) α   

  1 1 2 1 1 −2 − − z , × t + 3α 2α β 2α β         1 1 1 1 3 1 2 1 − − t+ α − − z. η= 2α β 2α β 2α β ξ =± −

M. Wang et al. / Physics Letters A 363 (2007) 96–101

6.2.3.

Acknowledgements

When 2γ = βδ, αβ > 0 and αγ < 0, which implies δ < 0. (6.4) Eq. (6.1) admits the bell type solitary wave solutions  1/p 1 u4 (z, t) = eiη , β 1 cosh p − (p+1)σ k 2 α ξ − σ3 3

where 1p+2 δ− σ3 = 4p+1 either ξ = kt +

101







1p+2 1+ δ 4p+1

 1 β − kz, (p + 2)σ3 3α

2 ,

η=

σ3 < −1,

1 2 z, t− 3α 27α 2

or

      1 1 2 1 β 1 −2 − − + kz, ξ = kt + 3α 2α β 2α β (p + 2)σ3        1 1 1 3 1 2 1 1 − − t+ α − − η= 2α β 2α β 2α β   β 1 3α 1 + − z. (p + 2)σ3 α 2 β Under the constraint condition (6.4), Eq. (6.1) also admits the kink type solitary wave solutions  2(p + 1) u5 (z, t) = − (p + 2)δ   1/p 1 1 2(p + 1) β p × − eiη , ± tanh ξ 2 2 2 (p + 2)2 δk 2 α where either   1 1 2 2(p + 1) − z, kz, η= t− ξ = kt + 2 3α (p + 2) δ 3α 27α 2 or       1 1 2 1 2(p + 1)β 1 −2 − − + ξ = kt + 3α kz, 2α β 2α β (p + 2)2 δ        1 1 1 1 3 1 2 1 − − t+ α − − η= 2α β 2α β 2α β   2(p + 1) β 1 3α + − z. β (p + 2)2 δ α 2

This work is supported in part by the Natural Science Foundation of Education Department of Henan Province of China (Grant No. 2006110002) and the Science Foundation of Henan University of Science and Technology (Grant Nos. 2004ZD002 and 2006ZY001). References [1] R. Radhakrishnan, A. Kundu, M. Lakshmanan, Phys. Rev. E 60 (1999) 3314. [2] D. Anderson, M. Lisak, Phys. Rev. A 27 (1983) 1393. [3] D. Mihalache, A. Truta, L.C. Crasovan, Phys. Rev. E 56 (1997) 1064. [4] M. Gedalin, T.C. Scott, Y.B. Band, Phys. Rev. Lett. 78 (1997) 448. [5] N. Akhmediev, A. Ankiewicz, Solitons: Nonlinear Pulses and Beams, Chapman & Hall, London, 1997. [6] R.H. Enns, S.S. Rangenekar, A.E. Kaplan, Phys. Rev. A 36 (1987) 1270. [7] A. Kumar, T. Kurz, W. Lauterbon, Phys. Lett. A 236 (1997) 367. [8] D.I. Pushkarov, S. Tanev, Opt. Commun. 124 (1996) 354. [9] S. Tanev, D.I. Pushkarov, Opt. Commun. 141 (1997) 322. [10] C. Zhou, X.T. He, S. Chen, Phys. Rev. A 46 (1992) 2277. [11] P. Honzatko, Opt. Commun. 127 (1996) 363. [12] W.P. Hong, Opt. Commun. 194 (2001) 217. [13] M.J. Ablowitz, P.A. Clarkson, Solitons, Nonlinear Evolution Equations and Inverse Scattering Transform, Cambridge Univ. Press, Cambridge, 1991. [14] R. Hirota, J. Math. Phys. 14 (1973) 810. [15] M.L. Wang, Phys. Lett. A 199 (1995) 169. [16] M.L. Wang, Phys. Lett. A 213 (1996) 279. [17] M.L. Wang, Y.B. Zhou, Z.B. Li, Phys. Lett. A 216 (1996) 67. [18] M.L. Wang, Y.M. Wang, Phys. Lett. A 287 (2001) 211. [19] S.K. Liu, Z.T. Fu, S.D. Liu, Q. Zhao, Phys. Lett. A 289 (2001) 69. [20] Z.Y. Yan, Chaos Solitons Fractals 18 (2003) 299. [21] E.J. Parkes, B.R. Duffy, Comput. Phys. Commun. 98 (1996) 288. [22] Z.B. Li, Y.P. Liu, Comput. Phys. Commun. 148 (2002) 256. [23] E.G. Fan, Phys. Lett. A 277 (2000) 212. [24] Z.Y. Yan, Phys. Lett. A 292 (2001) 100. [25] B. Li, Y. Chen, H.Q. Zhang, Chaos Solitons Fractals 15 (2003) 647. [26] E. Yomba, Chaos Solitons Fractals 20 (2004) 1135. [27] Sirendaoreji, S. Jiong, Phys. Lett. A 309 (2003) 387. [28] M.L. Wang, Y.B. Zhou, Phys. Lett. A 318 (2003) 84. [29] Y.B. Zhou, M.L. Wang, Y.M. Wang, Phys. Lett. A 308 (2003) 31. [30] M.L. Wang, X.Z. Li, Chaos Solitons Fractals 24 (2005) 1257. [31] X.Y. Li, S. Yang, M.L. Wang, Chaos Solitons Fractals 25 (2005) 629. [32] X.Z. Li, J.L. Zhang, Y.M. Wang, et al., Acta Phys. Sinica 53 (2004) 4045. [33] E. Yomba, Phys. Lett. A 336 (2005) 463.