A Padé compact high-order finite volume scheme for nonlinear Schrödinger equations

Applied Numerical Mathematics 85 (2014) 115–127

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A Padé compact high-order finite volume scheme for nonlinear Schrödinger equations Wei Gao ∗ , Hong Li, Yang Liu, XiaoXi Wei School of Mathematical Sciences, Inner Mongolia University, West Daxue Road 235, Huhhot, 010021, PR China

a r t i c l e

i n f o

Article history: Received 15 April 2013 Received in revised form 2 January 2014 Accepted 31 May 2014 Available online 11 July 2014 Keywords: Padé compact scheme Finite volume Schrödinger equation

a b s t r a c t In this work, a Padé compact high-order finite volume scheme is presented for the solution of one-dimensional nonlinear Schrödinger equations. The compact high-order finite volume schemes posses inherent conservation of the equations and high order accuracy within small stencils. Fourier error analysis demonstrates that the spectral resolution of the Padé compact finite volume scheme exceeds that of the standard finite volume schemes in terms of the same order of accuracy. Besides, the linear stability of the temporal discretization scheme is also performed by using the Fourier analysis. Numerical results are obtained for the nonlinear Schrödinger equations with various initial and boundary conditions, which manifests high accuracy and validity of the Padé compact finite volume scheme. © 2014 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Consider the following nonlinear Schrödinger equation

∂u ∂ 2u + ν 2 + q|u |2 u = 0, a ≤ x ≤ b, t > 0 (1) ∂t ∂x where ν = 0 and q ≤ 0 are real parameters and u = u (x, t ) is a complex-valued function. It is known as the cubic nonlinear i

Schrödinger equation (CNLS) which plays a key role in modeling some physical phenomena such as the evolution of waves in water and plasma, optical pulse and self-focusing effects in laser pulses [2,17,20]. The function u (x, t ) can describe the propagation of the weakly nonlinear, strongly dispersive and monochromatic wave [24]. Many numerical methods were presented in order to simulate the complicated soliton models of CNLS equations since few of analytical solutions of CNLS equations can be solved [1,21]. Finite element methods, finite difference methods and finite volume methods are major approximation approaches with the latter two methods prevailing in computation of fluid flows and nonlinear physics. Numerical solutions of the CNLS equation by the finite difference method can be found in references [3,6,10,22,27]. Additionally, several non-standard finite element methods such as the local discontinuous Galerkin method [13,26], the space–time finite element method [7,8], the B-spline finite element method [4] and so on were recently employed to numerically solve the Schrödinger equation. In contrast to these two methods, finite volume methods are preferred in computational fluid dynamics. The main reason is that they are easy to implement on uniform meshes and more handy to develop than the finite element methods do on nonuniform meshes. Moreover, finite volume methods are inherently conservative approaches as they result from an integral formulation of the conservation law.

*

Corresponding author. E-mail addresses: [email protected] (W. Gao), [email protected] (H. Li), [email protected] (Y. Liu), [email protected] (X. Wei).

http://dx.doi.org/10.1016/j.apnum.2014.05.016 0168-9274/© 2014 IMACS. Published by Elsevier B.V. All rights reserved.

116

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

Fig. 1. The cell-centered meshes.

Owing to the Lagrangian interpolation, larger stencils need to be supplied for classical numerical schemes to get higher order of accuracy, which result in considerable costs of computation. On the contrary, compact schemes by Lele [12] used smaller stencils to get higher order of accuracy than the standard numerical schemes do. They also have better spectral resolution than the standard numerical schemes with the same order of accuracy. Thus, compact schemes are a good option for simulation of wave dominated problems [12,28]. Some researchers employed compact finite difference methods for the numerical solution of the nonlinear Schrödinger equations [5,19,25]. In contrast, few of studies [9,18,11] were conducted on the solution of the CNLS equation by compact finite volume schemes. The aim of this paper is to present a 4th order Padé compact finite volume method for numerical solution of the CNLS equation. We compare the spectral resolution of the compact finite volume and standard finite volume central schemes. The linear stability of the temporal discretization scheme is analyzed by using the von Neumann stability analysis. Several test cases demonstrate high order accuracy and validity of the Padé compact finite volume scheme for computation of the CNSL equation. 2. The Padé compact finite volume method 2.1. Space discretization by the finite volume method The nonlinear Schrödinger equation can be rewritten as

∂u ∂ 2u = i ν 2 + iq|u |2 u ∂t ∂x

(2)

The space discretization in the finite volume framework is performed on a cell centered mesh as plotted in Fig. 1. By integrating on the interval [x j −1/2 , x j +1/2 ], Eq. (2) can be reformed as

du¯ (x j , t ) dt



 = i ν u (x j +1/2 , t ) − u (x j −1/2 , t ) + iq 





x j +1/2

|u |2 udx

(3)

x j −1/2

where

u¯ j (t ) =

1

x



x j +1/2

u (ξ, t ) dξ

(4)

x j −1/2

is the cell-averaged value of the function u (x, t ). We approximate Eq. (3) by the following numerical scheme

du¯ j (t ) dt











x j +1/2

= i ν u j +1/2 − u j −1/2 + iq

|u |2 udx

(5)

x j −1/2

where u¯ j (t ) is the numerical approximation to the cell average u¯ (x j , t ). It is the key point for the numerical solution of the nonlinear Schrödinger equation that how ones approximate the first derivative u  at the cell face and the definite integral on the right-hand side of (5). 2.2. Padé compact approximation Next, we will present the Padé compact approximation for u j +1/2 by using the cell average u¯ j (t ). A generic formulation of the Padécompact finite volume scheme [9,11] can be written as follows k 

am ε (−mx) u  + u  +

m =1

k  m =1

where k, l are natural numbers and

ε s f (x) = f (x + s)

am εmx u  =

1

x

 −

l  n =1

bn ε (−n+1/2)x u¯ +

l 

bn ε (n−1/2)x u¯

  + O x2(k+l)

(6)

n =1

ε s : C ∞ → C ∞ with s ∈ R is the shift operator defined as (7)

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

117

Table 1 Coefficients and leading terms of truncation error for Padé compact high-order finite volume schemes. Order

a1

4

1 10

a2 –

6

1 11

–

8

344 1179

23 2558

a3

b1

b2

b3

Leading terms of H.O.T.

–

6 5

–

–

x) − ( u (5) 200

–

23(x)6 55440

–

79(x)8 (9) − 2971080 u

–

51 44

3 44

–

265 262

155 786

4

u (7)

The unknown coefficients am , bn in (6) can be determined via a Taylor expansion of u about the point x j +1/2 . We make the leading error term of the Talyor expansion match computational precision, which results in the accuracy order. This yields a system of equations to be satisfied by the unknown coefficients. Due to the cell-averaged values employed on the right-hand-side of (6), the integration of u on the interval [x j −1/2 , x j +1/2 ] have to be done after the Taylor expansion of u about the point x j +1/2 . Consider the term u¯ j as an example. The Taylor expansion of the point-wise value of u about the point x j +1/2 reads as follows

2 

 ∂u 1 2 ∂ u u (x) = u j +1/2 + (x − x j +1/2 ) + (x − x j +1/2 ) + H.O.T. ∂ x j +1/2 2! ∂ x2 j +1/2

(8)

where H.O.T. means ‘high order terms’. Combining (8) with (4), one obtains

u¯ j =

=

1

x 1

x



x j +1/2

u (x) dx x j −1/2



xu j +1/2 +

(−x)2  (x)3  u j +1/2 + u j +1/2 + H.O.T. 1! 2 2! 3

(9)

According to Eqs. (6), (8) and (9) with k = l = 1, one obtains at the point x = x j +1/2

a1 ε (−x) u j +1/2 + u j +1/2 + a1 ε x u j +1/2 =

1

x

  (−b1 ε (−1/2)x u¯ j +1/2 + b1 ε (1/2)x u¯ j +1/2 ) + O x4

(10)

Then, the 4th-order Padé compact finite volume scheme for the value of the first derivative u  at the cell face x = x j +1/2 is presented as

1 10

u j +3/2 + u j +1/2 +

1 10

u j −1/2 =

6 5 x

(u¯ j +1 − u¯ j )

(11)

Although the Padé compact finite volume schemes can be presented with higher order than four (see Table 1), this is to be done on stencils larger than those used for the standard second order central scheme. The 4th order Padé compact finite volume scheme (11) uses the smallest stencil possible for central schemes (only two cells involved), which is an attractive feature from the computational point of view. Besides, the Padé type schemes show a better resolution than classical central schemes although accuracies of 4th order Padé compact and standard 4th order central schemes are of course the same. The integral on the right-hand side of (5) needs to be numerically computed due to its nonlinear integrand. The Simpson formula is employed to match the 4th order Padé compact finite volume scheme. It reads as



x j +1/2

f dx =

x 6

  [ f j −1/2 + 4 f j + f j +1/2 ] + O (x)4

(12)

x j −1/2

It is noted that terms f j −1/2 , f j and f j +1/2 in the Simpson formula are the point-wise values. Thus, they need to be reconstructed by the cell-averaged values u¯ j to at least 4th order of accuracy. A 4th-order compact finite volume scheme is also presented for the reconstruction of the cell face values f j −1/2 and f j +1/2 by using the cell-averaged values ¯f j . According to results in [9,11], the generic compact interpolation formula to approximate the cell face value satisfies

a b c β f j −3/2 + α f j −1/2 + f j +1/2 + α f j +1/2 + β f j +5/2 = ( ¯f j + ¯f j +1 ) + ( ¯f j +2 + ¯f j −1 ) + ( ¯f j +3 + ¯f j −2 ) 2

2

2

(13)

The parameters a, b, c , α , β in (13) can also be determined by a Taylor expansion at the point x j +1/2 . The resultant order of accuracy is derived by requiring terms with lower order vanishing up to some term. One obtains a 4th-order compact finite volume scheme as

1 4

f j −1/2 + f j +1/2 +

1 4

f j +3/2 =

3 ¯ ( f j + ¯f j +1 ) 4

(14)

118

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

Moreover, the term f j can be reconstructed from the point-wise values of f j −1/2 and f j +1/2 by using the 4th order cellcentered compact finite difference scheme by Lele [12] which reads as

1 22

1

f j −1 + f j +

22

f j +1 =

12 11

( f j −1/2 + f j −1/2 )

(15)

2.3. Fourier error analysis

∂u ∂ 2u =ν 2, ∂t ∂x

ν >0

u (x, t ) = u (x + L , t )

(16)

The one-dimensional unsteady linear diffusion equation with periodic boundary conditions is used as a model for analysis of the spectral resolution of the Padé compact finite volume schemes and the standard central schemes. The 2nd order and 4th order central finite volume schemes approximating u  j +1/2 read as follows

u  j +1/2 = u  j +1/2 =

1

x

(u¯ j +1 − u¯ j )

(17)



1 12x



15(u¯ j +1 − u¯ j ) − (u¯ j +2 − u¯ j −1 )

(18)

The accuracy of a numerical scheme is often expressed in terms of the order of the leading term of truncation error from the Taylor expansion. It provides the asymptotic rate of convergence of the numerical solution to the analytical solution. It presents no information on the spectrum of the truncation error of numerical schemes. The Fourier analysis is the classical technique to study the spectrum of the truncation error of the numerical scheme [23]. It was also used by Lele [12] and Mahesh [15] for studying compact finite difference schemes. Next, we are to use the Fourier analysis to characterize the spectrums of standard central and Padé compact finite volume schemes. An arbitrary periodic function can be decomposed into its Fourier components e i κ x , where κ is the wave number. We will examine how a given finite volume scheme approximates the first derivative of u (x) = e i κ x on uniform grids. The exact relation for the linear diffusion equation at the point x = x j +1/2 yields

du¯ dt

= −νκ 2 u¯

(19)

So, one obtains

u¯ (t ) = u¯ 0 e −νκ

2

t

(20)

where u¯ 0 = u¯ (0). This means that the cell averaged value u¯ decays at an exponential rate which is proportional to the square of the waver number κ . Substituting

u=e

iκ x

u¯ j =

,



x j +1/2

1

x

e i κ x dx

(21)

x j −1/2

into the standard 2nd order central finite volume scheme (17), one can obtain that

LHS = i κ e i κ (x j +x/2)

(22)

and

RHS = i



4

(κ x)x

sin

κ x

2

e i κ (x j +x/2)

2

= i κ ∗ e i κ (x j +x/2)

(23)

Here, κ ∗ = κ ∗ approximates the exact wave number κ is usually considered as a measure of the accuracy of the numerical scheme. The modified wave number of the 4th-order central finite volume scheme reads

[ (κ 4x)x (sin

κ∗ =



5

(κ x)x

κ x )2 ] is the modified wave number. The degree to which 2

sin

κ x 2

2 +

1 6(κ x)x





cos(2κ x) − cos(κ x)

(24)

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

119

Fig. 2. Decay ratios for different numerical schemes.

Similarly, the modified wave number of the 4th order Padé compact finite volume schemes satisfies

iκ ∗



1 10

e i κ (−x) + 1 +



1 10

e i κ (x)

= (−i )

 1  i κ (−x) e − 2 + e i κ (x) 5 x κ  x 6

(25)

which gives

κ∗ =

x 2 ) 4(sin κ  2

1

6

5 x κ  x 1 +

1 5

(26)

cos(κ x)

Generally, a similar approach result in the modified wave numbers of higher order compact Padé finite volume schemes represented as ∗

κ =

4

(κ x)x

sin

κ x 2

 l

n=1 bn sin((n − 1/2) k 1 + 2 m=1 am cos(m

κ x) κ x)

(27)

∗

Fig. 2 shows the decay ratio κκ as a function of the grid wave number κ x for standard 2nd and 4th order central finite volume schemes and 4th, 6th and 8th order Padé compact finite volume schemes. As seen in the figure, Padé compact finite volume schemes perform better spectral resolution than standard central schemes do in the high frequency range. 2.4. Treatment of boundary conditions When the control-cell faces intersect the domain boundary, the boundary conditions need to be implemented for the compact Padé finite volume scheme. The 4th order compact approximation for boundary conditions is indispensable for the schemes to keep global accuracy and stability as certified in [9,12]. In comparison with the compact finite difference schemes, the imposition of the boundary conditions is relatively simple for the compact finite volume schemes. The periodic and Dirichlet boundary conditions are imposed for the boundary face value u 1/2 in a straight way on the cell-centered mesh. The periodic and Neumann boundary conditions can also be directly set for the boundary-face derivative u 1/2 . For a diffusion equation with the Dirichlet boundary condition, the boundary face derivative u 1/2 could be in no conjunction with the prescribed value on the boundary. One can deal with this case by employing the biased 4th order compact finite volume formula [9]. It reads as

u 1/2 + 6u 3/2 = −

1



5

x 3

uB +

89 18

u¯ 1 −

127 18

u¯ 2 +

4 9

 u¯ 3

(28)

where u B denotes the assigned value on the left boundary. Note that the specular version of (28) can be used on the right boundary. 2.5. Time discretization and linear stability After the space discretization is fulfilled, we obtain a set of ordinary differential equations

du dt

= L (u )

(29)

120

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

where the operator L (u ) at the mesh point is approximated by the 4th order compact Padé finite volume scheme. This set of ordinary differential equations is discretized by the 4th order Runge–Kutta scheme (RK4) as follows



u (1 ) = L u n





u

(2 )

u

(4 )



t

n

(1 )

=L u + u 2

 t (2) u (3 ) = L u n + u 2



n

= L u + tu (3)

u n +1 = u n +

t  6



u (1) + 2u (2) + 2u (3) + u (4)



(30)

Next, the linear stability analysis for the 4th order Runge–Kutta scheme is effective to handle the restriction of the time step t. The linearized version of the nonlinear Schrödinger equation can be written as

ut = i ν u xx + icu

(31)

where c = max(q|u (x, 0)| ) means a reference constant. By integration on the interval [x j −1/2 , x j +1/2 ], the integral version of (31) can be written as 2

du¯ j dt

=

iν 

x

 (u x ) j +1/2 − (u x ) j −1/2 + ic u¯ j

(32)

For simplicity, the cell-averaged value u¯ j in (32) is subject to the periodic boundary condition. According to the superposition principal of the linear equation, we need only to consider the effect of a single Fourier wave

un = g n e i2mπ x/ L ,

m = 0, 1 , . . . , N − 1

(33)

where g denotes the amplification factor. It is noted that one must consider the cell-averaged value of u volume framework

u¯ nj =



x j +1/2

gn

x

e i2mπ x/ L dx = g n

sin( 2mNπ )

x j −1/2

mπ N

e i2mπ x j / L

n

in the finite

(34)

with N = Lx . Based on the space discretization of the 4th order Padé compact finite volume scheme, Eq. (32) can be reformed in a matrix form as

¯ du dt

=

iν

x

¯ Aux + ic u

(35)

where

⎛¯ ⎞ u1 ⎜ u¯ 2 ⎟ ⎟ ¯ =⎜ u ⎝ .. ⎠ ; . u¯ N

⎛

⎞ ⎛ (u ) x 3/2 ⎜ (u x )5/2 ⎟ ⎟; ux = ⎜ .. ⎠ ⎝ . (u x ) N +1/2

1 ⎜ −1

⎜ ⎝ ...

AN ×N = ⎜ . . .

...

0 1

... 0

... ... ...

... ... . . . −1 1 . . . 0 −1

⎞ −1 ... ⎟ ⎟ ... ⎟ ⎠ 0

(36)

1

and

Bux =

1

x

¯ Du

(37)

According to 4th order Padé compact finite volume scheme (11) with the periodic boundary condition, the coefficient matrices B and D are presented as

⎛

BN ×N

1

1 10

1 10

... ... ... ... ...

⎜ 1 ⎜ 10 ⎜... ⎜ =⎜ ⎜... ⎜... ⎜ ⎝...

1

... ... 1 ... 10 ... ... ... ... ... ... ... ... ... ...

... ... ... ... ...

... ... ... ... ...

1 10

1

...

1

1 10

⎞

...⎟ ⎟ ...⎟ ⎟ ...⎟ ⎟ ...⎟ ⎟ 1 ⎠ 10 1 10

(38)

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

⎛

DN ×N

⎜ ⎜ ⎜ ⎜ =⎜ ⎜ ⎜ ⎜ ⎝

− 65 0

... ... ... ...

6 5 − 65

6 5

... ... ... ... ...

... ... 6 ... 5 ... ... ... ... ... ... ... ... ... ...

... ... ... ... ... 0

... ... ... ... ... − 65

...

0

... ... ... ... ... 6 5 − 65

121

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(39)

Thus, Eq. (35) can be rewritten as

¯ du dt

=



iν

¯ = Hu¯ AB−1 D + icE u 2

(40)

(x)

A, B and D are tridiagonal and circulant matrices. E is the unit matrix. Eigenvalues of matrices A, B and D are presented as



     2mπ 2mπ − i ds − ds sin , λm = d + ds + ds cos N

N

m = 0, 1 , . . . , N − 1

(41)

where the parameters d, d s and d s mean diagonal, subdiagonal and superdiagonal values in matrices, respectively [14]. Consequently, one obtains that eigenvalues of matrices A, B and D are

λ(Am) = 1 − e i2mπ / N

(42)

1

(m)

λB = 1 + e i2mπ / N

(43)

5

(m)

λD =

 6  i2mπ / N e −1 5

(44)

and eigenvalues of the matrix AB−1 D are (m) (m)

m) λ(AB −1 D =

λ A λD

=−

(m)

λB

24

(sin(mπ / N ))2

5 1+

1 5

(45)

cos(2mπ / N )

Substituting (40) into (30) yields the matrix equation on the amplifier matrix G as

G = 1 + (tH) +

1 2

1

1

6

24

(tH)2 + (tH)3 +

(tH)4

(46)

According to von Neumann stability analysis, a necessary condition for the stability of numerical schemes is that every one of eigenvalues of the matrix G is to be less than one. Eigenvalues of the matrix H can be presented as

λ(Hm) = (−i )

24(sin(mπ / N ))2 + c 5 + cos(2mπ / N )

,

m = 0, 1 , . . . , N − 1

(47)

One obtains that eigenvalues of the matrix G satisfy

 1 2 1  3 4  1  t λ(Hm) , λ(Gm) = 1 + t λ(Hm) + t λ(Hm) + t λ(Hm) + 2

6

24

m = 0, 1 , . . . , N − 1

(48)

The modulus of the eigenvalue of the matrix G is

 (m) 2 (t |λ(Gm) |)8 (t |λ(Gm) |)6 λ  = − +1 G 576

72

(49)

(m)

So, |λG |2 < 1 results in the restriction condition of the time step for the 4th-order Runge–Kutta scheme

√

t <

2 2(x)2

ν [6 + c (x)2 ]

(50)

122

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

Table 2 Errors and the convergence order at t=1.0. Mesh h= h= h= h=

4 10 2 10 1 10 1 20

e L 1

Order

e L 2

Order

e L ∞

3.5579e −4

–

1.3120e −3

–

7.2751e −3

–

2.1152e −5

4.072

7.7352e −5

4.084

4.1862e −4

4.121

1.3286e −6

3.993

4.8238e −6

4.001

2.6212e −5

3.998

8.2779e −8

4.004

3.0120e −7

4.003

1.6281e −6

4.009

Order

Fig. 3. Numerical energy and charge by the 4th-order Padé compact finite volume scheme.

3. Numerical results and discussion Computational errors between the numerical and exact solutions are measured by the discrete L 2 norm and L ∞ norm defined by





Er2 = unum − uexa  L



2

  N    num 2  u  = x − u exa j j



Er∞ = unum − uexa  L

(51)

j =0

∞

   = max u num − u exa j j

(52)

0≤ j ≤ N

The order of convergence of the numerical solution are computed by

order =

ln(Erh /Er2h )

(53)

ln 2

Next, we select three test cases for verification of the 4th order Padé compact finite volume scheme. The time step t is specifically determined by ν = 1 and c = max(q|u (x, 0)|2 ) in Eq. (50). Moreover, the parameter q is set to 2 for the first and second case and 2N 2 ( N = 3, 4, 5) for the third case. For the first test case, we get c ≤ 2 due to max |u (x, 0)|2 ≤ 1. For the second case, we get c ≤ 8 due to max |u (x, 0)|2 ≤ 4. For the last test case, we get c ≤ q (q = 2N 2 ) due to max |u (x, 0)|2 ≤ 1. Consequently, the following inequality

√

t <

2 2(x)2 7

√

2

 0.4(x)

<

2 2(x)2



ν [6 + c (x)2 ]

(54)

holds for sufficiently small x. So, we take t = 0.4(x)2 for present computations. 3.1. One soliton The cubic nonlinear Schrödinger equation

i

∂ u ∂ 2u + 2 + q|u |2 u = 0 ∂t ∂x

(55)

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

123

Fig. 4. Modulus of the numerical and exact solutions.

Fig. 5. The interaction of two solitons.

provides a good mathematical model to simulate solitons. Its exact one-soliton solution reads

u (x, t ) = exp(2ix − 3it ) sech(x − 4t )

(56)

for q = 2. The Dirichlet boundary conditions are set as u (−20, t ) = u (20, t ) = 0. This case is to be considered to test the accuracy order of the Padé compact finite volume scheme. Table 2 shows the numerical solution could achieve the 4th order accuracy for L 2 and L ∞ norms. The nonlinear Schrödinger equation has two significant conservative laws on the charge and energy, e.g.

Q =

+∞ +∞     u (x, t )2 dx = u (x, 0)2 dx

−∞

−∞

(57)

124

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

Fig. 6. The bound state solutions of Eq. (60) with N = 3.

   +∞  2 +∞    ∂u   ∂ u (x, 0) 2 q  q 4    − u (x, 0)4 dx E= η  − |u | dx = η ∂x 2 ∂x  2 −∞

(58)

−∞

The computed energy E and charge Q are plotted in Fig. 3. It is seen that two conservative quantities could be kept constant during the computational time range. The numerical solution versus the analytical solution is plotted in Fig. 4, which shows that the curve of the numerical solution agrees well with the one of the analytical solution. 3.2. Interaction of two solitons Coupled with the initial condition and periodic boundary condition

u (x, 0) =

2 

exp

j =1

i 2

 c j (x − x j ) sech(x − x j )

(59)

the cubic nonlinear Schrödinger equation (55) produces a two-soliton solution [26] with c 1 = 4, c 2 = −4, x1 = −10, x2 = 10. The solution is computed in the interval [−20, 20] on a uniform mesh with 200 cells. We compared the present solution with the solution by the local discontinuous Galerkin method [26] at three time instants t = 0, 2.5, 5. Figs. 5(a)–5(c) shows that the numerical solutions by the 4th order Padé compact finite volume scheme are in good agreement with the ones by the local discontinuous Galerkin method. 3.3. Bound state solution This test case shows the bound state solutions of the Schrödinger equation

i

∂ u ∂ 2u + 2 + q|u |2 u = 0 ∂t ∂x

(60)

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

125

Fig. 7. The bound state solutions of Eq. (60) with N = 4.

with the initial condition

u (x, 0) = sech x

(61)

The coefficient q = 2N 2 results in a bound state of N solitons. The numerical solutions are computed with periodic boundary conditions in [−15, 15] on a uniform mesh with 750 cells. The solution for a bound state of solitons would be difficult to be computed due to its narrow structures when N ≥ 3 [26,16]. The numerical solutions of the bound state of solitons are computed for N = 3, 4, 5 at t = 0, 0.2, 0.4, 0.6, respectively. Figs. 6–8 plot the numerical solutions. It is shown in the figures that more sharp peaks appear in a small computational segment [−5, 5] with the coefficient q increasing, which is in quantitative agreement with the theoretical solutions [16]. 4. Conclusion A Padé compact finite volume scheme is presented for nonlinear Schrödinger equations. The compact finite volume schemes guarantee high order accuracy on small stencils and with inherent conservation. The spectral resolution of the present scheme is analyzed by using Fourier analysis. It shows that the compact finite volume schemes perform better spectral resolution than theirs counterparts for standard central schemes. Numerical results are obtained for the nonlinear Schrödinger equations with various initial and boundary conditions, which manifest high accuracy and efficiency of the Padé compact finite volume scheme. Acknowledgements The authors would like to thank the editor and the reviewers for their valuable comments and suggestions to improve the results of this paper. The paper is supported by National Natural Science Foundation of China (Grant No. 11361035, 11301258) and Key Project of Chinese Ministry of Education (12024), Natural Science Foundation of Inner Mongolia Autonomous Region (2011BS0102,

126

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

Fig. 8. The bound state solutions of Eq. (60) with N = 5.

2012MS0108) and the Scientific Research Projection of Higher Schools of Inner Mongolia Autonomous Region (NJZZ12011) and Program of Higher-Level Talents of Inner Mongolia University (SPH-IMU 30105-125133). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

M. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, SIAM, Philadelphia, 1981. G. Agrawal, Nonlinear Fiber Optics, Academic Press, San Diego, 2001. Q. Chang, E.S.E. Jia, Difference schemes for solving the generalized nonlinear Schrödinger equation, J. Comput. Phys. 148 (1999) 397–415. I. Da˘g, A quadratic B-spline finite element method for solving nonlinear Schrödinger equation, Comput. Methods Appl. Mech. Eng. 174 (1999) 247–258. M. Dehghan, A. Taleei, A compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficients, Comput. Phys. Commun. 181 (2010) 43–51. S. Iyengar, G. Jayaraman, V. Balasubramanian, Variable mesh difference schemes for solving a nonlinear Schrödinger equation with a linear damping term, Comput. Math. Appl. 40 (2000) 1375–1385. O. Karakashian, C. Makridakis, A space–time finite element method for the nonlinear Schrödinger equation: the discontinuous Galerkin method, Math. Comput. 67 (1998) 479–499. O. Karakashian, C. Makridakis, A space–time finite element method for the nonlinear Schrödinger equation: the continuous Galerkin method, SIAM J. Numer. Anal. 36 (1999) 1779–1807. M. Kobayashi, On a class of Padé finite volume methods, J. Comput. Phys. 156 (1999) 137–180. A. Kurkinaitis, F. Ivanauskas, Finite difference solution methods for a system of the nonlinear Schrödinger equations, Nonlinear Anal., Model. Control 19 (2004) 247–258. C. Lacor, S. Smirnov, M. Baelmans, A finite volume formulation of compact central schemes on arbitrary structured grids, J. Comput. Phys. 198 (2004) 533–566. S. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16–42. D. Levy, C.-W. Shu, Local discontinuous Galerkin methods for nonlinear dispersive equations, J. Comput. Phys. 196 (2004) 752–772. H. Lomax, T. Pulliam, D. Zingg, Fundamentals of Computational Fluid Dynamics, Springer, Berlin, 2001. K. Mahesh, A family of high order finite difference schemes with good spectral resolution, J. Comput. Phys. 145 (1998) 332–358. J. Miles, An envelope soliton problem, SIAM J. Appl. Math. 41 (1981) 227–230. C. Nore, A. Abid, M. Brachet, Small-Scale Structures in Three-Dimensional Hydrodynamics and Magnethydrodynamics Turbulence, Springer, Berlin, 1996. M. Pillar, E. Stalio, Finite-volume compact schemes on staggered grids, J. Comput. Phys. 197 (2004) 299–340. M. Smadi, D. Bahloul, A compact split step Padé scheme for high-order nonlinear Schrödinger equation with power law nonlinearity and fourth order dispersion, Comput. Phys. Commun. 182 (2011) 366–371. C. Sulem, P. Sulem, The Nonlinear Schrödinger Equation: Self-focusing and Wave Collapse, Springer, New York, 1999.

W. Gao et al. / Applied Numerical Mathematics 85 (2014) 115–127

127

[21] T. Taha, Inverse scattering transform numerical schemes for nonlinear evolution equations and the method of lines, Appl. Numer. Math. 20 (1996) 181–187. [22] E. Twizell, A. Bratsos, J. Newby, A finite-difference method for solving the cubic Schrödinger equation, Math. Comput. Simul. 43 (1997) 67–75. [23] R. Vichnevetsky, J. Bowles, Fourier Analysis of Numerical Approximation of Hyperbolic Equations, SIAM, Philadelphia, 1982. [24] G. Whithman, Linear and Nonlinear Waves, Wiley, New York, 1974. [25] S. Xie, G. Li, S. Yi, Compact finite difference schemes with high accuracy for one-dimensional nonlinear Schrödinger equation, Comput. Methods Appl. Mech. Eng. 198 (2009) 1052–1060. [26] Y. Xu, C.-W. Shu, Local discontinuous Galerkin method for nonlinear Schrödinger equations, J. Comput. Phys. 205 (2005) 72–97. [27] L. Zhang, A high accurate and conservative finite difference scheme for nonlinear Schrödinger equation, Acta Math. Appl. Sin. 28 (2005) 178–186. [28] S. Zhang, S. Jiang, C.-W. Shu, Development of nonlinear weighted compact schemes with increasingly higher order accuracy, J. Comput. Phys. 237 (2008) 7294–7321.