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

Authors: Qinghe Ming; Yanping Yang; Yonglei Fang
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 1. Introduction

In this paper, we are concerned with the numerical integration of the one-dimensional Schrödinger equation of the form(1.1)y′′(x)=(v(x)-E)y(x),
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–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)y′(x)=f(x,y),y(x0)=y0,
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:
(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),
where the coefficients aij,  ci,  bi,  i=1,…,s are constants, h is the step size, and the parameters γi,i=1,…,s are even functions of ν=hω. It is convenient to express the modified RK method (2.2) by the Butcher tableau as follows:
(2.3)c1γ1c2γ2a21⋮⋮⋮⋱csγsas1⋯ass-1b1⋯bs-1bs
or simply by (c,γ,A,b). The extra-frequency-depending parameters γi(ν),  ν=hω,  i=1,…,s are introduced to tune the traditional RK method to the special oscillatory structure of 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 aij,  ci,  bi,  i=1,…,s be functions of ν=hω (see [16, 18, 29]).Applying the modified RK method (2.2) to the test equation as follows:
(2.4)y′=iωy,ω>0
yields
(2.5)yn+1=R(iν)yn,ν=ωh.
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 thephase lag (or dispersion) and the amplification factor error (or dissipation error), respectively:
(2.6)P(ν)=ν-arg(R(iν)),D(ν)=1-|R(iν)|.
The method is said to be dispersive of orderq and dissipative of orderp if
(2.7)P(ν)=𝒪(νq+1),D(ν)=𝒪(νp+1).
If P(ν)=0 and D(ν)=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)R(iν)=U(ν2)+iνV(ν2),
where
(2.9)U(ν2)=1-t2ν2+t4ν4+⋯,V(ν2)=1-t3ν2+t5ν4+⋯
are polynomials in ν2, which are completely defined by the Runge-Kutta coefficients c,  A,  γ, and b. Therefore, we have
(2.10)P(ν)=ν-arctan(νV(ν2)U(ν2)),D(ν)=1-(U(ν2))2+ν2(V(ν2))2.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 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.


## 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:(3.1)01015γ2150310γ3340940045γ44445-56153290891193726561-253602187644486561-21272901190173168-3553346732524749176-51031865603538405001113125192-2187678411840
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.2)P(ν)=tan(ν)-MN,d(ν)=1-M2+N2,der·P(ν)=sec2(ν)-M′N-MN′N2,
where
(3.3)M=15ν(474651-688905γ5-96460ν2+17808γ2ν4-11130γ4(2ν2-125)-640γ3(371ν2-1500)),N=-7(225(1855+6400γ3+2650γ4-729γ5)ν2-4579200+21200(3γ2-8γ3-4)ν4+7632ν6).
Now, solving (3.2), we get γ values in terms of ν. Instead of giving the very complicated expressions for γi, for the purpose of practical computation, we present their Taylor expansions as follows:
(3.4)γ2=1-24179ν2698950-491109813ν4279160630000-3747105974663311ν68278634104933500000-5235134512534020593713ν850148377999575005150000000+⋯,γ3=1+901ν4998500-146822137ν68973020250000+390542419221781ν839422067166350000000-3826206167449731276763ν101074608099990892967500000000+⋯,γ4=1-1088ν41747375+1151652176ν63925696359375+5225984025866ν8239543810906640625+52283859929609197732ν109794605078041993193359375+⋯.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
(3.5)γi=1+γi(2)ν2+γi(4)ν4+γi(6)ν6+⋯,
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:(3.6)∑i‍biγi(2)=0;(ii)
order 4 requires in addition:(3.7)∑i‍biciγi(2)=0,∑ij‍biaijγj(2)=0;(iii)
order 5 requires in addition:(3.8)∑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.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 y′′=f(x,y). Van de Vyver [34] adapted this theory to RK methods. Following Van de Vyver’s approach, we consider the test equation as follows:
(4.1)y′=iλy,λ>0.
Applying the modified RK method (2.2) 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 (see [34]).
For the modified RK method (2.2) with stability function M(iθ,ν), the region in the θ-ν plane
(4.4)Ω:={(θ,ν):|M(iθ,ν)|≤1}
is called the region of imaginary stability. And any closed curve defined by |M(iθ,ν)|=1 is a stability boundary of the method.In Figure1 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 M(iθ,ν), the quantities
(4.5)P~(θ,ν)=θ-arg(M(iθ,ν)),D~(θ,ν)=1-|M(iθ,ν)|
are called the phase lag (dispersion) and amplification factor error (dissipation), respectively. If
(4.6)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.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)P~(θ,rθ)=-(r2-1)2(29955+11552r2)62905500θ7+𝒪(θ9),D~(θ,rθ)=(r2-1)(-13979+10200r2)50324400θ6+𝒪(θ8).
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
(5.1)v(x)=c0z(1-a(1-z)),
where z=(exp(a(x-b)+1))-1,c0=-50,a=5/3,b=7. The problem is solved in the interval [0,15]. Following [16, 36–38], we choose the fitting frequency
(5.2)ω={50+E,x∈[0,6.5],E,x∈[6.5,15].
In the numerical experiment we consider the resonance problem (E>0), the numerical results Ecalculated are compared with the analytical solution Eanalytical of the Woods-Saxon potential, rounded to six decimal places. In Figures 2, 3, 4, and 5, we plot the error -log10|Eanalytical-Ecalculated| versus N (with the integration step-size 1/2N) for Eanalytical=53.588872, 163.215341, 341.495874, and 989.701916, respectively.Figure 2
Efficiency curves forE=53.588872.Figure 3
Efficiency curves forE=163.215341.Figure 4
Efficiency curves forE=341.495874.Figure 5
Efficiency curves forE=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].



---

*Source: 867948-2012-10-10.xml*