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.
The numerical solution of second-order initial-value problems with periodical and/or oscillatory solutions as:
has attracted the attention of many authors in the last decades [1,2, 3,4,5, 6,7,8, 9,10,11, 12,13,14, 15,16,17, 18,19,20, 21,22,23, 24,25]. 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.  and Steifel and Bettis , 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 and a second condition corresponds to large values of , which can be adopted due to some physical considerations. Here, is a real number denoting the energy, is an integer representing the angular momentum, and is the potential function which satisfies as .
In the literature, a lot of research works have been done in developing numerical methods of different types for solving equation (1.2). Explicit multistep phase-fitted methods are developed by Lambert and Watson in , Anastassi and Simos in , Simos and Williams in , Alolyan and Simos , Simos in [35,36,37], and Obaidat and Mesloub in . An implicit multistep phase-fitted methods is constructed in . Predictor corrector methods are designed by Panopoulos et al. in , Stasinos and Simos in , and Simos in . Exponentially fitted multistep methods are provided by Simos in , Zhang et al. in , and G. Avdelas et al. in . Exponentially and trigonometrically fitted multistep methods are developed by Konguetsof and Simos designed in kong, and Simos constructed in . Runge-Kutta methods are presented by Dormand and Prince in , Yang et al. in , and Yang et al. in . Runge-Kutta-Nyström methods are constructed by Dormand et al. in . Methods of Obrechkof-type are derived by Krishnaiah in , Achar in , Van Daele and Berghe in , and Jain in .
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 , 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 phase-lag 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.
A numerical solution of a mathematical model of type (1.1) can be obtained using multi-step methods as:
by dividing the interval of definition of the initial value problem under consideration into subintervals each of length , using a finite set of equally spaced points , where , for . If is an even integer and and for , then such method is called symmetric multi-step numerical method.
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 ; so that equation (2.1) takes the form of a -step symmetric method as:
Now, applying the symmetric -step method (2.3) to a scalar test equation as:
gives the difference equation:
where , , , and is the step-size.
Equation (2.5) is accompanied with a characteristic equation as:
Lambert and Watson have defined the interval of periodicity of a symmetric linear multi-step method as follows:
for any for some real valued function .
The phase-lag of a symmetric multi-step method is defined as:
We say that the phase-lag of a given multi-step method has an order if is as . Moreover, the phase-lag of a given symmetric method with -steps can be found using the formula in the following result:
3 The method development
In this paper, we consider an implicit four-step symmetric method on the form:
where , .
Using in Theorem 2.4 implies the following proposition:
The implicit linear multi-step method defined by (3.1) has a phase-lag order and a phase-lag constant given as:
and , .
In the method (3.1), we choose the parameters , , and so that this method is consistent and has the highest algebraic order, while the coefficients , , and are considered as free parameters. Thus, we take , , and . Then, using these values in equation (3.1) yields the following numerical method:
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 , , and , we obtain:
As it appears in (3.6), the expressions of , , and involve the frequency . Moreover, the numerators and denominators in these expressions approach zero as . 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
4.1 Comparative error analysis
The local truncation error LTer associated with the new method, denoted as NewM, is obtained by using the free parameters , , and given by (3.7) into equation (3.3), then expanding both sides in Taylor series, which implies:
where the function can be set as ():
with , where is a potential function, is an estimation of the potential , , and 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 explicit four-step method in :
The implicit four-step method in :
The new implicit method (3.3):
Then, using the expressions of , 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 explicit four-step method in :
The method derived in :
The new method given by equation (3.3):
We see from these equations that the expressions of these truncation errors involve powers of the expression . Therefore, we should consider the following cases:
, i.e., when the values of the energy and the potential 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 . Hence, it follows that all these methods are of comparable accuracy, as these terms are identical in all these methods.
or , i.e., when the values of the energy and the potential are significantly different from each other, that is, 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 , whereas it increases as the first power of for the other methods. But, for the new method, 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 , it seems to be that, among these numerical techniques, the new derived one is the most efficient.
4.2 Stability analysis
To study the stability analysis of the new numerical method, we substitute the values of the parameters , , and given in equation (3.6) in the method (3.3), then apply the resulting scheme to the test equation , which produces a difference equation associated with the following characteristic equation:
where , , , and . Here, the frequency of the derived method (3.3) is different than the frequency of the test equation . The stability region of the constructed method is presented in 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).
5 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:
and another boundary condition will be adapted to certain physical considerations for some large value of . Here, in equation (5.1) is a constant that stands for the energy, is the angular momentum, and represents the potential which satisfies 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 , and choose as:
where is the energy, and use the Woods-Saxon potential:
where and . The behavior of is shown in Figure 5.
In Woods-Saxon potential, the frequency is not a function in , instead it is approximated by estimating the potential at the certain potential critical points .
For the numerical tests, we take as follows (see )
which has the two linearly independent solutions, and , where and are the Bessel and Neumann spherical functions, respectively. Hence, as , equation (5.1) has a solution in asymptotic form as:
where represents the phase shift that can be determined by the formula:
where and are two successive points in the asymptotic region, with and .
Now, given the energy , we will employ the new method (3.3) to solve the resonance problem; that is, to estimate the phase shift which has exact value that equals to . But, first we need to determine the start values for . Thus, in view of the condition (5.2) we get , and the remaining other values can be computed using Runge-Kutta-Nyström methods (see  and ). Starting with these values, we will estimate the value of at the point using the following numerical schemes:
The explicit 8-step eighth order linear symmetric method with constant coefficients given in ; denote it by .
The explicit 10-step tenth order linear symmetric method with constant coefficients developed in ; denote it as .
The explicit 12-step twelfth order linear symmetric method with constant coefficients developed in ; denote it by .
The explicit fourth order 4-step linear symmetric method with vanished phase-lag and its first derivative given in ; 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 ; denote it by SIM.
The explicit fourth order 4-step symmetric method with vanished phase-lag and its first and second derivatives given in ; 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 ; denote it by .
The eighth order three stage two-step linear symmetric method with vanished phase-lag and its first, second, third, and fourth derivatives developed in ; denote it as .
The three-stage fifth order two-derivative Runge-Kutta method, developed in ; denote it by .
The fourth order three derivative Runge-Kutta method with vanished phase-lag and its first derivative produced in ; denote it as .
The implicit sixth order 4-step linear symmetric method with vanished phase-lag and its first derivative in ; denote it by SHIM.
The explicit fourth order 4-step linear symmetric method with vanished phase-lag and its first derivative developed in ; denote it by EFSM.
The implicit fourth order 4-step symmetric classical method (with constant coefficients) given in Section 2.3; denote it as CLM.
The new developed implicit sixth-order 4-step linear symmetric method with vanished phase-lag and its first derivative given in equation (3.3) developed in Section 2; denote it by NewM.
The efficiency of the new developed method is tested against these methods via computing the CPU time (in seconds) required to get different accuracy digits in estimating the value of using the two energy values, and . Then, the absolute error :
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.
The authors extend their appreciation to the Deanship of Scientific Research at King Saud University for funding this work through research group no. (RG-1441-338).
Conflict of interest: Authors state no conflict of interest.
 A. D. Raptis and A. C. Allison, Exponential-fitting methods for the numerical solution of the Schrödinger equation, Comput. Phys. Commun. 14 (1978), 1–5. 10.1016/0010-4655(78)90047-4Search in Google Scholar
 G. D. Quinlan, Resonances and instabilities in symmetric multi-step methods, arXiv:astro-ph/9901136, (1999). Search in Google Scholar
 S. D. Capper, J. R. Cash, and D. R. Moore, Lobatto-Obrechkoff formulae for 2nd order two-point boundary value problems, J. Numer. Anal. Ind. Appl. Math. 1 (2006), no. 1, 13–25. Search in Google Scholar
 M. K. Jain, R. K. Jain, and U. A. Krishnaiah, Obrechkoff methods for periodic initial value problems of second order differential equations, J. Math. Phys. Sci. 15 (1981), no. 3, 239–250. Search in Google Scholar
 M. M. Chawla and P. S. Rao, A Noumerov-type method with minimal phase-lag for the integration of second order periodic initial value problems. II. Explicit method, J. Comput. Appl. Math. 15 (1986), no. 3, 329–337. 10.1016/0377-0427(86)90224-4Search in Google Scholar
 G. Dahlquist, On accuracy and unconditional stability of linear multistep methods for second order differential equations, BIT 18 (1978), no. 2, 133–136. 10.1007/BF01931689Search in Google Scholar
 J. M. Franco, An explicit hybrid method of Numerov type for second-order periodic initial-value problems, J. Comput. Appl. Math. 59 (1995), no. 1, 79–90. 10.1016/0377-0427(94)00011-OSearch in Google Scholar
 U. A. Krishnaiah, P-stable Obrechkoff methods with minimal phase-lag for periodic initial value problems, Math. Comp. 49 (1987), no. 180, 553–559. 10.1090/S0025-5718-1987-0906188-XSearch in Google Scholar
 G. Saldanha and S. D. Achar, Symmetric multistep methods with zero phase-lag for periodic initial value problems of second order differential equations, Appl. Math. Comput. 175 (2006), no. 1, 401–412. 10.1016/j.amc.2005.07.054Search in Google Scholar
 S. D. Achar, Symmetric multistep Obrechkoff methods with zero phase-lag for periodic initial value problems of second order differential equations, Appl. Math. Comput. 218 (2011), no. 5, 2237–2248. 10.1016/j.amc.2011.07.040Search in Google Scholar
 M. Van Daele and G. V. Berghe, P-stable exponentially-fitted Obrechkoff methods of arbitrary order for second-order differential equations, Numer. Algorithms 46 (2007), no. 4, 333–350. 10.1007/s11075-007-9142-ySearch in Google Scholar
 J. Vigo-Aguiar and H. Ramos, On the choice of the frequency in trigonometrically-fitted methods for periodic problems, J. Comput. Appl. Math. 277 (2015), 94–105. 10.1016/j.cam.2014.09.008Search in Google Scholar
 J. Vigo-Aguiar and H. Ramos, A strategy for selecting the frequency in trigonometrically-fitted methods based on the minimization of the local truncation errors and the total energy error, J. Math. Chem. 52 (2014), 1050–1058. 10.1007/s10910-013-0282-0Search in Google Scholar
 H. Ramos and J. Vigo-Aguiar, A trigonometrically-fitted method with two frequencies, one for the solution and another one for the derivative, Comput. Phys. Commun. 185 (2014), no. 4, 1230–1236. 10.1016/j.cpc.2013.12.021Search in Google Scholar
 H. Ramos and J. Vigo-Aguiar, Variable-stepsize Chebyshev-type methods for the integration of second-order I.V.P.’s, J. Comput. Appl. Math. 204 (2007), 102–113. 10.1016/j.cam.2006.04.032Search in Google Scholar
 H. Ramos and J. Vigo-Aguiar, A variable-step Numerov method for the numerical solution of the Schrödinger equation, J. Math. Chem. 37 (2005), 255–262. 10.1007/s10910-004-1467-3Search in Google Scholar
 Z. A. Anastassi and T. E. Simos, A parametric symmetric linear four-step method for the efficient integration of the Schrödinger equation and related oscillatory problems, J. Comput. Appl. Math. 236 (2012), no. 16, 3880–3889. 10.1016/j.cam.2012.03.016Search in Google Scholar
 I. Alolyan, Z. A. Anastassi, and T. E. Simos, A new family of symmetric linear four-step methods for the efficient integration of the Schrödinger equation and related oscillatory problems, Appl. Math. Comput. 218 (2012), no. 9, 5370–5382. 10.1016/j.amc.2011.11.020Search in Google Scholar
 Z. A. Anastassi, A new symmetric linear eight-step method with fifth trigonometric order for the efficient integration of the Schrödinger equation, Appl. Math. Lett. 24 (2011), no. 8, 1468–1472. 10.1016/j.aml.2011.03.035Search in Google Scholar
 T. E. Simos, Exponentially-fitted multiderivative methods for the numerical solution of the Schrödinger equation, J. Math. Chem. 36 (2004), 13–27. 10.1023/B:JOMC.0000034930.81720.47Search in Google Scholar
 P. I. Stasinos and T. E. Simos, Symmetric embedded predictor-predictor-corrector (EPPCM) methods with vanished phase-lag and its derivatives for second order problems, AIP Conf. Proc. 1906 (2017), no. 1, 200023, https://doi.org/10.1063/1.5012499. 10.1063/1.5012499Search in Google Scholar
 T. E. Simos, Predictor-corrector phase-fitted methods for Y″=F(X,Y) and an application to the Schrödinger equation, Int. J. Quantum Chem. 53 (1995), no. 5, 473–483. 10.1002/qua.560530504Search in Google Scholar
 T. E. Simos, Exponentially and trigonometrically fitted methods for the solution of the Schrödinger equation, Acta Appl. Math. 110 (2010), no. 3, 1331–1352. 10.1007/s10440-009-9513-6Search in Google Scholar
 A. Konguetsof and T. E. Simos, An exponentially-fitted and trigonometrically-fitted method for the numerical solution of periodic initial-value problems, Comput. Math. Appl. 45 (2003), no. 1–3, 547–554. 10.1016/S0898-1221(03)80036-6Search in Google Scholar
 Y. Zhang, X. You, and Y. Fang, Exponentially fitted multi-derivative linear methods for the resonant state of the Schrödinger equation, J. Math. Chem. 55 (2017), 223–237. 10.1007/s10910-016-0683-ySearch in Google Scholar
 T. E. Simos and P. S. Williams, A finite-difference method for the numerical solution of the Schrdinger equation, J. Comput. Appl. Math. 79 (1997), 189–205. 10.1016/S0377-0427(96)00156-2Search in Google Scholar
 I. Alolyan and T. A. Simos, Family of explicit linear six-step methods with vanished phase-lag and its first derivative, J. Math. Chem. 52 (2014), 2087–2118. 10.1007/s10910-014-0364-7Search in Google Scholar
 T. E. Simos, An explicit four-step method with vanished phase-lag and its first and second derivatives, J. Math. Chem. 52 (2014), no. 1, 833–855. 10.1007/s10910-013-0296-7Search in Google Scholar
 T. E. Simos, A new explicit four-step method with vanished phase-lag and its first and second derivatives, J. Math. Chem. 53 (2015), no. 1, 402–429. 10.1007/s10910-014-0431-0Search in Google Scholar
 A. Shokri and M. Tahmourasi, A new efficient implicit four-step method with vanished phase-lag and its first derivative for the numerical solution of the radial Schrödinger equation, J. Mod. Methods Numer. Math. 8 (2017), no. 1–2, 77–89. 10.20454/jmmnm.2017.1223Search in Google Scholar
 G. A. Panopoulos, Z. A. Anastassi, and T. Simos, A symmetric eight-step predictor-corrector method for the numerical solution of the radial Schrödinger equation and related IVPs with oscillating solutions, Comput. Phys. Commun. 182 (2011), no. 8, 1626–1637. 10.1016/j.cpc.2011.04.011Search in Google Scholar
 G. Avdelas, E. Kefalidis, and T. E. Simos, New P-stable eighth algebraic order exponentially-fitted methods for the numerical integration of the Schrödinger equation, J. Math. Chem. 31 (2002), no. 4, 371–404. 10.1023/A:1021020705327Search in Google Scholar
 Y. Yang, Y. Fang, K. Wang, and X. You, THDRK methods with vanished phase-lag and its first derivative for the Schrödinger equation, J. Math. Chem. 57 (2019), 1496–1507. 10.1007/s10910-019-01002-7Search in Google Scholar
 L. Gr. Ixaru and M. Rizea, A Numerov-like scheme for the numerical solution of the Schrödinger equation in the deep continuum spectrum of energies, Comput. Phys. Commun. 19 (1980), no. 1, 23–27. 10.1016/0010-4655(80)90062-4Search in Google Scholar
 Z. Zhou and T. E. Simos, A new two stage symmetric two-step method with vanished phase-lag and its first, second, third and fourth derivatives for the numerical solution of the radial Schrödinger equation, J. Math. Chem. 54 (2016), 442–465. 10.1007/s10910-015-0571-xSearch in Google Scholar
 Y. Lan and T. E. Simos, An efficient and economical high order method for the numerical approximation of the Schrödinger equation, J. Math. Chem. 55 (2017), 1755–1778. 10.1007/s10910-017-0757-5Search in Google Scholar
© 2021 Saleem Obaidat and Rizwan Butt, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.