A Method of Lines Approach in the Numerical Solution of 1-Dimensional Schrödinger ’ s Equation

This paper discusses the numerical solution of the 1-dimensional Schrödinger equation using method of lines approach (MOL) where the spatial dimension is discretize using some finite difference approximation leaving the time dimension to be the only independent variable in the resulting system of initial value problems. The effect of changing the discretization size on the stability of the solution procedure with respect to the absolute stability property of the numerical method used has been studied. In the study the use of Simpson’s rule in MATLAB has also been incorporated. The result indicates that there is some relationship between the discretization size and the convergence property. Further work will look into the modeling this equation to the application in Physics.


The method of Lines (MOL)
The method of lines (MOL, NUMOL, NMOL) (Schiesser, 1990;Hamdi et al., 2009) is a procedure for solving partial differential equations (PDEs) in which all but one dimension is discretized (Schiesser, 1990; Wikipedia, the free encyclopedia, 04 November 2009; Cutlip & Shacham, 2001;Pregla, 2008;Shakeri & Dehghan, 2008).The MOL is most of the time refers to the construction or analysis of numerical PDEs that proceeds by first discretizing the spatial derivatives only -normally using finite difference method (FD) (Ismail et al., 2007) -and leaving the time variable continuous.This will lead to ordinary differential equation (ODE) to which numerical method for initial value ordinary equation can be applied (Wikipedia, the free encyclopedia, 04 November 2009).The MOL also has been applied to both Physical and Chemical problems (Dubey, 2011) especially in the theoretical physics (Pregla, 2008;Shakeri & Dehghan, 2008;Demirkaya, 2005).The idea was first applied by German mathematician -Erich Rothe in 1930 to parabolic type of equations.

What the Paper Covers
In this paper, the use of the Method of Lines (MOL) approach to solve the Schrödinger's equation was employed.The method involved replacing the spatial (boundary-value) derivatives in the PDE with algebraic approximations and by doing that, the spatial derivatives became no longer stated explicitly in terms of spatial independent variable.Time which is the initial-value variable in the problem remained.The Mathematical software (MATLAB) and library routine specifically used by Schiesser and Griffith (2009)  in order to see the effect computational time and accuracy of the result in each case.

Application of Schrödinger's Equation
The latest achievement of Bose -Einstein condensation using ultra cold neutral bosonic gases is one of the areas of interest on this equation (de la Hoz & Vadillo, 2008).The field of application varies from optics, propagation of electric field in optical fibers, self -focusing and self -modulation of trains in monochromatic wave in nonlinear optics, collapse of Langmuir waves in plasma Physics to modeling deep water freak waves (so -called rogue waves) in ocean (Zisowsky & Ehrhardt, 2008;Jime'nez et al., 2003;Ramos & Villatoro, 1994) such as modeling of 2004 Indonesian Tsunami (Ismail et al., 2007).
It is pertinent to note that, this paper is aimed at determining the effect of computational effort with respect to accuracy of results and to identify the level of relationship between discretization size and the convergence property.

Method
Basically, to apply the MOL normally involves the following stages; partitioning the solution region into regions or layers, discretizing the differential equation in one coordinate direction (mostly x axis), transformation to obtain decoupled ODE, inverse transformation and introduction of boundary conditions and Solution of the equation.
The Schrödinger equation is given by and normalizing  and 2m to unity, the preceding normalized nonlinear or cubic 1 dimensional Schrödinger equation is obtained (Schiesser & Griffith, 2009) Where u is complex dependent variable x is a boundary-value (spatial) independent variable t is initial-value independent variable q is arbitrary parameter (constant) and i is an imaginary complex number 1  .The term "Cubic" in CSE initiate its name from the term 2 0, q u  ( u represent the absolute value of complex variable u) Let's consider equation (1.3) above; to solve that problem involves the steps that were mentioned above.First of all the spatial derivative u xx has to be replaced with an algebraic approximation such of FD.But before this, for numerical treatment of nonlinear Schrödinger equation the complex u(x,t) must be replaced by its real part u(x,t) and imaginary part ( , ) Substituting the new value of u(x,t) into equation (1.3), the two equations were obtain; Equating real and imaginary part to zero, we obtain, the following nonlinear system is obtained Solving equations (2.2a) and (2.2b) we get; For the (2.2a), (2.3b) Applying the three point central difference scheme to the equations (2.3a) and (2.3b) we get from (2.3a), and from (2.3b), h is the spacing between the discretize line and n is the number of grids.
Substituting equations (2.4a) and (2.4b) in equations (2.3a) and (2.3b) we get; ) ) Now for the next step is introducing the IC and BCs; equation (1.3) is first order in t and second order in x, it therefore requires 1 (one) IC and 2 (two) BCs respectively.The IC is taken from the analytical solution (with q=1) (Schiesser & Griffith, 2009).
(0.5 0.75 ) ( , ) 2 sec ( ) (3.6a) Since it has already been mentioned above that the complex value u(x,t) need to be analyze as a two real functions, we can therefore say from equation (3.6a) ( , ) 2 sec ( 1) u x t h x   (3.6b) Equation (3.6b) will be use for comparison with analytical solution.Then with t=0 in equation (3.6a), the IC is As for the introduction of the two BCs, this is not going to be possible since for a class of rapidly decreasing initial conditions i.e |u(x,t)| 0  as |x|   that, if PDE is analyzed over an essential infinite domain x     and if changes in the solution occur only over a finite interval in x, the BCs at infinite have no effect; In other words, we do not have to actually specify BCs (since they have no effect) as it will be clarified through PDE solution (Schiesser & Griffith, 2009;Samrout, 2007).
Accordingly, equations (1.3) and (2.7) constitute the coupled PDE problem, equations (2.6a) and (2.6b) is used in programming and analysis to evaluate the numerical MOL solution.Note that equation (2.6a) and (2.6b) are traveling wave solution since the argument of sech function x-t (Schiesser & Griffith, 2009).

Results and Discussion
The maximum values plotted indicated that, the movement of the soliton is from left to right.See Figure 3.1 -3.3 for the summary of the peak movement at t =0, 5, …, 30 which corresponds to x = 0,5,…,30.
The three invariants which have to do with the degree of accuracy of MOL also can be calculated and their maximum variations can be shown thus; In the first program, the maximum variation is given as; [( 13.8865 ( 12.7830)] 100 7.9% ( 13.8865) While in the last two program, their maximum variations are just 7.6% and 6.7% respectively.This variation shows how effective and efficient a MOL is, and how it is being affected when increment in grids point is made (the high the grid point the more the variation).Furthermore, there is very good agreement between the numerical and analytical solutions as been justified by the error differences in the plotted output of the Figure 3.1-3.3.In the previous paragraph it was shown what this program is aimed to achieve; the detail about the output of the program including efficiency and limitation of MOL.In the subsequent paragraph the discussion on the second and the last segment of the program output will put in place, including the behavior of the computing time when the discretization size and step size were varied so also the concept of the accuracy concurrently.
Starting from Table 1, where we try to reduced and maintain the step size ( 0 30 t   ) was reduced and then maintained, the number of grids were then varied (n = 301, 251 and 201 respectively); at a small x  (n = 301) there is more accuracy of the result (7.0X10 -4 ), even though the CPU time is high (51.22secs.).When the discretization size was increased (number of grids were reduced to n = 250), the accuracy of the result tend to reduce also but the CPU time reduced to about 48% of the first computing time.As x  was finally increase (number of grids reduced to n = 200), the accuracy drastically reduced to 2.9X10

Conclusion
The use of the Method of Lines (MOL) approach to solve the Schrödinger's equation was employed.The method involved replacing the spatial (boundary-value) derivatives in the PDE with algebraic approximations and by doing that, the spatial derivatives became no longer stated explicitly in terms of spatial independent variable.
The paper described the effect of computational effort with respect to accuracy of results and identified the level of relationship between discretization size and the convergence property in the numerical solution of 1-dimentional Schrödinger's equation.
From the study, it can be shown that, increasing the number of grids (decreasing the discretization size x  ) and maintaining the step size will lead to the high accuracy of the result, when solving the Schrödinger equation using MOL approach.While, increasing the discretization size x  (reducing the number of grids) and maintaining or reducing the step size can only leads to faster computing effort with less accuracy of result.The study also shows that, even though there exist an agreement between numerical and analytical solution, the more the number of grids (decrease in the discretization size) the lesser the error difference between them and vise verse.Further work will look into the modeling this equation to the application in Physics.

Figure 3 .
Figure 3.1 shows partial output of the first program solution of CSE when the discretization size x  =0.167 (MOL grid has 301), spatial domain is 10 40 x    and step size is on the interval 0 30 t   , Figure 3.2 Partial output of the second program solution of CSE when the discretization size x  =0.2 (MOL grid has 251), spatial domain is 10 40 x    and step size is 0 30 t   and finally, Figure 3.3 Partial output of the last program solution of CSE when the discretization size x  =0.25 (MOL grid has 201), spatial domain is 10 40 x    and step size is on the interval 0 30 t   .

- 3 Table 1 .
CPU time and accuracy when x  is being varied and t