Received 13 July 2016; accepted 22 August 2016; published 25 August 2016

1. Introduction
Consider the following nonlinear Schrödinger equation
(1)
With the boundary conditions
(2)
And the initial condition
(3)
where
,
is the complex-valued wave function.
and
are constant,
is a bounded real function. This equation plays important roles in nonlinear physics. It can describe many nonlinear phenomena including plasma physics [1] , hydrodynamics [1] [2] , self-focusing in laser pulses [3] , propagation of heat pulses in crystals, models of protein dynamics [4] , quantum mechanics [5] , models of energy transfer in molecular systems [6] and quantum mechanics and optical communication [7] - [9] and so on.
In the past few years a great deal of efforts has been expended to solve NLS equations. It is more difficult to find the analytical solutions of the NLS equation, so the study of the numerical solution of NLS equation in the theory and application is important. Its numerical solutions have been researched by many authors. For example, finite difference method [10] [11] , quasi-interpolation scheme [12] , quadratic B-spline finite element scheme [13] , compact split-step finite difference method and pseudo-spectral collocation method [14] [15] , exponential spline method [16] , spline methods [17] [18] , split-step orthogonal spline collocation method [19] , a high-order and accurate method [20] , linearly implicit conservative scheme [21] .
The aim of this paper is to give an exponential spline interpolation method for the NLS equation. The paper is organized as follows. In Section 2, construction of the method is presented. The stability analysis of the scheme is investigated in Section 3. In Section 4, the computation of conserved quantities and error norms are given. In Section 5, two numerical examples are presented to demonstrate our theoretical results. The last section is a brief conclusion.
2. Construction of Exponential Spline Interpolation Method
We set up a grid in the
plane with grid points
and uniform grid spacing h and k, where
and
.
In the interval
, a exponential spline function
is given by
(4)
where
are coefficients to be determined,
and
are the auxiliary functions which contain a stiffness parameter
which will be used to raise the accuracy of the method, on the support
and are given by
(5)
(6)
Since the Taylor series expansions of the hyperbolic functions are
(7)
(8)
We note that
and
tend to
and
in the limit of p tending to zero, and in the opposite limit of p tending to infinity the nonlinear terms in
and
vanish as
.
So the exponential spline defined above share a number of interesting properties:
(1) When
,
reduces to cubic spline; when
,
reduces to linear spline.
(2) A change of character of the exponential spline function is from linear to third order polynomial on adjacent support intervals.
(3) In the general case the stiffness parameters p are different on every interval which provides the extremely high flexibility of the exponential spline function.
We wish to find
in Equation (4),
, Letting
be the unknown second derivative of the exponential spline of interpolation at the grid points, we can obtain the following representation for
on
in terms of the known interpolation data
and the unknown spline second derivatives ![]()
(9)
The terms involving the values
and
represent the linear interpolation part of
. The terms involving the second derivatives
and
introduce the curvature.
The function
on the interval
is obtained with
replacing i in Equation (9).
The continuity requirement for the first derivative
at the point
yields the following equation:
(10)
where ![]()
Remark 1.
(1) By expanding Equation (10) in Taylor series, the truncation error for Equation (10) is of the form
(11)
where
.
For
,
, the truncation
error in space of the relation (10) is of
.
From Equation (10), we can obtain
(12)
Or
(13)
Further, when
, then
,
, the truncation error in space of the relation (10) is of
, Equation (2.7) can be rewritten as
(14)
(15)
In order to get the error estimates of Equation (10), we put
in Equation (12), where E and D are the shift and differential operators respectively, and expand them in powers of hD, we have
(16)
Or
(17)
At the grid point
, Equation (1) can be discretized by
(18)
From Equation (18), we have
(19)
(20)
(21)
Substituting Equation (19), Equation (20) and Equation (21) into Equation (15) and after some simplifications, we obtain
(22)
where ![]()
![]()
The local truncation error of the relation (22) is of
.
The boundary conditions (2) and the system given in the Equation (22) consists of
equations in
unknown. We can write this system in a matrix form as follows:
(23)
where
,
Once the vectors
are computed,
, unknown vectors can be found repeatedly by solving the recurrence relation (23).
3. Stability Analysis
Following the von Neumann technique, we first linearize the nonlinear term in Equation (18) by making the quantity
as locally constant
and assume that the numerical solution can be expressed by means of a Fourier series
(24)
where
,
is the amplitude at time level j,
is the wave number and h is the element size. Substituting Equation (24) into Equation (22), the amplification factor can be written as
(25)
Using Eulers formula, we have
![]()
where
,
Since
![]()
Thus this method is unconditionally stable.
4. Computation of Conserved Quantities and Error Norms
The nonlinear Schrödinger equation possesses two conservation quantities:
(1) Mass conservation:
(26)
Calculated by
(27)
(2) Energy conservation: If
and
are independent of t, then
(28)
Calculated by
(29)
where
and u are the approximate solution at n-th time step at j-th node and exact solution, respectively.
The maximum error norm
and discrete root mean square error norm
will be calculated
(30)
(31)
The relative error of numerical solution is defined as
(32)
5. Numerical Results
In the section, we present the results of our numerical experiments for the proposed scheme described in the previous section.
Example 1. Consider the one dimensional Gross-Pitaevskii equation
(33)
With the analytical solution
(34)
Conserved quantities and error norms at various times are recorded in Table 1. The real and imaginary parts of the numerical and exact solutions are tabulated in Table 2, the numerical results reveal the accuracy of the proposed method.
The absolute error at different space step sizes h at time
are shown in Figure 1, it can be seen that the absolute errors becomes smaller as decreasing h.
Example 2. Consider the equation (1) with ![]()
(35)
The exact solution of this problem is
(36)
![]()
Table 1. Conserved quantities and error norms at various times for example 1 with
.
![]()
![]()
Table 2.The real and imaginary parts of the numerical and exact solutions for Example 1 with
.
.
![]()
Figure 1. The absolute error at different h for example 1 with
.
![]()
Table 3. Conserved quantities and error norms at various times for example 2 with
.
![]()
![]()
Figure 2. The numerical solution at various times t = 1, 2, 3, 4 with
.
![]()
Figure 3. The numerical solutions and analytical solutions for k = 0.01, h = 0.1 at time t = 3.
![]()
Figure 4.The numerical solutions and analytical solutions for k = 0.01, h = 0.1 at time t = 4.
Conserved quantities and error norms at various times are presented in Table 3. The numerical results reveal that the values of
is almost constant while the values of
differ slightly and the errors are very small.
The numerical solutions at various times are given in Figure 2. The numerical solutions and analytical solutions at time
and
are shown in Figure 3 and Figure 4, respectively. The absolute error at time
and
are plotted in Figure 5 and Figure 6, respectively. It observed that (1) the propagation of solitary wave is rightward while preserving unchanged shape; (2) our method gives a good approximation compared with the exact solutions.
6. Conclusion
A numerical method based on exponential spline interpolation function is applied to study a class of nonlinear Schrödinger equation. We use exponential spline collocation method, which results in tri-diagonal systems of
![]()
Figure 5. The absolute error for k = 0.01, h = 0.1 at time t = 3.
![]()
Figure 6. The absolute error for k = 0.01, h = 0.1 at time t = 4.
equations that can be solved efficiently by the Thomas algorithm. The numerical simulations confirm and demonstrate the reliability and efficiency of the schemes and tell us that the method is applicable technique, relatively simple and approximates the exact solution very well.
Acknowledgements
The authors would like to thank the editor and the reviewers for their valuable comments. This work was supported by the Natural Science Foundation of Guangdong (2015A030313827).