Fourth derivative singularly P-stable method for the numerical solution of the Schrödinger equation

In this paper, we construct a method with eight steps that belongs to the family of Obrechkoff methods. Due to the explicit nature of the new method, not only does it not require another method as predictor, but it can also be considered as a suitable predictive technique to be used with implicit methods. Periodicity and error terms are studied when applied to solve the radial Schrödinger equation, considering different energy levels. We show its advantages in terms of accuracy, consistency, and convergence in comparison with other methods of the same order appearing in the literature.


Introduction
Our main goal is to discuss the one-dimensional radial Schrödinger equation of the form y (x) = l(l + 1) where V (x) and E are the potential and energy, respectively. The second additional condition is obtained from the physical meaning of the problem. The numerical solution of equations of type (1) is of special importance for applied scientists, because they appear in many areas of science when formulating natural problems such as in chemical physics and in quantum physics, among others. Since the analytical solution for problem (1) has a special and complex structure, its numerical approximation and the construction of accurate approximation methods are of great importance. In recent years, various methods have been designed by different authors to solve the Schrödinger equation numerically. We will mention a limited number of these methods, such as Runge-Kutta methods [1][2][3][4][5][6][7][8][9][10][11][12], hybrid methods [21][22][23], EF and TF methods [13][14][15][16][17][18][19][20].
Since the method we propose here belongs to the category of linear multistep multiderivative methods, we turn our attention to presenting some properties of these methods.
Consider the multiderivative symmetric multistep Obrechkoff method defined as In the case α j = α k-j and β ij = β i,k-j , for j = 0, 1, 2, . . . , k, method (2) will be symmetric. Also, its order will be q, provided the truncation error is TE = C q+2 h q+2 y (q+2) (η), x n-k+1 < η < x n+1 , in which C q+2 is a constant. To study the properties of stability of a method for solving second-order differential equations, Lambert and Watson [10] offered the scalar test equation If we apply (2) to (3), then we have the characteristic equation in which v = ωh and β ij ξ k-j , i = 1, 2, . . . , l.
Definition 1. 1 The interval (0, v 2 0 ) is called the interval of periodicity of method (2) provided for each v 2 in this interval (4) has complex roots so that at least two of the roots are on the unit circle and the rest are inside the unit circle. Also, if the interval of periodicity is (0, ∞), then method (2) has the property of P-stability.
In general, if we apply a 2k-step symmetric scheme to (3), its difference equation is obtained in the form in which A i (v), i = 0, 1, . . . , k, are polynomials depending on v. In addition, CE (the characteristic equation) corresponding to the difference equation (5) is Definition 1.2 For any scheme with CE (6), PL (phase-lag) is defined as the first term in the expansion of where θ (v) is a real function of v. If we can write t = O(v q+1 ) as v → ∞, then the order of PL is q. .
Achieving the P-stability in numerical schemes is of great importance because with Pstability property, we can say with confidence that the step-size constraint is eliminated. In other words, differential equations with extreme oscillations can be easily solved with P-stable methods, and the desired solution can be approximated accurately. Of course, the P-stability property is a general concept and is used for differential equations involving multiple frequencies. The concept of singular P-stability has been applied to problems with only one frequency. So when we exclusively talk about the interval and distinguish it from the periodicity interval, we are, in fact, talking about the same concepts, and there is no difference between the P-stability and the singular P-stability. Indeed, singular Pstability is the same as P-stability with one frequency. On the other hand, given that the Schrödinger equation discussed in this paper has only one frequency, we use the expression singular P-stability. We will show this property in the section devoted to the periodicity analysis.
The paper is organized as follows. Introduction and review of multiderivative methods are presented in Sect. 1. In Sect. 2, we propose the main method, and its coefficients will be produced with a new approach. Periodicity analysis and Schrödinger error coefficient analysis are provided in Sect. 3. Finally, the efficiency of the new scheme is demonstrated in Sect. 4 by implementing the new method on the Schrödinger equation with several energy levels.

Development and analysis 2.1 Development
Since the main purpose of this paper is to solve numerically (1) using symmetric multiderivative multistep methods, we consider an eight-step fourth-derivative method of the form where c 3 = 1 500 , c 2 = -1 500 , , and a j , b j , j = 0, 1, 2, 3, are real constant coefficients that should be computed. It is natural that all the properties and characteristics of any numerical method are deduced from its coefficients. However, the coefficients must be carefully generated in order to achieve superior properties. In this paper, several goals are pursued simultaneously so that we can produce a better and more accurate method. Not only do we want to control the interval of periodicity so that (0, +∞) is produced, but also we intend to have the Schrödinger error coefficient reduced. Applying the new method (7) to the scalar equation (3), CE of the method (7) will be as follows: where Now, to control the interval of periodicity, we will generate the coefficients so that CE of the scheme is To this end, the polynomials A 1 (v), A 2 (v), and A 3 (v) of the characteristic function must be equal to zero, which produces three equations. On the other hand, because the new method has eight free parameters, we need five more independent linear equations to obtain unique coefficients. Now, to control the Schrödinger error coefficient, suppose that PL and its derivatives up to four are zero. This leads to 8 equations with 8 unknowns. According to Theorem 1.3 with k = 4, PL is The local truncation error of (7), namely LTE New , is obtained using the usual Taylor series expansion and is given by LTE New = -326,687 5,670,000 h 10 ω 10 y n + 5ω 8 y (2) n + 10ω 6 y (4) n + 10ω 4 y (6) n + 5ω 2 y (8) n + y (10) Then the local truncation error of the classical form of (7) (the same structure with constant coefficients), namely LTE Class , will be LTE Class = -326,687 5,670,000 To study the efficiency of the new scheme for solving (1), we have to have its Schrödinger error. To this end, we can transfer (1) to Now, as per Ixaru and Rizea's paper [9], f (x) can be changed to and V c is the constant approximation of the potential and G = v 2 = V c -E (see [9]).

Theorem 2.1 The Schrödinger error of the new method increases as the second power of G.
Proof We know that, for any eighth algebraic-order linear multistep method, the general form of the LTEs is given by where, in the classical case, N i are constant numbers and in the frequency dependent cases, N i are functions of v and G, and j is the maximum power of G. Since G = V c -E, we can assume two cases according to the values of E: Obviously, the method with asymptotic form of LTE, which includes the minimum power of G, is the most accurate one. Now, (10) and (11) is -326,687 5,670,000 [ω 10 y n + 5ω 8 y (2) n + 10ω 6 y (4) n + 10ω 4 y (6) n + 5ω 2 y (8) n + y (10) n ] and -326,687 5,670,000 y (10) n , respectively. For the case G 0 or G 0, firstly, we should calculate higher derivatives of y. By simple calculation, we have and thus, for this method, the error increases as the second power of G.

Periodicity analysis
In order to confidently implement a numerical method on problems, we must have accurate information about how the method behaves. More specifically, we need to know the stability or instability, convergence or divergence of the method and even the maximum step length. For such studies, we must calculate the stability function of the method and then generate the stability region. This helps us to easily compare the proposed method with other methods. To generate the stability area of the scheme, we apply the main method to the problem and obtain its CE. Note that the frequency used in (12) is different from the frequency used in (3). If we assume s = φh and v = ωh, then we will produce y-axes with v and x-axis with s in the stability region. Therefore, the proposed method can be used for problems that have two frequencies. Applying (7) to (12), we will have in which v = ωh, s = φh, and A i (s, v) are functions of s and v. To save space, A i (s, v) are given in Appendix B.
As explained at the beginning of this section, the stability function of the new method has two parameters, s and v, each of which is derived from a frequency. We are looking for P-stability of the method to solve more problems with the proposed method. The colored parts in Fig. 1, which is the figure for the stability area of the method, show the stability parts, and the white parts show the instability region of the method. The method will be Pstable when the entire sv plane is colored. But this is when we talk about two-frequency issues. Since the equation discussed in this paper has one frequency, the concept of Pstability changes to the one of singular P-stability. It is enough that the bisector of the first quarter belongs to the colored part. Hence, the method is singular P-stable. Now we prove the property of singular P-stability in the next theorem algebraically. (7) is singularly P-stable.
If we assume λ 2k = exp(I 2kπ m ) exp(Iv) (see [22] So, from (15) and (16) Clearly, for these roots we have |λ 0 | = |λ 1 | = 1 and |λ i | ≤ 1, i = 2, 3, . . . , 7. Accordingly, when s = v, the periodicity interval of the new scheme is equal to (0, ∞), and the method is P-stable. Therefore, it is singularly P-stable. (14) is the same as the one obtained in [18,Theorem 5]. To explain this, since to prove the P-stability property of a linear multistep method, we have to show that the characteristic roots have modulus less (or equal) than one. In other words, they must lie in or on the unit circle. Although the structure of the proposed method in [18] is different from the method presented in this paper, we know the characteristic equation of an eight-step linear multistep method is equal to (8), and also its two variable characteristic equation is (13). Singularly P-stability means P-stability when s = v. To generate a system of equations for calculating the coefficients of the method, we have assumed that A 1 (v), A 2 (v), and A 3 (v) are equal to zero and the remaining equations are obtained from vanishing phase-lag and some of its derivatives. Now, by sub-

Remark 2.3 The characteristic equation
. Hence, the characteristic polynomial will be λ 8 + B(v)λ 4 + 1 = 0, and then all characteristic roots are equal or less than one.

Numerical results
This section is devoted to implementing the method built on the Schrödinger problem. In order to better judge the quality of the method, we have implemented it on two energy levels E = 341.495874 and E = 989.701916. We have compared the produced results with those of the other methods which are of the same order as the new one, and we showed the superiority of the new method. Since we need a value ω in numerical implementation, we need to specify this value. Different methods have been proposed in different papers to consider ω. We may mention ω = |V (x) -E| in which v(x) is the potential; in this article, we will use the Woods-Saxon potential function. The definition of the function V is (see Fig. 2) Here, in order to implement and produce numerical results of the new method (see [8]), we let This value is considered for the interval [0, 15].

Schrödinger equation-resonance problem
Consider the numerical solution of the radial time-independent Schrödinger equation (1) in the well-known case of the Woods-Saxon potential (17). To numerically solve this problem, we should approximate the true or infinite interval of integration by applying a finite interval. Since we need to illustrate our problem numerically, we take the domain of integration as x ∈ [0, 15]. We consider (1) in a relatively large domain of energies, i.e., E ∈ [1, 1000]. When it comes to positive energies, E = k 2 , the potential fades faster than the term l(l+1) x 2 , and the Schrödinger equation effectively reduces to for x greater than some value X. Equation (18) has two linearly independent solutions kxj l (kx) and kxn l (kx), where j l and n l are the spherical Bessel and Neumann functions, respectively. When x → ∞, the solution takes the asymptotic form where δ l is called scattering phase shift that may be calculated from the formula where x 1 and x 2 are distinct points in the asymptotic region (we choose x 1 as the righthand end point of the interval of integration and x 2 = x 1h) with S(x) = kxj l (kx) and C(x) = -kxn l (kx). The problem is dealt with as an initial value problem; thus, we have to have y 0 . We obtain y 0 from the initial condition. With these starting values, we evaluate at x 1 of the asymptotic region the phase shift δ l . For positive energies, we have resonance problem which is comprised either of finding the phase-shift δ l or finding those E, for E ∈ [1, 1000], at which δ l = π 2 . We solve the problem when the positive energies lie under potential barrier. The boundary conditions for this problem are The following methods have been used to compare the new method: • The twelve-step thirteenth algebraic-order method developed by Quinlan and Tremaine [15] which is indicated as I; • The eight-step method with PL and its first derivative equal to zero obtained in [5] which is indicated as II; • The ten-step method with PL and its first, second, and third derivatives equal to zero obtained in [7] which is indicated as III; • The exponentially-fitted four-step method developed by Raptis [16] which is indicated as IV; • The eight-step method with PL and its first and second derivative equal to zero obtained in [6] which is indicated as V; • The trigonometrically-fitted six-step method developed by Wang [22] which is indicated as VI; • The new explicit eight-step singularly P-stable multiderivative method developed in this paper which is indicated as New.
The computed eigenenergies are compared with the exact ones. In Fig. 3, we present the digits of accuracy given by log 10 (Err), where versus the CPU times for several methods used for calculating the eigenenergy E 2 = 341.495874. In Fig. 4, we present the maximum absolute error log 10 (Err) of the eigenenergy E 3 = 989.701916. In addition, it is pointed out in Fig. 5

Conclusion
By using higher-order derivatives in classical methods and combining them with the system of PL and its derivatives, we were able to create a new high-efficiency method (with lower CPU time) that can approximate different types of energy levels of the Schrödinger equation with high accuracy. In fact, by keeping the computation time low, we were able to produce better quality results, which is very important in numerical analysis.