Adaptive Time Steps Runge-Kutta Methods: Comparative Analysis of Simulation Time in Nonlinear and Harmonically Excited Pendulum and Duffing Oscillators

This work was carried out in collaboration between both authors. Author TAOS, the first and lead author conceived the research idea and discussed with author OOA, the second and corresponding author. Author OOA carried out extensive literature study on the subject matter. From the literature point of view, it was jointly agreed by both authors that the idea has potentials research gap to be filled if investigated and will further enrich the literature. Author TAOS sourced for relevant governing equations from the literature, modified where necessary and carried out computer simulation experiments. Results obtained were analyzed by the lead author and discussed with the second author. Author OOA developed the introduction section and harmonized it with sections (Methodology and Results) written by author TAOS to produce the first draft of the manuscript. Both authors jointly reviewed the draft manuscript and approved. ABSTRACT Time management without integrity compromise is an integral part of good engineering practice. The present study investigated for the required computation time in nonlinear and harmonically excited oscillators (Pendulum and Duffing). Adaptive time steps simulation of their governing equations with personal computer were performed by Runge-Kutta methods (RK2, RK3, RK4, RK5, RK5M) and one blend (RKB) comprising unsteady and steady solutions. The respective Pendulum and Duffing Poincare result at damp quality (4, 0.168), excitation amplitude (1.5, 0.21) and excitation frequency (2/3, 1.0) were used to validate the codes developed in FORTRAN environment. Actual simulation time was monitored at three different lengths of excitation periods (40000, 80000 and 120000) using the current time subroutine call command. Except for RK2, the validation Poincaré results compare well with the counterpart available in the literature for the oscillators. The actual computation time decrease rapidly with increasing order of Runge-Kutta method, but suffered relative increase for the blended method. The difference in computation time required between RK5 and RK5M is negligible for all studied cases. The Pendulum required longer actual computation time (4-115 seconds) than Duffing (3-52 seconds). The respective normalized range of time step for Pendulum and Duffing formed a simple average ratios {(1.0): (7.5): (13.3): (26.0): (24.0): (29.7)} and {(1.0): (3.7): (5.1): (9.1): (8.0): (11.2)} for RK2, RK3, RK4, RK5, RK5M and RKB. It is concluded that substantial time management can be achieved by the Runge-Kutta methods except RK2 that permitted strictly shorter simulation time steps leading to longer actual simulation time and consequently simulated largely an unacceptable Poincare result.


INTRODUCTION
Extensive studies have been carried out by seasoned researchers on pendulum and Duffing oscillator systems dynamics. The dynamics of self-excited Froude pendulum under the influence of external Froude forcing were examined by [1]. Chaotic vibrations and quasiperiodic rotation of the system was studied by means of Melnikov criterion, Lyapunov exponent, poincaré maps and bifurcation diagrams. In the self excited Froude pendulum where there is presence of additional characteristic self excitation, it was found that the type of regular motion is controlled by Rayleigh damping coefficients. Research article [2] was on global dynamics of an autoparametric four degree-of-freedom (DOF) spring-mass-pendulum system. A novel approach based on forward-time solutions for finite-time Lyapunov exponent was adopted in the study. The outcome of the work has shed some light on the significant contributions of other modes present in a multi-DOF system in the modifications of the entire system. Many other research efforts have also investigated the dynamics of unexcited pendulum (see [3][4][5][6][7]). Furthermore, some other researchers' are becoming more interested on studies which focus on the understanding of the dynamics of this pendulum system when it is excited. Research article [8] characterized the harmonically excited pendulum using fractal disk dimension. A fourth order Runge-Kutta method was used as a tool for simulating the system response over a 101 by 101 parameter points on plane defined by the forcing amplitude and the damping factor. The study revealed that the three adjustable parameters (Drive frequency, forcing amplitude and damping factor) considered exert pronounced effects on the steady state response of the non-dimensioned damped sinusoidally driven nonlinear pendulum. The authors' work has indeed enriched the available literature on the excited pendulum dynamics and its possible applications.
Chaotic responses of a harmonically excited spring-pendulum system were investigated by [9]. The original system was reduced to a second-order approximate system using the method of multiple scales. The results obtained by the authors showed that the secondorder approximations technique gives more satisfactory results. The results shared better qualitative agreement with the original system when compared to the first-order approximations. Chaotic dynamics of a spherical pendulum with a harmonically excited vibration suspension was studied by [10]. The governing equations for a lightly damped spherical pendulum were considered. The physical parameters leading to chaotic solutions in terms of spherical angles were derived by the vanishing Melnikov-Marsden(MHM) integral. The results obtained by the authors showed the existence of real zeros of the MHM integral can be interpreted as the possible chaotic motion of the harmonically forced spherical pendulum. This paper provides an interesting platform for understanding the possible dynamics of harmonically excited pendulum system.
The dynamics of a periodically driven Duffing resonator coupled elastically to a Van der Pol oscillator was explored by [11]. A two-time scales method was used to obtain the frequencyamplitude modulation. The study showed that the internal resonance leads to an anti-resonance response of the Duffing oscillator and a stagnant response of the Van der Pol oscillator. In addition, the couple oscillator system showed a hysteric response pattern and symmetry breaking facets. In [12] paper, the chaotic character of Duffing oscillator was analyzed and the Melnikov method of determining chaotic threshold of Duffing oscillators was explained. The authors' findings showed that the weak periodic signal can be easily and efficiently detected. Research article of [13] examined the nonlinear dynamics of plasma oscillations modelled by a forced modified Van der Pol-Duffing oscillator. The outcome of the study could be helpful to experimentalists who have keen interests in stabilizing an oscillating system with external forcing. Several other researchers have also investigated the dynamics of Duffing oscillator (selected samples included [14][15][16]). The growing interests for Duffing oscillator system dynamics when excited have motivated some researchers in this field in the recent time. An experimental investigation of the response of a harmonically excited Duffing oscillator has been carried out by [17]. The authors successfully fabricated a mechanical analogue of hard Duffing equation with strong nonlinearity. Findings obtained from the study of the forced response of the system showed complicated and chaotic dynamics at low frequency regime. The author of [18] studied the response of a lightly damped hard Duffing oscillator to harmonic excitations using numerical and experimental approaches. The outcome of both experimental and numerical study showed complex chaotic dynamics at low frequency regime with increase in excitation level. The influence of adding fast harmonic excitation on the entrainment area of the main parametric resonance for excited Van der Pol-Mathieu-Duffing oscillator has been examined by [19]. The method of averaging was adopted to derive a governing model for describing slow dynamic of the oscillator.
Bifurcation and frequency response curves were developed from the obtained equations. The outcome of the authors' study revealed that fast harmonic excitation is capable of altering the nonlinear characteristics spring behaviour. This paper is a good contribution to the understanding of Duffing oscillator dynamics. Authors of research article [20] considered the chaotic dynamics of a coupled Duffing oscillator which was excited by harmonic parametric excitations. It was shown in the authors' study that phase disorder by randomly choosing the initial phase excitations can suppress spatio-temporal chaos in the system coupled by chaotic Duffing oscillator. Results obtained by the authors are outstanding contributions to chaotic dynamics of harmonically excited Duffing oscillator. The authors of this present paper have equally carried out some studies on the dynamics of harmonically excited Duffing oscillators and found that the dynamical system have many interesting engineering applications (see [21][22][23]).
Runge-Kutta method has been found to be an effective tool for investigating nonlinear dynamics. In [24] paper, the solutions of the nonlinear Duffing oscillator with the damping effect are obtained using Fourth Order Runge-Kutta method. Adaptive time step has been demonstrated by many researchers to be highly useful in various science and engineering fields. Adaptive time stepping has been shown to be a resourceful tool for controlling the accuracy of simulations and on the enhancement of their efficiencies. Authors of [25] employed adaptive time Runge-Kutta algorithms as tool for seeking the chaotic solutions of a harmonically excited Duffing oscillator. Several other studies have also demonstrated the usefulness of adaptive time step technique (samples are [26][27][28]).
In the foregoing, concerted research efforts have been made on the investigation of harmonically excited pendulum and Duffing oscillator systems dynamics using Runge-Kutta techniques. The importance of time management without integrity compromise is enormous part of good engineering practices. This present study is motivated by the need to further enrich the available literature through computer time saving research efforts in nonlinear systems. This present work explored adaptive time steps Runge-Kutta method for comparative study of simulation time in selected nonlinear systems (Harmonically excited pendulum and Duffing oscillators).

Governing Equations and Runge-Kutta Methods
The governing equations (1) and (2) used for this study were obtained respectively from [29,30]. The respective drive parameters for Pendulum and Duffing oscillators are damp quality ( , q γ ), excitation amplitude ( , o g P ), excitation frequency ( , D ω ω ) and time (t).
To simulate equations (1) and (2) with any Runge-Kutta method require its first order pair transformation on assumptions that: x linear velocity = The transformed equations (3 & 4) and (5 & 6) are obtained respectively for the non-linear harmonically excited Pendulum and Duffing oscillators.
The focus of this present study is the same as [31] except that the under-listed popular Runge-Kutta methods operation was based on adaptive time steps and limited to a maximum of two halfsteps. When necessary increase and decrease in the time step ( t ∆ ) is made by equation (7) and (8) respectively at specified tolerance ( ε ) of 5 4 10 − × and computed error ( t ε ) during simulation of either Pendulum or Duffing oscillator by any of the method (codes developed in FORTRAN environment) from initial conditions (0, 0) through unsteady and steady simulation periods (120,000).  (6) can be sought using equation (9), with φ being an incremental weighting function. The general form for φ is given by equation (10) (11) to (17) with time step ( h t ← ∆ ) and Table 1 provided the corresponding Butcher's table. In addition, the details of coefficients for the five popular Runge-Kutta methods used in this study are provided in Table 2. i

Simulation Parameters
The present study simulated dynamics of Pendulum and Duffing oscillators at ( The time required to progress the simulation to 40000, 80000 and 12000 were obtained by system clock which responds to a subroutine call command-CALL SYSTEM_CLOCK (ITIME  The steady Poincare section simulated by RK2 is given by Fig. 2 which differ drastically in qualitative terms from the counterpart result obtained for other methods. This is a manifestation of the instability of Runge-Kutta of order two when driven by adaptive time step to predict the dynamics of harmonically excited nonlinear Pendulum. Interestingly, the second order Runge-Kutta method was able to simulate the Poincare section of harmonically excited Duffing oscillator to the same degree of quality as its counter methods as shown in Fig. 3. Hence, it can be argued that the stability of second order is system sensitive. In summary, the qualitative result in Figs. 1 and 3 are respectively the same as can be found in [29,30] hence validated computer codes developed for the present study.

RESULTS AND DISCUSSION
From Tables 3 to 5 it is shown that RK2 takes longest actual simulation time (38-115seconds) while RK5 and RK5M takes shortest actual simulation time (4-11seconds) for the specified three simulation length periods. Except for the RKB the normalized actual simulation time for the methods follows a rapidly decreasing trend as shown in Fig. 4 Fig. 8 compared well with the results reported by [25] each showing clearly two characteristic distinct peaks.

CONCLUSIONS
This study has shown that simulation time of quality Poincare section of harmonically excited nonlinear Pendulum and Duffing oscillators can generally be reduced drastically by use of adaptive Runge-Kutta methods. However exceptional observation was made in respect of the adaptive second order Runge-Kutta method. Relative to other methods, its simulation time steps selection range was narrowly constrained, total simulation time exponentially larger and worst. Still, the reliability of the simulated Poincare section cannot always be guaranteed. Therefore adaptive second order Runge-Kutta method has limited use for simulation of quality Poincare section of nonlinear oscillators and summarily of what use is time managed without reliable result?