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

Authors: Licheng Zhang; Theodore E. Simos
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

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.


---

## Body

## 1. 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:(1)d2dx2+ki2-lili+1x2-Viiyij=∑m=1NVimymj,where 1≤i≤N and m≠i and the boundary conditions are as follows:(2)yij=0atx=0(3)yij~kixjlikixδij+kikj1/2Kijkixnlikix,where jl(x) and nl(x) are the spherical Bessel and Neumann functions. We will examine the case in which all channels are open (see [1]).Defining a matrixK′ and diagonal matrices M, N by (see [1])(4)Kij′=kikj1/2Kij,Mij=kixjlikixδij,Nij=kixnlikixδijwe obtain a new form of the asymptotic condition (3): (5)y~M+NK′.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], the entrance channel by the quantum numbers (j,l), the exit channels by (j′,l′), and the total angular momentum by J=j+l=j′+l′, we find that (6)d2dx2+kj′j2-l′l′+1x2yj′l′Jjlx=2μħ2∑j′′∑l′′j′l′;JVj′′l′′;Jyj′′l′′Jjlx,where (7)kj′j=2μħ2E+ħ22Ijj+1-j′j′+1.E is the kinetic energy of the incident particle in the center-of-mass system, I 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.The multistep finite difference method has the general form(8)∑i=-mmciyn+i=h2∑i=-mmbifxn+i,yn+i,where xii=-mm are distinct points within the integration area and h given by h=|xi+1-xi|,i=1-m(1)m-1 isthe step size or step length of the integration.Remark 1.
A method (8) is called symmetric multistep method or symmetric 2m-step method if c-i=ci and b-i=bi,i=0(1)m.If we apply the symmetric2m-step method (i=-m(1)m), to the scalar test equation(9)y′′=-ϕ2y,we obtain the difference equation(10)Amvyn+m+⋯+A1vyn+1+A0vyn+A1vyn-1+⋯+Amvyn-m=0and the associated characteristic equation(11)Amvλm+⋯+A1vλ+A0v+A1vλ-1+⋯+Amvλ-m=0,where v=ϕh, h is the step length, and Aj(v)j=0(1)k are polynomials of v.We give some definitions.Definition 2 (see [2]).
For a symmetric2m-step method with characteristic equation given by (11) one will say that it has an interval of periodicity (0,v02) if, for all v∈(0,v02), the roots λi,i=1(1)2m, of (11) satisfy(12)λ1=eiθv,λ2=e-iθv,λi≤1,i=312m,where θ(v) is a real function of v.Definition 3 (see [2]).
One calls P-stable method a multistep method if its interval of periodicity is equal to(0,∞).Definition 4.
One calls singularly P-stable method a multistep method if its interval of periodicity is equal to(0,∞)-S (where S is a set of distinct points).Definition 5 (see [3, 4]).
For a symmetric2m-step method with the characteristic equation given by (11), one defines as the phase-lag the leading term in the expansion of (13)t=v-θv.Then, if the quantity t=O(vq+1) as v→∞, the order of the phase-lag is q.Definition 6 (see [5]).
One calls symmetric2m-step methodphase-fitted if its phase-lag is equal to zero.Theorem 7 (see [3]).
The symmetric2m-step method with characteristic equation given by (11) has phase-lag order q and phase-lag constant c given by(14)-cvq+2+Ovq+4=2Amvcos⁡mv+⋯+2Ajvcos⁡jv+⋯+A0v2m2Amv+⋯+2j2Ajv+⋯+2A1v.


## 2. 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:(15)y^n=yn-a0h2fn+1-2fn+fn-1,yˇn=yn-a1h2fn+1-2f^n+fn-1-2a2h2f^n,y^n+1/2=12yn+yn+1-h2a3fˇn+18-a3fn+1,y^n-1/2=12yn+yn-1-h2a3fˇn+18-a3fn-1,yn+1+a4yn+yn-1=h2b1fn+1+fn-1+b0fn+b2f^n+1/2+f^n-1/2,where a0=45469/862066800, a1=-2793/26878564, a2=-86919/13439282, a3=6719641/52720800, fn+i=y′′xn+i,yn+i,i=-1(1/2)1, f^n=y′′xn,y^n, fˇn=y′′xn,y˘n, and a4, bj,j=0(1)2, are free parameters.Application of the above-mentioned method (15) to the scalar test equation (9) leads to the difference equation (10) and the characteristic equation (11) with m=1 and(16)A1v=1+v2b1+12b2-v4b2a3-18+2v6b2a1a3+4v8b2a0a3a2-a1,A0v=a4+v2b0+b2+2v4b2a3+4v6b2a3a2-a1-8v8b2a0a3a2-a1.If we request the new method (15) to have vanished phase-lag and its first, second, and third derivatives, then the following system of equations is obtained:(17)Phase-Lag  (PL)=T0Tdenom=0,First  Derivative  of  the  Phase-Lag=T1Tdenom2=0,Second  Derivative  of  the  Phase-Lag=T2Tdenom3=0,Third  Derivative  of  the  Phase-Lag=T3Tdenom4=0,where Tj,j=0(1)3, are given in Appendix A.Now solving the system of (17) we determine the other coefficients of the new obtained three-stage two-step method:(18)b0=1136407T4Tdenom1,b1=1136407T5Tdenom1,b2=145469T6Tdenom1,a4=T7Tdenom2,where the formulae Tj,j=4(1)7, Tdenom1, and Tdenom2 are given in Appendix B.Additionally to the above formulae for the coefficientsa4, bj,j=0(1)2, we give also the following Taylor series expansions of these coefficients, for the case of heavy cancelations for some values of v in the formulae given by (18): (19)b0=7320563882+268231049v10262060331977892448-473946402833v12224533292438658249446400-841146961608850447v1426876635104907392458734080000-34711252155550520483v1623840650403457053406595478323200-450175574700292070703419v187016973932030008047188109768345600000+⋯,b1=51911383292+268231049v101572361991867354688-61738277963993v121347199754631949496678400-559868876269811659v14161259810629444354752404480000-602760578937670659833v163576097560518558010989321748480000-526438577422138492034251v1884203687184360096566257317220147200000+⋯,b2=-1997095823-268231049v10393090497966838672+15790029293123v12336799938657987374169600+181541151070183v1440314952657361088688101120+21664254981643052111v1689402439012963950274733043712000+8364955394669318397937v18842036871843600965662573172201472000+⋯,a4=-2+134317v1615752663892725760+1170023347v184217775757277322240000+⋯.The behavior of the coefficients is given in Figure 1.Figure 1
Behavior of the coefficients of the new obtained method given by (18) for several values of v=ϕh.The local truncation error of the new developed five-stage two-step method (15) with the coefficients given by (18)-(19), which is indicated as NM2S5S3DV, is given by(20)LTENM2S5S3DV=134317551343236245401600h16yn16+56ϕ10yn6+140ϕ12yn4+120ϕ14yn2+35ϕ16yn+Oh18.


## 3. Analysis of the Method

### 3.1. Error Analysis

The test equation used for the local truncation error (LTE) analysis is given by(21)y′′x=Vx-Vc+Gyx,where V(x) is a potential function, Vc is a constant value approximation of the potential for the specific x, and G=Vc-E and E is the energy. This is the radial time independent Schrödinger equation.We will study the following methods.


#### 3.1.1. Classical Method (i.e., Method (15) with Constant Coefficients)

(22) L T E C L = 134317 551343236245401600 h 16 y n 16 + O h 18 .


#### 3.1.2. The Five-Stage Two-Step Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives Developed in Section3

(23) L T E N M 2 S 5 S 3 D V = 134317 551343236245401600 h 16 y n 16 + 56 ϕ 10 y n 6 + 140 ϕ 12 y n 4 + 120 ϕ 14 y n 2 + 35 ϕ 16 y n + O h 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(24)LTE=hp∑i=0maiGi,where ai are constant numbers (classical methods) or formulae of ϕ (fitted methods) and p is the algebraic order of the specific method.Two cases of the parameterG are studied:(i)
The Energy and the Potential Are Closed to Each Other. Consequently, G=Vc-E≈0⇒Gi=0,i=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 Section3 are the same since the formulae of the LTE are free from G (i.e., LTE=hpa0 in (24)) and the free from G terms (i.e., the terms of the formulae which do not have the quantity G) 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 G, the methods are of comparable accuracy.(ii)
The Energy and the Potential Have Big Difference. Consequently, G≫0 or G≪0 and G is a large number. For this case the most accurate method is the method with the minimum power of G in the formulae of LTE (i.e., the method with the highest accuracy is the method with minimum i in (24)).We give now the asymptotic expansions of the local truncation errors.


#### 3.1.3. Classical Method

(25) L T E C L = 134317 551343236245401600 h 16 y x G 8 + ⋯ + O h 18 .


#### 3.1.4. The Four-Stage Two-Step P-Stable Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives Developed in Section3

(26) L T E N M 2 S 5 S 3 D V = 134317 9845414932953600 h 16 20 g x y x d 2 d x 2 g x + 15 d d x g x 2 y x + 10 d 3 d x 3 g x d d x y x + 31 d 4 d x 4 g x y x G 5 + ⋯ + O h 18 .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 G.
(ii) Fourteenth-algebraic order five-stage two-step method with vanished phase-lag and its first, second, and third derivatives developed in Section3: for this method the error increases as the fifth power of G.
Therefore, for large values of|G|=Vc-E, 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.


### 3.2. Stability and Interval of Periodicity Analysis

Let us define the scalar test equation for the stability and interval of periodicity analysis:(27)y′′=-ω2y.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 (28)A1s,vyn+1+yn-1+A0s,vyn=0,where(29)A1s,v=1+188b1+4b2s2+18b2-8a3+1s4+4a3b2v2a0a2+12a1s6-4a0a1a3b2s8,A0s,v=a4+b0+b2s2+4a3b2v2a2+12s4-8a3b2v2a0a2+12a1s6+8a0a1a3b2s8,s=ωh, and v=ϕh.Taking into account the coefficientsai,i=0,1, and bi,i=0(1)2, and their substitution into formulae (29), the new stability polynomials are produced: (30)A1s,v=T8Tdenom3,A0s,v=T9Tdenom3,where the formulae Tk,k=8,9, and Tdenom3 are given in Appendix D.s-v plane of the new produced method is shown in Figure 2.Figure 2
s-v plane of the new obtained five-stage symmetric two-step fourteenth-algebraic order method with vanished phase-lag and its first, second, and third derivatives.Remark 10.
Observings-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 thes-v plane, where s=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 Section3(i)
is of fourteenth algebraic order,(ii)
has vanished the phase-lag and its first, second, and third derivatives,(iii)
has an interval of periodicity equal to0,24 (when s=v).


## 3.1. Error Analysis

The test equation used for the local truncation error (LTE) analysis is given by(21)y′′x=Vx-Vc+Gyx,where V(x) is a potential function, Vc is a constant value approximation of the potential for the specific x, and G=Vc-E and E is the energy. This is the radial time independent Schrödinger equation.We will study the following methods.


### 3.1.1. Classical Method (i.e., Method (15) with Constant Coefficients)

(22) L T E C L = 134317 551343236245401600 h 16 y n 16 + O h 18 .


### 3.1.2. The Five-Stage Two-Step Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives Developed in Section3

(23) L T E N M 2 S 5 S 3 D V = 134317 551343236245401600 h 16 y n 16 + 56 ϕ 10 y n 6 + 140 ϕ 12 y n 4 + 120 ϕ 14 y n 2 + 35 ϕ 16 y n + O h 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(24)LTE=hp∑i=0maiGi,where ai are constant numbers (classical methods) or formulae of ϕ (fitted methods) and p is the algebraic order of the specific method.Two cases of the parameterG are studied:(i)
The Energy and the Potential Are Closed to Each Other. Consequently, G=Vc-E≈0⇒Gi=0,i=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 Section3 are the same since the formulae of the LTE are free from G (i.e., LTE=hpa0 in (24)) and the free from G terms (i.e., the terms of the formulae which do not have the quantity G) 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 G, the methods are of comparable accuracy.(ii)
The Energy and the Potential Have Big Difference. Consequently, G≫0 or G≪0 and G is a large number. For this case the most accurate method is the method with the minimum power of G in the formulae of LTE (i.e., the method with the highest accuracy is the method with minimum i in (24)).We give now the asymptotic expansions of the local truncation errors.


### 3.1.3. Classical Method

(25) L T E C L = 134317 551343236245401600 h 16 y x G 8 + ⋯ + O h 18 .


### 3.1.4. The Four-Stage Two-Step P-Stable Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives Developed in Section3

(26) L T E N M 2 S 5 S 3 D V = 134317 9845414932953600 h 16 20 g x y x d 2 d x 2 g x + 15 d d x g x 2 y x + 10 d 3 d x 3 g x d d x y x + 31 d 4 d x 4 g x y x G 5 + ⋯ + O h 18 .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 G.
(ii) Fourteenth-algebraic order five-stage two-step method with vanished phase-lag and its first, second, and third derivatives developed in Section3: for this method the error increases as the fifth power of G.
Therefore, for large values of|G|=Vc-E, 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.


## 3.1.1. Classical Method (i.e., Method (15) with Constant Coefficients)

(22) L T E C L = 134317 551343236245401600 h 16 y n 16 + O h 18 .


## 3.1.2. The Five-Stage Two-Step Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives Developed in Section3

(23) L T E N M 2 S 5 S 3 D V = 134317 551343236245401600 h 16 y n 16 + 56 ϕ 10 y n 6 + 140 ϕ 12 y n 4 + 120 ϕ 14 y n 2 + 35 ϕ 16 y n + O h 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(24)LTE=hp∑i=0maiGi,where ai are constant numbers (classical methods) or formulae of ϕ (fitted methods) and p is the algebraic order of the specific method.Two cases of the parameterG are studied:(i)
The Energy and the Potential Are Closed to Each Other. Consequently, G=Vc-E≈0⇒Gi=0,i=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 Section3 are the same since the formulae of the LTE are free from G (i.e., LTE=hpa0 in (24)) and the free from G terms (i.e., the terms of the formulae which do not have the quantity G) 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 G, the methods are of comparable accuracy.(ii)
The Energy and the Potential Have Big Difference. Consequently, G≫0 or G≪0 and G is a large number. For this case the most accurate method is the method with the minimum power of G in the formulae of LTE (i.e., the method with the highest accuracy is the method with minimum i in (24)).We give now the asymptotic expansions of the local truncation errors.


## 3.1.3. Classical Method

(25) L T E C L = 134317 551343236245401600 h 16 y x G 8 + ⋯ + O h 18 .


## 3.1.4. The Four-Stage Two-Step P-Stable Method with Vanished Phase-Lag and Its First, Second, and Third Derivatives Developed in Section3

(26) L T E N M 2 S 5 S 3 D V = 134317 9845414932953600 h 16 20 g x y x d 2 d x 2 g x + 15 d d x g x 2 y x + 10 d 3 d x 3 g x d d x y x + 31 d 4 d x 4 g x y x G 5 + ⋯ + O h 18 .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 G.
(ii) Fourteenth-algebraic order five-stage two-step method with vanished phase-lag and its first, second, and third derivatives developed in Section3: for this method the error increases as the fifth power of G.
Therefore, for large values of|G|=Vc-E, 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.


## 3.2. Stability and Interval of Periodicity Analysis

Let us define the scalar test equation for the stability and interval of periodicity analysis:(27)y′′=-ω2y.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 (28)A1s,vyn+1+yn-1+A0s,vyn=0,where(29)A1s,v=1+188b1+4b2s2+18b2-8a3+1s4+4a3b2v2a0a2+12a1s6-4a0a1a3b2s8,A0s,v=a4+b0+b2s2+4a3b2v2a2+12s4-8a3b2v2a0a2+12a1s6+8a0a1a3b2s8,s=ωh, and v=ϕh.Taking into account the coefficientsai,i=0,1, and bi,i=0(1)2, and their substitution into formulae (29), the new stability polynomials are produced: (30)A1s,v=T8Tdenom3,A0s,v=T9Tdenom3,where the formulae Tk,k=8,9, and Tdenom3 are given in Appendix D.s-v plane of the new produced method is shown in Figure 2.Figure 2
s-v plane of the new obtained five-stage symmetric two-step fourteenth-algebraic order method with vanished phase-lag and its first, second, and third derivatives.Remark 10.
Observings-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 thes-v plane, where s=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 Section3(i)
is of fourteenth algebraic order,(ii)
has vanished the phase-lag and its first, second, and third derivatives,(iii)
has an interval of periodicity equal to0,24 (when s=v).


## 4. Numerical Results

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


### 4.1. 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., [1–55]).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 inyn+1L is estimated by(31)LTE=yn+1H-yn+1L,where yn+1L gives the lower algebraic order solution which is obtained using the twelfth-algebraic order method developed in [54] and yn+1H 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)
IfLTE<acc then the step size is duplicated; that is, hn+1=2hn.(ii)
Ifacc≤LTE≤100acc then the step size remains stable; that is, hn+1=hn.(iii)
If100acc<LTE then the step size is halved and the step is repeated; that is, hn+1=1/2hn.In the above, hn is the step length used for the nth 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 thanacc we use the lower algebraic order solution yn+1L it is the higher algebraic order solution yn+1H which is accepted at each point as integration.


### 4.2. 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 V which can be presented as (see for details [1])(32)Vx,k^j′jk^jj=V0xP0k^j′jk^jj+V2xP2k^j′jk^jj.Based on (32), the coupling matrix element is given by(33)j′l′;JVj′′l′′;J=δj′j′′δl′l′′V0x+f2j′l′,j′′l′′;JV2x,where f2 coefficients are obtained from formulas given by Bernstein et al. [56] and k^j′j is a unit vector parallel to the wave vector kj′j and Pi,i=0,2, are Legendre polynomials (see for details [57]). The above leads to the new form of boundary conditions:(34)yj′l′Jjlx=0atx=0(35)yj′l′Jjlx~δjj′δll′exp⁡-ikjjx-12lπ,-kikj1/2SJjl;j′l′exp⁡ikj′jx-12l′π,where(36)S=I+iKI-iK-1.We will use the variable step method described in Section4.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:(37)2μħ2=1000.0,μI=2.351,E=1.1,V0x=1x12-21x6,V2x=0.2283V0x.For our test we takeJ=6 and we will consider excitation of the rotator from j=0 state to levels up to j′=2, 4, and 6 which is equivalent to sets offour, nine, and sixteen coupled differential equations, respectively. Using the methodology fully described by Bernstein [57] and Allison [1], we consider the potential infinite for x<x0. The wave functions then tend to zero and then the boundary conditions (34) are given by(38)yj′l′Jjlx0=0.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].(iii)
Method III. The embedded Runge-Kutta Dormand and Prince method 5(4) [49].(iv)
Method IV. The embedded Runge-Kutta method ERK4(2) developed in Simos [41].(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 two-step method developed in this paper.The results are presented in Table1. 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.Table 1
Coupled differential equations. Real time of computation (in seconds) (RTC) and maximum absolute error (MErr) to calculateS2 for the variable step methods, Method I–Method X. acc=10-6. We note that hmax is the maximum step size. N indicates the number of equations of the set of coupled differential equations.
Method N h m a x RTC MErr Method I 4 0.014 3.25 1.2 × 1 0 - 3 9 0.014 23.51 5.7 × 1 0 - 2 16 0.014 99.15 6.8 × 1 0 - 1  Method II 4 0.056 1.55 8.9 × 1 0 - 4 9 0.056 8.43 7.4 × 1 0 - 3 16 0.056 43.32 8.6 × 1 0 - 2  Method III 4 0.007 45.15 9.0 × 1 0 0 9 16  Method IV 4 0.112 0.39 1.1 × 1 0 - 5 9 0.112 3.48 2.8 × 1 0 - 4 16 0.112 19.31 1.3 × 1 0 - 3  Method V 4 0.448 0.14 3.4 × 1 0 - 7 9 0.448 1.37 5.8 × 1 0 - 7 16 0.448 9.58 8.2 × 1 0 - 7  Method VI 4 0.448 0.09 2.9 × 1 0 - 7 9 0.448 1.10 4.5 × 1 0 - 7 16 0.448 8.57 7.4 × 1 0 - 7  Method VII 4 0.448 0.06 1.3 × 1 0 - 7 9 0.448 1.04 1.7 × 1 0 - 7 16 0.448 7.58 2.9 × 1 0 - 7  Method VIII 4 0.448 0.04 8.8 × 1 0 - 8 9 0.448 1.02 9.2 × 1 0 - 8 16 0.448 7.48 8.9 × 1 0 - 8  Method IX 4 0.448 0.04 9.7 × 1 0 - 8 9 0.448 1.01 1.2 × 1 0 - 7 16 0.448 7.15 2.3 × 1 0 - 7  Method X 4 0.896 0.02 7.0 × 1 0 - 8 9 0.896 0.45 5.7 × 1 0 - 8 16 0.896 6.11 6.5 × 1 0 - 8


## 4.1. 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., [1–55]).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 inyn+1L is estimated by(31)LTE=yn+1H-yn+1L,where yn+1L gives the lower algebraic order solution which is obtained using the twelfth-algebraic order method developed in [54] and yn+1H 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)
IfLTE<acc then the step size is duplicated; that is, hn+1=2hn.(ii)
Ifacc≤LTE≤100acc then the step size remains stable; that is, hn+1=hn.(iii)
If100acc<LTE then the step size is halved and the step is repeated; that is, hn+1=1/2hn.In the above, hn is the step length used for the nth 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 thanacc we use the lower algebraic order solution yn+1L it is the higher algebraic order solution yn+1H which is accepted at each point as integration.


## 4.2. 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 V which can be presented as (see for details [1])(32)Vx,k^j′jk^jj=V0xP0k^j′jk^jj+V2xP2k^j′jk^jj.Based on (32), the coupling matrix element is given by(33)j′l′;JVj′′l′′;J=δj′j′′δl′l′′V0x+f2j′l′,j′′l′′;JV2x,where f2 coefficients are obtained from formulas given by Bernstein et al. [56] and k^j′j is a unit vector parallel to the wave vector kj′j and Pi,i=0,2, are Legendre polynomials (see for details [57]). The above leads to the new form of boundary conditions:(34)yj′l′Jjlx=0atx=0(35)yj′l′Jjlx~δjj′δll′exp⁡-ikjjx-12lπ,-kikj1/2SJjl;j′l′exp⁡ikj′jx-12l′π,where(36)S=I+iKI-iK-1.We will use the variable step method described in Section4.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:(37)2μħ2=1000.0,μI=2.351,E=1.1,V0x=1x12-21x6,V2x=0.2283V0x.For our test we takeJ=6 and we will consider excitation of the rotator from j=0 state to levels up to j′=2, 4, and 6 which is equivalent to sets offour, nine, and sixteen coupled differential equations, respectively. Using the methodology fully described by Bernstein [57] and Allison [1], we consider the potential infinite for x<x0. The wave functions then tend to zero and then the boundary conditions (34) are given by(38)yj′l′Jjlx0=0.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].(iii)
Method III. The embedded Runge-Kutta Dormand and Prince method 5(4) [49].(iv)
Method IV. The embedded Runge-Kutta method ERK4(2) developed in Simos [41].(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 two-step method developed in this paper.The results are presented in Table1. 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.Table 1
Coupled differential equations. Real time of computation (in seconds) (RTC) and maximum absolute error (MErr) to calculateS2 for the variable step methods, Method I–Method X. acc=10-6. We note that hmax is the maximum step size. N indicates the number of equations of the set of coupled differential equations.
Method N h m a x RTC MErr Method I 4 0.014 3.25 1.2 × 1 0 - 3 9 0.014 23.51 5.7 × 1 0 - 2 16 0.014 99.15 6.8 × 1 0 - 1  Method II 4 0.056 1.55 8.9 × 1 0 - 4 9 0.056 8.43 7.4 × 1 0 - 3 16 0.056 43.32 8.6 × 1 0 - 2  Method III 4 0.007 45.15 9.0 × 1 0 0 9 16  Method IV 4 0.112 0.39 1.1 × 1 0 - 5 9 0.112 3.48 2.8 × 1 0 - 4 16 0.112 19.31 1.3 × 1 0 - 3  Method V 4 0.448 0.14 3.4 × 1 0 - 7 9 0.448 1.37 5.8 × 1 0 - 7 16 0.448 9.58 8.2 × 1 0 - 7  Method VI 4 0.448 0.09 2.9 × 1 0 - 7 9 0.448 1.10 4.5 × 1 0 - 7 16 0.448 8.57 7.4 × 1 0 - 7  Method VII 4 0.448 0.06 1.3 × 1 0 - 7 9 0.448 1.04 1.7 × 1 0 - 7 16 0.448 7.58 2.9 × 1 0 - 7  Method VIII 4 0.448 0.04 8.8 × 1 0 - 8 9 0.448 1.02 9.2 × 1 0 - 8 16 0.448 7.48 8.9 × 1 0 - 8  Method IX 4 0.448 0.04 9.7 × 1 0 - 8 9 0.448 1.01 1.2 × 1 0 - 7 16 0.448 7.15 2.3 × 1 0 - 7  Method X 4 0.896 0.02 7.0 × 1 0 - 8 9 0.896 0.45 5.7 × 1 0 - 8 16 0.896 6.11 6.5 × 1 0 - 8


## 5. 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).



---

*Source: 8181927-2016-11-23.xml*