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

Abstract

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.

1. Introduction

In this paper, we are concerned with the numerical integration of the one-dimensional Schrödinger equation of the form 𝑦^′′(𝑥)=(𝑣(𝑥)−𝐸)𝑦(𝑥), (1.1) where the real number 𝐸 is the energy and the function 𝑣(𝑥) is the effective potential satisfying 𝑣(𝑥)→0 as 𝑥→∞. Two boundary conditions are associated with this equation: one is 𝑦(0)=0, and the other imposed at large 𝑥 is determined by physical considerations. The form of this second boundary condition depends crucially on the sign of the energy 𝐸. 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–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 fifth-order 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.

2. Preliminaries

We begin by considering the numerical integration of the initial value problem (IVP) of first-order differential equations in the following form: 𝑦^′(𝑥)=𝑓𝑥,𝑦,𝑦𝑥₀=𝑦₀, (2.1) 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 𝑠-stage-modified explicit Runge-Kutta (RK) method has the following scheme: 𝑌_𝑖=𝛾_𝑖𝑦_𝑛+ℎ^𝑖−1_𝑗=1𝑎_𝑖𝑗𝑓𝑥_𝑛+𝑐_𝑗ℎ,𝑌_𝑗𝑦,𝑖=1,…,𝑠,_𝑛+1=𝑦_𝑛+ℎ^𝑠_𝑖=1𝑏_𝑖𝑓𝑥_𝑛+𝑐_𝑖ℎ,𝑌_𝑖, (2.2) where the coefficients 𝑎_𝑖𝑗,  𝑐_𝑖,  𝑏_𝑖,  𝑖=1,…,𝑠 are constants, ℎ is the step size, and the parameters 𝛾_𝑖,𝑖=1,…,𝑠 are even functions of 𝜈=ℎ𝜔. It is convenient to express the modified RK method (2.2) by the Butcher tableau as follows: 𝑐₁𝛾₁𝑐₂𝛾₂𝑎₂₁⋮⋮𝑐⋮⋱_𝑠𝛾_𝑠𝑎_𝑠1⋯𝑎_𝑠𝑠−1𝑏₁⋯𝑏_𝑠−1𝑏_𝑠 (2.3) or simply by (𝑐,𝛾,𝐴,𝑏). The extra-frequency-depending parameters 𝛾_𝑖(𝜈),  𝜈=ℎ𝜔,  𝑖=1,…,𝑠 are introduced to tune the traditional RK method to the special oscillatory structure of the problem. We assume that lim_𝜈→0𝛾_𝑖(𝜈)=1,  𝑖=1,…,𝑠 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 𝑎_𝑖𝑗,  𝑐_𝑖,  𝑏_𝑖,  𝑖=1,…,𝑠 be functions of 𝜈=ℎ𝜔 (see [16, 18, 29]).

Applying the modified RK method (2.2) to the test equation as follows: 𝑦^′=𝑖𝜔𝑦,𝜔>0 (2.4) yields 𝑦_𝑛+1=𝑅(𝑖𝜈)𝑦_𝑛,𝜈=𝜔ℎ. (2.5) A comparison of the numerical solution with the exact solution leads to the notions of phase-lag 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: 𝑃(𝜈)=𝜈−arg(𝑅(𝑖𝜈)),𝐷(𝜈)=1−|𝑅(𝑖𝜈)|. (2.6) The method is said to be dispersive of order 𝑞 and dissipative of order 𝑝 if 𝑃(𝜈)=𝒪𝜈^𝑞+1,𝐷(𝜈)=𝒪𝜈^𝑝+1. (2.7) If 𝑃(𝜈)=0 and 𝐷(𝜈)=0, the method is called phase fitted (zero dispersive) and amplification-fitted (zero dissipative), respectively.
For modified RK method (2.2), we have 𝑅(𝑖𝜈)=𝑈𝜈²+𝑖𝜈𝑉𝜈², (2.8) where 𝑈𝜈²=1−𝑡₂𝜈²+𝑡₄𝜈⁴+⋯,𝑉𝜈²=1−𝑡₃𝜈²+𝑡₅𝜈⁴+⋯ (2.9) are polynomials in 𝜈², which are completely defined by the Runge-Kutta coefficients 𝑐,  𝐴,  𝛾, and 𝑏. Therefore, we have 𝑃(𝜈)𝜈𝑉𝜈=𝜈−arctan²𝑈𝜈²,𝐷(𝜈)=1−𝑈𝜈²²+𝜈²𝑉𝜈²². (2.10)

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 𝛾₂ (taking the orthers 𝛾₁=𝛾_𝑖=1 for 𝑖=3,…,6) determined by nullifying the quantity tan(𝜈)−𝜈((𝑉(𝜈²))/(𝑈(𝜈²))). 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.

3. Construction of the New Method

In this section, we are concerned with the following Runge-Kutta method given by the Butcher tableau as follows: 01015𝛾₂150310𝛾₃394004405𝛾₄44−455615329089119372−656125360218764448−65612120172919017−31683553346732524749−17651030186563503845001113125192−2187678411084 (3.1) If we choose 𝛾₂=𝛾₃=𝛾₄=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 𝛾₂,𝛾₃,𝛾₄ 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 𝜈,𝛾₂,𝛾₃,𝛾₄ as follows: 𝑃(𝜈)(𝜈)−𝑀=tan𝑁,𝑑(𝜈)=1−𝑀²+𝑁²,der⋅𝑃(𝜈)=sec²(𝜈)−𝑀^′𝑁−𝑀𝑁^′𝑁², (3.2) where 𝑀=15𝜈474651−688905𝛾₅−96460𝜈²+17808𝛾₂𝜈⁴−11130𝛾₄2𝜈²−125−640𝛾₃371𝜈²−1500,𝑁=−72251855+6400𝛾₃+2650𝛾₄−729𝛾₅𝜈²−4579200+212003𝛾₂−8𝛾₃−4𝜈⁴+7632𝜈⁶. (3.3) Now, solving (3.2), we get 𝛾 values in terms of 𝜈. Instead of giving the very complicated expressions for 𝛾_𝑖, for the purpose of practical computation, we present their Taylor expansions as follows: 𝛾₂=1−24179𝜈²−698950491109813𝜈⁴−2791606300003747105974663311𝜈⁶−82786341049335000005235134512534020593713𝜈⁸𝛾50148377999575005150000000+⋯,₃=1+901𝜈⁴−998500146822137𝜈⁶+8973020250000390542419221781𝜈⁸−394220671663500000003826206167449731276763𝜈¹⁰𝛾1074608099990892967500000000+⋯,₄=1−1088𝜈⁴+17473751151652176𝜈⁶+39256963593755225984025866𝜈⁸+23954381090664062552283859929609197732𝜈¹⁰9794605078041993193359375+⋯. (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 𝛾_𝑖=1+𝛾_𝑖⁽²⁾𝜈²+𝛾_𝑖⁽⁴⁾𝜈⁴+𝛾_𝑖⁽⁶⁾𝜈⁶+⋯, (3.5) we obtain the following additional conditions for the modified RK method (2.2) to be of up to order five (see [16]): (i)order 3 requires: _𝑖𝑏_𝑖𝛾_𝑖⁽²⁾=0; (3.6)(ii)order 4 requires in addition: _𝑖𝑏_𝑖𝑐_𝑖𝛾_𝑖⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑎_𝑖𝑗𝛾_𝑗⁽²⁾=0; (3.7)(iii)order 5 requires in addition: _𝑖𝑏_𝑖𝛾_𝑖⁽²⁾²=0,_𝑖𝑏_𝑖𝛾_𝑖⁽⁴⁾=0,_𝑖𝑏_𝑖𝑐²_𝑖𝛾_𝑖⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑐_𝑖𝑎_𝑖𝑗𝛾_𝑗⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑎_𝑖𝑗𝑐_𝑗𝛾_𝑗⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑎_𝑖𝑗𝑎_𝑗𝑘𝛾_𝑘⁽²⁾=0. (3.8)By simple calculation, it is verified that the new method is of algebraic order fifth. We denote the new method as MODRK5PLDPLAM.

4. Analysis of Stability and Phase Properties

In this section, we are interested in the stability and phase properties of the new method. Lambert and Watson’s stability theory [32] was reformulated by Coleman and Ixaru [33] for the periodicity of exponentially fitted symmetric methods for 𝑦^′′=𝑓(𝑥,𝑦). Van de Vyver [34] adapted this theory to RK methods. Following Van de Vyver’s approach, we consider the test equation as follows: 𝑦^′=𝑖𝜆𝑦,𝜆>0. (4.1) Applying the modified RK method (2.2) to test (4.1) yields the difference equation 𝑦_𝑛+1=𝑀(𝑖𝜃,𝜈)𝑦_𝑛,𝜃=𝜆ℎ, (4.2) where 𝑀()=(𝜈)𝑏𝑖𝜃,𝜈det𝐼−𝑖𝜃𝐴+𝑖𝜃𝛾^𝑇det(𝐼−𝑖𝜃𝐴) (4.3) with 𝐼 the 𝑠×𝑠 identity matrix.

Definition 4.1 (see [34]). For the modified RK method (2.2) with stability function 𝑀(𝑖𝜃,𝜈), the region in the 𝜃-𝜈 plane Ω∶={(𝜃,𝜈)∶|𝑀(𝑖𝜃,𝜈)|≤1} (4.4) is called the region of imaginary stability. And any closed curve defined by |𝑀(𝑖𝜃,𝜈)|=1 is a stability boundary of the method.

In Figure 1 we plot the region of imaginary stability for the method MODRK5PLDPLAM.

Figure 1 

Regions of imaginary stability of the MODDPHARK5 method.

Definition 4.2 (see [34]). For the modified RK method (2.2) with stability function 𝑀(𝑖𝜃,𝜈), the quantities 𝑃(𝜃,𝜈)=𝜃−arg(𝑀(𝑖𝜃,𝜈)),𝐷(𝜃,𝜈)=1−|𝑀(𝑖𝜃,𝜈)| (4.5) are called the phase lag (dispersion) and amplification factor error (dissipation), respectively. If 𝑃(𝜃,𝜈)=𝑐_𝜙𝜃^𝑞+1+𝒪𝜃^𝑞+3,𝐷(𝜃,𝜈)=𝑐_𝑑𝜃^𝑝+1+𝒪𝜃^𝑝+3, (4.6) the method is said to be of phase-lag order 𝑞 and dissipation order 𝑝, respectively, where the 𝑐_𝜙 and 𝑐_𝑑 are called the phase-lag constant and dissipation constant, respectively.

We note that, by definition, when 𝜈=𝜃 (𝜔=𝜆), it must be true that 𝑃(𝜃,𝜈)=0 and 𝐷(𝜃,𝜈)=0. In general, 𝜔≠𝜆 since the fitting frequency 𝜔 is just an estimate of the true frequency. Therefore the order of 𝑃(𝜃,𝜈)=0 and 𝐷(𝜃,𝜈)=0 in Definition 4.2 measure to what extent a modified RK method is accurate in phase and dissipation. Denoting the ratio 𝑟=𝜈/𝜃=𝜔/𝜆, we obtain the following expressions for the phase lag and the dissipation error of the new method MODRK5PLDPLAM: 𝑃()𝑟𝜃,𝑟𝜃=−²−1²29955+11552𝑟²𝜃62905500⁷+𝒪𝜃⁹,𝐷()=𝑟𝜃,𝑟𝜃²−1−13979+10200𝑟²𝜃50324400⁶+𝒪𝜃⁸. (4.7) Thus, the method MODRK5PLDPLAM has a phase lag of order six and a dissipation of order five.

5. 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, respectively. We compare the new method with some existing highly efficient methods in the literature.

The methods we choose for comparison are as follows: (i)PHARK5S: the phase-fitted fifth-order RK method given by Simos in [17],(ii)MODPHARK5S: the modified phase-fitted fifth-order RK method given by Simos and Aguiar in [18],(iii)MODPHARK5V: the higher-order method of the modified phase-fitted embedded RK5 (2.4) pair given by Van de Vyver in [16],(iv)ARK5: an adapted fifth-order RK method given by Fang et al. in [35],(v)PHADISRK5S: the phase-fitted and zero dissipation fifth-order RK method given by Tsitouras and Simos in [22],(vi)MODRK5PLDPLAM: the phase-fitted fifth-order method derived in this paper.

We consider the numerical integration of the Schrödinger equation (1.1) with the well-known Woods-Saxon potential 𝑣(𝑥)=𝑐₀𝑧(1−𝑎(1−𝑧)), (5.1) where 𝑧=(exp(𝑎(𝑥−𝑏)+1))⁻¹,𝑐₀=−50,𝑎=5/3,𝑏=7. The problem is solved in the interval [0,15]. Following [16, 36–38], we choose the fitting frequency √𝜔=[],√50+𝐸,𝑥∈0,6.5𝐸,𝑥∈[6.5,15]. (5.2) In the numerical experiment we consider the resonance problem (𝐸>0), the numerical results 𝐸_calculated are compared with the analytical solution 𝐸_analytical of the Woods-Saxon potential, rounded to six decimal places. In Figures 2, 3, 4, and 5, we plot the error −log₁₀|𝐸_analytical−𝐸_calculated| versus 𝑁 (with the integration step-size 1/2^𝑁) for 𝐸_analytical=53.588872, 163.215341, 341.495874, and 989.701916, respectively.

Figure 2 

Efficiency curves for 𝐸=53.588872.

Figure 3 

Efficiency curves for 𝐸=163.215341.

Figure 4 

Efficiency curves for 𝐸=341.495874.

Figure 5 

Efficiency curves for 𝐸=989.701916.

6. 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].

Acknowledgments

The authors are deeply grateful to the anonymous referees for their constructive comments and valuable suggestions. This research is partially supported by NSFC (no. 11101357), the foundation of Shandong Outstanding Young Scientists Award Project (no. BS2010SF031), the foundation of Scientific Research Project of Shandong Universities (no. J11LG69), and NSF of Shandong Province, China (no. ZR2011AL006).