Optimizing a Hybrid Two-Step Method for the Numerical Solution of the Schrödinger Equation and Related Problems with Respect to Phase-Lag

We use a methodology of optimization of the efficiency of a hybrid two-step 
method for the numerical solution of the radial Schrodinger equation and related problems with 
periodic or oscillating solutions. More specifically, we study how the vanishing of the phase-lag 
and its derivatives optimizes the efficiency of the hybrid two-step method.


Introduction
In this paper, we investigate the numerical solution of systems of second order differential equations of the form y x f x, y , 1.1 while the following initial conditions hold: where f is independent of y x .

Phase-Lag Analysis of Symmetric Multistep Methods
For the numerical solution of the initial value problem: consider a multistep method with m steps which can be used over the equally spaced intervals {x i } m i 0 ∈ a, b and h |x i 1 − x i |, i 0 1 m − 1.If the method is symmetric, then a i a m−i and b i b m−i , i 0 1 m/2 .When a symmetric 2m-step method, that is, for i −m 1 m, is applied to the scalar test equation: a difference equation of the form: is obtained, where H ωh, h is the step length, and A 0 H , A 1 H , . .., A m H are polynomials of H ωh.
The characteristic equation associated with 2.3 is given by: Theorem 2.1 see 3, 23 .The symmetric 2m-step method with characteristic equation given by 2.4 has phase-lag order q and phase-lag constant c given by The formula mentioned in the above theorem is a direct method for the computation of the phase-lag of any symmetric 2m-step method.

The General Family of Methods
Consider the following family of hybrid two-step methods see 24 :

3.1
The above-mentioned method belongs to the families of hybrid Runge-Kutta type symmetric two-step methods for the numerical solution of problems of the form u f x, u .In the above general form, the coefficient b 0 and b 1 are free parameters.In the above method, h Journal of Applied Mathematics is the step size of the integration and n is the number of steps, that is, y n is the approximation of the solution on the point x n , x n x 0 nh, and x 0 is the initial value point.

The Optimized Hybrid Method of the Family with Vanished Phase-Lag and Its First Derivative
Consider the method 3.1 .
If we apply the method 3.1 to the scalar test equation 2.2 , we obtain the difference equation 2.3 with m 1 and A j H , j 0, 1 given by Requiring the above method to have its phase-lag vanished and by using the formulae 2.5 for k 1 and 3.2 , we have the following equation: Demanding the method to have the first derivative of the phase-lag vanished as well, we have the equation: where DPL is the first derivative of the phase-lag.
Requiring now the coefficients of the new proposed method to satisfy the equations 3.3 -3.4 , we obtain the following coefficients of the new developed method:

3.5
For some values of |ω|, the formulae given by 3.5 are subject to heavy cancellations.In this case, the following Taylor series expansions should be used:

Error Analysis
We will study the following methods.

4.2
The error analysis is based on the following steps.
i The radial time independent Schr ödinger equation is of the form: ii Based on the paper of Ixaru and Rizea 2 , the function f x can be written in the form: where g x V x −V c g, where V c is the constant approximation of the potential and G v 2 V c − E.
iii We express the derivatives u i n , i 2, 3, 4, . .., which are terms of the local truncation error formulae, in terms of 4.4 .The expressions are presented as polynomials of G.
iv Finally, we substitute the expressions of the derivatives, produced in the previous step, into the local truncation error formulae.
We use the procedure mentioned above and the formulae: We consider two cases in terms of the value of E.
1 The energy is close to the potential, that is, G V c − E ≈ 0. Consequently, the free terms of the polynomials in G are considered only.Thus, for these values of G, the methods are of comparable accuracy.This is because the free terms of the polynomials in G are the same for the cases of the standard method and of the trigonometrically fitted methods.
2 G 0 or G 0. Then |G| is a large number.
Hence, we have the following asymptotic expansions of the Local Truncation Errors.

Standard Method
It holds that

New Method with Vanished Phase-Lag and Its First Derivative (Developed in Section 3.2)
It holds that

4.7
From the above equations, we have the following theorem.
Theorem 4.1.For the standard hybrid two-step method, the error increases as the fourth power of G.For the new method with vanished phase-lag and its first derivative (developed in Section 3.2), the error increases as the second power of G. So, for the numerical solution of the time independent radial Schrödinger equation, the new method with vanished phase-Lag and its first derivative is much more efficient, especially for large values of |G| |V c − E|.

Stability Analysis
Applying the new method to the scalar test equation: we obtain the following difference equation: where where The corresponding characteristic equation is given by.
Definition 5.1 see 25 .A symmetric 2m-step method with the characteristic equation given by 2.4 is said to have an interval of periodicity 0, v 2 0 if, for all v ∈ 0, v 2 0 , the roots z i , i 1, 2 satisfy where ζ v is a real function of zh and v zh.
Definition 5.2 see 25 .A method is called P-stable if its interval of periodicity is equal to 0, ∞ .
Definition 5.3.A method is called singularly almost P-stable if its interval of periodicity is equal to 0, ∞ − S where S is a set of distinct points only when the frequency of the phase fitting is the same as the frequency of the scalar test equation, that is, v H.
In Figure 2, we present the H-v plane for the method developed in this paper.A shadowed area denotes the H-v region where the method is stable, while a white area denotes the region where the method is unstable.Remark 5.4.For the solution of the Schr ödinger equation, the frequency of the exponential fitting is equal to the frequency of the scalar test equation.So, it is necessary to observe the surroundings of the first diagonal of the H-v plane.
In the case that the frequency of the scalar test equation is equal with the frequency of phase fitting, that is, in the case that v H i.e., see the surroundings of the first diagonal of the H-v plane , it is easy to see that the interval of periodicity of the new method developed in Section 3.2 is equal to 0, 39.47841760 .
From the above analysis, we have the following theorem.
Theorem 5.5.The method developed in Section 3.2 is of eighth algebraic order, has the phase-lag and its first derivative equal to zero, and has an interval of periodicity equals to 0, 39.47841760 .
Based on the analysis presented above, we studied the interval of periodicity of some well-known methods mentioned in the previous paragraph.The results presented in the Table 1.

Numerical Results
In order to study the efficiency of the new developed method, we apply it i to the radial time-independent Schr ödinger equation, and ii to a nonlinear orbital problem.

The Radial Schr ödinger Equation
The radial Schr ödinger equation can be presented as q r l l 1 r 2 V r − k 2 q r .6.1 The above equation presents the model for a particle in a central potential field where r is the radial variable see 26, 27 .
In 6.1 , we have the following terms.
i The function W r l l 1 /r 2 V r is called the effective potential.This satisfies W r → 0 as r → ∞.
ii The quantity k 2 is a real number denoting the energy.
iii The quantity l is a given integer representing the angular momentum.iv V is a given function which denotes the potential.
We note here that the models which are given via the radial Schr ödinger equation are boundary-value problems.In these cases, the boundary conditions are q 0 0 6.2 and a second boundary condition, for large values of r, determined by physical considerations.
In order to apply the new obtained method to the radial Schr ödinger equation, the value of parameter ω is needed.In 6.1 , the parameter ω is given by where V r is the potential and E is the energy.

Woods-Saxon Potential
We use as a potential the well-known Woods-Saxon potential which can be written as with y exp r − R 0 /a , u 0 −50, a 0.6, and X 0 7.0.The behaviour of Woods-Saxon potential is shown in Figure 3.
It is well known that for some potentials, such as the Woods-Saxon potential, the definition of parameter ω is given not as a function of r but as based on some critical points which have been defined from the investigation of the appropriate potential see for details 28 .For the purpose of obtaining our numerical results, it is appropriate to choose v as follows see for details 2, 26 : for r ∈ 0, 6.5 − 2h , √ −37.5 E, for r 6.5 − h √ −25 E, for r 6.5 √ −12.5 E, for r 6.5 h √ E, for r ∈ 6.5 2h, 15 .

6.5
For example, in the point of the integration region r 6.5, the value of ω is equal to √ −25 E. So, H ωh √ −25 Eh.In the point of the integration region r 6.5 − 3h, the value of ω is equal to √ −50 E, and so forth.

Radial Schrödinger Equation-The Resonance Problem
We consider the numerical solution of the radial Schr ödinger equation 6.1 in the wellknown case of the Woods-Saxon potential 6.4 .In order to solve this problem numerically, we must approximate the true infinite interval of integration by a finite interval.For the purpose of our numerical illustration, we take the domain of integration as r ∈ 0, 15 .We consider equation 6.1 in a rather large domain of energies, that is, E ∈ 1, 1000 .
In the case of positive energies, E k 2 , the potential decays faster than the term l l 1 /x 2 and the Schr ödinger equation effectively reduces to q r k 2 − l l 1 r 2 q r 0, 6.6 for r greater than some value R.
The above equation has linearly independent solutions krj l kr and krn l kr , where j l kr and n l kr are the spherical Bessel and Neumann functions, respectively.Thus, the solution of equation 6.1 when r → ∞ has the asymptotic form: where δ l is the phase shift that may be calculated from the formula: for r 1 and r 2 distinct points in the asymptotic region we choose r 1 as the right-hand end point of the interval of integration and r 2 r 1 − h with S r krj l kr and C r −krn l kr .Since the problem is treated as an initial-value problem, we need q j , j 0, 1 before starting a two-step method.From the initial condition, we obtain q 0 .The value q 1 is obtained by using high-order Runge-Kutta-Nystr öm methods see 29-31 .With these starting values, we evaluate at r 2 of the asymptotic region the phase shift δ l .
For positive energies, we have the so-called resonance problem.This problem consists either of finding the phase-shift δ l or finding those E, for E ∈ 1, 1000 , at which δ l π/2.We actually solve the latter problem, known as the resonance problem.
The boundary conditions for this problem are q 0 0, q r cos √ Er for large r.6.9 We compute the approximate positive eigenenergies of the Woods-Saxon resonance problem using the following.
i The eighth-order multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT8.
ii The tenth-order multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT10.
iii The twelfth-order multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT12.
iv The fourth-algebraic-order method of Chawla and Rao with minimal phase-lag 33 , which is indicated as Method MCR4.
v The hybrid sixth-algebraic-order method developed by Chawla and Rao with minimal phase-lag 4 , which is indicated as Method MCR6.
vi The standard form of the eighth-algebraic-order method developed in Section 3.2, which is indicated as Method NMCL.with the term standard we mean the method of Section 3.2 with constant coefficients.
vii The new developed hybrid two-step method with vanished phase-lag and its first derivative obtained in Section 3.2 , which is indicated as Method NM.The computed eigenenergies are compared with reference values the reference values are computed using the well-known two-step method of Chawla and Rao 4 with small step size for the integration.In Figures 4 and 5

A Nonlinear Orbital Problem
Consider the nonlinear system of equations:

6.11
The analytical solution of the problem is the following: The system of equations 6.11 has been solved for 0 ≤ x ≤ 10000 and W 10 using the methods mentioned in Section 6.1.
For this problem, we have ω 10.The numerical results obtained for the seven methods mentioned above were compared with the analytical solution.Figure 6 shows the absolute errors Err max defined by Err max log 10 max u x calculated − u x theoretical , v x calculated − v x theoretical , x ∈ 0, 10000 ,

6.13
for several values of the CPU time in seconds .The nonexistence of a value of accuracy digits indicates that for this value of CPU time, accuracy digits is less than 0.

Conclusions
The purpose of this paper was the optimization of the efficiency of a hybrid two-step method for the approximate solution of the radial Schr ödinger equation and related problems.We have described how the methodology of vanishing of the phase-lag and its first derivative optimize the behaviour of the specific numerical method.The results of the application of this methodology was a hybrid two-step method that is very efficient on any problem with oscillating solutions or problems with solutions contain the functions cos and sin or any combination of them.
From the results presented above, we can make the following remarks.
1 The standard form of the eighth-algebraic-order method developed in Section 3.2, which is indicated as Method NMCL, is more efficient than the fourth-algebraicorder method of Chawla and Rao with minimal phase-lag 33 , which is indicated as Method MCR4.
2 The tenth-order multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT10, is more efficient than the fourth-algebraic-order method of Chawla and Rao with minimal phase-lag 33 , which is indicated as Method MCR4.Method QT10 is also more efficient than the eighth-order multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT8.Finally, Method QT10 is also more efficient than the hybrid sixth-algebraicorder method developed by Chawla and Rao with minimal phase-lag 4 , which is indicated as Method MCR6.
3 The twelfth-order multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT12, is more efficient than the tenth-order, multistep method developed by Quinlan and Tremaine 32 , which is indicated as Method QT10.
4 Finally, the new developed hybrid two-step method with vanished phase-lag and its first derivative obtained in Section 3.2 , which is indicated as Method NM, is the most efficient one.
All computations were carried out on an IBM PC-AT compatible 80486 using doubleprecision arithmetic with 16 significant digits accuracy IEEE standard .

1 bFigure 1 :
Figure 1: Behaviour of the coefficients of the new proposed method given by 3.5 for several values of H ωh.

Figure 2 :
Figure 2: v-H plane of the new developed method.

Figure 4 :
Figure 4: Accuracy digits for several values of CPU Time in seconds for the eigenvalue E 2 341.495874.The nonexistence of a value of accuracy digits indicates that for this value of CPU, accuracy digits is less than 0.

Figure 5 :
Figure 5: Accuracy Digits for several values of CPU Time in Seconds for the eigenvalue E 3 989.701916.The nonexistence of a value of Accuracy Digits indicates that for this value of CPU, Accuracy Digits is less than 0.

Figure 6 :
Figure 6: Accuracy digits for several values of CPU time in seconds for the nonlinear orbital problem.The nonexistence of a value of accuracy digits indicates that for this value of CPU time, accuracy digits is less than 0.

Table 1 :
Comparative stability analysis for the methods mentioned in Section 5.
, we present the maximum absolute error Err max |log 10 Err |, where Err |E calculated − E accurate |, 6.10 of the eigenenergies E 2 341.495874 and E 3 989.701916,respectively, for several values of CPU time in seconds .We note that the CPU time in seconds counts the computational cost for each method.