Numerical Solution of Nonlinear Schrödinger Equation with Neumann Boundary Conditions using Quintic B-spline Galerkin Method

This paper is concerned with the numerical solution of the nonlinear Schrödinger (NLS) equation with Neumann boundary conditions by quintic B-spline Galerkin finite element method as the shape and weight functions over the finite domain. The Galerkin B-spline method is more efficient and simpler than the general Galerkin finite element method. For the Galerkin B-spline method, the Crank Nicolson and finite difference schemes are applied for nodal parameters and for time integration. Two numerical problems are discussed to demonstrate the accuracy and feasibility of the proposed method. The error norms , and conservation laws , are calculated to check the accuracy and feasibility of the method. The results of the scheme are compared with previously obtained approximate solutions and are found to be in good agreement.


Introduction and Governing Equation
In this article, quintic B-spline Galerkin finite element method is applied to find the numerical solution of nonlinear Schrödinger (NLS) equation: Equation ( 1) is called a self-focusing NLS equation (α > 0) and allows for bright soliton solutions, as well as the defocusing NLS equation (α < 0).u = u(x, t) is a complex-valued function over the real line, α is a positive number and i = √ −1.The initial and boundary conditions are as follows: u(a, t) = u(b, t) = u x (a, t) = u x (b, t) = u xx (a, t) = u xx (b, t) = 0. ( where r(x, t) and s(x, t) are real functions.Substituting Equation (4) into Equation (1), we obtain the coupled partial differential equations Finite element method is a powerful and established method used to approximate the solution of the partial differential equation.Besides finite element method, splines are also very useful to approximate the solution of partial differential equation with piecewise polynomial approximation.A finite element method with B-splines defines a new weighting approximate method and possesses computational advantages of B-splines and finite elements.Spline functions have been applied to develop numerical methods for the solution of nonlinear differential equations [1,2].The analytical solution of the NLS equation is solved by using the inverse scattering method by Zakharov and Shabat [3].In 1974, Zakharov and Manakov proved that the NLS equation is completely integrable [4].Many researchers have worked on the solution of the partial differential equations by using collocation finite element method based on splines.
Spline-based numerical methods have been proposed by many researchers to obtain the numerical solution of nonlinear evolutionary problems.Mittal [5] obtained the numerical solutions of the extended Fisher-Kolmogorov equation using the quintic B-spline collocation method.Saka and Dag [6] presented the Galerkin finite element method based on quartic B-spline functions to obtain the numerical solution of the regularized long-wave (RLW) equation.Gardner [7] proposed the Galerkin finite element method based on cubic B-spline to find the numerical solution of the RLW equation.Dogan [8] presented the Galerkin method based on linear space finite elements to the numerical solution of the RLW equation.Kutluay and Ucar [9] and Saka and Dag [10] found the numerical solution of the coupled Korteweg-de-Vries (KdV) and Korteweg-de-Vries-Burgers (KdVB) equations using the Galerkin finite element method based on quadratic and quartic B-spline functions, respectively.
Gorgulu et al. [11] used exponential B-splines Galerkin finite element method for solving the advection-diffusion equation.They developed a new algorithm by incorporating exponential B-spline functions with the Galerkin finite element method.This method gives satisfactory results.The exponential B-splines Galerkin method is also applied to solve the Burger's equation and the results are comparable with the quartic B-spline collocation method [2].
There are many non-spline numerical methods developed to solve the NLS equation; some of them are discussed here.Wang et al. [12] proposed a finite difference method using an artificial boundary conditions on an unbounded domain.In this scheme, extrapolation operator is applied to deal with the nonlinear term.Moreover, Barletti et al. [13] presented energy-conserving methods that can confer robustness on the numerical solution.Taleei and Dehghan [14] presented a time-splitting pseudo-spectral domain decomposition method, whereby the original equation is split into linear and nonlinear equations.The Chebyshev pseudo-spectral collocation method is used to solve the linear equation in the spatial dimension and Crank-Nicolson scheme in the temporal dimension.In this study, overlapping multi-domain scheme is chosen.Univariate multi-quadrics (MQ) quasi-interpolation method is developed where the spatial derivatives are calculated from the derivative of the quasi-interpolation.Some spline-based numerical methods are proposed to solve the NLS equation.Quartic spline approximation and semidiscretization were applied using finite difference method by Sheng et al. [15].Zlotnik and Zlotnik [16] were the first to implement the finite element method using non-discrete transparent boundary conditions.Naturally, higher-order finite element method is observed to converge faster.The exponential B-spline with collocation method was presented by Ersoy et al. [17].The Crank-Nicolson scheme is used for time integration and exponential cubic B spline functions for the space integration.Dag [18] proposed a Galerkin finite element method based on quadratic B-spline functions.
In this paper, numerical scheme for the Galerkin method with quintic B-spline is developed to solve the NLS equation.Numerical results are generated and compared with some of the afore-mentioned methods.
This paper is organized as follows.In Section 2, the fundamentals of quintic B-spline Galerkin method are introduced.In Section 3, the initial parameters to find the solutions of the system are calculated.Numerical results and test problems are discussed in Section 4 followed by the conclusion in Section 5.

Quintic B-spline Galerkin Method
We consider a mesh Π over the finite domain [a, b] divided uniformly by grid points x k with h = x k − x k−1 , k = 1, . . ., N. Quintic B-spline function is chosen as the weight and trial function.The quintic B-splines, B k (x), k = −2, . . ., N + 2 at the grid points x k forms a basis over the interval [a, b] as follows [19]: The global approximate solution, u N (x, t), for the NLS equation in Equation ( 1) is written in terms of the quintic B-spline function as where The functions ρ k (t) and ψ k (t) are time-dependent parameters that are determined from the boundary and residual conditions.We make use of a local coordinate transformation The quintic B-spline shape functions in Equation ( 6) can be defined in term of η, Since all other quintic B-spline functions are zero over the interval [x k , x k+1 ] except for B k−2 , B k−1 , . . . ,B k+3 , the approximation function in Equation (8) over the typical interval [x k , x k+1 ] can be written as Applying the quintic B-splines definition in Equations ( 6)- (8), the nodal values of s k , r k , and its first and second derivatives at the knots x k are found to be When using the Galerkin method on Equation ( 5) with weight function W(x), the weak form of Equation ( 5) over the finite interval where By using the weight function W(x) as quintic B-spline shape functions B k and inserting the quantities r N (x, t) and s N (x, t) in Equation ( 10) into the integral Equation ( 12) instead of r(x, t) and s(x, t) k + 3 and the dot "•" represents the derivative with respect to time t.The finite element in Equation ( 14) can be written in matrix form as where ) T are the element parameters.The element matrices A e ij , M e ij and C e ij are given by integrals where The element matrices in Equation ( 16) are calculated as The values of z L are obtained by writing s = s k +s k+1 2 and r = in Equation (13).Using the values of s N and r N at the grid points x k , we obtain Assembling all the elements in Equation ( 15) leads to the following system: where ρ = (ρ −2 , ρ −1 , . . ., ρ k , ρ k+1 ) T and ψ = (ψ −2 , ψ −1 , . . ., ψ k , ψ k+1 ) T are the global element parameters, and A, M and C(z L ) are the global matrices with generalized kth row given by: where Substituting the time derivatives of the parameters  into Equation (18) yields Equation (21).

The Initial Vectors
The initial parameters ρ 0 and ψ 0 can be determined using initial and boundary conditions to solve the system in Equation (21).The approximation can be re-written over the interval [a, b] at t = 0 as where all parameters ρ 0 and ψ 0 are determined.The functions s N and r N are required to satisfy the following relations at the nodal points x k : A pentadiagonal matrix system is obtained using the conditions in Equation ( 24).The initial vectors ρ 0 and ψ 0 can be calculated from the following matrix equations: The above system can be solved by inverting the matrices in MATLAB.The approximate solution of s N (x, t) and r N (x, t) can be calculated from ρ n and ψ n using Equation (21).

Numerical Experiments and Results
Two physical problems, single solitary wave and interaction of two solitary waves, were considered to assess the performance of the proposed method summarized in Equation (22).The performance and accuracy of the approach were tested using the L 2 and L ∞ norms defined as and where u exact and u N denote the exact and numerical solutions, respectively.Moreover, Equation (1) must satisfy the two conservation laws,

Problem 4.1 (Single Solitary Wave Solution)
A single solitary wave solution to the NLS equation is given as in [20]: The initial and boundary conditions are taken as u(−20, The values of the initial parameters from Equations ( 25) and ( 26) are calculated by using the initial and boundary conditions.The values of all parameters were chosen to be α = 2, S = 4, β = 1, h = 0.05, and ∆t = 0.002, 0.001.The parameter S represents the speed of the solitary wave whose magnitude depends on the real parameter β.The L 2 and L ∞ norms and conservation laws I 1 and I 2 are tabulated in Table 1 for ∆t = 0.002 and Table 2 for ∆t = 0.001.It is observed that the norms remain very small.The numerical results obtained by the present scheme are more accurate than the explicit, implicit, and split-step Fourier and other methods [19,21,22].The norms naturally decrease with the increase in number of partitions.We found a good result even for large step size, as displayed in Table 3.The L ∞ and L 2 norms converge to zero as the number of nodes increases.The numerical simulations are shown at different times over the region [−20, 20] in Figure 1.Furthermore, numerical results for different values of ∆ were generated to calculate the rate of convergence using the following formula [23]: where (2ℎ, 2∆) is either the  error norm or the  error norm in spatial and temporal directions.The error norms  ,  and order of convergence rate at time  = 1 is shown in Table 4.In Table 4, we see that this method is nearly of second-order convergence.Furthermore, numerical results for different values of ∆t were generated to calculate the rate of convergence using the following formula [23]: where E(2h, 2∆t) is either the L ∞ error norm or the L 2 error norm in spatial and temporal directions.The error norms L ∞ , L 2 and order of convergence rate at time t = 1 is shown in Table 4.In Table 4, we see that this method is nearly of second-order convergence.Table 5 presents comparison between our proposed method with other methods.The proposed scheme gave satisfactory results.Table 6 displays numerical results of present method compared to those of Taha et al. [21].In general, the present method generated more accurate results for the specific values of parameters.The simulation of the single soliton with amplitude equal to 2 is presented in Figure 2.   Table 6 displays numerical results of present method compared to those of Taha et al. [21].In general, the present method generated more accurate results for the specific values of parameters.The simulation of the single soliton with amplitude equal to 2 is presented in Figure 2.

Problem 4.2 (The Interaction of Two Solitary Waves)
In this problem we considered the behavior of two solitons moving in opposite directions using the following initial conditions [18,20,21]:

Problem 4.2 (The Interaction of Two Solitary Waves)
In this problem we considered the behavior of two solitons moving in opposite directions using the following initial conditions [18,20,21]: where β k , α and x k are constants.
The values of all parameters were chosen to be 05 and ∆t = 0.05.Two solitary waves were traveling in opposite direction with the same magnitude 1 and speed 10.One of the solitary waves placed at x = 10 was traveling to the left side with speed 4 and the second wave on the other side placed at x = −10 was traveling to the right side with speed 4. As we know, solitons move in opposite direction and they collide and separate.We noticed that the shape of the solitons did not change after the collision of both solitons as expected.The two solitary waves collided at times t = 1, 2, 3 and then separated at t = 4, 5, 6.We see that the solitons preserved the original shapes after the collision.The interactions of two solitons at different times t = 0, 1, 2, 2.8, 3, 4, 5 and 6 can be seen in Figure 3.
The L ∞ and L 2 norms and I 1 conservation law were calculated at various times with ∆t = 0.005, 0.002, as shown in Tables 7 and 8.  Furthermore, Table 9 presents numerical results of our proposed method compared with previous methods (e.g., [21]).The present method gives accurate results for the specific values of parameters.Furthermore, Table 9 presents numerical of our proposed method compared with previous methods (e.g., [21]).The present method gives accurate results for the specific values of parameters.

Conclusions
In this paper, the numerical solution of the NLS equation with Neumann boundary conditions is obtained by the Galerkin finite element method with quintic B-spline shape function.The accuracy and feasibility of the method was evaluated by two test problems related to single solitary wave and interaction of two solitary waves.The accuracy of numerical method was examined by showing reasonably small error norms L 2 and L ∞ .The interaction of two solitary waves was investigated and observed that the shape of the solitons did not change after the collision of both solitons as expected.The proposed method successfully simulated the soliton picture by choosing different parameters in the case of motion of single soliton and interaction of two soliton.The obtained results were compared with published results [17,19,21,22] and it was observed that the all results are acceptable and reflect the analytical solution.The rate of convergence was calculated and found to be almost of second-order convergence.In conclusion, the present Galerkin finite element scheme with quintic B-spline presents an acceptable soliton solution method for solving NLS equation.The simplicity of the quintic B-spline Galerkin finite element method is an advantage over the general Galerkin finite element method to obtain the numerical solution of NLS equation.The proposed scheme can easily be applied to solve various nonlinear differential equations.

Table 7 .
Norms and conservation laws for Problem 4.2 with h = 0.05 and ∆t = 0.005.

Table 8 .
Norms and conservation laws for Problem 4.2 with h = 0.05 and ∆t = 0.002.