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

Saleem Obaidat; Rizwan Butt

doi:10.1515/math-2021-0009

Open Access Published by De Gruyter Open Access May 3, 2021

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

Saleem Obaidat and Rizwan Butt

From the journal Open Mathematics

https://doi.org/10.1515/math-2021-0009

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.

Keywords: phase-lag; stability analysis; implicit symmetric method; Schrödinger equation; truncation error

MSC 2010: 65D15; 65D20; 34A50

1 Introduction

The numerical solution of second-order initial-value problems with periodical and/or oscillatory solutions as:

(1.1) q ″ ( x ) = f ( x , q ( x ) ) , q ( x 0 ) = q 0 and q ′ ( x 0 ) = q 0 ′ ,

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. [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:

(1.2) y ″ ( φ ) = l ( l + 1 ) φ 2 + V ( φ ) − τ 2 y ( φ ) .

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 φ → ∞ .

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 [32], Anastassi and Simos in [23], Simos and Williams in [33], Alolyan and Simos [34], Simos in [35,36,37], and Obaidat and Mesloub in [38]. An implicit multistep phase-fitted methods is constructed in [39]. Predictor corrector methods are designed by Panopoulos et al. in [40], Stasinos and Simos in [27], and Simos in [28]. Exponentially fitted multistep methods are provided by Simos in [26], Zhang et al. in [31], and G. Avdelas et al. in [41]. Exponentially and trigonometrically fitted multistep methods are developed by Konguetsof and Simos designed in kong, and Simos constructed in [29]. Runge-Kutta methods are presented by Dormand and Prince in [42], Yang et al. in [43], and Yang et al. in [44]. Runge-Kutta-Nyström methods are constructed by Dormand et al. in [45]. Methods of Obrechkof-type are derived by Krishnaiah in [14], Achar in [16], Van Daele and Berghe in [17], and Jain in [9].

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 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.

2 Preliminaries

A numerical solution of a mathematical model of type (1.1) can be obtained using multi-step methods as:

(2.1) ∑ i = 0 m a i q n + i = h 2 ∑ i = 0 m b i f ( x n + i , q n + i ) ,

by dividing the interval of definition of the initial value problem under consideration into m subintervals each of length h , using a finite set of equally spaced points { x i } i = 0 k , where h = ∣ x i + 1 − x i ∣ , for i = 0 , 1 , … , m − 1 . If m is an even integer and a j = a m − j and b j = b m − j for j = 0 , 1 , … , m 2 , 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:

(2.2) ℒ ( x ) = ∑ j = 0 m a j ψ ( x + j h ) − h 2 ∑ j = 0 m b j ψ ″ ( x + j h ) ,

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 m = 1 ; so that equation (2.1) takes the form of a 2 m -step symmetric method as:

(2.3) a 0 q n + ∑ k = 1 m a k ( q n + k + q n − k ) = h 2 b 0 f n + ∑ k = 1 m b k ( f n + k + f n − k ) ,

where f n ± k = f ( x n ± k , q n ± k ) .

Definition 2.1

[23] The multi-step method (2.1) is said to have an algebraic order r whenever the values of the corresponding operator ℒ vanish at any function p ( x ) = ∑ j = 0 r + 1 c j x j , for any real numbers c j , j = 0 , 1 , … , r + 1 .

Now, applying the symmetric 2 m -step method (2.3) to a scalar test equation as:

(2.4) q ″ = − ϕ 2 q ,

gives the difference equation:

(2.5) D m ( v ) q n + m + ⋯ + D 1 ( v ) q n + 1 + D 0 ( v ) q n + D 1 ( v ) q n − 1 + ⋯ + D m ( v ) q n − m = 0 ,

where D l ( v ) = a l + b l v 2 , l = 0 , … , m , v = ϕ h , and h is the step-size.

Equation (2.5) is accompanied with a characteristic equation as:

(2.6) D m ( v ) λ m + ⋯ + D 1 ( v ) λ + D 0 ( v ) + D 1 ( v ) λ − 1 + ⋯ + D m ( v ) λ − m = 0 .

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 2 k -step method has an interval of periodicity ( 0 , v 0 2 ) whenever the roots w ℓ , ℓ = 1 , 2 , … , 2 k of the characteristic equation (2.6) satisfy the conditions:

(2.7) w 1 = e I σ ( v ) , w 2 = e − I σ ( v ) , and ∣ w ℓ ∣ ≤ 1 , for ℓ = 3 , … , 2 k ,

for any v ∈ ( 0 , v 0 2 ) , for some real valued function σ ( v ) .

The phase-lag of a symmetric multi-step method is defined as:

Definition 2.3

[1,33] Each symmetric multi-step method with characteristic equation as (2.6) has a phase-lag that equals to the first non-vanishing term in the series expansion of:

(2.8) ζ = v − σ ( v ) .

We say that the phase-lag of a given multi-step method has an order κ if ζ is O ( v κ + 1 ) as v → ∞ . Moreover, the phase-lag of a given symmetric method with 2 k -steps can be found using the formula in the following result:

Theorem 2.4

[33] Any linear symmetric 2 k -step method with characteristic equation as (2.6) has a phase-lag order κ and a phase-lag constant ρ satisfying the relation:

(2.9) − ρ v κ + 2 + O ( v κ + 4 ) = D 0 ( v ) + ∑ j = 1 k 2 D j ( v ) cos ( j v ) ∑ j = 1 k 2 j 2 D j ( v ) .

3 The method development

In this paper, we consider an implicit four-step symmetric method on the form:

(3.1) y n + 2 + μ 1 ( y n + 1 + y n − 1 ) + μ 0 y n + y n − 2 = h 2 [ ν 2 ( f n + 2 + f n − 2 ) + ν 1 ( f n + 1 + f n − 1 ) + ν 0 f n ] ,

where f n + i = y ″ ( x n + i , y n + i ) , i = 0 , ± 1 , ± 2 .

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 phase-lag constant r given as:

(3.2) − r v s + 2 + O ( v s + 4 ) = ℋ 0 ( v ) ℋ 1 ( v ) ,

where

H 0 ( v ) = 2 D 2 ( v ) cos ( 2 v ) + 2 D 1 ( v ) cos ( v ) + D 0 ( v ) , H 1 ( v ) = 8 D 2 ( v ) + 2 D 1 ( v ) ,

and D i ( v ) = μ i + ν i v 2 , i = 0 , 1 , 2 .

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, we take μ 2 = 1 , μ 1 = − 15 16 , and μ 0 = − 1 8 . Then, using these values in equation (3.1) yields the following numerical method:

(3.3) y n + 2 − 15 16 ( y n + 1 + y n − 1 ) − 1 8 y 0 + y n − 2 = h 2 [ ν 2 ( f n + 2 + f n − 2 ) + ν 1 ( f n + 1 + f n − 1 ) + ν 0 f 0 ] .

Hence, in view of (3.2), the method (3.3) has a phase-lag given by:

(3.4) P L ( v ) = ϕ 1 ( v ) ϕ 2 ( v ) ,

where

ϕ 1 ( v ) ≔ − 1 8 + v 2 ν 0 + 2 cos v − 15 16 + v 2 ν 1 + 4 cos v ( 1 + v 2 ν 2 ) ,

and

ϕ 2 ( v ) ≔ 2 − 15 16 + v 2 ν 1 + 8 ( 1 + v 2 ν 2 ) .

The following conditions are imposed on the new method:

(3.5) P L ( v ) = 0 , P L ′ ( v ) = 0 , 16 ν 0 + 32 ν 1 + 32 ν 2 = 49 ,

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:

(3.6) ν 0 = ψ 1 ( v ) ψ 0 ( v ) , ν 1 = − ψ 2 ( v ) ψ 0 ( v ) , ν 2 = − ψ 3 ( v ) 2 ψ 0 ( v ) ,

where

ψ 1 ( v ) = 4 cos ( v ) + 60 cos ( v ) 2 − 4 cos ( 2 v ) − 124 cos ( v ) cos ( 2 v ) + 64 cos ( 2 v ) 2 − 2 v sin ( v ) + 2 v cos ( 2 v ) sin ( v ) + 49 v 3 cos ( 2 v ) sin ( v ) + 4 v sin ( 2 v ) − 4 v cos ( v ) sin ( 2 v ) − 98 v 3 cos ( v ) sin ( 2 v ) , ψ 2 ( v ) = 2 + 30 cos ( v ) − 34 cos ( 2 v ) − 30 cos ( v ) cos ( 2 v ) + 32 cos ( 2 v ) 2 + 15 v sin ( v ) − 15 v cos ( 2 v ) sin ( v ) − 30 v sin ( 2 v ) − 49 v 3 sin ( 2 v ) + 30 v cos ( v ) sin ( 2 v ) , ψ 3 ( v ) = − 4 − 56 cos ( v ) + 60 cos ( v ) 2 + 64 cos ( 2 v ) − 64 cos ( v ) cos ( 2 v ) − 32 v sin ( v ) + 49 v 3 sin ( v ) + 32 v cos ( 2 v ) sin ( v ) + 64 v sin ( 2 v ) − 64 v cos ( v ) sin ( 2 v ) , ψ 0 ( v ) = 16 v 3 ( − sin ( v ) + cos ( 2 v ) sin ( v ) + 2 sin ( 2 v ) − 2 cos ( v ) sin ( 2 v ) ) .

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:

(3.7) ν 0 = 1873 1920 + 359 v 2 11520 − 1059 v 4 358400 − 8209 v 6 212889600 − 2178023 v 8 498161664000 − 391187 v 10 996323328000 − ⋯ , ν 1 = 467 480 − 359 v 2 17280 + 4007 v 4 3628800 − 8747 v 6 958003200 − 77381 v 8 373621248000 − 107411 v 10 4483454976000 − ⋯ , ν 2 = 271 3840 + 359 v 2 69120 + 21667 v 4 58060800 + 108869 v 6 3832012800 + 7153117 v 8 2988969984000 + 3950327 v 10 17933819904000 + ⋯ .

Figures 1, 2, and 3 represent the behaviors of the parameters ν 0 , ν 1 , and ν 2 , respectively.

$Figure 1 Behavior of ν 0 {\nu }_{0} .$

Figure 1

Behavior of ν 0 .

$Figure 2 Behavior of ν 1 {\nu }_{1} .$

Figure 2

Behavior of ν 1 .

$Figure 3 Behavior of ν 2 {\nu }_{2} .$

Figure 3

Behavior of ν 2 .

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 ν 0 , ν 1 , and ν 2 given by (3.7) into equation (3.3), then expanding both sides in Taylor series, which implies:

(4.1) LTer NewM = − 359 h 8 138240 ( w 6 y ( 4 ) + 2 w 4 y ( 6 ) + w 2 y ( 8 ) ) + O ( h 10 ) .

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:

(4.2) y ″ ( x ) = p ( x ) y ( x ) ,

where the function p can be set as ([46]):

(4.3) p ( x ) = η ( x ) + G ,

with η ( x ) = Q ( x ) − Q c , where Q ( x ) is a potential function, Q c is an estimation of the potential Q ( x ) , G = Q c − E , 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 (3.3) by setting v = 0 in equation (3.7):
LTer CLM = − 359 h 8 138240 y n ( 8 ) + O ( h 10 ) .
The explicit four-step method in [34]:
LTer OLS2 = 787 h 8 12096 ( 4 w 6 y n ( 6 ) + 3 w 8 y n ( 4 ) + y n ( 8 ) ) + O ( h 10 ) .
The implicit four-step method in [39]:
LTer SHK2 = − 19 h 8 6048 ( w 4 y ( 4 ) + 2 w 2 y ( 6 ) + y ( 8 ) ) + O ( h 10 ) .
The new implicit method (3.3):
LTer NewM = − 359 h 8 138240 ( w 4 y ( 4 ) + 2 w 2 y ( 6 ) + y ( 8 ) ) + O ( h 10 ) .

To obtain the asymptotic expressions for the local truncation errors for these methods, first the values of the derivatives in their expressions when applied to the equation q n ″ = ( η ( x ) + G ) q ( x ) are computed, where q n is an estimation of q ( x n ) , which implies:

q n ( 2 ) = ( η ( x ) + G ) q ( x ) , q n ( 4 ) = η ″ ( x ) q ( x ) + 2 η ′ ( x ) q ′ ( x ) + ( η ( x ) + G ) 2 q ( x ) , q n ( 6 ) = η ( 4 ) ( x ) q ( x ) + 4 η ( 3 ) ( x ) q ′ ( x ) + 7 ( η ( x ) + G ) q ( x ) η ″ ( x ) + 4 ( η ′ ( x ) ) 2 q ( x ) + 6 ( η ( x ) + G ) q ′ ( x ) η ′ ( x ) + ( η ( x ) + G ) 3 q ( x ) , q n ( 8 ) = η ( 6 ) ( x ) q ( x ) + 6 η ( 5 ) ( x ) q ′ ( x ) + 16 ( η ( x ) + G ) q ( x ) η ( 4 ) ( x ) + 26 η ′ ( x ) η ( 3 ) ( x ) q ( x ) + 24 ( η ( x ) + G ) q ′ ( x ) ( η ′ ( x ) ) 3 + 15 ( η ″ ( x ) ) 2 q ( x ) + 48 η ′ ( x ) q ′ ( x ) η ″ ( x ) + 22 ( η ( x ) + G ) 2 q ( x ) η ″ ( x ) + 28 ( η ( x ) + G ) 2 q ( x ) ( η ′ ( x ) ) 2 + 12 ( η ( x ) + G ) 2 q ′ ( x ) η ′ ( x ) + ( η ( x ) + G ) 4 q ( x ) .

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):
LTEr CLM = h 8 − 359 138240 q ( x ) G 4 + ⋯ + O ( h 10 ) .
The explicit four-step method in [34]:
LTEr OSII = h 8 8657 6048 q ( x ) g ″ ( x ) + 787 1008 g ′ ( x ) q ′ ( x ) + 787 2016 y ( x ) g 2 ( x ) G 2 + ⋯ + O ( h 10 ) .
The method derived in [39]:
LTEr SHK2 = − 19 h 8 6048 [ ( 28 ( g ′ ( x ) ) 2 q ( x ) + 2 g ′ ( x ) q ′ ( x ) + 9 g ″ ( x ) q ( x ) ) G 2 + ⋯ ] + O ( h 10 ) .
The new method given by equation (3.3):
LTEr NewM = − 359 h 8 138240 [ ( 28 q ( x ) ( g ′ ( x ) ) 2 + 2 g ′ ( x ) q ′ ( x ) + 9 q ( x ) g ″ ( x ) ) G 2 + ⋯ ] + O ( h 10 ) .

We see from these equations that the expressions of these truncation errors involve powers of the expression G = V c − E . Therefore, we should consider the following cases:

Case 1

G ≈ 0 , i.e., 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.

Case 2

G ≪ 0 or G ≫ 0 , i.e., 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.

4.2 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 ″ = − θ 2 y , which produces a difference equation associated with the following characteristic equation:

D 2 λ 4 + D 1 λ 3 + D 0 λ 2 + D 1 λ + D 2 = 0 ,

where D 2 = 1 + ν 2 s 2 , D 1 = − 15 16 + ν 1 s 2 , D 0 = − 1 8 + ν 0 s 2 , and s = θ h . 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).

$Figure 4 New method stability region.$

Figure 4

New method stability region.

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:

(5.1) p ″ ( φ ) = l ( l + 1 ) φ 2 + V ( φ ) − γ 2 p ( φ ) ,

accompanied with two conditions:

(5.2) p ( 0 ) = 0 ,

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:

(5.3) ϕ = ∣ V ( φ ) − E ∣ ,

where E ≔ γ 2 is the energy, and use the Woods-Saxon potential:

(5.4) V ( φ ) = 50 u γ ( 1 + u ) 2 − 50 1 + u ,

where u = e φ − 7 γ and a = 0.6 . The behavior of V ( φ ) is shown in Figure 5.

$Figure 5 Behavior of Woods-Saxon potential.$

Figure 5

Behavior of Woods-Saxon potential.

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])

(5.5) ϕ = 50 + E , φ ∈ [ 0 , 6.5 ] , E , φ ∈ [ 6.5 , 15 ] .

For large values of φ , the potential V ( φ ) in equation (5.1) vanishes faster than l ( l + 1 ) φ 2 . Thus, equation (5.1) reduces to:

(5.6) p ″ + γ 2 − l ( l + 1 ) φ 2 = 0 ,

which has the two linearly independent solutions, γ φ j l ( γ x ) and γ φ N l ( γ φ ) , where j l ( γ φ ) and N l ( γ φ ) are the Bessel and Neumann spherical functions, respectively. Hence, as φ → ∞ , equation (5.1) has a solution in asymptotic form as:

p ≈ A γ φ j l ( γ x ) − B γ φ N l ( γ φ ) ≈ A C sin γ φ − l π 2 + tan ( δ l ) cos γ φ − l π 2 ,

where δ l represents the phase shift that can be determined by the formula:

(5.7) δ l = arctan y ( φ i ) S ( φ i + 1 ) − y ( φ i + 1 ) S ( φ i ) y ( φ i + 1 ) C ( φ i + 1 ) − y ( φ i ) C ( φ i ) ,

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 Q T ( 8 ) .
The explicit 10-step tenth order linear symmetric method with constant coefficients developed in [3]; denote it as Q T ( 10 ) .
The explicit 12-step twelfth order linear symmetric method with constant coefficients developed in [3]; denote it by Q T ( 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 2 S t .
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 3 S t .
The three-stage fifth order two-derivative Runge-Kutta method, developed in [43]; denote it by Y 2 D .
The fourth order three derivative Runge-Kutta method with vanished phase-lag and its first derivative produced in [44]; denote it as Y 3 D .
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.
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 δ l using the two energy values, E 1 = 341.495874 and E 2 = 989.701916 . Then, the absolute error M E R = ∣ log ( E R ) ∣ :

E R = π 2 − δ l , comp ,

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.

$Figure 6 Methods efficiency using E 1 = 341.495874 {E}_{1}=341.495874 .$

Figure 6

Methods efficiency using E 1 = 341.495874 .

$Figure 7 Methods efficiency using E 2 = 989.701916 {E}_{2}=989.701916 .$

Figure 7

Methods efficiency using E 2 = 989.701916 .

6 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.

Acknowledgments

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.

References

[1] R. M. Thomas, Phase properties of high order almost P-stable formulae, BIT 24 (1984), 225–238. 10.1007/BF01937488Search in Google Scholar

[2] 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

[3] D. G. Quinlan and S. Tremaine, Symmetric multistep methods for the numerical integration of planetary orbits, Astron. J. 100 (1990), no. 5, 1694–1700. 10.1086/115629Search in Google Scholar

[4] J. Coleman and L. Ixaru, P-stability and exponential-fitting methods for y″=f(x,y), IMA J. Numer. Anal. 16 (1996), 179–199. 10.1093/imanum/16.2.179Search in Google Scholar

[5] G. D. Quinlan, Resonances and instabilities in symmetric multi-step methods, arXiv:astro-ph/9901136, (1999). Search in Google Scholar

[6] Z. Wang, P-stable linear symmetric multistep methods for periodic initial-value problems, Comput. Phys. Commun. 171 (2005), 162–174. 10.1016/j.cpc.2005.05.004Search in Google Scholar

[7] 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

[8] H. Van de Vyver, Phase-fitted and amplification-fitted two-step hybrid methods for y″=f(x,y), J. Comput. Appl. Math. 209 (2007), no. 1, 33–53. 10.1016/j.cam.2006.10.025Search in Google Scholar

[9] 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

[10] E. Stiefel and D. G. Bettis, Stabilization of Cowell’s method, Numer. Math. 13 (1969), 154–175. 10.1007/BF02163234Search in Google Scholar

[11] 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

[12] 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

[13] 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

[14] 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

[15] 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

[16] 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

[17] 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

[18] 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

[19] 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

[20] 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

[21] 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

[22] 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

[23] 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

[24] 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

[25] 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

[26] 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

[27] 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

[28] 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

[29] 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

[30] 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

[31] 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

[32] J. D. Lambert and I. A. Watson, Symmetric multistep methods for periodic initial value problems, IMA J. Appl. Math. 18 (1976), no. 2, 189–202. 10.1093/imamat/18.2.189Search in Google Scholar

[33] 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

[34] 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

[35] T. E. Simos, On the explicit four-step methods with vanished phase-lag and its first derivative, Appl. Math. Inf. Sci. 8 (2014), no. 2, 447–458. 10.12785/amis/080201Search in Google Scholar

[36] 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

[37] 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

[38] S. Obaidat and S. Mesloub, A new explicit four-step symmetric method for solving Schrödingeras equation, Mathematics 7 (2019), no. 11, 1124. 10.3390/math7111124Search in Google Scholar

[39] 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

[40] 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

[41] 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

[42] J. R. Dormand and P. J. Prince, A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 6 (1980), 19–26. 10.1016/0771-050X(80)90013-3Search in Google Scholar

[43] Y. Yang, Y. Fang, and X. You, Modified two-derivative Runge-Kutta methods for the Schrödinger equation, J. Math. Chem. 56 (2017), 799–812. 10.1007/s10910-017-0830-0Search in Google Scholar

[44] 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

[45] J. R. Dormand, M. E. El-Mikkawy, and P. J. Prince, Families of Runge-Kutta-Nyström formulae, IMA J. Numer. Anal. 7 (1987), 235–250. 10.1093/imanum/7.2.235Search in Google Scholar

[46] 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

[47] 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

[48] 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

Received: 2020-07-23

Revised: 2020-11-27

Accepted: 2021-01-09

Published Online: 2021-05-03

This work is licensed under the Creative Commons Attribution 4.0 International License.

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

Abstract

1 Introduction

2 Preliminaries

Definition 2.1

Definition 2.2

Definition 2.3

Theorem 2.4

3 The method development

Proposition 3.1

4 The method error and stability analysis

4.1 Comparative error analysis

Case 1

Case 2

4.2 Stability analysis

5 Numerical results and discussion

6 Conclusions

Acknowledgments

References

Journal and Issue

Articles in the same Issue