An Optimized Runge-Kutta Method for the Numerical Solution of the Radial Schrödinger Equation

An optimized explicit modified Runge-Kutta RK method for the numerical integration of the radial Schrödinger equation is presented in this paper. This method has frequency-depending coefficients with vanishing dispersion, dissipation, and the first derivative of dispersion. Stability and phase analysis of the new method are examined. The numerical results in the integration of the radial Schrödinger equation with the Woods-Saxon potential are reported to show the high efficiency of the new method.


Introduction
In this paper, we are concerned with the numerical integration of the one-dimensional Schrödinger equation of the form where the real number E is the energy and the function v x is the effective potential satisfying v x → 0 as x → ∞. Two boundary conditions are associated with this equation: one is y 0 0, and the other imposed at large x is determined by physical considerations. The form of this second boundary condition depends crucially on the sign of the energy E. Such problems are frequently encountered in a variety of scientific fields and engineering applications 1-9 . Concerning the oscillatory character of the solution to the Schrödinger equation 1.1 , there have appeared a lot of numerical integrators of adapted type, a pronounced class of which is based on important properties such as the phase lag and the amplification see 10-18 . These are actually two different kinds of truncation errors.
The first is the angle between the analytical solution and the numerical solution, and the second is the distance from a standard cyclic solution. If a good frequency is estimated in advance, then it is a good choice to construct numerical methods with zero dispersion or/and zero dissipation. These techniques are called phase fitted or/and zero dissipation. Related work can be founded in [19][20][21] . For Runge-Kutta methods, Simos and Aguiar 18 constructed a modified Runge-Kutta method for the numerical integration of the Schrödinger equation by phase fitting based on the fifth-order RK method. Recently, Van de Vyver 16 gave an embedded pair of modified RK methods by nullifying the phase-lags of the fifthorder method and the fourth-order method. And in 22 , Tsitouras and Simos constructed phase-fitted and zero dissipation fifth-order Runge-Kutta method for the numerical solution of oscillatory problems.
In this paper, inspired by the ideas in 23-28 , we construct a new kind of modified fifth-order Runge-Kutta method by nullifying the dispersion, the dissipation, and the first derivative of the dispersion. In Section 2, the preliminaries of the phase properties of explicit modified Runge-Kutta methods are introduced. In Section 3, the coefficients of a new kind of optimized modified RK method are obtained. Section 4 examines the stability and phase properties of the new method. In Section 5, the numerical experiments are reported.

Preliminaries
We begin by considering the numerical integration of the initial value problem IVP of firstorder differential equations in the following form: whose solution shares an oscillatory character. We follow the convention to assume that the frequency is known to be ω in advance or can be accurately estimated. An s-stage-modified explicit Runge-Kutta RK method has the following scheme: the problem. We assume that lim ν → 0 γ i ν 1, i 1, . . . , s so that as ν → 0, the modified RK method 2.2 reduces to a traditional RK method. An alternative approach adopted by, for example, exponential/trigonometric fitting techniques, is to let some of the coefficients a ij , c i , b i , i 1, . . . , s be functions of ν hω see 16, 18, 29 . Applying the modified RK method 2.2 to the test equation as follows: y iωy, ω > 0 2.4 yields y n 1 R iν y n , ν ωh.

2.5
A comparison of the numerical solution with the exact solution leads to the notions of phaselag and dissipation error defined as follows.
Definition 2.1. The following two quantities are called the phase lag or dispersion and the amplification factor error or dissipation error , respectively: The method is said to be dispersive of order q and dissipative of order p if If P ν 0 and D ν 0, the method is called phase fitted (zero dispersive) and amplificationfitted (zero dissipative), respectively.
For modified RK method 2.2 , we have are polynomials in ν 2 , which are completely defined by the Runge-Kutta coefficients c, A, γ, and b. Therefore, we have Based on the fifth algebraic order six-stage Dormand and Prince Runge-Kutta method, Simos and Aguiar 18 obtained an explicit modified RK method with one parameter γ 2 taking the orthers γ 1 γ i 1 for i 3, . . . , 6 determined by nullifying the quantity 4 Mathematical Problems in Engineering tan ν − ν V ν 2 / U ν 2 . In 22 , Tsitouras and Simos presented an optimized Runge-Kutta method by nullifying the dispersion and the dissipation. In this paper, we construct a new optimized Runge-Kutta method by nullifying the dispersion, the dissipation, and the first derivative of the dispersion.

Construction of the New Method
In this section, we are concerned with the following Runge-Kutta method given by the Butcher tableau as follows: If we choose γ 2 γ 3 γ 4 1, the classical Runge-Kutta method with order fifth derived by Dormand and Prince 30 is recovered. In order to construct the new embedded RK pair, we set γ 2 , γ 3 , γ 4 free and keep the rest of the coefficients. Motivated by the ideas in 23-28 , we obtain the dispersion, the dissipation, and the first derivative of the dispersion of this method, which depend on ν, γ 2 , γ 3 , γ 4 as follows: 3.4 In order to check the algebraic order of the newly obtained modified RK method, we note that the order conditions listed in 31 for traditional RK methods are not sufficient for the modified RK method 2.2 . Writing

Analysis of Stability and Phase Properties
In this section, we are interested in the stability and phase properties of the new method. Lambert  is called the region of imaginary stability. And any closed curve defined by |M iθ, ν | 1 is a stability boundary of the method.
In Figure 1 we plot the region of imaginary stability for the method MODRK5PLDPLAM. are called the phase lag dispersion and amplification factor error dissipation , respectively. If the method is said to be of phase-lag order q and dissipation order p, respectively, where the c φ and c d are called the phase-lag constant and dissipation constant, respectively.
We note that, by definition, when ν θ ω λ , it must be true that P θ, ν 0 and D θ, ν 0. In general, ω / λ since the fitting frequency ω is just an estimate of the true frequency. Therefore the order of P θ, ν 0 and D θ, ν 0 in Definition 4.2 measure to what extent a modified RK method is accurate in phase and dissipation. Denoting the ratio r ν/θ ω/λ, we obtain the following expressions for the phase lag and the dissipation error of the new method MODRK5PLDPLAM:

4.7
Thus, the method MODRK5PLDPLAM has a phase lag of order six and a dissipation of order five.

Numerical Experiments
In this section, we test the numerical performance of the new fifth-order method in the integration of the radial Schrödinger equation with the well-known Woods-Saxon potential,  In the numerical experiment we consider the resonance problem E > 0 , the numerical results E calculated are compared with the analytical solution E analytical of the Woods-Saxon potential, rounded to six decimal places. In Figures 2, 3

Conclusions and Discussions
Based on the classical fifth RK method of Dormand and Prince 30 , a new optimized explicit modified RK method with modifying parameters is obtained by nullifying the dispersion, the dissipation, and the first derivative of the dispersion. The numerical results stated in Figures  2-5 illustrate the higher efficiency of the new method compared to some highly efficient methods in the recent literature 16-18, 22, 35 .