# A New Trigonometrically Fitted Two-Derivative Runge-Kutta Method for the Numerical Solution of the  Schrödinger Equation and Related Problems

Authors: Yanwei Zhang; Haitao Che; Yonglei Fang; Xiong You
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

A new trigonometrically fitted fifth-order two-derivative Runge-Kutta method with
variable nodes is developed for the numerical solution of the radial Schrödinger equation
and related oscillatory problems. Linear stability and phase properties of the new
method are examined. Numerical results are reported to show the robustness and competence
of the new method compared with some highly efficient methods in the recent
literature.


---

## Body

## 1. Introduction

We are interested in initial value problems related to systems of first-order ordinary differential equations (ODEs) of the form(1)y′=f(x,y),y(x0)=y0,
whose solutions have a pronounced oscillatory character. Such problems often occur in many applied fields such as celestial mechanics, physical chemistry, quantum mechanics, and electronics. The investigation of the numerical solution of (1) has been the subject of extensive research activity during the last few decades [1–16]. Such special optimized methods fall into two classes. The first class consists of the methods with constant coefficients. These methods can be applied to any problem with periodic solution. Methods of the second class have coefficients depending on the frequency of problem. When a good estimate of the frequency is known in advance one can use procedures such as exponential/trigonometric fitting or phase fitting. For example, Chen et al. [17] considered phase-fitted and amplification-fitted RK methods whose updates are also phase fitted and amplification fitted. Their methods are shown to be more efficient than the codes in the literature for some typical test problems and for the Lotka-Volterra system and a two-gene regulatory network in biology as well. You et al. [18] investigated the trigonometrically fitted Scheifele methods for oscillatory problems. A good theoretical foundation of the exponential fitting techniques for multistep methods was presented by Gautschi [19] and Lyche [20].Recently, Simos et al. [21] introduced an exponentially fitted explicit Runge-Kutta method which integrates exactly the model problemy′=iωy. Van de Vyver [22] constructed a new fourth-order explicit Runge-Kutta method based on Simos’ approach [21] for the numerical integration of the Schrödinger equation. And another exponentially fitted Runge-Kutta method with four stages was constructed by Vanden Berghe et al. [23] which exactly integrates differential initial value problems whose solutions are linear combinations of functions of the form{exp(ωx),exp(-ωx)}. More recently, Chan and Tsai [24] presented a family of two-derivative Runge-Kutta (TDRK) methods. An advantage of the two-derivative Runge-Kutta methods over classical Runge-Kutta methods is that they can reach higher order with fewer function evaluations.Inspired by the approaches of Simos [21] and of Chan and Tsai [24], we will construct in this paper a new two-derivative Runge-Kutta method for the numerical integration of the Schrödinger equation and related oscillatory problems. The remainder of this paper is organized as follows: in Section 2, we introduce the scheme of two-derivative Runge-Kutta methods and present (up to fifth) order conditions for TDRK methods. In Section 3, we construct a new trigonometrically fitted TDRK method and analyze its linear stability and phase properties. In Section 4, some numerical examples are given to show the effectiveness and competence of our new method compared to the selected methods in the recent literature. Section 5 is devoted to conclusive remarks.


## 2. Two-Derivative Runge-Kutta Methods

We start with a special form of explicit two-derivative Runge-Kutta (TDRK) methods proposed by Chan and Tsai [24]:
(2)Yi=yn+cihf(xn,yn)+h2∑j=1i-1‍aijg(xn+cjh,Yj),11111111111111111111111111i=2,…,s,yn+1=yn+hf(xn,yn)+h2∑i=1s‍big(xn+cih,Yi),
whereg(x,y)=y′′(x):=(∂f(x,y)/∂x)+(∂f(x,y)/∂y)f(x,y). The coefficients of the scheme can be expressed by the Butcher tableau
(3)cAbT=0c2a21⋮⋮⋱csas1⋯ass-1b1⋯bs-1bs
or simplify by(c,A,b). This explicit TDRK method involves only one evaluation of the functionfandsevaluations of the functiongper step.According to Chan and Tsai [24], the conditions for order up to five are listed as follows.(i)
Order 2 requires(4)∑i=1s‍bi=12.(ii)
Order 3 requires, in addition,(5)∑i=2s‍bici=16.(iii)
Order 4 requires, in addition,(6)∑i=2s‍bici2=112.(iv)
Order 5 requires, in addition,(7)∑i=2s‍bici3=120,∑i=3s‍∑j=2i-1‍biaijcj=1120.In practice, the following simplifying assumption is useful:
(8)∑i=1s‍aij=12ci2,fori=2,…,s.
Choosingc3=4/5and solving (4)–(8) yield a three-stage TDRK method of Chan and Tsai [24] which is given by the following tableau:
(9)0000131180045-212542125054892825336


## 3. A Fifth-Order TDRK Method with Frequency-Dependent Coefficients

In order to adapt TDRK methods to the problem (1) whose solutions share an oscillatory feature withωbeing an accurate estimate of the principal frequency, we allow the coefficients to depend onω. In particular, we consider in this paper the three-stage explicit TDRK method given by the following Butcher tableau:
(10)0000c2(ν)a21(ν)00c3(ν)a31(ν)a32(ν)0b1(ν)b2(ν)b3(ν)


### 3.1. Construction of the New Method

If we require the TDRK method (10) to integrate exactlyexp(iωx); that is,
(11)exp(iν)=cos(ν)+isin(ν)=1+iν-k1ν2-ik2ν31.+k3ν4+ik4ν5-k5ν6,ν=ωh,
then we have
(12)cos(ν)=1-k1ν2+k3ν4-k5ν6,sin(ν)=ν-k2ν3+k4ν5,
where
(13)k1=b1(ν)+b2(ν)+b3(ν),k2=b2(ν)c2(ν)+b3(ν)c3(ν),k3=b2(ν)a21(ν)+b3(ν)(a31(ν)+a32(ν)),k4=b3(ν)a32(ν)c2(ν),k5=b3(ν)a32(ν)a21(ν).
Solving (4)–(6) fors=3, the simplifying conditions (8) and (12), we obtain the coefficients expressed in terms ofc3(ν)(14)b1(ν)=(6(ν2-4)+c3(ν)(72-24ν2+ν4)1+(24-72c3(ν))cos(ν)(ν2-4)1+6(1-2c3(ν))νsin(ν))×(6c3(ν)M)-1,a32(ν)=c3(ν)(-6ν+ν3+6sin(ν))2N(3ν4M(ν2-4+4cos(ν)+νsin(ν))),b2(ν)=(2c3(ν)-1)ν2(ν3-6ν+6sin(ν))2(3MN),b3(ν)=(-4+ν2+4cos(ν)+νsin(ν))(c3(ν)N),c2(ν)=(24-12ν2+ν4-24cos(ν))(-12ν2+2ν4+12νsin(ν)),a21(ν)=c2(ν)22,a31(ν)=c3(ν)22-a32(ν),
in which
(15)M=24-12ν2+ν4-24cos(ν),N=-24-12(c3(ν)-1)ν2+(2c3(ν)-1)ν4+24cos(ν)+12c3(ν)νsin(ν).
Assuming thatc3(ν)is a constant and letting
(16)b2(ν)c23(ν)+b3(ν)c3(ν)3=120+𝒪(ν2),
we obtainc3=4/5. It is easy to check that
(17)b3(ν)a32(ν)c2(ν)=1120-ν25040+⋯.
The limit of the above equation asν→0verifies the second equation of (7). Thus we obtain a new TDRK method of order five given by (14) withc3=4/5. We denote this method by TFTDRKV5.For small values of|ν|the above formula is subject to heavy cancelations and in that case the following series expansions should be used:
(18)b1(ν)=548+ν2896-11ν42257920-37ν61390878720+4759ν845565186867200+⋯,b2(ν)=928-3ν25488+529ν464538880+4507ν6208718737920-74446597ν8191453523919257600+⋯,b3(ν)=25336-25ν243904-515ν4154893312+1115ν6222633320448+2904025ν810210854609027072+⋯,a21(ν)=118+ν21512+23ν47620480-ν648898080-7153ν819222813209600+⋯,a31(ν)=-2125+13ν21750-277ν42205000+713ν6582120000-1269227ν8133491758400000+⋯,a32(ν)=42125-13ν21750+277ν42205000-713ν6582120000+1269227ν8133491758400000+⋯,c2(ν)=13+ν2504+ν4317520-47ν6586776960-629ν8961140660480+⋯.It can be seen that whenν→0, this new method reduces to the fifth-order TDRK method (9) of Chan and Tsai [24]. The local truncation error of the above method is given by
(19)L.T.E11=-h621600{y(6)(xn)-gy(xn,yn)y(4)(xn)111111111111+20gyy(xn,yn)f(xn,yn)y(3)(xn)111111111111+20gxy(xn,yn)y(3)(xn)}+𝒪(h7).


### 3.2. Stability and Phase Analysis of the New Method

In this section, we investigate the linear stability of the new method. We consider the following test equation:(20)y′=iλy,λ>0.
When applied to (20) the TDRK method (2) produces the difference equation
(21)yn+1=M(iθ,ν)yn,θ=λh,i2=-1,
where
(22)M(iθ,ν)=(1-θ2b(ν)T(I+θ2A(ν))-1e)+i(×(I+θ2A(ν))-1c(ν)θ(1-θ2)b(ν)Thhhh×(I+θ2A(ν))-1c(ν))
withIthe identity matrix ande=(1,…,1)T.Definition 1.
For the TDRK method (2) with stability functionM(iθ,ν), the region in theθ-νplane
(23)Ω:={(θ,ν):|M(iθ,ν)|≤1}
is called the imaginary stability region of the method. And any closed curve defined by|M(iθ,ν)|=1is a stability boundary of the method.The imaginary stability region of the TFTDRKV5 is plotted in Figure1.Figure 1
Stability region of method TFTDRKV5.Definition 2 (see [25]).
For the TDRK method (2) with stability functionM(iθ,ν), the quantities
(24)P~(θ,ν)=θ-arg(M(iθ,ν)),D(θ,ν)=1-|M(iθ,ν)|
are called thephase-lag (dispersion) and amplification factor error (dissipation), respectively. If
(25)P~(θ,ν)=cϕθq+1+𝒪(θq+3),D~(θ,ν)=cdθp+1+𝒪(θp+3),
the method is said to have phase-lag order  q  and dissipation order  p, respectively.Denoting the ratior=ν/θ=ω/λ, we obtain the following expressions for the phase-lag and the amplification of the method TFTDRKV5:
(26)P~(θ,rθ)=-r2-15040θ7+𝒪(θ9),D~(θ,rθ)=r2-15760θ8+𝒪(θ10).
Thus, the method TFTDRKV5 has a phase-lag of order six and a dissipation of order seven.


## 3.1. Construction of the New Method

If we require the TDRK method (10) to integrate exactlyexp(iωx); that is,
(11)exp(iν)=cos(ν)+isin(ν)=1+iν-k1ν2-ik2ν31.+k3ν4+ik4ν5-k5ν6,ν=ωh,
then we have
(12)cos(ν)=1-k1ν2+k3ν4-k5ν6,sin(ν)=ν-k2ν3+k4ν5,
where
(13)k1=b1(ν)+b2(ν)+b3(ν),k2=b2(ν)c2(ν)+b3(ν)c3(ν),k3=b2(ν)a21(ν)+b3(ν)(a31(ν)+a32(ν)),k4=b3(ν)a32(ν)c2(ν),k5=b3(ν)a32(ν)a21(ν).
Solving (4)–(6) fors=3, the simplifying conditions (8) and (12), we obtain the coefficients expressed in terms ofc3(ν)(14)b1(ν)=(6(ν2-4)+c3(ν)(72-24ν2+ν4)1+(24-72c3(ν))cos(ν)(ν2-4)1+6(1-2c3(ν))νsin(ν))×(6c3(ν)M)-1,a32(ν)=c3(ν)(-6ν+ν3+6sin(ν))2N(3ν4M(ν2-4+4cos(ν)+νsin(ν))),b2(ν)=(2c3(ν)-1)ν2(ν3-6ν+6sin(ν))2(3MN),b3(ν)=(-4+ν2+4cos(ν)+νsin(ν))(c3(ν)N),c2(ν)=(24-12ν2+ν4-24cos(ν))(-12ν2+2ν4+12νsin(ν)),a21(ν)=c2(ν)22,a31(ν)=c3(ν)22-a32(ν),
in which
(15)M=24-12ν2+ν4-24cos(ν),N=-24-12(c3(ν)-1)ν2+(2c3(ν)-1)ν4+24cos(ν)+12c3(ν)νsin(ν).
Assuming thatc3(ν)is a constant and letting
(16)b2(ν)c23(ν)+b3(ν)c3(ν)3=120+𝒪(ν2),
we obtainc3=4/5. It is easy to check that
(17)b3(ν)a32(ν)c2(ν)=1120-ν25040+⋯.
The limit of the above equation asν→0verifies the second equation of (7). Thus we obtain a new TDRK method of order five given by (14) withc3=4/5. We denote this method by TFTDRKV5.For small values of|ν|the above formula is subject to heavy cancelations and in that case the following series expansions should be used:
(18)b1(ν)=548+ν2896-11ν42257920-37ν61390878720+4759ν845565186867200+⋯,b2(ν)=928-3ν25488+529ν464538880+4507ν6208718737920-74446597ν8191453523919257600+⋯,b3(ν)=25336-25ν243904-515ν4154893312+1115ν6222633320448+2904025ν810210854609027072+⋯,a21(ν)=118+ν21512+23ν47620480-ν648898080-7153ν819222813209600+⋯,a31(ν)=-2125+13ν21750-277ν42205000+713ν6582120000-1269227ν8133491758400000+⋯,a32(ν)=42125-13ν21750+277ν42205000-713ν6582120000+1269227ν8133491758400000+⋯,c2(ν)=13+ν2504+ν4317520-47ν6586776960-629ν8961140660480+⋯.It can be seen that whenν→0, this new method reduces to the fifth-order TDRK method (9) of Chan and Tsai [24]. The local truncation error of the above method is given by
(19)L.T.E11=-h621600{y(6)(xn)-gy(xn,yn)y(4)(xn)111111111111+20gyy(xn,yn)f(xn,yn)y(3)(xn)111111111111+20gxy(xn,yn)y(3)(xn)}+𝒪(h7).


## 3.2. Stability and Phase Analysis of the New Method

In this section, we investigate the linear stability of the new method. We consider the following test equation:(20)y′=iλy,λ>0.
When applied to (20) the TDRK method (2) produces the difference equation
(21)yn+1=M(iθ,ν)yn,θ=λh,i2=-1,
where
(22)M(iθ,ν)=(1-θ2b(ν)T(I+θ2A(ν))-1e)+i(×(I+θ2A(ν))-1c(ν)θ(1-θ2)b(ν)Thhhh×(I+θ2A(ν))-1c(ν))
withIthe identity matrix ande=(1,…,1)T.Definition 1.
For the TDRK method (2) with stability functionM(iθ,ν), the region in theθ-νplane
(23)Ω:={(θ,ν):|M(iθ,ν)|≤1}
is called the imaginary stability region of the method. And any closed curve defined by|M(iθ,ν)|=1is a stability boundary of the method.The imaginary stability region of the TFTDRKV5 is plotted in Figure1.Figure 1
Stability region of method TFTDRKV5.Definition 2 (see [25]).
For the TDRK method (2) with stability functionM(iθ,ν), the quantities
(24)P~(θ,ν)=θ-arg(M(iθ,ν)),D(θ,ν)=1-|M(iθ,ν)|
are called thephase-lag (dispersion) and amplification factor error (dissipation), respectively. If
(25)P~(θ,ν)=cϕθq+1+𝒪(θq+3),D~(θ,ν)=cdθp+1+𝒪(θp+3),
the method is said to have phase-lag order  q  and dissipation order  p, respectively.Denoting the ratior=ν/θ=ω/λ, we obtain the following expressions for the phase-lag and the amplification of the method TFTDRKV5:
(26)P~(θ,rθ)=-r2-15040θ7+𝒪(θ9),D~(θ,rθ)=r2-15760θ8+𝒪(θ10).
Thus, the method TFTDRKV5 has a phase-lag of order six and a dissipation of order seven.


## 4. Numerical Results

In this section we carry out six numerical results to illustrate the performance of the new method. The criterion used in the numerical comparison is the decimal logarithm of the maximum global error (LOG10(ERROR)) versus the computational effort measured by the CPU seconds (CPU SECONDS) required by each problem. The methods used for comparison are listed in the following.(i)
EFRK4: the fourth-order exponentially fitted Runge-Kutta method given by Vanden Berghe et al. in [23].(ii)
RK4V: the fourth-order optimized Runge-Kutta method given by Van de Vyver in [22].(iii)
RK5S: the fifth-order trigonometrically fitted Runge-Kutta method derived by Anastassi and Simos in [26].(iv)
RK4S: the fourth-order exponentially fitted Runge-Kutta method derived by Simos in [21].(v)
EFRK5: the fifth-order trigonometrically fitted Runge-Kutta method derived by Sakas and Simos in [27].(vi)
TFTDRKV5: the fifth-order TFTDRK method with one evaluation of functionfand three evaluations of functiongper step derived in Section 3.1 of this paper.Problem 1.
Consider the numerical integration of the Schrödinger equation(27)y′′(x)=(v(x)-E)y(x),
with the well-known Woods-Saxon potential
(28)v(x)=c0z(1-a(1-z)),
wherez=(exp(a(x-b)+1))-1,c0=-50,a=5/3, andb=7.The domain of numerical integration is[0,15]. It is appropriate to chooseωas follows (see [10]):
(29)ω={50+E,x∈[0,6.5],E,x∈[6.5,15].
In this experiment we consider the resonance problem (E>0). The numerical results are compared with the analytical solution of the Woods-Saxon potential, rounded to six decimal places. In Figures 2, 3, 4, and 5, we plot the error of|Eanalytical-Ecalculated|versus the computational effort measured by CPU seconds required by each method forEanalytical=53.588872,163.215341,341.495874,989.701916, respectively.Figure 2
Efficiency curves for Problem1:E=53.588872.Figure 3
Efficiency curves for Problem1:E=163.215341.Figure 4
Efficiency curves for Problem1:E=341.495874.Figure 5
Efficiency curves for Problem1:E=989.701916.Problem 2.
We consider the inhomogeneous equation(30)y′′(x)+100y(x)=99sin(x),y(0)=1,y′(0)=11,
whose exact solution is
(31)y(x)=cos(10x)+sin(10x)+sin(x).
In this test we chooseω=10. The numerical results given in Figure 6 are computed with the stepsizesh=2-i-4,i=1,…,4,  for all the methods on the interval[0,1000].Figure 6
Efficiency curves for Problem2.Problem 3.
We consider the linear problem studied in [9]
(32)y′′+10000y=(10000-4x2)cos(x2)-2sin(x2),y(0)=1,y′(0)=100,
whose exact solution is
(33)y(x)=sin(100x)+cos(x2).
In this test we takeω=100. The numerical results given in Figure 7 are computed with the stepsizesh=2-9-i,i=1,…,4,  for all the methods on the interval[0,100].Figure 7
Efficiency curves for Problem3.Problem 4.
We consider the following “almost periodic” orbit problem studied by Franco and Palacios [28]:
(34)z′′+z=ϵeiψx,z(0)=1,z′(0)=i,z∈ℂ,
or equivalently
(35)u′′+u=ϵcos(ψx),u(0)=1,u′(0)=0,v′′+v=ϵsin(ψx),v(0)=0,v′(0)=1,
whereϵ=0.001andψ=0.01.The analytic solution to the problem is given by
(36)u(x)=1-ϵ-ψ21-ψ2cos(x)+ϵ1-ψ2cos(ψx),v(x)=1-ϵψ-ψ21-ψ2sin(x)+ϵ1-ψ2sin(ψx).
In this test we takeω=1. The numerical results given in Figure 8 are computed with the stepsizesh=2-i,i=1,…,4,  for all the methods on the interval[0,1000].Figure 8
Efficiency curves for Problem4.Problem 5.
We consider the two-dimensional problem studied in [29]
(37)y′′(t)+(13-12-1213)y(t)=(f1(t)f2(t)),y(0)=(10),y′(0)=(-48),
where
(38)f1(t)=9cos(2t)-12sin(2t),f2(t)=-12cos(2t)+9sin(2t).
The analytical solution is
(39)y(t)=(sin(t)-sin(5t)+cos(2t)sin(t)-sin(5t)+sin(2t)).
In this test we chooseω=5. The numerical results given in Figure 9 are computed with the stepsizesh=2-2-i,i=1,…,5,  for all the methods on the interval[0,100].Figure 9
Efficiency curves for Problem5.Problem 6.
We consider the Fermi-Pasta-Ulam problem studied in [30]. This problem consists of a chain of springs, where nonlinear springs alternate with stiff harmonic springs. The problem can be described by a Hamiltonian system with Hamiltonian function
(40)H=12∑i=1N‍(ui′2+vi′2)+ω22∑i=1N‍vi2+14∑i=1N‍(ui+1-vi+1-ui-vi)4,
whereu0=v0=uN+1=vN+1=0. This is equivalent to the nonlinear oscillatory problem of the form
(41)ui′′=-∂Hp∂ui,i=1,…,N,vi′′+ω2vi=-∂Hp∂vi,i=1,…,N,
where
(42)Hp=14∑i=1N‍(ui+1-vi+1-ui-vi)4.
Since the analytic solution is not available, we study the numerical total energyHon the interval[0,tend]. In this test we chooseN=3withω=100andtend=10and we take the initial values
(43)u1(0)=1,u1′(0)=1,v1=1ω,v1′(0)=1,
the remaining initial values being zero.
The numerical results given in Figure10 are computed with the stepsizes h=2-8-i,i=1,…,4, for all the methods.
In Figures2–10, we can see that the new method TFTDTRK5V is the most efficient among the six methods we select for comparison.Figure 10
Efficiency curves for Problem6.


## 5. Conclusions

A new trigonometrically fitted two-derivative Runge-Kutta method of algebraic order five is derived in this paper. We also analyze the linear stability and phase properties of the new method. Numerical results are reported to show the robustness and competence of the new method compared with some highly efficient methods in the recent literature. To explain the superiority of the new method, we observe that, compared with RK methods of a specific number of stages, TDRK methods have the possibility of gaining a higher algebraic order than traditional RK methods with the same number of stages. We conclude that, for problems with oscillatory solutions, trigonometrically fitted TDRK methods are more accurate and more efficient than the adapted (exponentially or trigonometrically fitted) RK methods.



---

*Source: 937858-2013-12-03.xml*