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

Authors: Yonglei Fang; Qinghong Li; Qinghe Ming; Kaimin Wang
Journal: Abstract and Applied Analysis
 (2012)

---

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


---

## Body

## 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)y′′(x)=(v(x)-E)y(x),
where E is a real number denoting the energy, the function v(x) is the effective potential satisfying v(x)→0 as x→∞. 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)y′(x)=f(x,y),y(x0)=y0,
we consider the s-stage modified explicit Runge-Kutta method of the form
(2.2)Yi=γiyn+h∑j=1i-1‍aijf(xn+cjh,Yj),i=1,…,s,yn+1=yn+h∑i=1s‍bif(xn+cih,Yi),
which can be expressed in Butcher tableau as
(2.3)cγAbT=01c2γ2a21⋮⋮⋮⋱csγsas1⋯ass-1b1⋯bs
or equivalently by the triplet (c,γ,A,b) with c=(0,c2,…,cs)T,  γ=(γ1,…,γs), A=(aij)s×s, b=(b1,…,cs). Here, following the approach of exponential fitting and/or phase fitting in [5, 7, 14], the frequency-depending parameters γi=γi(ν),  ν=hω,  i=1,…,s 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
(2.4)γi=1+γi(2)ν2+γi(4)ν4+γi(6)ν6+⋯,
where γi(2j)=(d2jγi/dν2j)(0), j=1,2,…, the third-to-fifth order conditions are listed as follows (see Vyver [5]).(i)
Order 3 requires(2.5)∑i‍biγi(2)=0;(ii)
order 4 requires, in addition(2.6)∑i‍biciγi(2)=0,∑i,j‍biaijγj(2)=0;(iii)
order 5 requires, in addition(2.7)∑i‍bi(γi(2))2=0,∑i‍biγi(4)=0,∑i‍bici2γi(2)=0,∑ij‍biciaijγj(2)=0,∑ij‍biaijcjγj(2)=0,∑ij‍biaijajkγk(2)=0.


### 2.2. Dispersion and Dissipation of Modified Runge-Kutta Methods

Applying the modified RK method (2.2) to the test differential equation
(2.8)y′=iωy,ω>0,
yields
(2.9)yn+1=R(iν)yn,ν=hω,
where
(2.10)R(iν)=1-iνbT(I-iνA)-1γ=1-iνbTγ+ν2bTAγ+iν3bTA2γ-⋯+(iν)sbTAs-1γ
with I the s×s identity matrix.Definition 2.1.
The quantities(2.11)P(ν)=ν-arg(R(iν)),D(ν)=1-|R(iν)|,
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 q and dissipative of order r if
(2.12)P(ν)=𝒪(νq+1),D(ν)=𝒪(νr+1),
respectively.Writing(2.13)R(ν)=U(ν2)+iνV(ν2),
we have
(2.14)P(ν)=ν-arctan(νV(ν2)U(ν2)),D(ν)=1-U2(ν2)+ν2V2(ν2).
From the formula (2.10), U(ν2) and V(ν2) are polynomials in ν2:
(2.15)U(ν2)=1-ν2bTAγ+ν4bTA3γ-⋯,V(ν2)=bTγ-ν2bTA2γ+ν4bTA4γ-⋯,
which are completely determined by c,  A,  γ, and b.


## 2.1.  Modified Runge-Kutta Methods

For the numerical integration of the initial-value problem of first-order differential equations(2.1)y′(x)=f(x,y),y(x0)=y0,
we consider the s-stage modified explicit Runge-Kutta method of the form
(2.2)Yi=γiyn+h∑j=1i-1‍aijf(xn+cjh,Yj),i=1,…,s,yn+1=yn+h∑i=1s‍bif(xn+cih,Yi),
which can be expressed in Butcher tableau as
(2.3)cγAbT=01c2γ2a21⋮⋮⋮⋱csγsas1⋯ass-1b1⋯bs
or equivalently by the triplet (c,γ,A,b) with c=(0,c2,…,cs)T,  γ=(γ1,…,γs), A=(aij)s×s, b=(b1,…,cs). Here, following the approach of exponential fitting and/or phase fitting in [5, 7, 14], the frequency-depending parameters γi=γi(ν),  ν=hω,  i=1,…,s 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
(2.4)γi=1+γi(2)ν2+γi(4)ν4+γi(6)ν6+⋯,
where γi(2j)=(d2jγi/dν2j)(0), j=1,2,…, the third-to-fifth order conditions are listed as follows (see Vyver [5]).(i)
Order 3 requires(2.5)∑i‍biγi(2)=0;(ii)
order 4 requires, in addition(2.6)∑i‍biciγi(2)=0,∑i,j‍biaijγj(2)=0;(iii)
order 5 requires, in addition(2.7)∑i‍bi(γi(2))2=0,∑i‍biγi(4)=0,∑i‍bici2γi(2)=0,∑ij‍biciaijγj(2)=0,∑ij‍biaijcjγj(2)=0,∑ij‍biaijajkγk(2)=0.


## 2.2. Dispersion and Dissipation of Modified Runge-Kutta Methods

Applying the modified RK method (2.2) to the test differential equation
(2.8)y′=iωy,ω>0,
yields
(2.9)yn+1=R(iν)yn,ν=hω,
where
(2.10)R(iν)=1-iνbT(I-iνA)-1γ=1-iνbTγ+ν2bTAγ+iν3bTA2γ-⋯+(iν)sbTAs-1γ
with I the s×s identity matrix.Definition 2.1.
The quantities(2.11)P(ν)=ν-arg(R(iν)),D(ν)=1-|R(iν)|,
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 q and dissipative of order r if
(2.12)P(ν)=𝒪(νq+1),D(ν)=𝒪(νr+1),
respectively.Writing(2.13)R(ν)=U(ν2)+iνV(ν2),
we have
(2.14)P(ν)=ν-arctan(νV(ν2)U(ν2)),D(ν)=1-U2(ν2)+ν2V2(ν2).
From the formula (2.10), U(ν2) and V(ν2) are polynomials in ν2:
(2.15)U(ν2)=1-ν2bTAγ+ν4bTA3γ-⋯,V(ν2)=bTγ-ν2bTA2γ+ν4bTA4γ-⋯,
which are completely determined by c,  A,  γ, and b.


## 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(3.1)cAbTb*T=01015γ2150310γ3340940045γ44445-56153290891193726561-253602187644486561-21272901190173168-3553346732524749176-5103186560113538405001113125192-21876784118403538405001113125192-21876784118405179576000757116695393640-920973392001872100140
Choosing γ2=γ3=γ4=1, the classical embedded RK5(4) pair derived by Dormand and Prince [29] is recovered, where the method (c,A,b) is of order five and (c,A,b*) is of order four. It should be noted that the first method (c,A,b) 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 γ2. 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:
(3.2)PH(ν)=tan(ν)-M1N1,PL(ν)=tan(ν)-M2N2,Der.PH(ν)=sec2(ν)-M1′N1-N1′M1N12,
in which
(3.3)M1=15ν(107145+48230ν2-8904γ2ν4+5565γ4(2ν2-125)+320γ3(371ν2-1500)),N1=7(225(563+3200γ3+1325γ4)ν2+10600(3γ2-4-8γ3)ν4+3816ν6-2289600),M2=ν(216532500+173461225ν2-113080800γ2ν2-9283904ν4-17051160γ2ν4+133560ν6+83475γ4(397ν2-23580)-8000γ3(181704-55090ν2+371ν4)),N2=-3205440000+75(2536615+13703168γ3+5129817γ4)ν2-3710(9263+2352γ2+24160γ3+225γ4)ν4+26712(136+25γ2)ν6.
Solving the systems of (3.2) we obtain
(3.4)γ2=(-25ν(102879000-16468980ν2+5693624ν6+52050ν8)cos(ν)+4473000000cos(ν)2sin(ν)+ν(-25(76041000-26581020ν23569296ν4+469645ν6)cos(3ν)-2ν(1225233750+380478750ν2+59883450ν4-19537525ν6-233289ν8+3710ν10)×sin(ν)-50ν(46352550-1042218ν2-804802ν4+97573ν6)sin(3ν)))/(5ν4M),γ3=(ν(∑211041495000+ν2(-11689730450+ν2(-1264487200+7ν2(58586875+212ν2(45ν2-8753))))∑)×cos(ν)-212026500000cos(ν)2sin(ν)+5ν(∑-5(-39400200-181538138ν2+7037128ν4+878157ν6)cos(ν)-2ν(7258856775+ν2(-1213296080+7ν2(17587275+53ν2(480ν2-56077))))×sin(ν)-10ν(149294235-22496216ν2+2498685ν4)sin(3ν)∑))/(200ν2M),(3.5)γ4=(∑-ν(14108160000-5363145725ν2+258455175ν4+22475450ν6-6105436ν8+89040ν10)×cos(ν)+18288000000cos(ν)2sin(ν)+ν(∑-25(167193600+12945389ν2-2421645ν4+68704ν6)cos(3ν)-2ν(2084367750+552861425ν2-100100925ν4+15041455ν6-1192856ν8+14840ν10)×sin(ν)+50ν(-6590070-3738885ν2+176272ν4)sin(3ν)∑))/(25ν2M),
where
(3.6)M=(ffν(2013125+23057425ν2-1145578ν4+7420ν6)cos(ν)+200850000cos(ν)2sin(ν)-25ν(ff(1059125-97573ν2)cos(3ν)+2ν(2650105-310653ν2+1855ν4)sin(ν)+389770νsin(3ν)ff)).
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
(3.7)γ2=1-22051ν21175900-23696602511ν410453520523600-1936263085085455799ν6921922108777593000000-132294694246155651595747ν8112745173602003447304800000+⋯,γ3=1-28991ν45879500+2740771225019ν61866700093500000-633801949482499571ν81317031583967990000000-708641114757728596912360427ν102113972005037564636965000000000+⋯,γ4=1+35008ν410289125-77729377736ν6408340645453125-1838102360373886ν88002796083138828125-5030427774325408868897437ν1019267974004248636014003906250+⋯.
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
(3.8)DH(ν)=-215377ν6197551200+637415987609ν89408168471240000-4487714698722250553ν1066378391831986696000000+⋯,DL(ν)=-10822793ν69260212500+7211815545941ν8261338013090000000-154473286018569080117ν101382883163166389500000000+⋯,
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 y′′=f(x,y). Vyver [44] formulated this theory to RK methods. Following Van de Vyver’s idea, we consider the test equation
(4.1)y′=iλy,λ>0.
Applying the modified RK method (3.1) to test (4.1) yields the difference equation
(4.2)yn+1=M(iθ,ν)yn,θ=λh,
where
(4.3)M(iθ,ν)=det(I-iθA+iθγ(ν)bT)det(I-iθA)
with I the s×s identity matrix.Definition 4.1.
For the modified RK method (3.1) with stability function M(iθ,ν), the quantities
(4.4)P~(θ,ν)=θ-arg(M(iθ,ν)),D~(θ,ν)=1-|M(iθ,ν)|
are called the phase lag (dispersion) and amplification factor error (dissipation), respectively. If
(4.5)P~(θ,ν)=cϕθq+1+𝒪(θq+3),D~(θ,ν)=cdθp+1+𝒪(θp+3),
the method is said to be of phase lag orderq and dissipation orderp, respectively, where the cϕ and cd are called the phase lag constant and dissipation constant, respectively.For the convenience of analyzing the phase lag and the dissipation, we denote the ratior=ν/θ. Then the phase lags and the dissipations of the higher-order method and the lower order method are
(4.6)P~H(θ,rθ)=-(r2-1)2(176385+198413r2)370408500θ7+𝒪(θ9),D~H(θ,rθ)=(164626-154357r2-656400r4)592653600θ6+𝒪(θ8),P~L(θ,rθ)=-(r2-1)(3421869+621392r2)4233240000θ5+𝒪(θ7),D~L(θ,rθ)=(-31690505+26395047r2-37995714r4)37040850000θ6+𝒪(θ8),
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
(5.1)v(x)=c0z(1-a(1-z)),
in which z=(exp(a(x-b)+1))-1,  c0=-50,  a=5/3, and  b=7. The domain of numerical integration is [0,15]. It is appropriate to choose ω as follows [5, 46]:
(5.2)ω={50+E,x∈[0,6.5],E,x∈[6.5,15].
In the numerical experiments we consider the resonance problem (E>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 -log10|Eanalytical-Ecalculated| versus the integration step-size (1/2N) for Eanalytical=53.588872,163.215341,341.495874,and989.701916, respectively.Figure 1
Efficiency curves forE=53.588872.Figure 2
Efficiency curves forE=163.215341.Figure 3
Efficiency curves forE=341.495874.Figure 4
Efficiency curves forE=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])
(5.3)v(x)=l(l+1)x2+500(1x12-1x6),
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 k2=2500 and k2=10000. We choose the fitting frequency ω=k. For the calculation of the phase shifts, we show the number of function evaluations as a function of l=0,1,…,10 in Figures 5-6.Figure 5
Efficiency curves fork=50.Figure 6
Efficiency curves fork=100.Figures1–6 show that our new methods are more efficient than the other methods we select for comparison.


## 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
(5.1)v(x)=c0z(1-a(1-z)),
in which z=(exp(a(x-b)+1))-1,  c0=-50,  a=5/3, and  b=7. The domain of numerical integration is [0,15]. It is appropriate to choose ω as follows [5, 46]:
(5.2)ω={50+E,x∈[0,6.5],E,x∈[6.5,15].
In the numerical experiments we consider the resonance problem (E>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 -log10|Eanalytical-Ecalculated| versus the integration step-size (1/2N) for Eanalytical=53.588872,163.215341,341.495874,and989.701916, respectively.Figure 1
Efficiency curves forE=53.588872.Figure 2
Efficiency curves forE=163.215341.Figure 3
Efficiency curves forE=341.495874.Figure 4
Efficiency curves forE=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])
(5.3)v(x)=l(l+1)x2+500(1x12-1x6),
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 k2=2500 and k2=10000. We choose the fitting frequency ω=k. For the calculation of the phase shifts, we show the number of function evaluations as a function of l=0,1,…,10 in Figures 5-6.Figure 5
Efficiency curves fork=50.Figure 6
Efficiency curves fork=100.Figures1–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.



---

*Source: 641236-2012-11-14.xml*