An Efficient Numerical Method for the Solution of the Schrödinger Equation

The development of a new five-stage symmetric two-step fourteenth-algebraic order method with vanished phase-lag and its first, second, and third derivatives is presented in this paper for the first time in the literature. More specifically we will study (1) the development of the new method, (2) the determination of the local truncation error (LTE) of the new method, (3) the local truncation error analysis which will be based on test equation which is the radial time independent Schrödinger equation, (4) the stability and the interval of periodicity analysis of the new developed method which will be based on a scalar test equation with frequency different than the frequency of the scalar test equation used for the phase-lag analysis, and (5) the efficiency of the new obtained method based on its application to the coupled Schrödinger equations.


Introduction
The approximate solution of the close-coupling differential equations of the Schrödinger type is studied in this paper.The above-mentioned problem has the following form: where 1 ≤  ≤  and  ̸ =  and the boundary conditions are as follows: ∼      (  )   + (     ) where   () and   () are the spherical Bessel and Neumann functions.We will examine the case in which all channels are open (see [1]).
Defining a matrix   and diagonal matrices ,  by (see [1]) we obtain a new form of the asymptotic condition (3): In several scientific areas (e.g., quantum chemistry, theoretical physics, material science, atomic physics, and molecular physics) there exists a real problem which is the rotational excitation of a diatomic molecule by neutral particle impact.The mathematical model of this problem can be expressed with close-coupling differential equations arising from the Schrödinger equation.Denoting, as in [1], 2 Advances in Mathematical Physics the entrance channel by the quantum numbers (, ), the exit channels by (  ,   ), and the total angular momentum by  =  +  =   +   , we find that where is the kinetic energy of the incident particle in the centerof-mass system,  is the moment of inertia of the rotator, and  is the reduced mass of the system.The above-described problem will be solved numerically via finite difference method of the form of special multistep method.
If we apply the symmetric 2-step method ( = −(1)), to the scalar test equation we obtain the difference equation and the associated characteristic equation where V = ℎ, ℎ is the step length, and   (V)  = 0(1) are polynomials of V.
We give some definitions.
Definition 4. One calls singularly P-stable method a multistep method if its interval of periodicity is equal to (0, ∞) −  (where  is a set of distinct points).
Definition 5 (see [3,4]).For a symmetric 2-step method with the characteristic equation given by (11), one defines as the phase-lag the leading term in the expansion of Then, if the quantity  = (V +1 ) as V → ∞, the order of the phase-lag is .
Theorem 7 (see [3]).The symmetric 2-step method with characteristic equation given by (11) has phase-lag order  and phase-lag constant  given by

The New Five-Stage Fourteenth-Algebraic Order P-Stable Two-Step Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives
We consider the following family of five-stage symmetric two-step methods: Advances in Mathematical Physics

Analysis of the Method
where () is a potential function,   is a constant value approximation of the potential for the specific , and  =   −  and  is the energy.This is the radial time independent Schrödinger equation.We will study the following methods.ℎ 16 ( (16)    + 56 10  (6)   + 140 12  (4)  + 120 14  (2)   + 35 16   ) +  (ℎ 18 ) . ( If we substitute the higher order derivatives requested in the LTE formulae, which are obtained using the test problem (21), into the LTE expressions, we produce new formulae of LTE which have the general form where   are constant numbers (classical methods) or formulae of  (fitted methods) and  is the algebraic order of the specific method.Two cases of the parameter  are studied: (i) The Energy and the Potential Are Closed to Each Other.Consequently,  =   −  ≈ 0 ⇒   = 0,  = 1, 2, . ... Therefore, the local truncation error for the classical method (constant coefficients) and the local truncation error for the five-stage two-step method with eliminated phase-lag and its first, second, and third derivatives developed in Section 3 are the same since the formulae of the LTE are free from  (i.e., LTE = ℎ   0 in (24)) and the free from  terms (i.e., the terms of the formulae which do not have the quantity ) in the local truncation errors in this case are the same and are given in Appendix C. From the above it is easy to see that, for these values of , the methods are of comparable accuracy.
(ii) The Energy and the Potential Have Big Difference.Consequently,  ≫ 0 or  ≪ 0 and || is a large number.For this case the most accurate method is the method with the minimum power of  in the formulae of LTE (i.e., the method with the highest accuracy is the method with minimum  in ( 24)).
We give now the asymptotic expansions of the local truncation errors.

The Four-Stage Two-Step P-Stable Method with Vanished
Phase-Lag and Its First, Second, and Third Derivatives Developed in Section 3 The above leads us to the following theorem.Theorem 8. (i) Classical method (i.e., method (15) with constant coefficients): for this method the error increases as the eighth power of .
(ii) Fourteenth-algebraic order five-stage two-step method with vanished phase-lag and its first, second, and third derivatives developed in Section 3: for this method the error increases as the fifth power of .
Therefore, for large values of || = |  − |, the new developed fourteenth-algebraic order five-stage two-step method with vanished phase-lag and its first, second, and third derivatives developed in Section 3 is the most accurate method for the numerical solution of the radial Schrödinger equation.

Stability and Interval of Periodicity Analysis.
Let us define the scalar test equation for the stability and interval of periodicity analysis: Remark 9. Comparing the test equations ( 9) and ( 27), we have that the frequency  is not equal to the frequency ; that is,  ̸ = .
The difference equation which is produced after application of the new method (15) with the coefficients given by ( 18)-( 19) to the scalar test equation ( 27) is given by where s (test problem) = ℎ, and V = ℎ.
Taking into account the coefficients   ,  = 0, 1, and   ,  = 0(1)2, and their substitution into formulae (29), the new stability polynomials are produced: where the formulae   ,  = 8, 9, and  denom3 are given in Appendix D. -V plane of the new produced method is shown in Figure 2.
Remark 10.Observing -V region we can define two areas: (i) The shadowed area which defines the area where the method is stable.
(ii) The white area which defines the area where the method is unstable.
For problems of the form of the coupled equations arising from the Schrödinger equation and related problems (in quantum chemistry, material science, theoretical physics, atomic physics, astronomy, astrophysics, physical chemistry, and chemical physics), the area of the region which plays critical role in the stability of the numerical methods is the surroundings of the first diagonal of the -V plane, where  = V.Studying this case, we found that the interval of periodicity is equal to (0, 24).
The above development leads to the following theorem.
Theorem 11.The five-stage symmetric two-step method developed in Section 3 (ii) has vanished the phase-lag and its first, second, and third derivatives, (iii) has an interval of periodicity equal to (0, 24) (when  = V).

Numerical Results
We will apply the new produced method to the approximate solution for coupled differential equations of the Schrödinger type.

Error Estimation.
For the numerical solution of the previously referred problem we will use variable step methods.An algorithm is called variable step when it is based on the change of the step size of the integration using a local truncation error estimation (LTEE) procedure.In the past decades several methods of constant or variable steps have been developed for the approximation of coupled differential equations of the Schrödinger type and related problems (see, e.g., ).
In our numerical tests we use local estimation procedure which is based on an embedded pair and on the fact that we have better approximation for the problems with oscillatory and/or periodical solutions having maximal algebraic order of the methods.
The local truncation error in   +1 is estimated by where   +1 gives the lower algebraic order solution which is obtained using the twelfth-algebraic order method developed in [54] and   +1 gives the higher order solution which is obtained using the five-stage symmetric two-step method of fourteenth algebraic order with vanished phase-lag and its first, second, and third derivatives developed in Section 3.
In our numerical tests the changes of the step sizes are reduced on duplication of step sizes.We use the following procedure: (i) If LTE < acc then the step size is duplicated; that is, ℎ +1 = 2ℎ  .
(iii) If 100 acc < LTE then the step size is halved and the step is repeated; that is, In the above, ℎ  is the step length used for the th step of the integration and acc is the requested accuracy of the local truncation error (LTE).
Remark 12.The local extrapolation technique is also used; that is, while for a local truncation error estimation less than acc we use the lower algebraic order solution   +1 it is the higher algebraic order solution   +1 which is accepted at each point as integration.

Coupled Differential
Equations.We will present the numerical solution of the coupled Schnrödinger equations (1) using the boundary conditions ( 2) and (3).Mathematical models of this form are observed in many real problems in quantum chemistry, material science, theoretical physics, atomic physics, physical chemistry and chemical physics, and so forth.The methodology fully described in [1] will be used for the approximate solution for this problem.
Equation (1) contains the potential  which can be presented as (see for details [1]) Based on (32), the coupling matrix element is given by where  2 coefficients are obtained from formulas given by Bernstein et al. [56] and k   is a unit vector parallel to the wave vector k    and   ,  = 0, 2, are Legendre polynomials (see for details [57]).The above leads to the new form of boundary conditions: where We will use the variable step method described in Section 4.1 in order to compute the cross sections for rotational excitation of molecular hydrogen by impact of various heavy particles.
We use the following parameters in our example: For our test we take  = 6 and we will consider excitation of the rotator from  = 0 state to levels up to   = 2, 4, and 6 which is equivalent to sets of four, nine, and sixteen coupled differential equations, respectively.Using the methodology fully described by Bernstein [57] and Allison [1], we consider the potential infinite for  <  0 .The wave functions then tend to zero and then the boundary conditions (34) are given by For comparison purposes the following variable step methods are used: (i) Method I.The iterative Numerov method of Allison [1].
(ii) Method II.The variable step method of Raptis and Cash [47].
Advances in Mathematical Physics (v) Method V.The embedded symmetric two-step method developed in [50].
(vi) Method VI.The embedded symmetric two-step method developed in [52].
(vii) Method VII.The embedded symmetric two-step method developed in [52].
(viii) Method VIII.The embedded symmetric two-step method developed in [53].
(ix) Method IX.The embedded symmetric two-step method developed in [54].
(x) Method X.The developed embedded symmetric twostep method developed in this paper.
The results are presented in Table 1.More specifically, we present (1) the requested real time of computation by the variable step methods mentioned above in order to calculate the square of the modulus of S matrix for sets of 4, 9, and 6 coupled differential equations and (2) the maximum error in the computation of the square of the modulus of S matrix.

Conclusions
In this paper we introduce, for the first time in the literature, a new five-stage symmetric two-step fourteenth-algebraic order family of methods and we produced a method of the family with vanished phase-lag and its first, second, and third derivatives.In this paper, (1) we introduced the new family of methods, (2) we developed the new method of the new family with vanished phase-lag and its first, second, and third derivatives, (3) we analyzed the local truncation error, (4) we analyzed the interval of periodicity and the stability of the new developed method, (5) finally, we analyzed the efficiency of the new obtained method on the numerical solution for coupled Schrödinger equations.
From the analysis presented above, we conclude that the new developed method is much more efficient than the known ones for the numerical solution for the coupled Schrödinger equations.
All computations were carried out on IBM PC-AT compatible 80486 using double precision arithmetic with 16 significant digits' accuracy (IEEE standard).

8 Figure 1 :
Figure1: Behavior of the coefficients of the new obtained method given by(18) for several values of V = ℎ.

Figure 2 :
Figure2: -V plane of the new obtained five-stage symmetric twostep fourteenth-algebraic order method with vanished phase-lag and its first, second, and third derivatives.

Table 1 :
Coupled differential equations.Real time of computation (in seconds) (RTC) and maximum absolute error (MErr) to calculate || 2 for the variable step methods, Method I-Method X. acc = 10 −6 .We note that ℎ max is the maximum step size. indicates the number of equations of the set of coupled differential equations.