A New Optimized Runge-Kutta Pair for the Numerical Solution of the Radial Schrödinger Equation

Abstract

A new embedded pair of explicit modified Runge-Kutta (RK) methods for the numerical integration of the radial Schrödinger equation is presented. The two RK methods in the pair have algebraic orders five and four, respectively. The two methods of the embedded pair are derived by nullifying the phase lag, the first derivative of the phase lag of the fifth-order method, and the phase lag of the fourth-order method. Nu merical experiments show the efficiency and robustness of our new methods compared with some well-known integrators in the literature.

1. Introduction

In molecular dynamics, quantum physics, and chemistry, no other equation has been studied more profoundly than the Schrödinger equation [1–3]. The one-dimensional Schrödinger equation has the form 𝑦^′′(𝑥)=(𝑣(𝑥)−𝐸)𝑦(𝑥), (1.1) where 𝐸 is a real number denoting the energy, the function 𝑣(𝑥) is the effective potential satisfying 𝑣(𝑥)→0 as 𝑥→∞. There have been a lot of numerical methods, such as exponentially fitted and phase fitted integrators based on the oscillatory property of the solution of the Schrödinger equation (1.1) (see, e.g., [4–13]). In [7], Simos and Aguiar 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, Vyver improved this method and gave an embedded pair of modified RK methods by nullifying the dispersion (phase lag) of the fifth-order method and the fourth-order method [5].

In this paper, we derive a new kind of phase fitting RK embedded pair by nullifying the phase lag, the first derivative of phase lag of the fifth-order method, and the phase lag of the fourth-order method.

2. Order Conditions and Phase Properties of Modified Runge-Kutta Methods

2.1. Modified Runge-Kutta Methods

For the numerical integration of the initial-value problem of first-order differential equations 𝑦^′(𝑥)=𝑓𝑥,𝑦,𝑦𝑥₀=𝑦₀, (2.1) we consider the 𝑠-stage modified explicit Runge-Kutta method of the form 𝑌_𝑖=𝛾_𝑖𝑦_𝑛+ℎ^𝑖−1_𝑗=1𝑎_𝑖𝑗𝑓𝑥_𝑛+𝑐_𝑗ℎ,𝑌_𝑗𝑦,𝑖=1,…,𝑠,_𝑛+1=𝑦_𝑛+ℎ^𝑠_𝑖=1𝑏_𝑖𝑓𝑥_𝑛+𝑐_𝑖ℎ,𝑌_𝑖, (2.2) which can be expressed in Butcher tableau as 𝑐𝛾𝐴𝑏^𝑇=01𝑐₂𝛾₂𝑎₂₁⋮⋮𝑐⋮⋱_𝑠𝛾_𝑠𝑎_𝑠1⋯𝑎_𝑠𝑠−1𝑏₁⋯𝑏_𝑠 (2.3) or equivalently by the triplet (𝑐,𝛾,𝐴,𝑏) with 𝑐=(0,𝑐₂,…,𝑐_𝑠)^𝑇,  𝛾=(𝛾₁,…,𝛾_𝑠), 𝐴=(𝑎_𝑖𝑗)_𝑠×𝑠, 𝑏=(𝑏₁,…,𝑐_𝑠). Here, following the approach of exponential fitting and/or phase fitting in [5, 7, 14], the frequency-depending parameters 𝛾_𝑖=𝛾_𝑖(𝜈),  𝜈=ℎ𝜔,  𝑖=1,…,𝑠 are introduced to adapt the traditional RK method to the oscillatory feature of the solution to the problem [15–27], for example, in this paper, to minimize the dispersion and/or dissipation (see next section).

The algebraic order conditions presented in [28] are not fit for the modified RK method (2.2). Writing 𝛾_𝑖=1+𝛾_𝑖⁽²⁾𝜈²+𝛾_𝑖⁽⁴⁾𝜈⁴+𝛾_𝑖⁽⁶⁾𝜈⁶+⋯, (2.4) where 𝛾_𝑖^(2𝑗)=(𝑑^2𝑗𝛾_𝑖/𝑑𝜈^2𝑗)(0), 𝑗=1,2,…, the third-to-fifth order conditions are listed as follows (see Vyver [5]).(i)Order 3 requires _𝑖𝑏_𝑖𝛾_𝑖⁽²⁾=0; (2.5)(ii) order 4 requires, in addition _𝑖𝑏_𝑖𝑐_𝑖𝛾_𝑖⁽²⁾=0,_𝑖,𝑗𝑏_𝑖𝑎_𝑖𝑗𝛾_𝑗⁽²⁾=0; (2.6)(iii)order 5 requires, in addition _𝑖𝑏_𝑖𝛾_𝑖⁽²⁾²=0,_𝑖𝑏_𝑖𝛾_𝑖⁽⁴⁾=0,_𝑖𝑏_𝑖𝑐²_𝑖𝛾_𝑖⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑐_𝑖𝑎_𝑖𝑗𝛾_𝑗⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑎_𝑖𝑗𝑐_𝑗𝛾_𝑗⁽²⁾=0,_𝑖𝑗𝑏_𝑖𝑎_𝑖𝑗𝑎_𝑗𝑘𝛾_𝑘⁽²⁾=0. (2.7)

2.2. Dispersion and Dissipation of Modified Runge-Kutta Methods

Applying the modified RK method (2.2) to the test differential equation 𝑦^′=i𝜔𝑦,𝜔>0, (2.8) yields 𝑦_𝑛+1=𝑅(i𝜈)𝑦_𝑛,𝜈=ℎ𝜔, (2.9) where 𝑅(i𝜈)=1−i𝜈𝑏^𝑇(𝐼−i𝜈𝐴)⁻¹𝛾=1−i𝜈𝑏^𝑇𝛾+𝜈²𝑏^𝑇𝐴𝛾+i𝜈³𝑏^𝑇𝐴²𝛾−⋯+(i𝜈)^𝑠𝑏^𝑇𝐴^𝑠−1𝛾 (2.10) with 𝐼 the 𝑠×𝑠 identity matrix.

Definition 2.1. The quantities 𝑃(𝜈)=𝜈−arg(𝑅(i𝜈)),𝐷(𝜈)=1−|𝑅(i𝜈)|, (2.11) are called the dispersion (or phase lag) and the dissipation (or amplification factor error) of the method (2.2), respectively. And the method is said to be dispersive of order 𝑞 and dissipative of order 𝑟 if 𝑃(𝜈)=𝒪𝜈^𝑞+1,𝐷(𝜈)=𝒪𝜈^𝑟+1, (2.12) respectively.

Writing 𝑅(𝜈)=𝑈𝜈²+i𝜈𝑉𝜈², (2.13) we have 𝑃(𝜈)𝜈𝑉𝜈=𝜈−arctan²𝑈𝜈²,𝐷(𝜈)=1−𝑈²𝜈²+𝜈²𝑉²𝜈². (2.14) From the formula (2.10), 𝑈(𝜈²) and 𝑉(𝜈²) are polynomials in 𝜈²: 𝑈𝜈²=1−𝜈²𝑏^𝑇𝐴𝛾+𝜈⁴𝑏^𝑇𝐴³𝑉𝜈𝛾−⋯,₂=𝑏^𝑇𝛾−𝜈²𝑏^𝑇𝐴²𝛾+𝜈⁴𝑏^𝑇𝐴⁴𝛾−⋯, (2.15) which are completely determined by 𝑐,  𝐴,  𝛾, and 𝑏.

3. Construction of the New Embedded Pair

In this section, we are concerned with the embedded modified pair which is simplify denoted by the Butcher tableau 𝑐𝐴𝑏^𝑇𝑏^∗𝑇=01015𝛾₂150310𝛾₃394004405𝛾₄44−455615329089119372−656125360218764448−65612120172919017−31683553346732524749−1765103011865613503845001113125192−21876784110843503845001113125192−21876784110845179057600757116695393640−920973392001871210040 (3.1) Choosing 𝛾₂=𝛾₃=𝛾₄=1, the classical embedded RK5(4) pair derived by Dormand and Prince [29] is recovered, where the method (𝑐,𝐴,𝑏) is of order five and (𝑐,𝐴,𝑏^∗) is of order four. It should be noted that the first method (𝑐,𝐴,𝑏) in this pair shares the FSAL property in the sense that it uses only six function evaluations at each step with seven stages. Simos and Aguiar [7] presented a modified phase fitted RK method by determining the one-parameter 𝛾₂. Following this approach, Vyver [5] constructed a phase fitted embedded RK5(4) pair. Our main aim in this section is to derive a more efficient embedded RK pair.

Inspired by the ideas in [30–40], with a variant expression of dispersion (see Simos [41]), we compute the dispersion of the higher-order method and the dispersion of the lower-order method and the first derivative of dispersion of the higher-order method of the pair (3.1) as follows: 𝑃_𝐻(𝜈)(𝜈)−𝑀=tan¹𝑁₁,𝑃_𝐿(𝜈)(𝜈)−𝑀=tan²𝑁₂,Der.𝑃_𝐻(𝜈)=sec²(𝜈)−𝑀^′₁𝑁₁−𝑁^′₁𝑀₁𝑁²₁, (3.2) in which 𝑀₁=15𝜈107145+48230𝜈²−8904𝛾₂𝜈⁴+5565𝛾₄2𝜈²−125+320𝛾₃371𝜈²,𝑁−1500₁=7225563+3200𝛾₃+1325𝛾₄𝜈²+106003𝛾₂−4−8𝛾₃𝜈⁴+3816𝜈⁶,𝑀−2289600₂=𝜈216532500+173461225𝜈²−113080800𝛾₂𝜈²−9283904𝜈⁴−17051160𝛾₂𝜈⁴+133560𝜈⁶+83475𝛾₄397𝜈²−23580−8000𝛾₃181704−55090𝜈²+371𝜈⁴,𝑁₂=−3205440000+752536615+13703168𝛾₃+5129817𝛾₄𝜈²−37109263+2352𝛾₂+24160𝛾₃+225𝛾₄𝜈⁴[phone]+25𝛾₂𝜈⁶. (3.3) Solving the systems of (3.2) we obtain 𝛾₂=−25𝜈102879000−16468980𝜈²+5693624𝜈⁶+52050𝜈⁸cos(𝜈)[phone]cos(𝜈)²sin(𝜈)+𝜈−2576041000−26581020𝜈²3569296𝜈⁴+469645𝜈⁶cos(3𝜈)−2𝜈1225233750+380478750𝜈²[phone]𝜈⁴−19537525𝜈⁶−233289𝜈⁸+3710𝜈¹⁰×sin(𝜈)−50𝜈46352550−1042218𝜈²−804802𝜈⁴+97573𝜈⁶sin(3𝜈)/5𝜈⁴𝑀,𝛾₃=𝜈211041495000+𝜈²−11689730450+𝜈²−1264487200+7𝜈²58586875+212𝜈²45𝜈²(𝜈)(𝜈)−8753×cos−212026500000cos²sin(𝜈)+5𝜈−5−39400200−181538138𝜈²+7037128𝜈⁴+878157𝜈⁶cos(𝜈)−2𝜈7258856775+𝜈²−1213296080+7𝜈²17587275+53𝜈²480𝜈²(𝜈)−56077×sin−10𝜈149294235−22496216𝜈²+2498685𝜈⁴sin(3𝜈)/200𝜈²𝑀,𝛾₄=−𝜈14108160000−5363145725𝜈²[phone]𝜈⁴[phone]𝜈⁶−6105436𝜈⁸+89040𝜈¹⁰×cos(𝜈)[phone]cos(𝜈)²sin(𝜈)+𝜈−25167193600+12945389𝜈²−2421645𝜈⁴+68704𝜈⁶cos(3𝜈)−2𝜈2084367750+552861425𝜈²−100100925𝜈⁴[phone]𝜈⁶−1192856𝜈⁸+14840𝜈¹⁰×sin(𝜈)+50𝜈−6590070−3738885𝜈²+176272𝜈⁴sin(3𝜈)/25𝜈²𝑀, (3.4)(3.5) where 𝑀=𝜈2013125+23057425𝜈²−1145578𝜈⁴+7420𝜈⁶cos(𝜈)[phone]cos(𝜈)²sin(𝜈)−25𝜈1059125−97573𝜈²cos(3𝜈)+2𝜈2650105−310653𝜈²+1855𝜈⁴(𝜈)(.sin+389770𝜈sin3𝜈) (3.6) For small values of 𝜈, say |𝜈|<0.04, the above formulae are subject to heavy cancelations and in that case the following Taylor series expansions must be used 𝛾₂=1−22051𝜈²−117590023696602511𝜈⁴−104535205236001936263085085455799𝜈⁶−921922108777593000000132294694246155651595747𝜈⁸𝛾112745173602003447304800000+⋯,₃=1−28991𝜈⁴+58795002740771225019𝜈⁶−1866700093500000633801949482499571𝜈⁸−1317031583967990000000708641114757728596912360427𝜈¹⁰𝛾2113972005037564636965000000000+⋯,₄=1+35008𝜈⁴−1028912577729377736𝜈⁶−4083406454531251838102360373886𝜈⁸−80027960831388281255030427774325408868897437𝜈¹⁰19267974004248636014003906250+⋯. (3.7) It is easy to check that the two schemes in the new pair (3.1) with 𝛾-values (3.5) are of algebraic orders five and four, respectively. The Taylor series of the dissipations of the new fifth-order method and the fourth-order method are given by 𝐷_𝐻(𝜈)=−215377𝜈⁶+197551200637415987609𝜈⁸−94081684712400004487714698722250553𝜈¹⁰𝐷66378391831986696000000+⋯,_𝐿(𝜈)=−10822793𝜈⁶+92602125007211815545941𝜈⁸−261338013090000000154473286018569080117𝜈¹⁰1382883163166389500000000+⋯, (3.8) respectively.

4. Stability Analysis

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

Definition 4.1. For the modified RK method (3.1) with stability function 𝑀(i𝜃,𝜈), the quantities 𝑃(𝜃,𝜈)=𝜃−arg(𝑀(i𝜃,𝜈)),𝐷(𝜃,𝜈)=1−|𝑀(i𝜃,𝜈)| (4.4) are called the phase lag (dispersion) and amplification factor error (dissipation), respectively. If 𝑃(𝜃,𝜈)=𝑐_𝜙𝜃^𝑞+1+𝒪𝜃^𝑞+3,𝐷(𝜃,𝜈)=𝑐_𝑑𝜃^𝑝+1+𝒪𝜃^𝑝+3, (4.5) 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.

For the convenience of analyzing the phase lag and the dissipation, we denote the ratio 𝑟=𝜈/𝜃. Then the phase lags and the dissipations of the higher-order method and the lower order method are 𝑃_𝐻()𝑟𝜃,𝑟𝜃=−²−1²176385+198413𝑟²𝜃370408500⁷+𝒪𝜃⁹,𝐷_𝐻(𝜃,𝑟𝜃)=164626−154357𝑟²−656400𝑟⁴𝜃592653600⁶+𝒪𝜃⁸,𝑃_𝐿()𝑟𝜃,𝑟𝜃=−²−13421869+621392𝑟²𝜃4233240000⁵+𝒪𝜃⁷,𝐷_𝐿(𝜃,𝑟𝜃)=−31690505+26395047𝑟²−37995714𝑟⁴𝜃37040850000⁶+𝒪𝜃⁸, (4.6) respectively. Thus, the higher-order method has a phase lag of order six and a dissipation of order five and the low-order method is of phase lag order four and dissipation order five.

5. Numerical Experiments

In this section, we will compare the numerical performance of the new pair with some existing well-known methods proposed in the scientific literature.

5.1. Comparison with Fixed Step-Size Methods

The following fixed step-size methods are selected for comparison: (i)PHARK5S: the phase fitted fifth-order RK method derived by Simos [6]; (ii)MODPHARK5S: the modified phase fitted fifth-order RK method obtained by Simos and Aguiar in [7]; (iii)MODPHARK5V: the higher-order method of the modified phase fitted embedded RK5(4) pair obtained by Vyver presented in [5]; (iv)ARK5: an adapted fifth-order RK method given by Fang et al. in [45];(v)MODDPHARK5: the higher-order method of the new pair.

We consider the numerical integration of the Schrödinger equation (1.1) with the well-known Woods-Saxon potential 𝑣(𝑥)=𝑐₀𝑧(1−𝑎(1−𝑧)), (5.1) in which 𝑧=(exp(𝑎(𝑥−𝑏)+1))⁻¹,  𝑐₀=−50,  𝑎=5/3, and  𝑏=7. The domain of numerical integration is [0,15]. It is appropriate to choose 𝜔 as follows [5, 46]: √𝜔=[],√50+𝐸,𝑥∈0,6.5𝐸,𝑥∈[6.5,15]. (5.2) In the numerical experiments we consider the resonance problem (𝐸>0), the numerical results were compared with the analytical solution of the Woods-Saxon potential, rounded to six decimal places. In Figures 1, 2, 3, and 4, we plot the error of −log₁₀|𝐸_analytical−𝐸_calculated| versus the integration step-size (1/2^𝑁) for 𝐸_analytical=53.588872,163.215341,341.495874,and989.701916, respectively.

Figure 1 

Efficiency curves for 𝐸=53.588872.

Figure 2 

Efficiency curves for 𝐸=163.215341.

Figure 3 

Efficiency curves for 𝐸=341.495874.

Figure 4 

Efficiency curves for 𝐸=989.701916.

5.2. Comparison with Variable Step-Size Methods

Next we select the following embedded RK5(4) pairs: (i)PHARK5(4)S: the phase fitted embedded RK5(4) pair derived by Simos [6];(ii)MODPHARK5(4)V: the modified phase fitted embedded RK5(4) pair obtained by Vyver presented in [5]; (iii)MODDPHARK5(4): the new embedded RK5(4) pair given in this paper.

We consider the numerical integration of the Schrödinger equation (1.1) with the well-known Lennard-Jones potential (see [5]) 𝑣(𝑥)=𝑙(𝑙+1)𝑥²1+500𝑥¹²−1𝑥⁶, (5.3) and we compute some phase shifts for this potential. In the numerical results, we compute the phase shifts correct to four decimal places for the energies 𝑘²=2500 and 𝑘²=10000. We choose the fitting frequency 𝜔=𝑘. For the calculation of the phase shifts, we show the number of function evaluations as a function of 𝑙=0,1,…,10 in Figures 5-6.

Figure 5 

Efficiency curves for 𝑘=50.

Figure 6 

Efficiency curves for 𝑘=100.

Figures 1–6 show that our new methods are more efficient than the other methods we select for comparison.

6. Conclusions and Discussions

A new kind of modified phase fitted explicit embedded RK pair for the numerical integration of the radial Schrödinger equation is presented in this paper. This new pair is based on the classical RK5(4) pair obtained by Dormand and Prince [29]. The phase fitted technique can be regarded as an improvement of the ideas from [5, 7, 30, 31]. The two schemes in this pair are of orders five and four, respectively. Numerical experiments show the effectiveness and competence of the new pair.

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 Shangdong Universities (no. J11LG69), NSF of Shandong Province, China (no. ZR2011AL006), NSF of Universities of Anhui Province, China (no. KJ2010A248), and the Scientific Research Start-up Fund of Chuzhou University, China (no. 2010qd03).