A new implicit symmetric method of sixth algebraic order with vanished phase-lag and its first derivative for solving Schrödinger's equation

Abstract: In this article, we have developed an implicit symmetric four-step method of sixth algebraic order with vanished phase-lag and its first derivative. The error and stability analysis of this method are investigated, and its efficiency is tested by solving efficiently the one-dimensional time-independent Schrödinger’s equation. The method performance is compared with other methods in the literature. It is found that for this problem the new method performs better than the compared methods.


Introduction
The numerical solution of second-order initial-value problems with periodical and/or oscillatory solutions as: , and , has attracted the attention of many authors in the last decades . The aim of these studies was the production of efficient, fast, and reliable algorithms for solving this model. These algorithms are of two main types, namely, algorithms with constant coefficients as given in Jain et al. [9] and Steifel and Bettis [10], and other ones with variable coefficients depending on the problem frequency as presented in [18][19][20][21][22][23][24][25][26][27][28][29][30][31]. In practice, equation (1.1) has been used to present mathematical models in several disciplines, such as, chemistry, quantum chemistry, physics, quantum mechanics, etc. The one-dimensional time-independent Schrödinger equation: is an example. Usually, the solution of equation (1.2) is studied under two boundary conditions, the first one is ( ) = y 0 0, and a second condition corresponds to large values of φ, which can be adopted due to some physical considerations. Here, τ 2 is a real number denoting the energy, l is an integer representing the angular momentum, and V is the potential function which satisfies ( ) → V φ 0 as → ∞ φ .
The aim of this paper is to develop a new efficient numerical method of the second type for solving the one-dimensional time-independent Schrödinger's equation. Our approach is based on the symmetric multi-step methods introduced by Quinlan and Tremaine [3], as well as the recent methodology for developing numerical methods, for approximating the periodical solutions of certain initial-value problems, which requires nullifying the method phase-lag and some of its successive derivatives. Following this methodology, we will develop an implicit four-step symmetric numerical method with vanished phaselag and its first derivative for solving equation (1.2). The rest of this article is organized as the following. In Section 2, we present the basic theory from the literature. In Section 3, we present the method development. In Section 4, we provide the error and the stability analysis of the new method. In Section 5, we present the application of the method and the numerical results. Finally, in Section 6, we give some conclusions.

Preliminaries
A numerical solution of a mathematical model of type (1.1) can be obtained using multi-step methods as: Usually, an operator is associated with the linear multi-step method (2.1); it is given by: where ψ is a two times continuously differentiable function.
Without loss of generality, as scaling equation (2.1) by a nonzero constant will not affect the forthcoming results, we may take = a 1 m ; so that equation (2.1) takes the form of a m 2 -step symmetric method as: gives the difference equation: Lambert and Watson have defined the interval of periodicity of a symmetric linear multi-step method as follows: Definition 2.2. [32] We say that a symmetric k 2 -step method has an interval of periodicity ( ) v 0, 0 2 whenever the roots ℓ w , ℓ = … k 1, 2, , 2 of the characteristic equation (2.6) satisfy the conditions: , a n d 1 , f o r 3 , , 2, The phase-lag of a symmetric multi-step method is defined as: 33] Each symmetric multi-step method with characteristic equation as (2.6) has a phaselag that equals to the first non-vanishing term in the series expansion of: We say that the phase-lag of a given multi-step method has an order κ if ζ is ( ) Moreover, the phase-lag of a given symmetric method with k 2 -steps can be found using the formula in the following result: Any linear symmetric k 2 -step method with characteristic equation as (2.6) has a phase-lag order κ and a phase-lag constant ρ satisfying the relation:

The method development
In this paper, we consider an implicit four-step symmetric method on the form: Using = k 2 in Theorem 2.4 implies the following proposition: Proposition 3.1. The implicit linear multi-step method defined by (3.1) has a phase-lag order s and a phaselag constant r given as: In the method (3.1), we choose the parameters μ 0 , μ 1 , and μ 2 so that this method is consistent and has the highest algebraic order, while the coefficients ν 0 , ν 1 , and ν 2 are considered as free parameters. Thus, 15 16 , and = − μ 0 1 8 . Then, using these values in equation (3.1) yields the following numerical method: Hence, in view of (3.2), the method (3.3) has a phase-lag given by: The following conditions are imposed on the new method: which will be used to find the values of the free parameters in equation (3.3). Thus, solving the system (3.5) for the parameters ν 0 , ν 1 , and ν 2 , we obtain: where  As it appears in (3.6), the expressions of ν 0 , ν 1 , and ν 2 involve the frequency = v ϕh. Moreover, the numerators and denominators in these expressions approach zero as → v 0. To avoid the arising huge round-off errors due to this situation, the following series expansions of these expressions are used:    4 The method error and stability analysis

Comparative error analysis
The local truncation error LTer associated with the new method, denoted as NewM, is obtained by using the free parameters ν 0 , ν 1 , and ν 2 given by (3.7) into equation (3.3), then expanding both sides in Taylor series, which implies: To investigate the LTer, presented in equation (4.1), we apply our method to the radial time-independent Schrödinger's equation (1.2) which is in the form: where the function p can be set as ( [46]): , and E is the energy.
Next, the truncation error formula (4.1) is applied to the above test equation to compare the behavior of the error LTer given in (4.1) with the behavior of other error formulas corresponding to truncation errors associated with similar methods in the literature. Our comparison will be based on the following methods: • The classical method; that is, the standard method with constant coefficients obtained from     Then, using the expressions of ( ) q n j , = … j 2, 4, , 8 in the above error formulas of the above methods gives the following asymptotic expansions for the LTer corresponding to each one of these numerical methods: • The classical scheme (3.1), i.e., the numerical method with constant coefficients obtained by setting = v 0 in equation (3.3): • The explicit four-step method in [34]:  We see from these equations that the expressions of these truncation errors involve powers of the expression = − G V E c . Therefore, we should consider the following cases: , when the values of the energy E and the potential V c are close to each other. In this case, in comparing the expressions of the truncation errors in these methods, we need to check only the terms free of G. Hence, it follows that all these methods are of comparable accuracy, as these terms are identical in all these methods.
, when the values of the energy E and the potential V c are significantly different from each other, that is, | | G is large.
In this case, we see from the above asymptotic error expressions that for the classical method the asymptotic error increases as the third power of G, whereas it increases as the first power of G for the other methods. But, for the new method, G has the lowest coefficients among all the above methods. Thus, for the numerical solution of the radial time-independent Schrödinger's equation with large values of | | G , it seems to be that, among these numerical techniques, the new derived one is the most efficient.

Stability analysis
To study the stability analysis of the new numerical method, we substitute the values of the parameters ν 0 , ν 1 , and ν 2 given in equation (3.6) in the method (3.3), then apply the resulting scheme to the test equation ″ = − y θ y 2 , which produces a difference equation associated with the following characteristic equation:  Figure 4, in which the shaded parts represent the region on which the method is stable. Hence, in view of this region it follows that the periodicity interval of this new numerical scheme is (0, 10.5177).

Numerical results and discussion
To examine the efficiency of the constructed method, we apply it to solve the radial time-independent Schrödinger equation: accompanied with two conditions:

2)
and another boundary condition will be adapted to certain physical considerations for some large value of φ. Here, γ 2 in equation (5.1) is a constant that stands for the energy, l is the angular momentum, and V represents the potential which satisfies ( ) → V φ 0 whenever → ∞ φ . The new constructed method involves non-constant parameters, as they depend on the frequency of the studied model. Thus, in order to apply this method to solve equation (5.1), first of all we need to specify a value for the frequency ϕ. Therefore, we consider (5.1) with = l 0, and choose ϕ as: where ≔ E γ 2 is the energy, and use the Woods-Saxon potential: and = a 0.6. The behavior of ( ) V φ is shown in Figure 5. In Woods-Saxon potential, the frequency ϕ is not a function in φ, instead it is approximated by estimating the potential ( ) V φ at the certain potential critical points [46]. For the numerical tests, we take ϕ as follows (see [46]) , 0, 6.5 ,  where δ l represents the phase shift that can be determined by the formula: where + φ i 1 and φ i are two successive points in the asymptotic region, with ( ) = ( ) S φ γφj γφ l and ( ) = ( ) C φ γφN γφ l . Now, given the energy E, we will employ the new method (3.3) to solve the resonance problem; that is, to estimate the phase shift δ l which has exact value that equals to π 2 . But, first we need to determine the start values p i for = … i 0, 1, , 3. Thus, in view of the condition (5.2) we get = p 0 0 , and the remaining other values can be computed using Runge-Kutta-Nyström methods (see [45] and [42]). Starting with these values, we will estimate the value of δ l at the point + φ i 1 using the following numerical schemes: • The explicit 8-step eighth order linear symmetric method with constant coefficients given in [3]; denote it by ( ) QT 8 .
• The explicit 10-step tenth order linear symmetric method with constant coefficients developed in [3]; denote it as ( ) QT 10 .
• The explicit 12-step twelfth order linear symmetric method with constant coefficients developed in [3]; denote it by ( ) QT 12 .
• The explicit fourth order 4-step linear symmetric method with vanished phase-lag and its first derivative given in [35]; denote it as SIMI. • The explicit fourth order 4-step linear symmetric method with vanished phase-lag and its first and second derivatives constructed in [36]; denote it by SIM. • The explicit fourth order 4-step symmetric method with vanished phase-lag and its first and second derivatives given in [37]; denote it as SIMN. • The eighth order two-stage two-step linear symmetric method with vanished phase-lag and its first, second, third, and fourth derivatives obtained in [47]; denote it by S St 2 . • The eighth order three stage two-step linear symmetric method with vanished phase-lag and its first, second, third, and fourth derivatives developed in [48]; denote it as S St 3 . • The three-stage fifth order two-derivative Runge-Kutta method, developed in [43]; denote it by Y D 2 . • The fourth order three derivative Runge-Kutta method with vanished phase-lag and its first derivative produced in [44]; denote it as Y D 3 . • The implicit sixth order 4-step linear symmetric method with vanished phase-lag and its first derivative in [39]; denote it by SHIM. • The explicit fourth order 4-step linear symmetric method with vanished phase-lag and its first derivative developed in [38]; denote it by EFSM. is graphed versus the required consumed CPU time (in seconds) corresponding to the energies E 1 and E 2 , as shown in Figures 6 and 7, respectively.

Conclusions
The results in Figures 6 and 7 show that the efficiency of the new developed implicit method (NewM) in solving Schrödinger equation with different energy values is confirmed, as it gives the highest accuracy among the other compared methods in less CPU time.