1 Introduction

The nonlinear Schrödinger (NLS) equation is one of the most important mathematical models of modern science, which describes many physical phenomena, and has wide applications in many different fields, e.g., nonlinear optics^[1], plasma physics^[2], semiconductor industry^[3], and fluid dynamics^[4-5].

The initial-boundary value problem of the one-dimensional generalized NLS (GNLS) equation defined on an arbitrary interval x ∈ [L₀, L₁] with the coordinate transformation x=(x-L₀)/(L₁ -L₀) can be written as follows^[1-30]:

(1)

where is the complex unit, ψ (x, t) is a complex-valued wave function or order parameter, and the known real functions r_k (x) and s_k (x) (k=0, 1, 2) describe the initial and boundary conditions. The real parameter α (t) is often taken as a constant in most physics models^[1-5], and we consider it as a real function of the time t in this study. The external potential ν (x, t) is a given real function, whose specific form depends on specific applications. In Bose-Einstein condensation, it is usually chosen as either a harmonic confining potential ν (x, t)=-x²/2 or an optical lattice potential ν (x, t)=A cos(Lx), where A and L are constants^[6-7]. In Eq. (1), the nonlinearities f(ρ) and g(ρ) are real-valued smooth functions with respect to the density ρ=|ψ|², the specific forms of which depend on the practical applications. In Eq. (1), the most popular nonlinear term f(ρ)=λρ, in which λ is a constant^[8-23], is extended to the general form f(ρ), which can be chosen as an arbitrary function as needed, e.g.,

where c₀, s, and c₁ are real constants. The general damping term g(ρ), which is usually ignored in applications^{[6, 8-23]}, is also considered, since the damping effect may play an important role and should not be ignored in some physical processes, e.g., the inelastic collisions in Bose-Einstein condensation^{[7, 24-25, 31-32]}.

Due to the broad applications of Eq. (1), finding accurate and efficient solutions for such equations has attracted considerable research attention in the past few years. Several analytical or semi-analytical methods, e.g., the homotopy analysis method (HAM)^[20-21], the Adomian decomposition method^[20], the differential transform method^[22], and the symbolic expansion method^[23], have been proposed to study NLS equations. However, the applicable scope of these analytical techniques seems to be limited to a few special problems^[20-23]. As a result, numerous numerical methods have been developed to understand the physical behavior of NLS equations^{[7-19, 25-30]}. Gao et al.^[13] used a finite volume method to investigate the NLS equations. Liao et al.^[14] used a symmetric fourth-order compact difference scheme to solve Eq. (1) without the damping term. Wang and Zhang^[19] developed a split-step orthogonal spline collocation method^[19] to solve the NLS equations. Dag^[28] proposed a quadratic B-spline finite element method to study the one-dimensional cubic NLS equations. Mocz and Succi^[29] studied the NLS equations with the smoothed-particle hydrodynamics method. Moreover, many other numerical procedures for solving the NLS equations can be found in Refs. [7]-[12], [15]-[18], [25]-[27], and [30]. These existing numerical methods are effectively applied to the solution of the NLS equations under certain conditions. However, the convergence rates in space of almost all of these methods seem to be confined to 4^{[8-9, 12-18, 30]}, and many of them are only 2^{[15-18, 30]}. Since a high-order accurate method allows for using lower computational cost with respect to a low-order accurate method to ensure a similar accuracy, researchers have been trying to find a high-order accurate method to numerically solve Eq. (1). Bao et al.^{[7, 25-26]} proposed a time-splitting pseudospectral method with high-order accuracy. However, in such a procedure, there is a need to obtain the density ρ by analytically solving a nonlinear differential equation^{[7, 25]}. It is an extremely difficult task for the general damping term g(ρ), although Bao et al.^[31-32] have proposed an approach to find the solution of such an equation and provided the exact solutions for several special forms of the damping term.

Recently, Wang^[33] and Liu et al.^[34-38] developed two wavelet methods for solving nonlinear boundary value problems and Burgers' equations. The proposed wavelet algorithms show a much better accuracy and a much faster convergence rate than many other existing numerical methods^[33-38]. The convergence order in space of such a wavelet algorithm for Burgers' equation can exceed 5. Following these previous studies^[36-38], we combine a wavelet Galerkin technique with the Runge-Kutta method to uniformly solve Eq. (1). Firstly, we propose an approximation scheme for the bounded functions based on Coiflet-type wavelet expansion and a new boundary extension method. Then, a Galerkin procedure based on such a wavelet approximation is used to transform Eq. (1) into a system of nonlinear ordinary differential equations, which can be solved by the classical fourth-order explicit Runge-Kutta method. Finally, a systematic investigation on the efficiency and accuracy of the proposed wavelet formulation for solving Eq. (1) is conducted by use of eight widely considered examples. A comparison between the present wavelet solutions and those achieved from many other existing numerical methods is made to demonstrate the effectiveness of the proposed high-order accurate wavelet algorithm.

2 Sampling approximation of an interval-bounded L²-function

A wavelet multiresolution analysis is an increasing sequence of the closed subspaces of , which can verify the following properties:

(ⅰ) , and is dense in ;

(ⅱ) f(x) ∈ V_j ⇔ f(2x) ∈ V_j+1 for all and all ;

(ⅲ) f(x) ∈ V₀ ⇔ f(x-k) ∈ V₀ for all and all ;

(ⅳ) There exists a function φ(x) ∈ V₀, i.e., the scaling function, such that the sequence is a Riesz basis of the space V₀^[39-40]. Specially, if φ(x) is the scaling function of Coiflets, which is similar to the symmetric wavelets designed by Daubechies, at the request of Coifman, to have scaling functions with vanishing moments^[39], all bounded functions f(x) can be approximated by^{[33, 36-37]}

(2)

where M₁ =∫xφ(x)dx is the first-order moment of the scaling function φ(x), and j is the decomposition level. When , the accuracy of Eq. (2) can be estimated as follows^{[33, 37]}:

(3)

where N is the number of the vanishing moment of the wavelet function corresponding to the scaling function φ(x), the constant C depends only on the smoothness of the function f(x), and the non-negative integers n < N. In this study, the Coiflet-type orthogonal scaling function φ(x), in which M₁=7 and N=6, is adopted.

It can be seen from Eq. (2) that the orthogonal scaling functions originally form a function basis on the whole real line. When one uses them in the approximation of a function defined on a bounded interval by simply taking the function values outside the interval as zero, some instability problems will arise^{[36, 39]}. Therefore, we need to do extra treatments to avoid this drawback.

Considering a bounded function g(x) defined on the interval [0, 1], by applying the power series expansion at x=0 and 1, we have

(4)

in which the parameters Δ₀ ≥ 2^-j(3N-M₁ -2) and Δ₁ ≥ 2^-j(M₁ -1). In order to determine the coefficients d_{0, i} and d_{1, i}, we assign, respectively, x=k/2^j (k=0, 1, …, M) for the first equation in Eq. (4) and x=1-k/2^j for the last one. Then, we have

(5)

where

and the subscripts k, l=0, 1, …, M. From Eq. (5), we have

(6)

in which the coefficients p_{0, i, k} and p_{1, i, k} are determined by

Substituting Eq. (6) into Eq. (4), we have

(7)

where

(8)

In the present study, we choose M=7 for Eqs. (5)-(8).

When a bounded function f(x) is approximated on the interval [0, 1] by use of the Coiflet-type scaling function with the compact support [0, 3N-1], Eq. (2) can be rewritten as follows:

(9)

Substituting Eq. (7) into Eq. (9) yields

(10)

where

(11)

Equation (10) provides a wavelet approximation for a function defined on a bounded interval, in which the expansion coefficients are just the function samplings at each nodal point. Since Eq. (10) is valid for all the bounded functions g(x) defined on the interval [0, 1], for the arbitrary nonlinear functional N satisfying that N[g(x)] is bounded on [0, 1], we have

(12)

whose accuracy can also be determined by Eq. (3).

3 Solution of the NLS equation

At first, let

where u₁ (x, t) and u₂ (x, t) are real-valued functions. Then, the one-dimensional generalized NLS equation (see Eq. (1)) can be rewritten as a system of two coupled real differential equations as follows:

(13)

where n=1, 2, and m=3-n.

Following the wavelet approximations (10) and (12) and using the boundary conditions in Eq. (1), we can write the unknown function u_n(x, t), nonlinear terms f(u_n² +u_m²)u_m and g(u_n² +u_m²)u_n, and the term ν (x, t)u_m in Eq. (13) as follows:

(14)

(15)

(16)

(17)

Substituting Eqs. (14)-(17) into Eq. (13) yields

(18)

Multiplying both sides of Eq. (18) by ϕ_{j, l}(x) (l=1, 2, … 2^j-1) and performing integration over the interval [0, 1], we have

(19)

(20)

where

and the vectors C_n(t) and D_n(t) depending on the boundary conditions can be expressed by

(21)

(22)

In Eq. (20), the operational rules

will always hold for all the vectors

The generalized connection coefficient can be exactly and readily obtained without the numerical integral but based on the database independent of the specific problems (see Ref. [38]).

Equation (19) can be solved directly by use of various time integration schemes, e.g., the Runge-Kutta method, the Adams-Bashforth scheme, and the Crank-Nicolson method. In this study, the classical fourth-order explicit Runge-Kutta method is used. With the method, we have

(23)

where t_n=nΔt≥ 0, and Δt is the time step.

Iteratively using Eq. (23) and directly using the initial condition given in Eq. (1), we can obtain the unknown vectors U(t_n)={U₁(t_n), U₂(t_n)}^T at each time step, which can be used to reconstruct the unknown function ψ (x, t) in terms of Eq. (14).

For the computational complexity of the present method, we first discuss the characteristics of the matrices A and B defined by Eq. (20). Considering the Coiflet scaling function φ(x) with the compact support [0, 17], we have . Therefore, A and B are banded matrices with the bandwidth b=33 for the decomposition level j>5^[38]. Moreover, because the Coiflet scaling function φ(x) holds the orthogonality, based on the definition (11), we have

for 8≤l≤2^j-10 or 8≤k≤2^j-10^[38]. Therefore, when the decomposition level j is larger than 5, the matrix A can be expressed as follows:

(24)

where the matrices P and Q are 7 × 7 and 9 × 9 dimensional square matrices, respectively, and I denotes the (2^j-17) × (2^j-17) dimensional identity matrix. It can be seen from Eq. (24) that the computational cost for estimating the inverse matrix A^-1 of the matrix A is extremely low, and the inverse matrix A^-1 is almost a diagonal matrix (see Eq. (24)). Following the above analysis and the procedure for calculating the value of Γ_{l, k}^{j, m[38]}, we can conclude that the calculated amount for calculating the matrix A^-1B used in Eq. (20) is about 200N_s simple arithmetic operations, where N_s=2^j-1 is the number of the grid points in space. Besides, we can also find that the matrix A^-1B is nearly a banded matrix with the bandwidth b=33 for the decomposition level j>5.

Considering the above characteristics of the matrices A^-1 and A^-1B, we can obtain readily that about 100N_s simple arithmetic operations (about 50N_s multiplications and 50N_s subtractions or additions) are needed for calculating the vector H_n(U₁(t), U₂(t)) defined by Eq. (20), and about 2N_s values of functions (about N_s values are needed to be calculated for each function of f(ρ) and g(ρ)). Therefore, about 200N_s simple arithmetic operations and 2N_s functions are needed for estimating the vector H(U(t)) defined by Eq. (19). Finally, we can obtain from Eq. (23) that the computational complexity of the present method at each time step is about 800N_s simple arithmetic operations, and about 4N_s values of each function of f(ρ) and g(ρ) are needed to be calculated for the decomposition level j>5.

4 Numerical examples

In the following, we solve numerically eight widely considered examples to demonstrate that the proposed wavelet method is very efficient and accurate for the solution of one-dimensional generalized NLS equations. To effectively evaluate the accuracy of the numerical solutions, we consider the error norms E_max and E₂ and the corresponding convergence rates R_max and R₂, which are, respectively, specified as follows^{[9, 13]}:

(25)

(26)

(27)

where ψ_num and ψ_exact are the numerical and exact solutions, respectively, N_s =2^j-1 is the number of the grid points in space, and ν = max or 2. We note that when the time step Δt is small enough to guarantee that the errors mainly come from the spatial discretization, the convergence rates R_max and R₂ defined by Eq. (27) are exactly the spatial convergence rates of the numerical method^{[9, 13]}. In the following, all the presented wavelet solutions are obtained by use of such an enough small time step Δt. Accordingly, the convergence rates R_max and R₂ estimated by Eq. (27) are the spatial convergence rates of the proposed wavelet method.

Example 1  We consider the undamped linear Schrödinger equation as follows^[8]:

(28)

whose exact solution is

Table 1 shows the error norm E_max and convergence rate R_max of the numerical solutions at the time T=0.5, which are obtained, respectively, by use of the present wavelet method and the compact finite difference scheme (CFDS)^[8]. Figure 1 displays the conservation errors of the mass E_Q =|Q(t)-Q(0)| and the energy E_E =|E(t)-E(0)|, where Q(t) and E(t) are, respectively, the total mass and the energy estimated by the numerical solution^[12].

Example 2  Consider the Schrödinger equation with cubic nonlinearity as follows^[8-18]:

(29)

with the Dirichlet boundary conditions extracted from the exact solution

(11)

In Table 2, we list the error norm E_max of the numerical solutions at different time T, which are achieved, respectively, from the proposed wavelet formulation, CFDS^[8], exponential spline solution (ESS)^[9], Crank-Nicolson methods (CNM)^{[8, 10]}, and the conservative finite difference method (CFDM)^[11]. Figures 2 and 3 plot the error norms of the present wavelet solutions with Δt=0.000 1 as a function of the number of the grid points in space N_s and those obtained with the ESS^[9], compact split-step finite difference method (CSSFDM)^[12], finite volume method (FVM)^[13], fourth-order compact difference scheme (FOCDS)^[14], split-step finite difference method (SSFDM)^[15], energy conservative difference scheme (ECDS)^[16], Crank-Nicolson difference scheme (CNDS)^[17], and linearly implicit conservative difference method (LICDM)^[18].

Example 3  Consider the following NLS equation^{[9, 12]}:

(30)

subject to the Dirichlet boundary conditions extracted from the exact solution

(31)

Table 3 gives a comparison of the error norms between the present wavelet solutions and those obtained, respectively, by use of the ESS^[9] and CSSFDM^[12] at different time T. The relation between the error norms of the wavelet solution with the time step Δt=0.000 1 and the number of the grid points in space N_s is shown in Fig. 4.

Example 4  Consider the following one-dimensional Gross-Pitaevskii equation^{[9, 12, 15, 21-22]}:

(32)

The exact solution is ψ (x, t)=sin(x)exp(-3it/2). Table 4 shows the error norms and convergence rates of the numerical solutions at time T=1, which are achieved, respectively, from the proposed wavelet method, ESS^[9], and CSSFDM^[12]. Figure 5 shows the error norm E_max of the present wavelet solution and that given by the SSFDM^[15] for various numbers of grid points in space N_s. Moreover, the conservation errors of mass and energy are presented in Fig. 6.

Example 5  Consider the following NLS equation with variable coefficients^{[9, 12]}:

(33)

with the initial and Dirichlet boundary conditions extracted from the exact solution

(34)

In Table 5, we make a comparison of error norms between the present wavelet solutions and those achieved, respectively, from the ESS^[9] and CSSFDM^[12] for various T. The relation between the error norms of the wavelet solution at various T with the time step Δt=0.000 1 and the number of grid points in space N_s is presented in Fig. 7.

Example 6  Consider the following NLS equation^[19]:

(35)

subject to the Dirichlet boundary conditions extracted from the solitary wave solution ψ (x, t)=exp(i(x-t)-(x-8-t)²). The error norms of the present wavelet solution with Δt=0.001 at different T for various N_s are shown in Fig. 8. The error norms achieved, respectively, from the SSFDM^{[15, 19]}, time-splitting spectral method (TSSM)^{[19, 25]}, and split-step orthogonal spline collocation method (SSOSCM)^[19] at T=1, N_s=128, and Δt=0.001 are 1.90×10^-2, 6.02×10^-4, and 2.55×10^-4, respectively.

Example 7  We consider the following cubic NLS equation with a periodic potential^[14]:

(36)

where k=4. Its exact solution is ψ (x, t)=exp(sin (2kx)-2ik²t). Table 6 shows the error norms E_max of the numerical solutions at T=1, which are obtained, respectively, by use of the proposed wavelet method, FOCDS^[14], and time-splitting Fourier pseudospectral method (TSFPM)^[14].

Example 8  Finally, we consider the damped NLS equation with cubic and quintic nonlinearities as follows:

(37)

The exact solution of Eq. (37) is ψ (x, t)=exp(2iπx-t). In Fig. 9, we plot the error norms of the wavelet solution at different T with the time step Δt=0.000 01 as a function of the number of the grid points in space N_s.

Examples 1-8 indicate that the present wavelet procedure has a good capability in solving generalized NLS equations. The results in Tables 1-6 and Figs. 2-9 demonstrate that the solutions achieved from the proposed wavelet method with coarser mesh have a much better numerical accuracy than those obtained by other existing numerical methods, e.g., the ESS^[9], FVM^[13], TSSM^{[19, 25]}, SSOSCM^[19], and various finite difference schemes^{[8, 10-12, 14-18]}. We also can find out that the present wavelet method has a high-order of convergence in space for generalized NLS equations, which obviously exceeds those of almost all existing numerical methods^{[8-19, 25]}. The order of convergence in space of the proposed wavelet formulation is about 9.0 for Eq. (29), and about 8.0 for Eqs. (30), (32), and (35). Moreover, it can be seen from Figs. 1 and 6 that the proposed wavelet method with a sufficiently fine time-space mesh can effectively conserve the total energy and the total mass of the Schrödinger systems (28) and (32), which hold the conservation laws of energy and mass.

5 Conclusions

In this paper, a sampling approximation is introduced for a function defined on a bounded interval by combining the Coiflet-type wavelet expansion and the boundary extension technique. Then, based on such a wavelet approximation, a high-order accurate wavelet formulation is proposed to solve uniformly the one-dimensional generalized NLS equation with damping. By studying eight widely considered test problems, the present wavelet algorithm is improved to have considerably high precision and fast convergence rate in space compared with many other existing numerical methods.