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

Abstract

An explicit Runge-Kutta-Nyström method is developed for solving second-order differential equations of the form 𝑞^′′=𝑓(𝑡,𝑞) 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.

1. Introduction

This paper deals with numerical method for second-order ODEs, in which the derivative does not appear explicitly,𝑞^′′=𝑓(𝑡,𝑞),𝑞(𝑡₀)=𝑞₀,𝑞^′(𝑡₀)=𝑞^′₀, (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 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 𝑞_𝑛+1=𝑞_𝑛+ℎ𝑞^′_𝑛+ℎ^2𝑚_𝑖=1𝑏_𝑖𝑓(𝑡_𝑛+𝑐_𝑖ℎ,𝑄_𝑖𝑞),_′𝑛+1=𝑞^′_𝑛+ℎ^𝑚_𝑖=1𝑏^′_𝑖𝑓(𝑡_𝑛+𝑐_𝑖ℎ,𝑄_𝑖), (2.1) where 𝑄_𝑖=𝑞_𝑛+𝑐_𝑖ℎ𝑞^′_𝑖+ℎ^2𝑖−1_𝑗=1𝑎_𝑖𝑗𝑓(𝑡_𝑛+𝑐_𝑗ℎ,𝑄_𝑗). (2.2)

The RKN parameters 𝑎_𝑖𝑗,𝑏_𝑗,𝑏^′_𝑗, and 𝑐_𝑗 are assumed to be real and 𝑚 is the number of stages of the method. Introduce the 𝑚-dimensional vectors 𝑐, 𝑏, and 𝑏^′ and 𝑚×𝑚 matrix 𝐴, where 𝑐=[𝑐₁,𝑐₂,…,𝑐_𝑚]^𝑇, 𝑏=[𝑏₁,𝑏₂,…,𝑏_𝑚]^𝑇, 𝑏^′=[𝑏^′₁,𝑏^′₂,…,𝑏^′_𝑚]^𝑇, and 𝐴=[𝑎_𝑖𝑗], respectively. An explicit of RKN methods of order 𝑟 can be expressed in Butcher notation by the table of coefficients 𝑐𝐴𝑏^𝑇𝑏^′𝑇. (2.3)

Now, the phase-lag analysis of the method (2.1) is investigated using the homogeneous test equation (see [6])𝑞^′′=(𝑖𝜆)²𝑞(𝑡),𝑖=√−1. (2.4)

By applying the general method (2.1) to the test equation (2.4) we obtain the following recursive relation:𝑞_𝑛+1ℎ𝑞^′_𝑛+1𝑞=𝐷^𝑛ℎ𝑞^′_𝑛,𝐷=𝐴(𝑧²)𝐵(𝑧²)𝐴^′(𝑧²)𝐵^′(𝑧²),𝑧=𝜆ℎ, (2.5) where 𝐴,𝐴^′,𝐵, and 𝐵^′ are polynomials in 𝑧², completely determined by the parameters of the method (2.1). The characteristic equation for (2.5) is given by𝜉²−tr(𝐷(𝑧²))𝜉+det(𝐷(𝑧²))=0. (2.6) The exact solution of (2.4) is given by𝑞(𝑡_𝑛)=𝜎₁[exp(𝑖𝜆)]^𝑛+𝜎₂[exp(−𝑖𝜆)]^𝑛, (2.7) where 𝜎_1,2=(1/2)[𝑞₀±(𝑖𝑞^′₀)/𝜆] or 𝜎_1,2=|𝜎|exp(±𝑖𝜒) where |𝜎|=(Re(𝜎_1,2))²+(Im(𝜎_1,2))² is the length of the vector 𝜎_1,2. Substituting in (2.7), we have𝑞(𝑡_𝑛)=2|𝜎|cos(𝜒+𝑛𝑧). (2.8)

Now let us assume that the eigenvalues of 𝐷 are 𝜌₁, 𝜌₂ and the corresponding eigenvectors are [1,𝜈₁]^𝑇, [1,𝜈₂]^𝑇, 𝜈_𝑖=𝐴^′/(𝜌_𝑖−𝐵^′), 𝑖=1,2. The numerical solution of (2.5) is𝑞_𝑛=𝑐₁𝜌^𝑛₁+𝑐₂𝜌^𝑛₂, (2.9) where 𝑐₁𝜈=−²𝑞₀−ℎ𝑞^′₀𝜈₁−𝜈₂,𝑐₂=𝜈₁𝑞₀−ℎ𝑞^′₀𝜈₁−𝜈₂. (2.10)

If 𝜌₁,𝜌₂ are complex conjugates, then 𝑐_1,2=|𝑐|exp(±𝑖𝜔) and 𝜌_1,2=|𝜌|exp(±𝑖𝜑) where |𝑐| and |𝜌| is the length of the vector 𝑐_1,2 and 𝜌_1,2, respectively. By substituting in (2.9), we have𝑞_𝑛=2|𝑐||𝜌|^𝑛cos(𝜔+𝑛𝜑). (2.11)

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 𝜙(𝑧)=𝑧−𝜑,𝛼(𝑧)=1−|𝜌| (2.12) are the phase-lag (or dispersion) and dissipation (or amplification error), respectively. If 𝜙(𝑧)=𝑂(𝑧^𝑢+1), then the RKN method is said to have phase-lag order 𝑢 and if 𝛼(𝑧)=𝑂(𝑧^𝑣+1), then the RKN method is said to have dissipation order 𝑣. If at a point 𝑧, 𝛼(𝑧)=0, then the RKN method has zero dissipation.

From Definition 2.1 it follows that𝜙(𝑧)=𝑧−cos⁻¹(𝑅(𝑧²)2𝑆(𝑧²)),𝛼(𝑧)=1−𝑆(𝑧²), (2.13) where 𝑅(𝑧²) and 𝑆(𝑧²) defined by𝑅(𝑧²)=tr(𝐷)=2+^𝑚_𝑖=1(−1)^𝑖𝛼_𝑖𝑧^2𝑖,𝑆(𝑧²)=det(𝐷)=1+^𝑚_𝑖=1(−1)^𝑖𝛽_𝑖𝑧^2𝑖. (2.14)

Following [6], for a zero-dissipative RKN method and since 𝑅(𝑧²) is at most degree 2𝑚 then the maximal attainable order of dispersion is 𝑢=2𝑚. 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𝛼₃=1.360 (2.15)

We next discuss the stability properties of method for solving (1.1) by considering once again the expression (2.5). The matrix 𝐷 in (2.5) can be written as𝐷(𝐻)=1−𝐻𝑏^𝑇(𝐼+𝐻𝐴)⁻¹𝑒1−𝐻𝑏^𝑇(𝐼+𝐻𝐴)⁻¹𝑐−𝐻𝑏^′𝑇(𝐼+𝐻𝐴)⁻¹𝑒1−𝐻𝑏^′𝑇(𝐼+𝐻𝐴)⁻¹𝑐, (2.16) where 𝐻=𝑧², 𝑒=[1⋯1]^𝑇, 𝑐=[𝑐₁⋯𝑐_𝑚]^𝑇. The 𝐷(𝐻) 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,𝐻_𝑝) if |𝑅(𝐻)|<2 and 𝑆(𝐻)≡1 for all 𝐻∈(0,𝐻_𝑝). Notice that 𝑆(𝐻)≡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,𝐻_𝑝) is called interval of periodicity of the method (2.1) if, for all 𝐻∈(0,𝐻_𝑝), |𝜉_1,2|=1 and 𝜉₁≠𝜉₂.

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]:order1∶𝑏^′𝑇𝑒=1,order2∶𝑏^𝑇1𝑒=2,𝑏^′𝑇1𝑐=2,order3∶𝑏^𝑇1𝑐=6,12𝑏^′𝑇𝑐²=16. (3.1)

All indexes are from 1 to 𝑚. Also the Nyström row sum conditions need to be satisfied12𝑐²_𝑖=^𝑚_𝑗=1𝑎_𝑖𝑗(𝑖=1,…,𝑚). (3.2)

In addition, we try to minimize the following norms of truncation error which is discussed by Dormand [9]:‖𝜏⁽⁴⁾‖₂=(²_𝑖=1(𝜏_𝑖⁽⁴⁾)²)^1/2,‖𝜏^′(4)‖₂=(³_𝑖=1(𝜏^′_𝑖⁽⁴⁾)²)^1/2, (3.3) where 𝜏_𝑖⁽⁴⁾ and 𝜏^′_𝑖⁽⁴⁾ are defined as in [2],𝜏₁⁽⁴⁾=12𝑏^𝑇𝑐²−124,𝜏₂⁽⁴⁾=𝑏^𝑇1𝐴𝑒−,𝜏24_′1(4)=16𝑏^𝑇𝑐³−124,𝜏^′₂⁽⁴⁾=𝑏^′𝑇1(𝑐.𝐴𝑒)−8,𝜏^′₃⁽⁴⁾=𝑏^′𝑇1𝐴𝑐−.24 (3.4)

The algebraic conditions (3.1)-(3.2) together with dispersion relation (2.15) and zero-dissipation condition 𝑆(𝐻)≡1 (i.e., 𝛽₃=0) formed a system of nonlinear equations. Letting 𝑐₃ as the free parameter and solving all nonlinear equations simultaneously yield the following one-parameter family of solution in term of 𝑐₃:𝑐₂=1√30(−𝐹+3(5+5)),𝑎₂₁=1√1800(−15−35+𝐹)²,𝑎₃₁1=−𝑐160³(𝐹(√5+3)(1−√5+2𝑐₃)+4√5−20−20𝑐₃−12𝑐₃√𝑎5),₃₂=1𝑐160³(𝐹(√5+3)(1−√5+2𝑐₃)+60𝑐₃−12𝑐₃√5+4√𝑏5−20),₁=−[𝐹(√5−5)(10+2√5−45𝑐₃−9𝑐₃√5+30𝑐₃²)+1800𝑐₃²−900𝑐₃−180𝑐₃√5]/[240𝑐₃(15+3√5−𝐹−30𝑐₃𝑏)],₂√=−𝐹(5−5)(−5−√5+10𝑐₃)+60√5+300−600𝑐₃80(−15−3√5+𝐹+30𝑐₃),𝑏₃𝐹=−6𝑐₃(15+3√5−𝐹−30𝑐₃),𝑏^′₁=[𝐹(√5−5)(5+√5−2𝑐₃√5−20𝑐₃+15𝑐₃²)+300𝑐₃²+60√5−100−120𝑐₃√5]/[40𝑐₃(−15−3√5+𝐹+30𝑐₃𝑏)],_′2√=−𝐹(5−5)(−2+3𝑐₃)−180𝑐₃+1208(−15−3√5+𝐹+30𝑐₃),𝑏^′₃=−5+3√5−𝐹2𝑐₃(15+3√5−𝐹−30𝑐₃), (3.5) where 𝐹=30(−1+√5).

The above set of solution will give ‖𝜏⁽⁴⁾‖₂=2.635231381×10⁻² and ‖𝜏^′(4)‖₂ can be written in terms of 𝑐₃. By using minimize command in MAPLE, ‖𝜏^′(4)‖₂ has a minimum value zero at 𝑐₃=1.047071398. Since RKN is a one-step method, we only consider the case 𝑐₃∈[0,1] where ‖𝜏^′(4)‖₂ lies between [1.872175321×10⁻²,5.796540889×10⁻²]. We developed a equidistant nodes code where 𝑐₃=0.0+0.01𝑖,𝑖=0,1,2,…. The method gives an accurate results for most of the problems tested when 𝑐₃=1, which will also gives ‖𝜏^′(4)‖₂=1.872175321×10⁻² and the dispersion constant is −(1/40320)𝑧⁷+𝑂(𝑧⁹) 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 Table 1 where 𝑐₂=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.

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

Figure 1 

(Energy Conservation). The logarithm error of energy (ERR) at each integration point when solving the harmonic oscillator for $\omega =1$ with initial condition $q_0=1\text{,}$$q_0^{\prime }=0$ and $\mathrm{\Delta }t=1/5$.

Figure 2 

The error at each integration point when solving the harmonic oscillator with $\omega =8$ with initial condition $q_0=1\text{,}$$q_0^{\prime }=-2$ and $\mathrm{\Delta }t=1/20$.

Figure 3 

The global error at each integration point when solving the “almost” periodic problem (Problem 2) with $\mathrm{\Delta }t=1/40$.

Figure 4 

The global error at each integration point when solving the “almost” periodic problem (Problem 3) with $\mathrm{\Delta }t=1/2$.

Figure 5 

The global error at each integration point when solving the inhomogeneous problem with $\mathrm{\Delta }t=1/10$.

Figure 6 

The error at each integration point when solving the inhomogeneous system with $\mathrm{\Delta }t=1/20$.

Figure 7 

The error at each integration point when solving the Duffing's equation with $\mathrm{\Delta }t=1/2$.

Figure 8 

The efficiency curve for Problem 1 with $\omega =5\text{,}$$t_{\text{end}}=10^5$ and time step $\mathrm{\Delta }t=0.1/2^i\text{,}$$i=0,\dots ,4$.

Figure 9 

The efficiency curve for Problem 2 with $t_{\text{end}}=2·10^5$ and time step $\mathrm{\Delta }t=1/(5·2^i)\text{,}$$i=0,\dots ,4$.

Figure 10 

The efficiency curve for Problem 3 with $t_{\text{end}}=10^5$ and time step $\mathrm{\Delta }t=0.8/2^i\text{,}$$i=0,\dots ,4$.

Figure 11 

The efficiency curve for Problem 4 with $t_{\text{end}}=10^5$ and time step $\mathrm{\Delta }t=0.1/2^i,$$i=1,\dots ,5$.

Figure 12 

The efficiency curve for Problem 5 with $t_{\text{end}}=5·10^4$ and time step $\mathrm{\Delta }t=0.03-0.05i\text{,}$$i=0,\dots ,4$.

Figure 13 

The efficiency curve for Duffing's equation with $t_{\text{end}}=10^5$ and time step $\mathrm{\Delta }t=1/2^i,$$i=1,\dots ,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 𝑞^′′(𝑡)=−𝜔²𝑞(𝑡),𝑞(0)=𝑞₀,𝑞^′(0)=𝑞^′₀,𝑡∈[0,𝑡_end]. (4.1)
Exact solution: 𝑞(𝑡)=𝑐₁sin(𝜔𝑡)+𝑐₂cos(𝜔𝑡)
Total energy as given in [1],𝐸(𝑞,𝑞^′𝑞)=²2+𝑞^′22=𝑎²2, (4.2) where 𝑎 depends on the initial conditions.

Problem 2 (An “almost” Periodic Orbit problem studied by Stiefel and Bettis [12]). One has 𝑧^′′+𝑧=0.001𝑒^𝑖𝑡,𝑧(0)=1,𝑧^′(0)=0.9995𝑖,𝑡∈[0,𝑡_end]. (4.3)
Exact solution:𝑧(𝑡)=(1−0.0005𝑖𝑡)𝑒^𝑖𝑡. (4.4)
We write in the equivalent form𝑞₁^′′+𝑞₁(𝑡)=0.001cos(𝑡),𝑞₁(0)=1,𝑞^′₁𝑞(0)=0,_2′′+𝑞₂(𝑡)=0.001sin(𝑡),𝑞₂(0)=0,𝑞^′₂(0)=0.9995. (4.5) Exact solution 𝑞₁(𝑡)=cos(𝑡)+0.0005𝑡sin(𝑡), 𝑞₂(𝑡)=sin(𝑡)−0.0005𝑡cos(𝑡)

Problem 3 (An “almost” Periodic Orbit problem studied by Franco and Palacios [13]). One has 𝑧^′′+𝑧=𝜖𝑒^𝑖𝜓𝑡,𝑧(0)=1,𝑧^′(0)=𝑖,𝑡∈[0,𝑡_end], (4.6) where 𝜖=0.001 and 𝜓=0.01. The analytical solution 𝑧(𝑡)=𝑢(𝑡)+𝑖𝑣(𝑡) is given by 𝑢(𝑡)=1−𝜖−𝜓²1−𝜓²𝜖cos(𝑡)+1−𝜓²cos(𝜓𝑡),𝑣(𝑡)=1−𝜖𝜓−𝜓²1−𝜓²𝜖sin(𝑡)+1−𝜓²sin(𝜓𝑡). (4.7)

Problem 4 (Problem studied by van der Houwen and Sommeijer [6]). One has 𝑞^′′(𝑡)=−𝑣²𝑞(𝑡)+(𝑣²−1)sin(𝑡),𝑞(0)=1,𝑞^′(0)=𝑣+1, (4.8) where 𝑣≫1,𝑡∈[0,𝑡_end].
Exact solution is 𝑞(𝑡)=cos(𝑣𝑡)+sin(𝑣𝑡)+sin(𝑡). Numerical result is for the case 𝑣=10.

Problem 5 (Inhomogeneous system studied by Franco [5]). One has 𝑞^′′(𝑡)+1012−992−9921012𝑞(𝑡)=𝜀932cos(2𝑡)−992sin(2𝑡)932sin(2𝑡)−992,1cos(2𝑡)𝑞(0)=−1+𝜀,𝑞^′(0)=−1010+2𝜀,𝑡∈[0,𝑡_end]. (4.9)
Exact solution:𝑞(𝑡)=−cos(10𝑡)−sin(10𝑡)+𝜀cos(2𝑡)cos(10𝑡)+sin(10𝑡)+𝜀sin(2𝑡). (4.10)

Problem 6 (The undamped Duffing's equations as given in [14]). One has 𝑞^′′(𝑡)+𝑞(𝑡)+(𝑞(𝑡))³=0.002cos(1.01𝑡),𝑞(0)=0.200426728067,𝑞^′(0)=0,𝑡∈[0,𝑡_end]. (4.11) The exact solution computed by the Galerkin method and given by 𝑞(𝑡)=⁴_𝑖=0𝑎_2𝑖+1cos[1.01(2𝑖+1)𝑡], (4.12) with 𝑎₁=0.200179477536, 𝑎₃=0.246946143×10⁻³, 𝑎₅=0.304014×10⁻⁶, 𝑎₇=0.374×10⁻⁹, and 𝑎₉<10⁻¹².

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, Δ𝑡 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].

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.

Acknowledgment

This work is partially supported by IPTA Fundamental Research Grant, Universiti Putra Malaysia (Project no. 05-10-07-385FR).