The Use of Phase Lag and Amplification Error Derivatives for the Construction of a Modified Runge-Kutta-Nyström Method

Abstract

A new modified Runge-Kutta-Nyström method of fourth algebraic order is developed. The new modified RKN method is based on the fitting of the coefficients, due to the nullification not only of the phase lag and of the amplification error, but also of their derivatives. Numerical results indicate that the new modified method is much more efficient than other methods derived for solving numerically the Schrödinger equation.

1. Introduction

Among the most commonly used methods in the numerical integration of second order differential equations, in which the first derivative terms are omitted, are Runge-Kutta (-Nyström) methods. These methods have been used widely due to their simplicity and their accuracy and mainly because they are one-step methods and thus they require no additional starting values.

In the last decades many researchers developed optimized Runge-Kutta (-Nyström) methods, based mostly on the exponential fitting and the phase lag properties [1–14].

In the last years, a new methodology has been developed, which is based on the phase lag derivatives. Researchers that have been using the mentioned methodology achieve higher accuracy at their methods [6, 8, 10, 11, 13, 14].

In the present paper a new modified Runge-Kutta-Nyström method is constructed. The new method contains four additional variable coefficients (in comparison with the classical RKN method), which depend on 𝑧=𝑤ℎ, where 𝑤 is the dominant frequency of the problem and ℎ is the step length of integration. In order to evaluate the above coefficients, the new method combines the nullification of phase lag and amplification factor with the nullification of their derivatives.

The new modified RKN method that has been obtained will be used for the numerical solution of second order differential equations and more specifically for the numerical integration of the radial Schrödinger equation.

2. Modified Runge-Kutta-Nyström Method

The 𝑚-stage modified RKN method for the equation 𝑑²𝑦(𝑡)𝑑𝑡²=𝑓𝑡,𝑦(𝑡) (1) is of the form 𝑦_𝑛=𝑦_𝑛−1𝑔_𝑚+ℎ𝑦^′_𝑛−1+ℎ^2𝑚_𝑖=1𝑏_𝑖𝑓𝑡_𝑛−1+𝑐_𝑖ℎ,𝑓_𝑖,𝑦^′_𝑛=𝑦^′_𝑛−1+ℎ^𝑚_𝑖=1𝑏^′_𝑖𝑓𝑡_𝑛−1+𝑐_𝑖ℎ,𝑓_𝑖, (2) where 𝑓_𝑖=𝑦_𝑛−1𝑔_𝑖+ℎ𝑐_𝑖𝑦^′_𝑛−1+ℎ^2𝑖−1_𝑗=1𝛼_𝑖𝑗𝑓𝑡_𝑛−1+𝑐_𝑗ℎ,𝑓_𝑗,𝑖=1,…,𝑚. (3)

For the classical method 𝑔_𝑖=1, 𝑖=1(1)𝑚. In the present paper and based on the requirement of the development of the new method, values 𝑔_𝑖, 𝑖=1(1)𝑚 are variable and depend on 𝑧 (which is the product of the frequency 𝑤 and the step size ℎ). In Section 4 we will present a development of a four-stage modified Runge-Kutta-Nyström method of algebraic order four.

3. Phase Lag and Amplification Error Analysis of the Modified Runge-Kutta-Nyström Methods

Consider the problem 𝑑²𝑦(𝑡)𝑑𝑡²=(𝑖𝑤)²𝑦(𝑡)⟹𝑦^′′(𝑡)=−𝑤²𝑦(𝑡),𝑤∈𝑅, (4) with exact solution 𝑦(𝑡)_theoretical=𝑒^{±𝑖𝑤𝑡}. (5)

By applying a numerical method for the solution of (4), we obtain a numerical solution of the form 𝑦(𝑡)_approximate=𝑎(𝑤)𝑒^{±𝑖𝜃(𝑤)𝑡}. (6) Comparing the theoretical (5) and the numerical solution (6), we obtain the following two definitions (see [2] for details).

Definition 1. The difference 1−𝑎(𝑤) (7) is called dissipation error.

Definition 2. The difference 𝜃(𝑤)−𝑤 (8) is called dispersion error or phase lag.

Based on the above definitions, it is easy for one to see that we have two new errors, the dissipation and the dispersion errors. These errors can be expressed via Taylor series around an initial value 𝑤₀ of the frequency. Vanishing the appropriate derivatives of the Taylor series expansion (first, second, etc.) we optimize the specific error.

In order to develop the new method we use the test equation (4). By applying the modified RKN method (2) and (3) to the test equation (4) we obtain the numerical solution 𝑦_𝑛ℎ𝑦^′_𝑛=𝐷^𝑛𝑦₀ℎ𝑦^′₀,𝐴𝑧𝐷=²𝐵𝑧²𝐴̇𝑧²𝐵̇𝑧²,𝑧=𝑤ℎ, (9) where 𝐴, 𝐵, 𝐴̇, 𝐵̇ are polynomials in 𝑧², completely determined by the parameters of the method (2) and (3).

The eigenvalues of the amplification matrix 𝐷(𝑧²) are the roots of the characteristic equation 𝑟²−𝑡𝑟𝐷𝑧²𝑟+det𝐷𝑧²=0. (10)

In phase analysis one compares the phases of exp(𝑧) with the phases of the roots of the characteristic equation (10). The following definition is originally formulated by van der Houwen and Sommeijer [1].

Definition 3 (phase lag). Apply the RKN method (2) and (3) to the test equation (4). Then we define the phase lag Φ(𝑧)=𝑧−arccos(𝑡𝑟(𝐷)/2det(𝐷)). If Φ(𝑧)=𝑂(𝑧^𝑞+1), then the RKN method is said to have phase-lag order 𝑞. In addition, the quantity 𝑎(𝑧)=1−det(𝐷) is called amplification error. If 𝑎(𝑧)=𝑂(𝑧^𝑟+1), then the RKN method is said to be dissipative of order 𝑟.

From Definition 3, it follows that Φ(𝑧)𝑅𝑧=𝑧−arccos²2𝑄𝑧²,𝑎(𝑧)=1−𝑄𝑧², (11) where 𝑄(𝑧²) and 𝑅(𝑧²) are the determinant and the trace of the amplification matrix 𝐷(𝑧²), respectively. Also it is already known from [1], that we can write the functions 𝑅(𝑧²) and 𝑄(𝑧²) in the following form: 𝑅𝑧²=2−𝑟₁𝑧²+𝑟₂𝑧⁴−𝑟₃𝑧⁴+⋯+−𝑟_𝑖𝑧^2𝑖,𝑄𝑧²=1−𝑞₁𝑧²+𝑞₂𝑧⁴−𝑞₃𝑧⁴+⋯+−𝑞_𝑖𝑧^2𝑖. (12)

If at a point 𝑧, 𝑎(𝑧)=0, then the Runge-Kutta-Nyström method has zero dissipation at this point.

We can also put forward an alternative definition for the case of infinite order of phase lag.

Definition 4 (phase lag of order infinity). To obtain phase lag of order infinity the relation Φ(𝑧)=𝑧−arccos(𝑅(𝑧²)/2𝑄(𝑧²))=0 must hold.

From Definition 4 we have the following theorem.

Theorem 5. If one has phase lag of order infinity and at a point 𝑧, 𝑎(𝑧)=0, then 𝑅𝑧𝑧−arccos²2𝑄𝑧²=0,1−𝑄𝑧²⟹𝑅𝑧=0,²(𝑧),𝑄𝑧=2cos₂=1, (13) for more details see [7].

Lemma 6. For the construction of a method with nullification of phase lag, amplification error, and their derivatives, one must satisfy the conditions 𝑅(𝑧²)=2cos(𝑧), 𝑄(𝑧²)=1, 𝑅^′(𝑧²)=−2sin(𝑧), 𝑄^′(𝑧²)=0.

4. Derivation of the New Modified RKN Method

The new method that we are going to develop in this section, is a four-stage explicit Runge-Kutta-Nyström method with the FSAL technique (first stage as last), so the method actually uses three stages at each step for the function evaluations. From (2) and (3), the four-stage explicit modified RKN method can be written in the following form: 𝑦_𝑛=𝑦_𝑛−1𝑔₄+ℎ𝑦^′_𝑛−1+ℎ²𝑏₁𝑓₁+𝑏₂𝑓₂+𝑏₃𝑓₃+𝑏₄𝑓₄,𝑦^′_𝑛=𝑦^′_𝑛−1+ℎ𝑏^′₁𝑓₁+𝑏^′₂𝑓₂+𝑏^′₃𝑓₃+𝑏^′₄𝑓₄, (14) where 𝑓₁=𝑓𝑡_𝑛−1,𝑦_𝑛−1𝑔₁,𝑓₂=𝑓𝑡_𝑛−1+𝑐₂ℎ,𝑦_𝑛−1𝑔₂+𝑐₂ℎ𝑦^′_𝑛−1+ℎ²𝑎₂₁𝑓₁,𝑓₃=𝑓𝑡_𝑛−1+𝑐₃ℎ,𝑦_𝑛−1𝑔₃+𝑐₃ℎ𝑦^′_𝑛−1+ℎ²𝑎₃₁𝑓₁+𝑎₃₂𝑓₂𝑓₄=𝑓𝑡_𝑛−1+𝑐₄ℎ,𝑦_𝑛−1𝑔₄+𝑐₄ℎ𝑦^′_𝑛−1+ℎ²𝑎₄₁𝑓₁+𝑎₄₂𝑓₂+𝑎₄₃𝑓₃. (15)

By substituting the coefficients that have been used by the DEP algorithm in [15], (14) takes the following form: 𝑦_𝑛=𝑦_𝑛−1𝑔₄+ℎ𝑦^′_𝑛−1+ℎ²1𝑓14¹+8𝑓27²+25𝑓189³,𝑦^′_𝑛=𝑦^′_𝑛−11+ℎ𝑓14¹+32𝑓81²+250𝑓567³+5𝑓54⁴, (16) where 𝑓₁=𝑓𝑡_𝑛−1,𝑦_𝑛−1𝑔₁,𝑓₂=𝑓𝑡_𝑛−1+14ℎ,𝑦_𝑛−1𝑔₂+14ℎ𝑦^′_𝑛−1+1ℎ32²𝑓₁,𝑓₃=𝑓𝑡_𝑛−1+710ℎ,𝑦_𝑛−1𝑔₃+710ℎ𝑦^′_𝑛−1+ℎ²7𝑓1000¹+119𝑓500²,𝑓₄=𝑓𝑡_𝑛−1+ℎ,𝑦_𝑛−1𝑔₄+ℎ𝑦^′_𝑛−1+ℎ²1𝑓14¹+8𝑓27²+25𝑓189³. (17)

By applying numerical method (16) to the test equation (4), we compute the polynomials 𝐴, 𝐴̇, 𝐵, 𝐵̇ in terms of the modified Runge-Kutta-Nyström parameters. From these polynomials we obtain the expressions of 𝑅(𝑧²) and 𝑄(𝑧²). Then, according to Lemma 6 we solve the system of the four equations(𝑅(𝑧²)=2cos(𝑧), 𝑄(𝑧²)=1, 𝑅^′(𝑧²)=−2sin(𝑧), and 𝑄^′(𝑧²)=0) and thus we obtain the expressions of the coefficients 𝑔_𝑖, 𝑖=1(1)4, which are fully dependent on the product of the step length ℎ and the frequency 𝑤: 𝑔₁=5×6571982880𝑧⁷sin(𝑧)[phone]𝑧⁶−10298016𝑧⁸−5238722304𝑧⁴+61200𝑧¹⁰−1445𝑧¹²[phone]𝑧⁵sin(𝑧)[phone]𝑧²−7003998720sin(𝑧)𝑧³(𝑧)𝑧(𝑧)𝑧−87071293440+43535646720sin+7931520cos⁶[phone]cos(𝑧)𝑧⁴−29023764480cos(𝑧)𝑧²(𝑧)×𝑧+87071293440cos₄289𝑧⁸−12240𝑧⁶+203040𝑧⁴−1555200𝑧²+4665600⁻¹,𝑔₂5=−×31536−103538248704𝑧⁴−1175462461440+5554512𝑧¹²−120046752𝑧¹⁰−341029232640sin(𝑧)𝑧³−653034700800cos(𝑧)𝑧²−80053𝑧¹⁴[phone]𝑧⁶[phone]𝑧⁸[phone]𝑧⁹sin(𝑧)−3013231104𝑧⁷sin(𝑧)[phone]𝑧⁵sin(𝑧)[phone]cos(𝑧)𝑧⁸−9456238080cos(𝑧)𝑧⁶[phone]cos(𝑧)𝑧⁴[phone]sin(𝑧)𝑧+653034700800𝑧²(𝑧)×𝑧+1175462461440cos₄289𝑧⁸−12240𝑧⁶+203040𝑧⁴−1555200𝑧²+4665600⁻¹,𝑔₃1=−×6307200−9526307𝑧¹⁶[phone]𝑧¹⁴−4411486944𝑧¹²[phone]𝑧¹¹sin(𝑧)[phone]𝑧¹⁰cos(𝑧)−140829169536𝑧¹⁰−659696244480𝑧⁹sin(𝑧)[phone]𝑧⁸−2654019841536cos(𝑧)𝑧⁸[phone]𝑧⁷sin(𝑧)−36326761721856𝑧⁶[phone]cos(𝑧)𝑧⁶−97876195983360𝑧⁵sin(𝑧)[phone]𝑧⁴−260162575073280cos(𝑧)𝑧⁴[phone]sin(𝑧)𝑧³[phone]cos(𝑧)𝑧²−626390885007360𝑧²(𝑧)𝑧(𝑧)×𝑧−188073993830400sin−376147987660800cos+376147987660800₄289𝑧⁸−12240𝑧⁶+203040𝑧⁴−1555200𝑧²+4665600⁻¹,𝑔₄1=−×70956−80053𝑧¹²+3390480𝑧¹⁰+2298780𝑧⁸[phone]𝑧⁷sin(𝑧)−1744296768𝑧⁶[phone]cos(𝑧)𝑧⁶−5178046176𝑧⁵sin(𝑧)[phone]𝑧⁴−21763204416cos(𝑧)𝑧⁴[phone]sin(𝑧)𝑧³−131211601920𝑧²[phone]cos(𝑧)𝑧²(𝑧)(𝑧)×−362343559680sin𝑧+393634805760−724687119360cos289𝑧⁸−12240𝑧⁶+203040𝑧⁴−1555200𝑧²+4665600⁻¹. (18)

For small values of 𝑧 the following Taylor series expansions are used: 𝑔₁=1+86𝑧365²+45119𝑧1655640⁴+180461𝑧74503800⁶+3464911𝑧23602803840⁸+1124771𝑧1471133664000¹⁰,𝑔₂=1−387𝑧5840²+36731𝑧2207520⁴+1554263𝑧1192060800⁶+2028793𝑧17483558400⁸[phone]𝑧286380686592000¹⁰,𝑔₃=1+387𝑧18250²−25481237𝑧1103760000⁴+2106899𝑧1862595000⁶[phone]𝑧4370889600000⁸[phone]𝑧5369637873600000¹⁰,𝑔₄=1+52027𝑧21286800⁶+675821𝑧3576182400⁸+767177𝑧53642736000¹⁰[phone]𝑧50982056294400¹². (19)

5. Algebraic Order

In this section we study the algebraic order of the new modified RKN method. We require that 𝑦(𝑡_𝑛)=𝑦_𝑛 and 𝑦^′(𝑡_𝑛)=𝑦^′_𝑛. By using the chain rule and 𝑦^′′(𝑡)=𝑓(𝑡,𝑦), the Taylor expansion of the exact solution 𝑦(𝑡_𝑛+ℎ) is 𝑦𝑡_𝑛+1(𝑡)+𝑑=𝑦𝑦(𝑡)1𝑑𝑡ℎ+2𝐹ℎ²+16𝐹_𝑡+16𝐹_𝑦𝑑𝑦(𝑡)ℎ𝑑𝑥³+1𝐹24^𝑡,𝑡+1𝑑12𝑦(𝑡)𝐹𝑑𝑡^𝑡,𝑦+1𝑑24𝑦(𝑡)𝑑𝑡²𝐹_𝑦,𝑦+1𝐹24^𝑦𝐹ℎ⁴+𝑂ℎ⁵. (20)

By expanding in Taylor series the corresponding numerical solution 𝑦_𝑛+1 of a classical explicit Runge-Kutta-Nyström method, where 𝑐₁=0, we have 𝑦𝑡_𝑛(𝑡)+𝑑+ℎ=𝑦𝑦(𝑡)ℎ+𝑏𝑑𝑡₁𝐹+𝑏₂𝐹+𝑏₃𝐹+𝑏₄𝐹ℎ²+𝑏₂𝐹_𝑡𝑐₂+𝐹_𝑦𝑐₂𝑑𝑦(𝑡)𝑑𝑡+𝑏₃𝐹_𝑡𝑐₃+𝐹_𝑦𝑐₃𝑑𝑦(𝑡)𝑑𝑡+𝑏₄𝐹_𝑡𝑐₄+𝐹_𝑦𝑐₄𝑑𝑦(𝑡)ℎ𝑑𝑡³+𝑏₂12𝐹_𝑡,𝑡𝑐²₂+𝐹_𝑡,𝑦𝑐²₂𝑑𝑦(𝑡)+1𝑑𝑡2𝑐²₂𝑑𝑦(𝑡𝑑𝑡)²𝐹_𝑦,𝑦+𝐹_𝑦𝑎₂₁𝐹+𝑏₃𝐹_𝑦𝑎₃₁𝐹+𝐹_𝑦𝑎₃₂𝐹+𝐹_𝑡,𝑦𝑐²₃𝑑𝑦(𝑡)+1𝑑𝑡2𝑐²₃𝑑𝑦(𝑡)𝑑𝑡²𝐹_𝑦,𝑦+12𝐹_𝑡,𝑡𝑐²₃+𝑏₄12𝐹_𝑡,𝑡𝑐²₄+𝐹_𝑦𝑎₄₂𝐹+𝐹_𝑦𝑎₄₃𝐹+12𝑐²₄𝑑𝑦(𝑡)𝑑𝑡²𝐹_𝑦,𝑦+𝐹_𝑦𝑎₄₁𝐹+𝐹_𝑡,𝑦𝑐²₄𝑑𝑦(𝑡)ℎ𝑑𝑡⁴+𝑂ℎ⁵. (21)

By equating (20) and (21) with respect to ℎ, we obtain the algebraic order conditions for a fourth algebraic order Runge-Kutta-Nyström method. So the equations obtained for 𝑦 are.

Order 2: ⁴_𝑖=1𝑏_𝑖=12, (22)

order 3: ⁴_𝑖=1𝑏_𝑖𝑐_𝑖=16, (23)

order 4: ⁴_𝑖=1𝑏_𝑖𝑐²_𝑖=1,12⁴_𝑖=1𝑏_𝑖𝑎_𝑖𝑗=1.24 (24)

By following the same procedure for 𝑦^′, we obtain the corresponding algebraic order conditions, which are.

Order 1: ⁴_𝑖=1𝑏^′_𝑖=1, (25)

order 2: ⁴_𝑖=1𝑏^′_𝑖𝑐_𝑖=12, (26)

order 3: ⁴_𝑖=1𝑏^′_𝑖𝑐²_𝑖=13,⁴_𝑖=1𝑏^′_𝑖𝑎_𝑖𝑗=16, (27)

order 4: ⁴_𝑖=1𝑏^′_𝑖𝑐³_𝑖=14,⁴_𝑖=1𝑏^′_𝑖𝑐_𝑖𝑎_𝑖𝑗=18,⁴_𝑖=1𝑏^′_𝑖𝑎_𝑖𝑗𝑐_𝑗1.=.24 (28)

For the classical RKN method [15] that we use in the present paper, the above conditions are satisfied, as this method is of fourth algebraic order.

In order to verify the algebraic order of the the new modified explicit Runge-Kutta-Nyström method, first we will produce the algebraic order conditions. To do that, we apply the above methodology for the new method. At first, we want to extract the conditions for 𝑦 and so we expand the modified method (14) in Taylor series. The result from the computations is given below: 𝑦𝑡_𝑛(𝑡)𝑔+ℎ=𝑦₄+𝑑𝑦(𝑡)ℎ+𝑏𝑑𝑡₁𝐺1+𝑏₂𝐺2+𝑏₃𝐺3+𝑏₄𝐺4ℎ²+𝑏₂𝐺2_𝑡𝑐₂+𝐺2_𝑦𝑐₂𝑑𝑦(𝑡)𝑑𝑡+𝑏₃𝐺3_𝑡𝑐₃+𝐺3_𝑦𝑐₃𝑑𝑦(𝑡)𝑑𝑡+𝑏₄𝐺4_𝑡𝑐₄+𝐺4_𝑦𝑐₄𝑑𝑦(𝑡ℎ𝑑𝑡)³+𝑏₂12𝐺2_𝑡,𝑡𝑐²₂+𝐺2_𝑡,𝑦𝑐²₂𝑑𝑦(𝑡)+1𝑑𝑡2𝑐²₂𝑑𝑦(𝑡)𝑑𝑡²𝐺2_𝑦,𝑦+𝐺2_𝑦𝑎₂₁𝐺1+𝑏₃𝐺3_𝑦𝑎₃₁𝐺1+𝐺3_𝑦𝑎₃₂𝐺2+𝐺3_𝑡,𝑦𝑐²₃𝑑𝑦(𝑡)+1𝑑𝑡2𝑐²₃𝑑𝑦(𝑡)𝑑𝑡²𝐺3_𝑦,𝑦+12𝐺3_𝑡,𝑡𝑐²₃+𝑏₄𝐺4_𝑦𝑎₄₂𝐺2+𝐺4_𝑦𝑎₄₃𝐺3+𝐺4_𝑡,𝑦𝑐²₄𝑑𝑦(𝑡)𝑑𝑡+𝐺4_𝑦𝑎₄₁+1𝐺12𝑐²₄𝑑𝑦(𝑡)𝑑𝑡²𝐺4_𝑦,𝑦+12𝐺4_𝑡,𝑡𝑐²₄ℎ⁴+𝑂ℎ⁵. (29)

By equating (20) and (29) with respect to ℎ, we obtain the algebraic order conditions of a fourth algebraic order modified Runge-Kutta-Nyström method of the form (14). The equations obtained for 𝑦 are.

Order 2: ⁴_𝑖=1𝑏_𝑖1𝐺𝑖=2𝐹, (30)

order 3: ⁴_𝑖=1𝑏_𝑖𝑐_𝑖𝐺𝑖_𝑡+𝐺𝑖_𝑦𝑑𝑦(𝑡)=1𝑑𝑡6𝐹_𝑡+𝐹_𝑦𝑑𝑦(𝑡),𝑑𝑡 (31)

order 4: ⁴_𝑖=1𝑏_𝑖𝑐²_𝑖𝐺𝑖_𝑡,𝑡+𝐺𝑖_𝑡,𝑦+𝑑𝑦(𝑡)𝑑𝑡²𝐺𝑖_𝑦,𝑦=1𝐹12^𝑡,𝑡+𝐹_𝑡,𝑦+𝑑𝑦(𝑡)𝑑𝑡²𝐹_𝑦,𝑦,_4𝑖=1𝑏_𝑖𝑎_𝑖𝑗𝐺𝑖_𝑦1𝐺𝑗=𝐹24^𝑦𝐹. (32)

By following the same procedure for the approximate solution of 𝑦^′, we obtain the algebraic order conditions for 𝑦^′, which are.

Order 1: ⁴_𝑖=1𝑏^′_𝑖𝐺𝑖=𝐹, (33)

order 2: ⁴_𝑖=1𝑏^′_𝑖𝑐_𝑖𝐺𝑖_𝑡+𝐺𝑖_𝑦𝑑𝑦(𝑡)=1𝑑𝑡2𝐹_𝑡+𝐹_𝑦𝑑𝑦(𝑡),𝑑𝑡 (34)

order 3: ⁴_𝑖=1𝑏^′_𝑖𝑐²_𝑖𝐺𝑖_𝑡,𝑡+2𝐺𝑖_𝑡,𝑦𝑑𝑦(𝑡)+𝑑𝑑𝑡𝑦(𝑡)𝑑𝑡²𝐺𝑖_𝑦,𝑦=13𝐹_𝑡,𝑡+2𝐹_𝑡,𝑦𝑑𝑦(𝑡)+𝑑𝑑𝑡𝑦(𝑡)𝑑𝑡²𝐹_𝑦,𝑦,_4𝑖=1𝑏^′_𝑖𝑎_𝑖𝑗𝐺𝑖_𝑦1𝐺𝑗=6𝐹_𝑦𝐹, (35)

order 4: ⁴_𝑖=1𝑏^′_𝑖𝑐³_𝑖𝐺𝑖_{𝑡,𝑡,𝑡}+𝐺𝑖_{𝑡,𝑡,𝑦}𝑑𝑦(𝑡)𝑑𝑡+𝐺𝑖_{𝑡,𝑦,𝑦}𝑑𝑦(𝑡)𝑑𝑡²+𝐺𝑖_{𝑦,𝑦,𝑦}𝑑𝑦(𝑡)𝑑𝑡³=14𝐹_{𝑡,𝑡,𝑡}+𝐹_{𝑡,𝑡,𝑦}𝑑𝑦(𝑡)𝑑𝑡+𝐹_{𝑡,𝑦,𝑦}𝑑𝑦(𝑡)𝑑𝑡²+𝐹_{𝑦,𝑦,𝑦}𝑑𝑦(𝑡)𝑑𝑡³_4𝑖=1,𝑏^′_𝑖𝑐_𝑖𝑎_𝑖𝑗𝐺𝑗𝐺𝑖_𝑡,𝑦+𝑑𝑦(𝑡)𝑑𝑡𝐺𝑖_𝑦,𝑦=18𝐹𝐹_𝑡,𝑦+𝐹_𝑦,𝑦,_4𝑖=1𝑏^′_𝑖𝑎_𝑖𝑗𝑐_𝑗𝐺𝑖_𝑦𝐺𝑗_𝑡+𝐺𝑖_𝑦𝐺𝑗_𝑦𝑑𝑦(𝑡)=1𝑑𝑡𝐹24^𝑦𝐹_𝑡+𝐹²_𝑦𝑑𝑦(𝑡).𝑑𝑡 (36)

As it is proved, in order to construct an explicit modified Runge-Kutta-Nyström method of fourth algebraic order, (30)–(36) must be satisfied.

For the proposed modified RKN method (16) developed in Section 4, we observe that the Taylor expansions for 𝑔_𝑖 parameters that were obtained are dependent on 𝑧. What is important for the new modified method is to maintain the fourth algebraic order as the step size ℎ takes small values. As already mentioned, 𝑧 is the product of the frequency 𝑤 and the step size ℎ; thus for 𝑧→0 (30)–(36) should verify the fourth algebraic order of the new modified method.

For 𝑧→0 it is easy to prove that lim_𝑧→0𝑔_𝑖=1, 𝑖=1(1)4. From the last relation, it is extracted that lim_𝑧→0𝐺𝑖=𝐹, 𝑖=1(1)4, and thus (30)–(36) reduced to (22)–(28) which are satisfied for the constant coefficients of DEP algorithm [15].

So it is proved that the new RKN method is of fourth algebraic order. Moreover the local truncation error in 𝑦 is given from the following equation: LTE=𝑦_𝑛+1−𝑦𝑡_𝑛+ℎ, (37) accordingly the local truncation error in its first derivative (𝑦^′) is given from the following equation: LTE_der=𝑦^′_𝑛+1−𝑦^′𝑡_𝑛+ℎ. (38)

At this point, we have already compute the Taylor expansions of(a)the exact solution 𝑦(𝑡_𝑛+ℎ) and the corresponding numerical solution 𝑦_𝑛+1,(b)the first derivative 𝑦(𝑡_𝑛+ℎ) of the exact solution and the first derivative 𝑦^′_𝑛+1 of the corresponding numerical solution.

The LTE verifies the fourth algebraic order of the new modified method. From the above procedure the local truncation error in 𝑦 is 1LTE=𝑓2160^𝑦𝑓_𝑡+𝑓_𝑦𝑑𝑦(𝑡)ℎ𝑑𝑡⁵+𝑂ℎ⁶. (39)

Respectively, the local truncation error in 𝑦^′ is LTE_der=12592031𝑓²_𝑦𝑓−12𝑓_𝑡,𝑦𝑓_𝑡+31𝑓_𝑦𝑓_𝑡,𝑡𝑑+19𝑦(𝑡𝑑𝑡)²𝑓_𝑦,𝑦𝑓_𝑦+50𝑓_𝑦𝑑𝑦(𝑡)𝑓𝑑𝑡^𝑡,𝑦𝑑−12𝑦(𝑡)𝑓𝑑𝑡^𝑦,𝑦𝑓_𝑡ℎ⁵+𝑂ℎ⁶. (40)

From (30)–(36) and for 𝑧→0 we prove that the new method is of fourth algebraic order. Moreover the forms (39) and (40) show that the new modified RKN method is of fourth algebraic order, because all the coefficients up to the power of ℎ⁴ vanish.

6. Numerical Results

6.1. Schrödinger Equation

The one-dimensional or radial Schrödinger equation has the form 𝑦^′′(𝑥)+𝑙()𝐸−𝑙+1𝑥²−𝑉(𝑥)𝑦(𝑥)=0,where0≤𝑥<∞. (41)

We call the term 𝑙(𝑙+1)/𝑥² the centrifugal potential, and the function 𝑉(𝑥) the electric potential. In (41), 𝐸 is a real number denoting the energy, and 𝑙 is a quantum number. The function 𝑊(𝑥)=𝑙(𝑙+1)/𝑥²+𝑉(𝑥) denotes the effective potential, where lim_𝑥→∞𝑉(𝑥)=0 and so lim_𝑥→∞𝑊(𝑥)=0. The boundary condition 𝑦(0)=0, together with a second boundary condition, for large values of 𝑥, is determined by the physical considerations.

6.1.1. Woods-Saxon Potential with Positive Energies

For the purpose of our numerical illustration we will take the domain of integration as 𝑥∈[0,15], using the Woods-Saxon potential: 𝑉(𝑥)=𝑢₀+𝑢1+𝑞¹𝑞1+𝑞²,𝑞=exp𝑥−𝑥₀𝛼,where𝑢₀=−50,𝛼=0.6,𝑥₀=7,𝑢₁𝑢=−⁰𝛼. (42)

In the case of positive energies (𝐸=𝑘²), the potential (𝑉(𝑥)) dies away faster than the centrifugal potential (𝑙(𝑙+1)/𝑥²), so for a large number for 𝑥, Schrödinger equation effectively reduces to 𝑦^′′(𝑥)+𝑘²−𝑙(𝑙+1)𝑥²𝑦(𝑥)=0. (43)

(43) has two linearly independent solutions, 𝑘𝑥𝑗_𝑙(𝑘𝑥) and 𝑘𝑥𝑛_𝑙(𝑘𝑥), where 𝑗_𝑙 and 𝑛_𝑙 are the spherical Bessel and Neumann functions, respectively. When 𝑥→∞, the solution of (41) takes the following asymptotic form: 𝑦(𝑥)≃𝐴𝑘𝑥𝑗_𝑙(𝑘𝑥)−𝐵𝑘𝑥𝑛_𝑙(𝑘𝑥)≃𝐷sin𝑘𝑥−𝑙𝜋2+tan𝛿_𝑙cos𝑘𝑥−𝑙𝜋2, (44) where 𝛿_𝑙 is the scattering phase shift that may be calculated from the below formula: tan𝛿_𝑙=𝑦𝑥_𝑖𝑆𝑥_𝑖+1−𝑦𝑥_𝑖+1𝑆𝑥_𝑖𝑦𝑥_𝑖+1𝐶𝑥_𝑖−𝑦𝑥_𝑖𝐶𝑥_𝑖+1, (45) where 𝑆(𝑥)=𝑘𝑥𝑗_𝑙(𝑘𝑥) and 𝐶(𝑥)=𝑘𝑥𝑛_𝑙(𝑘𝑥) and 𝑥_𝑖<𝑥_𝑖+1 both exist in the asymptotic region.

For positive energies and for 𝑙=0, we calculate the phase shift (𝛿_𝑙) and then we compare it with the accurate value which is 𝜋/2. The boundary conditions for this eigenvalue problem are 𝑦(0)=0 and 𝑦(𝑥)=cos(√𝐸𝑥) for large 𝑥.

We use the following eigenenergies: 𝐸₁𝐸=53.588872,₂𝐸=163.215341,₃𝐸=341.495874,₄=989.701916. (46)

6.1.2. Woods-Saxon Potential with Negative Energies

In the case of negative energies (𝐸<0), we consider the eigenvalue problem with boundary conditions: 𝑦(0)𝑦(𝑥)−√=0,=exp−𝐸𝑥forlarge𝑥. (47)

In order to solve this problem numerically, by a chosen eigenvalue, we integrate forward from the point 𝑥=0 and backward from the point 𝑥=15, matching up the solution at some internal point in the range of integration.

For this problem we use the following eigenenergies: 𝐸₁=−49.457788728,𝐸₂𝐸=−38.122785096,₃=−22.588602257,𝐸₄=−3.908232481, (48)

6.1.3. Lennard-Jones Potential

In order to investigate the case of 𝑙≠0, we examine the Lennard-Jones potential which is given by the equation 𝑉(𝑥)1=𝑚𝑥¹²−1𝑥⁶, (49) where 𝑚=500, for 𝑙={0,1,2,…,10},𝐸={25,100}, and step length of integration ℎ=0.1. We compare the phase shifts with the values found in [16] and present the decimal digits that the approximate solutions that agree with the reference solution.

6.2. Numerical Results

In this section four Runge-Kutta-Nyström methods are compared (including the new method). These methods have four algebraic orders with four stages; also all of them are using the FSAL properties. The methods used in the comparison have been denoted by(i)MRKNDPAF4: the new fourth-order RKN method with four stages (three effective stages with FSAL property), derived in Section 3,(ii)RKNPAF4: the fourth-order RKN method with four stages (three effective stages with FSAL property), phase lag, and amplification error of order infinity of Papadopoulos et al. [7],(iii)DEPRKN4: the high-order method of pair RKN 4(3)4 method of Dormand et al. [15],(iv)EFRKN4: the fourth-order exponential fitted RKN method with four stages (three effective stages with FSAL property) of Franco [4].

One way to measure the efficiency of the method is to compute the accuracy in the decimal digits, that is, −log₁₀(𝑒𝑟𝑟𝑜𝑟𝑎𝑡𝑡ℎ𝑒𝑒𝑛𝑑𝑝𝑜𝑖𝑛𝑡), when comparing the phase shift to the actual value 𝜋/2 versus the computational effort measured by the log₁₀(𝑛𝑢𝑚𝑏𝑒𝑟𝑜𝑓𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛𝑒𝑣𝑎𝑙𝑢𝑎𝑡𝑖𝑜𝑛𝑠𝑟𝑒𝑞𝑢𝑖𝑟𝑒𝑑). The range of integration steps that have been used is ℎ=1/2^𝑁,𝑁=3(1)8 for the Woods-Saxon potential with positive energies and ℎ=1/2^𝑁,𝑁=2(1)6 for the Woods-Saxon potential with negative energies.

In the case of Lennard-Jones potential, we compare the phase shifts with the values found in [16] and present the decimal digits that the approximate solutions that agree with the reference solution. The results are given in Tables 1 and 2. The step length of integration is ℎ=0.1.

Table 1 

Lennard-Jones potential with 𝐸=25.

Table 2 

Lennard-Jones potential with 𝐸=100.

The frequency is given by the suggestion of Ixaru and Rizea [17]: √𝑤=[],√𝐸+50,𝑥∈0,6.5𝐸,𝑥∈[6.5,15]. (50)

The numerical results were obtained by using the high-level language MATLAB. In Figures 1, 2, 3, 4, 5, 6, 7 and 8 we display the efficiency curves, that is, the accuracy versus the computational cost measured by the number of function evaluations required by each method.

Figure 1 

Efficiency for the Schrödinger equation using E = 53.588872.

Figure 2 

Efficiency for the Schrödinger equation using E = 163.215341.

Figure 3 

Efficiency for the Schrödinger equation using E = 341.495874.

Figure 4 

Efficiency for the Schrödinger equation using E = 989.701916.

Figure 5 

Efficiency for the Schrödinger equation using E = − 49.457788728.

Figure 6 

Efficiency for the Schrödinger equation using E = − 38.122785096.

Figure 7 

Efficiency for the Schrödinger equation using E = − 22.588602257.

Figure 8 

Efficiency for the Schrödinger equation using E = − 3.908232481.

Numerical results indicate that the new method derived in Section 4 is very efficient in solving numerically the Schrödinger equation and more accurate than the other methods at all the eigenvalues.

More specifically in the case of Woods-Saxon potential with positive energies, the new method remains more accurate than the RKNPAF4 and EFRKN4 methods up to two decimals at all eigenvalues. Also our method is more accurate than the classical DEPRKN4 method, by two decimals for the eigenvalue 𝐸=53.588872, by three decimals for the eigenvalue 𝐸=163.215341, and by four decimals for the next two eigenvalues 𝐸={341.495874,989.701916}.

In the case of Woods-Saxon potential with negative energies, the new method has almost the same accuracy with the RKNPAF4 and EFRKN4 methods, but it remains (even a little) more accurate at the majority of the negative eigenvalues. In comparison with the DEPRKN4 method, the new method is much more accurate and specifically by one digit for the eigenvalue 𝐸=−49.457788728 and by two digits for the rest of the negative eigenvalues.

At last, for the Lennard-Jones potential the new method is more accurate than all the other methods for eigenvalues 𝐸={25,100} and 𝑙={0,1,2,…,10}.

7. Conclusions

The modified RKN method, developed in this paper, is much more efficient than the classical one, in any case. The new method remained more efficient for all the eigenvalues and in some cases was more accurate than the other methods up to two decimals. Moreover we observe that the accuracy difference between the new method and the other methods increased as the eigenvalue increased (for details about the original proof see [17]).

Disclosure

T. E. Simos is an active member of the European Academy of Sciences and the European Academy of Sciences and Arts and a corresponding member of the European Academy of Arts, Sciences and Humanities.