A New Explicit Four-Step Symmetric Method for Solving Schrödinger’s Equation

In this article we have developed a new explicit four-step linear method of fourth algebraic order with vanished phase-lag and its first derivative. The efficiency of the method is tested by solving effectively the one-dimensional time independent Schrödinger’s equation. The error and stability analysis are studied. Also, the new method is compared with other methods in the literature. It is found that this method is more efficient than these methods.


Introduction
In this article the numerical solution of certain second-order initial-value problems with periodical and/or oscillatory solutions as η (x) = f (x, η), η(x 0 ) = η 0 and η (x 0 ) = η 0 , is investigated. As it is seen in Equation (1), this model does not involve the first derivative explicitly. In practice, there are several mathematical models in the form of Equation (1), characterizing certain scientific problems in the fields of chemistry, quantum chemistry, physics, quantum mechanics, etc., [1]. In the literature, there are various algorithms addressing the numerical solution of boundary and initial value problems of this type. These methods can be classified in two main categories. The first category involves methods with constant coefficients, and the second one involves methods using variable coefficients. Methods of the first type are developed in Gautschi [2], Jain and et al. [3], and Steifel and Bettis in [4], while methods of the second kind are presented in Chawla and et al. [5,6], Dahlquist [7], Franco [8], Lambert and Watson [9], Krishnaiah [10], Simos and et al. [11,12], Saldanha and Achar [13], Achar [14], and Daele and Vanden Berghe [15].
The purpose of this article is to develop a new efficient method, with variable coefficients, for solving the time independent Schrödinger's equation. The methodology in developing this method is based on the symmetric multistep methods introduced by Quinlan and Tremaine [22]. Its efficiency in handling numerically initial value problems with oscillatory solutions comes from the requirements imposed on its phase-lag and some of its derivatives, as well as its high algebraic order.
The article is organized as follows: Section 2 involves the basic concepts from the literature required for developing the new method. In Section 3 we present the development methodology of this method. The analysis of the error and stability of the new method is considered in Section 4. Finally, in Section 5 we present the numerical results.

Preliminaries
Numerical solution of the initial value problem Equation (1), can be derived by using multistep methods with k steps as: using a set of equally spaced points {x i } k i=0 discretizing the interval [a, b] into k subintervals, with a step size h, where h = |x i+1 − x i |, for i = 0, 1, ..., k − 1. When k is an even integer, α ν = α k−ν and β ν = β k−ν for ν = 0, 1, ..., k 2 , then the method Equation (2) is said to be symmetric method. The linear multistep method given in Equation (2) is associated with an operator as where u is any two times continuously differentiable function.
Without loss of generality we may take α k = 1, hence Equation (2) can be written as Definition 1. [23] The multistep method Equation (2) is said to have an algebraic order q when the associated linear operator L vanishes for any algebraic polynomial in the linear span of {1, x, x 2 , ..., x q+1 }.
Definition 3. [24,25] The phase-lag, of any symmetric multistep method corresponding to the characteristic Equation (7), is defined as the first non vanishing term in the power series expansion of The phase-lag is said to have an order p whenever ξ is O(v p+1 ) as v → ∞. The phase-lag of any linear symmetric 2m−step method can be computed using the following theorem proved in [24]: Theorem 1. Any linear symmetric 2m-step method with characteristic equation given by Equation (7) has a phase-lag order q and a phase-lag constant c that satisfies

Method Development
In this article we consider an explicit linear symmetric four-step method as: where f n+i = y (x n+i , y n+i ), i = 0, ±1.
The following proposition follows radially from Theorem 1 with m = 2.

Proposition 1.
The linear symmetric four-step method Equation (11) has a phase-lag order q and a phase-lag constant c given by: where The parameters are chosen so that the method is consistent. Thus, setting α 1 = − 3 40 and α 0 = − 37 20 in Equation (11) we get the following method: Hence, in view of Equation (12), the method Equation (13) has a phase-lag given by: where The method in Equation (13) is subject to the following two conditions: Solving the linear system Equation (15) for the coefficients β 0 and β 1 we obtain As it is seen in Equation (16), the coefficients β i for i = 0, 1, depend on v = φ h. Moreover, these coefficients have indeterminate forms as v approaches zero, as their numerators and denominators vanish as v → 0. Hence, the division of such expressions will cause huge round-off errors. To overcome this situation, these coefficients are replaced by their Taylor series expansions, which are given by: The behavior of the coefficients β 0 and β 1 , given by Equation (16), is shown in Figures 1 and 2.

Comparative Error Analysis
The local truncation error, denoted by LTE, for the new developed method, named as NL4SM, is derived through substituting the values of the coefficients given by Equation (17) in Equation (13), then expanding both sides in Taylor series, to give: To investigate the behavior of the LTE, given by Equation (18), we will use radial time independent Schrödinger's equation in which the function η can be expressed as (see [17]) with is a potential function, V c is a constant approximating the potential at a certain point x, and G = V c − E, where E represents the energy. Next, we will compare the behavior of the LTE in Equation (18) when it is applied to the above test equation with other known methods of the same type in the literature. Namely, we consider the following methods:

•
The classical method, that is the method corresponding to Equation (13) when the coefficients are constants • The methods derived in [26]: n + 2 w 2 y • The newly developed method given by Equation (13) LTE NL4SM = 643 h 6 9600 y To determine the asymptotic expressions of the LTEs corresponding to these methods, we have computed the derivatives occurred in the expressions of these errors, using the equation y n = (ζ(x) + G) y(x), where y n represents an approximation of y(x) at x n , and we obtain Thus, substituting these derivatives in the above error formulas, lead to the following asymptotic expansions of the LTEs for the above mentioned methods: • The four-step method developed in [26] LTE L4SMSI = h 6 161 2400 (ζ(x)) 2 y(x) + 161 1200 • The four-step method developed in [26] • The new method given by Equation (11) As it is seen from the above equations the LTEs involve powers of G = V c − E. Thus, we will consider the following two cases: Case 1. G ≈ 0, i.e., the energy, E, is close to the potential V c . Therefore, in comparing the LTEs in these methods only the terms free of G are considered. Hence, for these terms all the methods are of comparable accuracy, as all their terms free of G are identical. Case 2. G 0 or G 0, that is |G| is considerably large.
As it is seen from the above equations, the asymptotic error in the case of the classical method increases as the third power of G, while in the last three methods, it increases as the first power of G, but in the new constructed method it has the lowest coefficients. Thus, for the numerical integration of the radial time independent Schrödinger's equation, it seems that the newly developed method is the most efficient, especially when |G| is large.

Stability Analysis
Substituting the values of the coefficients given by Equation (16) in the method Equation (13), then applying this method to the test problem y = −θ 2 y, leads to a difference equation associated with the following characteristic equation Let us mention that the frequency φ in the developed method Equation (13) is not equal to θ, which indicates the frequency of the above test problem. The stability region of the new method is presented in Figure 3. The shaded region represents the region where the method is stable, from which it follows that the method has a periodicity interval equal to (0, 32.1489).

Numerical Results and Discussion
To test the efficiency of the new method it is applied to solve numerically the radial time independent Schrödinger's equation: subject to two conditions: and one more boundary condition, imposed for large values of ξ, under certain physical considerations, where γ 2 is a real number denotes the energy, l is an integer represents the angular momentum, and V is the potential function which satisfies V(ξ) → 0 as ξ → ∞.
The coefficients in the new developed method are variables depending on the frequency of the problem under consideration. Therefore, we need to specify the value of the frequency φ. We will apply the new method to solve Equation (21) with l = 0, thus we take φ as: where E := γ 2 is the energy, and we will use the Woods-Saxon potential: where q = e ξ−7 a and a = 0.6. The behavior of V(ξ) is shown in Figure 4. In Woods-Saxon potential, the parameter φ is not presented as a function of ξ, instead it is estimated through approximating the potential, V(ξ), at certain critical points of this potential [17]. For our numerical tests we choose φ as follows (see [17]) For large values of ξ the potential V(ξ), in Equation (21), decays faster than l(l+1) ξ 2 , thus the radial Schrödinger equation is reduced to which has two linearly independent solutions, namely, γξ j l (γx) and γξ N l (γξ), where j l (γξ) and N l (γξ) are the well know Bessel and Neumann spherical functions, respectively. Hence, as ξ → ∞, Equation (21) has a solution in asymptotic form which is given by: where δ l indicates the phase shift, and it can be computed using the equation where ξ i+1 and ξ i are two different points in the asymptotic region, with S(ξ) = γ ξ j l (γ ξ) and C(ξ) = γ ξ N l (γ ξ). Given the energy, E, we will use the developed method to solve the resonance problem, that is to approximate the phase shift, δ l , which has the exact value π 2 . In order to start the method Equation (13) we need to determined the values of y j for j = 0, 1, ..., 3. From the condition Equation (22) we have y 0 = 0. The remaining values are computed using Runge-Kutta-Nyström methods [27,28]. Based on these values, we will compute the value of δ l at x i+1 in the asymptotic region by using the following methods:

•
The method of order eight produced in [22], which is denoted by QT(8).

•
The method of order ten developed in [22], which is denoted by QT(10).

•
The method of order eight developed in [22], which is denoted by QT(12).

•
The four-step method produced in [26], which is denoted by SI MI.

•
The four-step method developed in [29], which is denoted by SI M.

•
The method presented in [30], which is denoted by SI MN.

•
The method developed in [31], which is denoted by S2St.

•
The method produced in [20], which is denoted by S3St.

•
The method developed in [32], which is denoted by Y2D.

•
The method produced in [21], which is denoted by Y3D.

•
The four-step classical method described in Section 2.3, which is denoted by CLM.

•
The new four-step method presented in Section 2, as given in Equation (13), which is denoted by NewM.

Conclusions
As it is seen from the results in Figures 5-7, the newly developed explicit four-step method (NewM) is more efficient than the other methods in solving the resonance problem for the different energy values; E 1 = 53.588872, E 2 = 341.495874, and E 3 = 989.701916. It is found that among all the compared methods the new method (NewM) provides the highest accuracy with a lower cost of consumed CPU time.