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

An explicit Runge-Kutta-Nyström method is developed for solving second-order differential equations of the form q′′ 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.


Introduction
This paper deals with numerical method for second-order ODEs, in which the derivative does not appear explicitly, q f t, q , q t 0 q 0 , q t 0 q 0 , 1.1 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 zerodissipation 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.

General Theory
An explicit Runge-Kutta-Nystr öm RKN method for the numerical integration of the initial value problem IVP is given by where The .

2.3
Now, the phase-lag analysis of the method 2.1 is investigated using the homogeneous test equation see 6 By applying the general method 2.1 to the test equation 2.4 we obtain the following recursive relation: where A, A , B, and B are polynomials in z 2 , completely determined by the parameters of the method 2.1 .The characteristic equation for 2.5 is given by The exact solution of 2.4 is given by where is the length of the vector σ 1,2 .Substituting in 2.7 , we have q t n 2|σ| cos χ nz .

2.8
Now let us assume that the eigenvalues of D are ρ 1 , ρ 2 and the corresponding eigenvectors are 1, The numerical solution of 2.5 is where

2.14
Following 6 , for a zero-dissipative RKN method and since R z 2 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 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 where

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
All indexes are from 1 to m.Also the Nyström row sum conditions need to be satisfied In addition, we try to minimize the following norms of truncation error which is discussed by Dormand 9 : where

3.4
The algebraic conditions 3.1 -3.2 together with dispersion relation 2.15 and zerodissipation condition S H ≡ 1 i.e., β 3 0 formed a system of nonlinear equations.
Letting c 3 as the free parameter and solving all nonlinear equations simultaneously yield the following one-parameter family of solution in term of c 3 : where The above set of solution will give τ 4 2 2.635231381 × 10 −2 and τ 4 2 can be written in terms of c 3 .By using minimize command in MAPLE, τ 4 2 has a minimum value zero at c 3 1.047071398.Since RKN is a one-step method, we only consider the case c 3 ∈ 0, 1 where τ 4 2 lies between 1.872175321 × 10 −2 , 5.796540889 × 10 −2 .We developed a equidistant nodes code where c 3 0.0 0.01i, i 0, 1, 2, . . . .The method gives an accurate results for most of the problems tested when c 3  1, which will also gives τ 4 2 1.872175321 × 10 −2 and the dispersion constant is − 1/40320 z 7 O z 9 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 .
In Table 1 where c 2 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

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.Figures 1, 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 .
To illustrate some properties of zero-dissipative with minimal phase-lag integrator, the following physical problems are tested.
Problem 1 Harmonic oscillator .One has q t −ω 2 q t , q 0 q 0 , q 0 q 0 , t ∈ 0, t end .4.1 The logarithm error of energy ERR at each integration point when solving the harmonic oscillator for ω 1 with initial condition q 0 1, q 0 0 and Δt 1/5.The error at each integration point when solving the harmonic oscillator with ω 8 with initial condition q 0 1, q 0 −2 and Δt 1/20.
Exact solution: q t c 1 sin ωt c 2 cos ωt Total energy as given in 1 , where a depends on the initial conditions.Exact solution:

4.5
Exact solution q 1 t cos t 0.0005t sin t , q 2 t sin t − 0.0005t cos t where 0.001 and ψ 0.01.The analytical solution z t u t iv t is given by Problem 4 Problem studied by van der Houwen and Sommeijer 6 .One has where v 1, t ∈ 0, t end .Exact solution is q 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 The exact solution computed by the Galerkin method and given by q t 4 i 0 a 2i 1 cos 1.01 2i 1 t , 4.12 with a 1 0.200179477536, a 3 0.246946143 × 10 −3 , a 5 0.304014 × 10 −6 , a 7 0.374 × 10 −9 , and a 9 < 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.Figures 8-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 .

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.
RKN parameters a ij , b j , b j , and c j 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 c 1 , c 2 , . . ., c m T , b b 1 , b 2 , . . ., b m T , b b 1 , b 2 , . . ., b m T , and A a ij ,respectively.An explicit of RKN methods of order r can be expressed in Butcher notation by the table of coefficients

Figure 1 :
Figure1: Energy Conservation .The logarithm error of energy ERR at each integration point when solving the harmonic oscillator for ω 1 with initial condition q 0 1, q 0 0 and Δt 1/5.

Figure 2 :
Figure 2:The error at each integration point when solving the harmonic oscillator with ω 8 with initial condition q 0 1, q 0 −2 and Δt 1/20.

Figure 3 :Figure 4 :
Figure 3: The global error at each integration point when solving the "almost" periodic problem Problem 2 with Δt 1/40.

Figure 5 :Figure 6 :Problem 3
Figure 5: The global error at each integration point when solving the inhomogeneous problem with Δt 1/10.
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, H p if |R H | < 2 and S H ≡ 1 for all H ∈ 0, H p .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, H p is called interval of periodicity of the method 2.1 if, for all H ∈ 0, H p , |ξ 1,2 | 1 and ξ 1 / ξ 2 .

Table 2 :
Summary of the characteristics of the RKN methods.
Note: u means dispersion order, v means dissipation order, D.C means dissipation constant, S.I means stability interval, P.I means periodicity interval.