# A Zero-Dissipative Runge-Kutta-Nyström Method with Minimal Phase-Lag

Authors: Norazak Senu; Mohamed Suleiman; Fudziah Ismail; Mohamed Othman
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

An explicit Runge-Kutta-Nyström method is developed for solving second-order
differential equations of the formq′′=f(t,q) where the solutions are oscillatory. The method has zero-dissipation with minimal phase-lag at a cost of three-function
evaluations per step of integration. Numerical comparisons with RKN3HS, RKN3V, RKN4G, and RKN4C methods show the preciseness and effectiveness of the method developed.


---

## Body

## 1. Introduction

This paper deals with numerical method for second-order ODEs, in which the derivative does not appear explicitly,(1.1)q′′=f(t,q),q(t0)=q0,q′(t0)=q0′,
for which it is known in advance that their solutions are oscillating. Such problems often arise in different areas of engineering and applied sciences such as celestial mechanics, quantum mechanics, elastodynamics, theoretical physics, chemistry, and electronics (see, e.g., [1–3]).For ODEs of type (1.1), it is preferable to use a direct numerical method, instead of reducing it into first-order system. One such method is Runge-Kutta-Nyström (RKN) method. Oscillatory problems (1.1) are usually considered as difficult to handle. Many methods with constant coefficients have already been derived for solving (1.1), see [2, 4–6]. Explicit RKN methods which relate to dispersion (or phase-lag) and dissipation (or amplification error) was first introduced by van der Houwen and Sommeijer [6], these properties are useful when dealing with periodic problems rather than just minimizing the truncation error of the methods. The objective of the study here is to construct RKN method with zero-dissipation and minimal phase-lag to ensure that the problem will be integrated as precisely as possible and the approximate solutions remain in phase especially for the harmonic oscillatory problem. To our knowledge, method of algebraic order greater than three with zero-dissipation and minimal phase-lag has not been done yet so far. Therefore, this motivates us to start the investigation with method of algebraic order three with zero-dissipation and minimal phase-lag.When solving (1.1) numerically, attention has to be given to the algebraic order of the method used, since this is the main criterion for achieving high accuracy for long-range integration. Therefore, it is desirable to have a lower stage RKN method with maximal order. This will reduced the computational cost. Moreover, if it is initially known that the solution of (1.1) is periodic in nature then it is essential to consider phase-lag and dissipation. These are actually two types of truncation errors beside the truncation error due to the algebraic order. The first is the angle between the true and the approximate solutions, while the second is the distance from a standard cyclic solution. Essentially, the method derived should ideally have the properties of zero dissipation, minimal phase-lag and small truncation error coefficients. In this paper, we developed an efficient a zero-dissipative explicit RKN method with minimal phase-lag in constant time step mode.


## 2. General Theory

An explicit Runge-Kutta-Nyström (RKN) method for the numerical integration of the initial value problem (IVP) is given by(2.1)qn+1=qn+hqn′+h2∑i=1mbif(tn+cih,Qi),qn+1′=qn′+h∑i=1mbi′f(tn+cih,Qi),
where
(2.2)Qi=qn+cihqi′+h2∑j=1i-1aijf(tn+cjh,Qj).The RKN parametersaij,bj,bj′, and cj are assumed to be real and m is the number of stages of the method. Introduce the m-dimensional vectors c, b, and b′ and m×m matrix A, where c=[c1,c2,…,cm]T, b=[b1,b2,…,bm]T, b′=[b1′,b2′,…,bm′]T, and A=[aij], respectively. An explicit of RKN methods of order r can be expressed in Butcher notation by the table of coefficients   (2.3)cAbTb′T.Now, the phase-lag analysis of the method (2.1) is investigated using the homogeneous test equation (see [6])(2.4)q′′=(iλ)2q(t),i=-1.By applying the general method (2.1) to the test equation (2.4) we obtain the following recursive relation:(2.5)[qn+1hqn+1′]=D[qnhqn′],D=[A(z2)B(z2)A′(z2)B′(z2)],z=λh,
where A,A′,B, and B′ are polynomials in z2, completely determined by the parameters of the method (2.1). The characteristic equation for (2.5) is given by(2.6)ξ2-tr(D(z2))ξ+det(D(z2))=0.
The exact solution of (2.4) is given by(2.7)q(tn)=σ1[exp(iλ)]n+σ2[exp(-iλ)]n,
where σ1,2=(1/2)[q0±(iq0′)/λ] or σ1,2=|σ|exp(±iχ) where |σ|=(Re(σ1,2))2+(Im(σ1,2))2 is the length of the vector σ1,2. Substituting in (2.7), we have(2.8)q(tn)=2|σ|cos(χ+nz).Now let us assume that the eigenvalues ofD are ρ1, ρ2 and the corresponding eigenvectors are [1,ν1]T, [1,ν2]T, νi=A′/(ρi-B′), i=1,2. The numerical solution of (2.5) is(2.9)qn=c1ρ1n+c2ρ2n,
where
(2.10)c1=-ν2q0-hq0′ν1-ν2,c2=ν1q0-hq0′ν1-ν2.Ifρ1,ρ2 are complex conjugates, then c1,2=|c|exp(±iω) and ρ1,2=|ρ|exp(±iφ) where |c| and |ρ| is the length of the vector c1,2 and ρ1,2, respectively. By substituting in (2.9), we have(2.11)qn=2|c||ρ|ncos(ω+nφ).Equations (2.8) and (2.11) led us to the following definition.Definition 2.1.
For the RKN method corresponding to the characteristic equation (2.6) the quantities
(2.12)ϕ(z)=z-φ,α(z)=1-|ρ|
are the phase-lag (or dispersion) and dissipation (or amplification error), respectively. If ϕ(z)=O(zu+1), then the RKN method is said to have phase-lag order u and if α(z)=O(zv+1), then the RKN method is said to have dissipation order v. If at a point z, α(z)=0, then the RKN method has zero dissipation.From Definition2.1 it follows that(2.13)ϕ(z)=z-cos-1(R(z2)2S(z2)),α(z)=1-S(z2),
where R(z2) and S(z2) defined by(2.14)R(z2)=tr(D)=2+∑i=1m(-1)iαiz2i,S(z2)=det(D)=1+∑i=1m(-1)iβiz2i.Following [6], for a zero-dissipative RKN method and since R(z2) is at most degree 2m then the maximal attainable order of dispersion is u=2m. The method with maximum order of the phase-lag is said to be a method with minimal phase-lag. The dispersion relations up to order four are already satisfied due to consistency condition of the method. Altogether, the condition to have phase-lag order six, the highest possible for a three-stage zero-dissipative explicit method, is(2.15)α3=1360.We next discuss the stability properties of method for solving (1.1) by considering once again the expression (2.5). The matrix D in (2.5) can be written as(2.16)D(H)=(1-HbT(I+HA)-1e1-HbT(I+HA)-1c-Hb'T(I+HA)-1e1-Hb'T(I+HA)-1c),
where H=z2, e=[1⋯1]T, c=[c1⋯cm]T. The D(H) is called the stability matrix. Following [7], we say that the numerical method (2.1) has interval of periodicity or interval of zero dissipation (0,Hp) if |R(H)|<2 and S(H)≡1 for all H∈(0,Hp). Notice that S(H)≡1 is a necessary condition for the existence of a nonempty periodicity interval. Method with zero-dissipation is also known as method with dissipation of order infinity.Definition 2.2.
An interval(0,Hp) is called interval of periodicity of the method (2.1) if, for all H∈(0,Hp), |ξ1,2|=1 and ξ1≠ξ2.


## 3. Construction of the Method

In the following we shall derive zero-dissipative minimal phase-lag RKN method. The RKN parameters must satisfy the following algebraic conditions as given in [8, 9]:(3.1)order1:b'Te=1,order2:bTe=12,b'Tc=12,order3:bTc=16,12b'Tc2=16.All indexes are from 1 tom. Also the Nyström row sum conditions need to be satisfied(3.2)12ci2=∑j=1maij(i=1,…,m).In addition, we try to minimize the following norms of truncation error which is discussed by Dormand [9]:(3.3)∥τ(4)∥2=(∑i=12(τi(4))2)1/2,∥τ'(4)∥2=(∑i=13(τi'(4))2)1/2,
where τi(4) and τi'(4) are defined as in [2],(3.4)τ1(4)=12bTc2-124,τ2(4)=bTAe-124,τ1'(4)=16bTc3-124,τ2'(4)=b'T(c.Ae)-18,τ3'(4)=b'TAc-124.The algebraic conditions (3.1)-(3.2) together with dispersion relation (2.15) and zero-dissipation condition S(H)≡1 (i.e., β3=0) formed a system of nonlinear equations. Letting c3 as the free parameter and solving all nonlinear equations simultaneously yield the following one-parameter family of solution in term of c3:(3.5)c2=130(-F+3(5+5)),a21=11800(-15-35+F)2,a31=-1160c3(F(5+3)(1-5+2c3)+45-20-20c3-12c35),a32=1160c3(F(5+3)(1-5+2c3)+60c3-12c35+45-20),b1=-[F(5-5)(10+25-45c3-9c35+30c32)+1800c32-900c3-180c35]/[240c3(15+35-F-30c3)],b2=-F(5-5)(-5-5+10c3)+605+300-600c380(-15-35+F+30c3),b3=-F6c3(15+35-F-30c3),b1′=[F(5-5)(5+5-2c35-20c3+15c32)+300c32+605-100-120c35]/[40c3(-15-35+F+30c3)],b2′=-F(5-5)(-2+3c3)-180c3+1208(-15-35+F+30c3),b3′=-5+35-F2c3(15+35-F-30c3),
where F=30(-1+5).The above set of solution will give∥τ(4)∥2=2.635231381×10-2 and ∥τ′(4)∥2 can be written in terms of c3. By using minimize command in MAPLE, ∥τ′(4)∥2 has a minimum value zero at c3=1.047071398. Since RKN is a one-step method, we only consider the case c3∈[0,1] where ∥τ′(4)∥2 lies between [1.872175321×10-2,5.796540889×10-2]. We developed a equidistant nodes code where c3=0.0+0.01i,i=0,1,2,…. The method gives an accurate results for most of the problems tested when c3=1, which will also gives ∥τ′(4)∥2=1.872175321×10-2 and the dispersion constant is -(1/40320)z7+O(z9) with a periodicity interval of approximately (0,7.571916). The new method which is denoted by RKN3NEW can be represented by the following array (see Table 1).Table 1
The RKN3NEW method.In Table1 where c2=0.520623384839812. The coefficients of RKN3NEW is generated using computer algebra package MAPLE whose environment variable Digits which control the number of significant digits which is set to 15. Table 2 shows the comparisons of the characteristics of the method derived against the methods due to Van de Vyver [2], RKN3V Calvo and Sanz-Serna [10], RKN4C van der Houwen and Sommeijer [6], RKN3HS and García et al. [11], RKN4G.Table 2
Summary of the characteristics of the RKN methods.
Methoduv∥τ(r+1)∥2∥τ′(r+1)∥2D.CS.I/P.IRKN3NEW6∞2.64×10-21.87×10-2-(0,7.57)RKN3V6∞2.70×10-39.01×10-4-(0,7.57)RKN4C4∞4.01×10-43.10×10-47.37×10-5(0,9.33)RKN3HS1231.12×10-21.59×10-21.83×10-2(0,9.42)RKN4G456.94×10-43.13×10-31/1440(0,8.77)Note:u means dispersion order, v means dissipation order, D.C means dissipation constant, S.I means stability interval, P.I means periodicity interval.Please note that the RKN3NEW, RKN3V, RKN4G, and RKN3HS methods have three function evaluations per step. The RKN4C method has five stages with FSAL technique applied meaning that it has four effective stages per step.


## 4. Problem Tested

In order to evaluate the preciseness and effectiveness of the new RKN3NEW method, we solved several physical problems which have oscillatory solutions using constant time step. The RKN3NEW algorithm is coded and executed on Intel Pentium IV processor with double precision arithmetic. Figures1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13 show the numerical results for all methods used. These codes have been denoted by the following.(i)
RKN3NEW: The new zero-dissipative three-stage third-order RKN method with minimal phase-lag derived in this paper.(ii)
RKN3V: The zero-dissipative three-stage third-order RKN method given in Van de Vyver [2].(iii)
RKN4C: The zero-dissipative five-stage fourth-order RKN method with FSAL technique given in Calvo and Sanz-Serna [10].(iv)
RKN3HS: The dissipative three-stage third-order RKN method given in van der Houwen and Sommeijer [6].(v)
RKN4G: The dissipative three-stage fourth-order RKN method derived by García et al. [11].Figure 1
(Energy Conservation). The logarithm error of energy (ERR) at each integration point when solving the harmonic oscillator forω=1 with initial condition q0=1, q0′=0 and Δt=1/5.Figure 2
The error at each integration point when solving the harmonic oscillator withω=8 with initial condition q0=1, q0′=-2 and Δt=1/20.Figure 3
The global error at each integration point when solving the “almost” periodic problem (Problem2) with Δt=1/40.Figure 4
The global error at each integration point when solving the “almost” periodic problem (Problem3) with Δt=1/2.Figure 5
The global error at each integration point when solving the inhomogeneous problem withΔt=1/10.Figure 6
The error at each integration point when solving the inhomogeneous system withΔt=1/20.Figure 7
The error at each integration point when solving the Duffing's equation withΔt=1/2.Figure 8
The efficiency curve for Problem1 with ω=5, tend=105 and time step Δt=0.1/2i, i=0,…,4.Figure 9
The efficiency curve for Problem2 with tend=2·105 and time step Δt=1/(5·2i), i=0,…,4.Figure 10
The efficiency curve for Problem3 with tend=105 and time step Δt=0.8/2i, i=0,…,4.Figure 11
The efficiency curve for Problem4 with tend=105 and time step Δt=0.1/2i,i=1,…,5.Figure 12
The efficiency curve for Problem5 with tend=5·104 and time step Δt=0.03-0.05i, i=0,…,4.Figure 13
The efficiency curve for Duffing's equation withtend=105 and time step Δt=1/2i,i=1,…,5.To illustrate some properties of zero-dissipative with minimal phase-lag integrator, the following physical problems are tested.Problem 1 (Harmonic oscillator).
One has(4.1)q′′(t)=-ω2q(t),q(0)=q0,q′(0)=q0′,t∈[0,tend].
Exact solution:q(t)=c¯1sin(ωt)+c¯2cos(ωt)
Total energy as given in [1],(4.2)E(q,q′)=q22+q'22=a22,
where a depends on the initial conditions.Problem 2 (An “almost” Periodic Orbit problem studied by Stiefel and Bettis [12]).
One has(4.3)z′′+z=0.001eit,z(0)=1,z′(0)=0.9995i,t∈[0,tend].
Exact solution:(4.4)z(t)=(1-0.0005it)eit.
We write in the equivalent form(4.5)q1′′+q1(t)=0.001cos(t),q1(0)=1,q1′(0)=0,q2′′+q2(t)=0.001sin(t),q2(0)=0,q2′(0)=0.9995.
Exact solution q1(t)=cos(t)+0.0005tsin(t),  q2(t)=sin(t)-0.0005tcos(t)Problem 3 (An “almost” Periodic Orbit problem studied by Franco and Palacios [13]).
One has(4.6)z′′+z=ϵeiψt,z(0)=1,z′(0)=i,t∈[0,tend],
where ϵ=0.001 and ψ=0.01. The analytical solution z(t)=u(t)+iv(t) is given by
(4.7)u(t)=1-ϵ-ψ21-ψ2cos(t)+ϵ1-ψ2cos(ψt),v(t)=1-ϵψ-ψ21-ψ2sin(t)+ϵ1-ψ2sin(ψt).Problem 4 (Problem studied by van der Houwen and Sommeijer [6]).
One has(4.8)q′′(t)=-v2q(t)+(v2-1)sin(t),q(0)=1,q′(0)=v+1,
where v≫1,t∈[0,tend].
Exact solution isq(t)=cos(vt)+sin(vt)+sin(t). Numerical result is for the case v=10.Problem 5 (Inhomogeneous system studied by Franco [5]).
One has(4.9)q′′(t)+(1012-992-9921012)q(t)=ɛ(932cos(2t)-992sin(2t)932sin(2t)-992cos(2t)),q(0)=(-1+ɛ1),q′(0)=(-1010+2ɛ),t∈[0,tend].
Exact solution:(4.10)q(t)=(-cos(10t)-sin(10t)+ɛcos(2t)cos(10t)+sin(10t)+ɛsin(2t)).Problem 6 (The undamped Duffing's equations as given in [14]).
One has(4.11)q′′(t)+q(t)+(q(t))3=0.002cos(1.01t),q(0)=0.200426728067,q′(0)=0,t∈[0,tend].
The exact solution computed by the Galerkin method and given by
(4.12)q(t)=∑i=04a2i+1cos[1.01(2i+1)t],
with a1=0.200179477536, a3=0.246946143×10-3, a5=0.304014×10-6, a7=0.374×10-9, and a9<10-12.The results show the typical properties of zero-dissipative with minimal phase-lag integrator, RKN3NEW which we have derived. We compare the method with the dissipative method of high order of dispersion, RKN3HS [6]. Figure 1 shows the error of energy at each integration point. It can be seen that the zero-dissipative with minimal phase-lag algorithm conserved the energy with energy error one order magnitude smaller than the energy error for dissipative algorithm. In Figures 2–7, we plotted the global error versus the time of integration for different time step, Δt for various physical problems. From the figures we observed that the global error produced by the RKN3NEW method do not increased over time. This means that the approximate solutions remain in phase over a long-range of integration. Clearly, the global error propagated faster over the time for dissipative RKN3HS method. Next we study the global error growth and the efficiency of the method over a long period of integration.Figures8–13 presented the efficiency of the method developed by plotting the graph of the decimal logarithm of the maximum global error against the logarithm number of function evaluations for long periods of computations. The RKN3NEW is significantly more efficient than the RKN4C method and the dissipative RKN3HS and RKN4G methods for the homogeneous and nonhomogeneous problems. It is also found that the RKN3NEW is competitive when compared with the symplectic RKN3V method. This validate the fact that the zero-dissipation with minimal phase-lag is an important element when solving ODEs in which the solutions are in oscillating form. The dissipative methods RKN3HS and RKN4G are less accurate and hence the global error accumulated faster when the integration is done over a longer interval. For the nonlinear oscillator problem, it is necessary to compare the method with another method of the same algebraic order because the global error is dominated by the truncation error due to algebraic order rather than the dispersion error, see [15].


## 5. Conclusion

In this paper we have derived an explicit zero-dissipative RKN3NEW method with minimal phase-lag which is suitable for problems with oscillating solutions especially for long-range integration. From the numerical results, we can conclude that the new RKN3NEW method is more promising compared to dissipative method RKN3HS, RKN4G and with symplectic method, RKN4C and as competitive when compared with symplectic RKN3V method. One can obtain higher-order accuracy by extending the idea given in this paper.



---

*Source: 591341-2010-06-27.xml*