EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6042 ISSN 1307-5543 – ejpam.com Published by New York Business Global A Novel Explicit Two-Derivative Runge-Kutta-Nyström Method with Energy Conservation for the Integration of Second-Order Periodic ODEs Zhuoyu Sun1,∗, K.C. Lee1,∗, N.H.A. Aziz3, I. Hashim1, M.A Alias1, N. Senu2,4 1 Department of Mathematical Sciences, Universiti Kebangsaan Malaysia, 43600 UKM Bangi, Selangor, Malaysia 2 Institute for Mathematical Research, Universiti Putra Malaysia, 43400 UPM, Serdang, Malaysia 3 Department of Mathematical Sciences, Faculty of Intelligent Computing, Universiti Malaysia Perlis (UniMAP), Kampus Alam UniMAP Pauh Putra, 02600 Arau, Perlis, Malaysia 4 Department of Mathematics and Statistics, Universiti Putra Malaysia, 43400 UPM, Serdang, Malaysia Abstract. A novel trigonometrically-fitted explicit two-derivative Runge-Kutta-Nyström (TFET- DRKN(5)) method with three-stage and fifth-order for solving a class of special second-order (sys- tem) ODEs in the form of u′′ = f (t, u) with periodic solutions is proposed. Order conditions of the new explicit two-derivative Runge-Kutta-Nyström (ETDRKN(5)) method are derived using Taylor expansion and comparison of step size, h over Taylor method and general formula of the ETDRKN(5) method. Trigonometrically-fitting technique is implemented into the ETDRKN(5) method to form the TFETDRKN(5) method. Stability analysis of the new proposed method is thoroughly investigated and discussed. Algebraic order of the ETDRKN(5) method is investigated. Numerical experiments for the TFETDRKN(5) method are conducted versus error accuracy, num- ber of function evaluations and computational time. Numerical tables and graphs demonstrate that the TFETDRKN(5) method has higher effectiveness and accuracy compared to selected existing methods. Further study for one typical real-word experiment is conducted. Besides, Hamiltonian energy, Lagrangian energy and momentum conservation of proposed method are investigated to outlook the energy conservation property. The related energy exchange of the above three energies during the tested real-world experiment is illustrated. 2020 Mathematics Subject Classifications: 65L05, 65L06, 65L20, 70H12 Key Words and Phrases: Second-order ordinary differential equations, trigonometrical inte- gration, explicit two-derivative Runge-Kutta-Nyström method, stability analysis, algebraic order analysis, energy conservation ∗Corresponding author. ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6042 Email addresses: a2408527882@163.com (Z. Sun), kclee 1017@ukm.edu.my (K. C. Lee) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 2 of 45 1. Introduction In this study, We consider the special second-order ordinary differential equations (ODEs) u′′ (t) = f (t, u (t)) , u (t0) = u0, u ′ (t0) = u′0, u : R → Rk, f : R× Rk → Rk, t ∈ [t0, tend] . (1) Second-order ODEs with periodic solutions are common in many scientific and engineering fields, modeling systems with inherent periodic behavior. These equations are used in ap- plications such as mechanical oscillators, electrical circuits with inductors and capacitors, and molecular vibrations, capturing the essential dynamics of systems that exhibit regular oscillations (see [1],[2] and so on). In the beginning, Runge-Kutta (RK) methods were designed to solve the first-order ODEs [3, 4]. With the development, Nyström [5] extended RK methods to second-order ODEs, now known as Runge-Kutta-Nyström (RKN) methods, building upon the origi- nal ideas of his predecessors. Compared with the general RK methods solving first order ODEs, RKN methods are designed specifically for second-order ODEs, which makes them more efficient and accurate for such types of problems. They use both the first and second derivatives of the solution, which helps to achieve higher accuracy. However, many stan- dard RKN methods do not account for the property of periodic solutions, often leading to unsatisfactory numerical results. To address this, many researchers have attempted to modify RKN methods by incorporating trigonometrically fitting techniques. Among these efforts, some have proven highly successful. The main idea behind trigonometrically fitted methods is to seek an exact integration for differential equations whose solutions can be expressed as linear combinations of certain functions {cos(λt), sin (λt) , λ > 0}. Paternoster [6] considered the construction of RKN methods for ODEs with oscillatory solutions and derived RK and RKN methods which integrate trigonometric polynomials exactly by using the linear stage representation of a RK method given in Albrecht’s ap- proach [7]. Franco [8] constructed new explicit RKN methods up to order 5, specially adapted to the numerical integration of perturbed oscillators. Li [9] developed fifth and sixth-order trigonometrically fitted three-derivative Runge-Kutta (TFTHDRK) by using rooted trees theory and B-series. At the same year, Li et al. [10] Extended explicit pseudo two-step Runge–Kutta–Nyström (EEPTSRKN) methods have been proposed for the nu- merical integration of oscillatory systems, inheriting the framework of the explicit pseudo two-step Runge-Kutta-Nyström methods. These methods are designed to integrate ex- actly the unperturbed problem y′′ +My = 0, achieve a maximum step order of s+ 2 and stage order of s+ 1. And very recently, Demba et al. successively developed explicit, embedded explicit, phase- and amplification-fitted 5(4) diagonally explicit, phase- and amplification-fitted 5(4) diagonally implicit Runge-Kutta-Nyström methods, making a significant contribu- tion to the progress of trigonometrical integration research (see[11–14] for further de- tails). These advancements in RKN methods hava significantly improved the efficiency and accuracy of solving oscillatory and periodic initial value problems. Several develop- Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 3 of 45 ing methods, including the fourth-order four-stage explicit trigonometrically-fitted RKN methods, demonstrate smaller global errors compared to existing methods [11]. Embed- ded pairs, such as the 4(3) explicit trigonometrically-fitted RKN methods, offer reduced computational costs with fewer function evaluations per step [12]. Additionally, phase- and amplification-fitted methods have been developed, maintaining high-order conver- gence and providing improved stability intervals [13]. Notably, the PFAFRKN5(3) and PFAF-DIRKN5(4)4 pairs show enhanced accuracy and efficiency over their counterparts, offering promising results for solving complex oscillatory problems [14]. Motivated by the goal of enhancing the order and numerical accuracy of RKN meth- ods, many researchers have modified the classical RKNmethods into two-derivative Runge- Kutta-Nyström (TDRKN) methods by incorporating the second derivatives of f -evaluation into the formulation. Chen et al. [15] extended the traditional RKN methods for general second-order ODEs to TDRKN methods involving the third derivative. The order crite- ria for TDRKN methods were derived using a novel version of Nyström tree theory and the associated B-series theory. They developed a two-stage explicit TDRKN method of fourth-order and a three-stage explicit TDRKN method of fifth-order. Jator [16] intro- duced a trigonometrically-fitted implicit third derivative Runge-Kutta-Nyström method (TTRKNM), with coefficients dependent on the frequency and step size, for periodic ini- tial value problems (IVPs). The TTRKNM consists of a pair of methods derived from its continuous version, which are used to produce simultaneous approximations of the solution and its first derivative at each point in the interval of interest. The stability properties of the method were discussed, and numerical experiments were conducted to demonstrate its accuracy and efficiency. Chen et al. [17] proposed a new family of modified TDRKN methods for solving second-order oscillatory ODEs. Order conditions were derived us- ing Nyström tree theory and B-series theory. They established trigonometrically-fitted conditions and constructed two practical explicit trigonometrically-fitted TDRKN (TFT- DRKN) methods. The phase properties of the new integrator were examined, and their periodicity regions were determined. Numerical experiments demonstrated the efficiency and competence of the new methods. For a comprehensive discussion on the construction and analysis of trigonometrically-fitted or exponentially-fitted TDRKN methods, readers are directed to the works in [18–20] and so on. Among articles about solving second-order ODEs using TDRKN methods, those ad- dressing general second-order ODEs account for the vast majority. There are numerous RKN methods for integrating a general class of second-order IVPs. However, there is a lack of research applying TDRKN methods to solving special classes of second-order ODEs in the form of u′′ = f (t, u). Additionally, there is a shortage of studies on con- structing trigonometrically-fitted TDRKN methods for directly solving special classes of second-order ODEs with periodic solutions. Hence, we focus on developing a new effi- cient trigonometrically-fitted TDRKN (TFETDRKN(5)) method with high effectiveness and accuracy and energy conservation for solving special classes of second-order periodic ODEs like (1). In sections 2 and 3 we construct a new ETDRKN(5) method with three-stage and fifth- order. In section 4 we apply the trigonometrically-fitting technique to the ETDRKN(5) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 4 of 45 method to derive the TFETDRKN(5) method, which constitutes the primary focus of this study. In section 5 we conduct a thorough stability analysis of the TFETDRKN(5) method involving stability matrix, dispersion and dissipation error analysis, and stability proper- ties and regions. In section 6 we prove the algebraic order of the ETDRKN(5) method is 5. In sections 7 and 8 we conduct numerical experiments for the TFETDRKN(5) method against error accuracy, number of function evaluations and computational time compared with some existing similar methods to demonstrate that the TFETDRKN(5) method has higher accuracy and effectiveness. In section 9 we demonstrate the enhanced accuracy and energy conservation of the ETDRKN(5) method, which is achieved by employing the trigonometrical integration technique. Section 10 is dedicated to the concluding remarks. 2. The Formulation of the Two-Derivative Runge-Kutta-Nyström Method The TDRKN method is a numerical technique designed to solve second-order ODEs. These methods are extensions of standard Runge-Kutta (Runge-Kutta-Nyström) meth- ods tailored for second-order ODEs (system), making them highly effective for problems involving dynamics like mechanical systems, wave motion, or electrical circuits. The key advantage of TDRKN methods lies in their ability to directly handle second-order ODEs by using both position and velocity (or their equivalents) in the system and incorporate higher-order accuracy while maintaining computational efficiency, leading to greater effi- ciency and accuracy compared to traditional methods that require transforming the system into a set of first-order equations. Letting u′′′ (t) = g (t, u (t) , u′ (t)) and combining u′′ (t) = f (t, u (t)) in the problem (1), we can convert this problem to the following form( u′′ (t) f ′ (t, u (t)) ) = ( f (t, u (t)) g ( t, u (t) , u′ (t) )) , (u (t0) f (t0, u (t0)) ) = ( u0 u′′0 ) . (2) The two-derivative Runge-Kutta-Nyström (TDRKN(5)) method with s-stage, which is developed by incorporating the third derivative, u ′′′ (t) into the formulation, is illustrated as un+1 = un + hu′n + h2 2 f (tn, un) + h3 s∑ i=1 d̄ig ( tn + cih, Ui, U ′ i ) , u′n+1 = u′n + hf (tn, un) + h2 s∑ i=1 d̃ig ( tn + cih, Ui, U ′ i ) , Ui = un + cihu ′ n + 1 2 (cih) 2f (tn, un) + h3 s∑ i=1 Āijg ( tn + cih, Uj , U ′ j ) , U ′ i = u′n + cihf (tn, un) + h2 s∑ i=1 Ãijg ( tn + cih, Uj , U ′ j ) , (3) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 5 of 45 where ci, d̄i, d̃i, Āij , Ãij , i = 1, · · · , s are real numbers. The scheme (3) can be expressed as in Kronecker’s block product notation: un+1 = un + hu′n + h2 2 f (tn, un) + h3 ( d̄⊗ Ik×k ) G ( U,U ′) , u′n+1 = u′n + hf (tn, un) + h2 ( d̃⊗ Ik×k ) G ( U,U ′) , U = un + h ( c⊗ u′n ) + 1 2 h2 ( ccT ⊗ f (tn, un) ) + h3 ( Ā⊗ Ik×k ) G ( U,U ′) , U ′ = e⊗ u′n + h ( cT ⊗ f (tn, un) ) + h2 ( Ã⊗ Ik×k ) G ( U,U ′) , (4) where e = (1, · · · , 1)T , c = (c1, · · · , cs)T , d̄ = ( d̄1, · · · , d̄s ) , d̃ = ( d̃1, · · · , d̃s ) are s- dimensional vectors, Ā = ( Āij ) s×s , Ã = ( Ãij ) s×s are s × s matrices and Ik×k is k × k identity matrix. The block vectors in Rs×k are U = (U1, · · · , Us) , U ′ = (U ′ 1, · · · , U ′ s) , G (U,U ′) = (g (t0 + c1h, U1, U ′ 1) , · · · , g (t0 + csh, Us, U ′ s)) . (5) An alternative expression of the scheme (3) is given as follows: un+1 = un + hu′n + h2 2 f (tn, un) + h3 s∑ i=1 d̄ili, u′n+1 = u′n + hf (tn, un) + h2 s∑ i=1 d̃ili, (6) where li = g ( tn + cih, un + cihu ′ n + 1 2 (cih) 2 f (tn, un) + h3 s∑ i=1 Āij lj , u ′ n + cihf (tn, un) + h2 s∑ i=1 Ãij lj ) . (7) It is convenient to represent the scheme (3) using the following Butcher tableau: Table 1: Butcher Tableau for the TDRKN(5) method. c Ā Ã d̄ d̃i The TDRKN(5) method is explicit if Āij = 0, Ãij = 0 for i ⩽ j and implicit if Āij ̸= 0, Ãij ̸= 0 for i ⩽ j and involves only one evaluation of f and many g evaluations of per step. We select to construct a new explicit method, named as the ETDRKN(5) method. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 6 of 45 3. Construction of the Fifth-Order Efficient Two-Derivative Runge-Kutta-Nyström Method First, we utilize the Taylor series expansion to determine the coefficients of the ET- DRKN(5) method of three-stage and fifth-order. By equating this expansion to the theo- retical solution, which is also represented by a Taylor series, and performing some simpli- fying assumptions, we derive a system of non-linear equations via Maple. These equations are known as the order conditions of the ETDRKN(5) method. The order conditions for u: Third-order: 3∑ i=1 d̄i = 1 6 . (8) Fourth-order: 3∑ i=1 d̄ici = 1 24 . (9) Fifth-order: 3∑ i=1 d̄ici 2 = 1 60 . (10) Sixth-order: 3∑ i=1 d̄ici 3 = 1 120 , (11a) 3∑ i=1  i−1∑ j=1 d̄iĀij  = 1 720 , (11b) 3∑ i=1  i−1∑ j=2 d̄iĀijcj  = 1 720 . (11c) The order conditions for u′: Second-order: 3∑ i=1 d̃i = 1 2 . (12) Third-order: 3∑ i=1 d̃ici = 1 6 . (13) Fourth-order: 3∑ i=1 d̃ici 2 = 1 12 . (14) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 7 of 45 Firth-order: 3∑ i=1 d̃ici 3 = 1 20 , (15a) 3∑ i=1  i−1∑ j=2 d̃iÃijcj  = 1 120 , (15b) 3∑ i=1  i−1∑ j=1 d̃iĀij  = 1 120 . (15c) Sixth-order: 3∑ i=1 d̃ici 4 = 1 30 , (16a) 3∑ i=1  i−1∑ j=1 d̃iciĀij  = 1 180 , (16b) 3∑ i=1  i−1∑ j=2 d̃iĀijcj  = 1 720 , (16c) 3∑ i=1  i−1∑ j=2 d̃iÃijcj 2  = 1 360 , (16d) 3∑ i=1  i−1∑ j=2 d̃iciÃijcj  = 1 180 . (16e) Referring to idea of reducing the number of N-trees required for order conditions obtained in [15], we propose the following relations (t0 + cih) 2 = t0 2 + 2t0 · (cih) + (cih) 2 = t0 2 + ∫ t0+cih t0 2t0dτ+ ∫ t0+cih t0 ∫ τ t0 2dσdτ, (t0 + cih) 3 = t0 3 + 3t0 2 · (cih) + 3t0 · (cih)2 + (cih) 3 = t0 3 + ∫ t0+cih t0 3t0 2dτ + ∫ t0+cih t0 ∫ τ t0 6t0dσdτ+ ∫ t0+cih t0 ∫ τ t0 ∫ σ t0 6dςdσdτ, (17) which yield the simplifying assumptions Ãe = 1 2 c2, Āe = 1 6 c3, i.e., s∑ j=1 Ãij = 1 2 ci 2, s∑ j=1 Āij = 1 6 ci 3. (18) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 8 of 45 We utilize the order conditions (8-15c,12-15c,11a) and simplifying assumptions (18) (i = 2, j = 1) to generate the coefficients of the ETDRKN(5) method which are presented in the following Butcher tableau: Table 2: Butcher Tableau for the ETDRKN(5) method with three-stage and fifth-order. 0 0 0 0 0 0 0 c2 Ā21 0 0 Ã21 0 0 c3 Ā31 Ā32 0 Ã31 Ã32 0 d̄1 d̄2 d̄3 d̃1 d̃2 d̃3 After solving the order conditions mentioned above, we obtains the solution with one free coefficient Ā31 as follow: Ā32 = −Ā31 + 1 30 − √ 5 75 . (19) Then, we generate all coefficients of the ETDRKN(5) method: c1 = 0, c2 = 1 2 + √ 5 10 , c3 = 1 2 − √ 5 10 , Ā21 = 1 30 + √ 5 75 , Ā32 = Ā31 + 1 30 − √ 5 75 , Ã21 = 3 20 + √ 5 20 , Ã31 = 0, Ã32 = 3 20 − √ 5 20 , d̄1 = 1 24 , d̄2 = 1 16 − √ 5 48 , d̄3 = 1 16 + √ 5 48 , d̃1 = 1 12 , d̃2 = 5 24 − √ 5 24 , d̃3 = 5 24 + √ 5 24 . (20) Euclidean norms (2-norms) [21] of the error terms for the sixth-order approximations of u and u′ of the ETDRKN(5) method are defined as ∥∥∥τ (6)∥∥∥ 2 = √√√√ K∑ i=1 ( τ (6) i )2 , ∥∥∥τ ′(6) ∥∥∥ 2 = √√√√ L∑ i=1 ( τi′ (6) )2 , ∥∥∥τ (6)g ∥∥∥ 2 = √√√√K+L∑ i=1 [ (τ (6) i )2 + (τ ′(6) i )2 ] , (21) where K and L are the total number of local truncation errors for u and u′, τ (6) and τ ′(6) are the sixth-order local truncation error norms for u and u′ respectively, and ∥∥∥τ (6)g ∥∥∥ 2 is sixth-order global truncation error norm used to decide the values of Ā31 and Ā32. We determine Ā31 by minimizing ∥∥∥τ (6)g ∥∥∥ 2 , as a result, it generates a minimum value 8.487 × 10−3 at Ā31 = −1288 452405 , which produces Ā32 = 98209 2714430 − √ 5 75 . In the meanwhile, at Ā31 = − 1288 452405 , we have ∥∥τ (6)∥∥ 2 ≈ 1.626 × 10−3 and ∥∥∥τ ′(6)∥∥∥ 2 ≈ 4.599 × 10−3. So far, we have preliminarily completed the construction of the ETDRKN(5) method. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 9 of 45 4. Implementation of the Trigonometrically-Fitting Technique Firstly, we look back to the scheme (3), we add χ̄i in front of un in U and χ̃i in front of f (tn, un) in U ′ i and set χ̄ = (χ̄1, · · · , χ̄s) T and χ̃ = (χ̃1, · · · , χ̃s) T . The introduction of these two new coefficent vectors is the first step for our trigonometrically-fitting technique. Rewriting the scheme (3), we get un+1 = un + hu′n + h2 2 f (tn, un) + h3 s∑ i=1 d̄ig ( tn + cih, Ui, U ′ i ) , u′n+1 = u′n + hf (tn, un) + h2 s∑ i=1 d̃ig ( tn + cih, Ui, U ′ i ) , Ui = χ̄iun + cihu ′ n + 1 2 (cih) 2f (tn, un) + h3 s∑ i=1 Āijg ( tn + cih, Uj , U ′ j ) , U ′ i = u′n + χ̃icihf (tn, un) + h2 s∑ i=1 Ãijg ( tn + cih, Uj , U ′ j ) , (22) The Butcher Tableau for the renewable ETDRKN(5) method is Table 3: Butcher Tableau for the renewable ETDRKN(5) method. c χ̄ χ̃ Ā Ã d̄ d̃i Then, we integrate the exponential terms, eiλt and e−iλt at every stage and w = λh with λ ∈ R, we obtain e±iciw = χ̄i ± iciw − 1 2 (ciw) 2 ∓ iw3 3∑ i=1 Āije ±icjw, e±iciw = 1± iχ̃iciw − w2 3∑ i=1 Ãije ±icjw. (23) The equations corresponding to u and u′ are as follows e±iw = 1± iw − 1 2 w2 ∓ iw3 3∑ i=1 d̄ie ±iciw, e±iw = 1± iw − w2 3∑ i=1 d̃ie ±iciw. (24) The relation cosw = 1 2 ( eiw + e−iw ) , sinw = 1 2i ( eiw − e−iw ) , cos (ciw) = 1 2 ( eiciw + e−iciw ) , sin (ciw) = 1 2i ( eiciw − e−iciw ) (25) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 10 of 45 are substituted in the equation (25), then we obtain the trigonometric functions in terms of w cosw = 1− 1 2 w2 + w3 3∑ i=1 d̄i sin (ciw), sinw = w − w3 3∑ i=1 d̄i cos (ciw), cosw = 1− w2 3∑ i=1 d̃i cos (ciw), sinw = w − w2 3∑ i=1 d̃i sin (ciw), cos (ciw) = 1− 1 2 (ciw) 2 + w3 3∑ i=1 Āij sin (cjw), sin (ciw) = ciw − w3 3∑ i=1 Āij cos (cjw), cos (ciw) = 1− w2 3∑ i=1 Ãij cos (cjw), sin (ciw) = ciw − w2 3∑ i=1 Ãij sin (cjw). (26) The trigonometrically-fitted coefficients Āij (w) and Ãij (w) can be obtained by Āi,i−1 (w) = cos (cjw)− 1 + 1 2(ciw) 2 − w3 i−2∑ j=1 Āij sin (cjw) w3 sin (ci−1w) , Ãi,i−1 (w) = 1− cos (cjw)− w2 i−2∑ j=1 Ãij cos (cjw) w2 cos (ci−1w) . (27) Then χ̄i (w) and χ̃i (w) are determined on the coefficients Āi,i−1 (w) and Ãi,i−1 (w) through χ̄i (w) = cos (ciw) + 1 2 (ciw) 2 − w3 3∑ i=1 Āij sin (cjw), Ãi,i−1 (w) = sin (ciw) + w2 3∑ i=1 Ãij sin (cjw) ciw . (28) Modifying the equations (27,28), we obtain the trigonometrically-fitted coefficients Āij (w), Ãij (w), χ̄i (w) and χ̃i (w) Ā21 (w) = c2w − sin (c2w) w3 , Ā32 (w) = c3w − sin (c3w)− Ā31w 3 w3 cos (c2w) , Ã21 (w) = 1− cos (c2w) w2 , Ã32 (w) = 1− cos (c3w)− Ã31w 2 w2 , χ̄2 (w) = cos (c2w) + 1 2 (c2w) 2, χ̄3 (w) = cos (c3w) + 1 2 (c3w) 2 − Ā32 sin (c2w) , χ̃2 (w) = sin (c3w) c3w , χ̃3 (w) = sin (c3w) + w2Ã32 sin (c2w) c3w . (29) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 11 of 45 Then, by substituting all the coefficients (20) of the ETDRKN(5) method into the equations (29) and employing the eighth-order Taylor series expansion, we compute the frequency-dependent coefficients for Āij (w), Ãij (w), χ̄i (w) and χ̃i (w) (see eqs 59 in Appendices). As a result, we conduct a novel trigonometrically-fitted explicit two-derivative Runge- Kutta-Nyström with the three-stage and fifth-order (TFETDRKN(5)) method. The Butcher Tableau for the TFETDRKN(5) method is Table 4: Butcher Tableau for the TFETDRKN(5) method. c1 (w) 0 0 0 0 0 0 0 0 c2 (w) χ̄2 (w) χ̃2 (w) Ā21 (w) 0 0 Ã21 (w) 0 0 c3 (w) χ̄3 (w) χ̃3 (w) Ā31 (w) Ā32 (w) 0 Ã31 (w) Ã32 (w) 0 d̄1 (w) d̄2 (w) d̄3 (w) d̃1 (w) d̃2 (w) d̃3 (w) Observe that as w approaches zero, all trigonometrically-fitted coefficients derived above for the TFETDRKN(5) method converge to the original constant coefficients for the ET- DRKN(5) method. Remark 1. The addition of vectors χ̄ and χ̃, and the generation of trigonometrically-fitted coefficients is to make the proposed method more accurate in solving second-order ODEs with periodic solutions, which is also the key point of this paper. 5. Stability Analysis 5.1. Stability matrix For testing the stability properties of the TFETDRKN(5) method, we use the homo- geneous test equation ([17]) u′′ = −µ2u, µ > 0 (30) where µ is the natural frequency. By applying the scheme (3) to the test equation (30), we obtain ( un+1 hu′n+1 ) = N ( z2;w )( un hu′n ) , (31) where N ( z2;w ) =  1− z2 2 + z4d̄(w) T ( I + z2Ã (w) )−1 χ̃ (w) ce −z2 + z4d̃(w) T ( I + z2Ã (w) )−1 χ̃ (w) ce 1− z2d̄(w) T ( I + z2Ã (w) )−1 e 1− z2d̃(w) T ( I + z2Ã (w) )−1 e  , z = µh (32) is called the stability matrix of the TFETDRKN(5) method. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 12 of 45 Remark 2. λ in w = λh is called fitted frequency, while µ in z = µh is called natural frequecny. For the trigonometrically-fitted methods applied to problems like (1), the fitted frequency of methods generally differs from the natural frequency of the theoretical solution. However, for the linear oscillator u′′ (t)+λ2u (t) = 0, the method’s fitted frequency matches the natural frequency of the theoretical solution [22]. So we will conduct the stability analysis considering two cases, λ = µ and λ ̸= µ. 5.2. Dispersion and dissipation error analysis Stability behavior of the numerical solution depends on eigenvalues ξi(i = 1, 2) of the stability matrix (32). Eliminating u′n+1 and u′n from the equations (31) by replacing the subscript 0 by 1 and 1 by 2 yields the difference equation u2 − tr (N)u1 + det (N)u0 = 0. (33) Accordingly, we obtain the characteristic polynomial σ (z, w; ξ) σ (z, w; ξ) = ξ2 − tr (N) ξ + det (N) = 0, (34) where tr (N) and det (N) are the trace and the determinant of the stability matrix (32), ξ1 and ξ2 are ξ1 = tr (N)− √ tr2 (N)− 4 det (N) 2 , ξ2 = tr (N) + √ tr2 (N)− 4 det (N) 2 . (35) Definition 1 ([17]). For the characteristic polynomial (34), we set the ratio κ = λ µ , then introduce two quantities ϕ (z;κ) = z − arccos ( tr (N) 2 √ det (N) ) , ψ (z;κ) = 1− √ det (N), (36) which are denoted as the dispersion error and the dissipation error respectively. TFET- DRKN(5) method is dispersive and dissipative of order s and order k respectively, if ϕ (z;κ) = O ( zs+1 ) , ψ (z;κ) = O ( zk+1 ) . (37) And the necessary and sufficient condition for the TFETDRKN(5) method is dispersive and dissipative of order s and order k respectively, is that the expression γ (z;κ) and ϑ (z;κ) satisfies [23] γ (z;κ) = tr (N)− 2 √ det (N) cos (z) = O ( zs+2 ) , ϑ (z;κ) = det (N)− 1 = O ( zk+1 ) . (38) By calculating via Maple, for λ = µ, we have γ (w) = 1 4 w4+ ( − 4 45 + √ 5 120 ) w6+O ( w8 ) , ϑ (w) = −1 6 w4+ ( − √ 5 160 + 1 160 ) w6+O ( w8 ) , (39) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 13 of 45 for λ ̸= µ, we have γ (z;κ) = 1 4 z4 + 1 7200 ((12 √ 5 + 20)κ2 + 33 √ 5− 645)z6 +O ( z8 ) , ϑ (z;κ) = −1 6 z4 − √ 5 7200 ((4 √ 5 + 12)κ2 − 13 √ 5 + 33)z6 +O ( z8 ) . (40) So the TFETDRKN(5) method is dispersive of order two and dissipative of order three from conditions (38). 5.3. Stability properties and regions In this subsection, we study the stability properties and regions of the TFETDRKN(5) method comprehensively to predict its efficiency region and range corresponding to differ- ent step sizes for two cases, λ = µ and λ ̸= µ. 5.3.1. For the case of same frequency Definition 2 ([14]). We define the stability region for the TFETDRKN(5) method as follows: RS = {w : |ξi| < 1, i = 1, 2} , where ξi (i = 1, 2) are the eigenvalues of N ( w2 ) . The stability regions for different ranges of the TFETDRKN(5) method with λ = µ are shown in Figs. 1-6. The stability region for w is the intersection region of stability region for w1 and w2. Note that the stability region still exists when w is relatively larger, but the detailed research in such case is not very meaningful because in numerical simulations, the step size h usually takes at least one decimal place and the frequency λ will not be very large, so w will not be so lager. Fig. 7 and Figs. 8, 9 show the distribution relationship of ξ, ξ1 and ξ2 versus w. More- over, the curves of ξi, abs(ξi)(i = 1, 2) and zeros of denominator are also shown intuitively. According to such three figures, the approximate w range that meets Here is your expres- sion with quotation marks added (using proper LaTeX quoting): latex“|ξi| < 1, i = 1, 2” can be concluded. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 14 of 45 Figure 1: Stability region for the real part of w ∈ [−6, 6]. Figure 2: Partial stability region for the real part of w1 ∈ [−6, 6]. Figure 3: Partial stability region for the real part of w2 ∈ [−6, 6]. Figure 4: Stability region for the real part of w ∈ [−20, 20] Figure 5: Partial stability region for the real part of w1 ∈ [−20, 20]. Figure 6: Partial stability region for the real part of w2 ∈ [−20, 20]. Figure 7: Stability property of characteristic polynomial for ξ against w. Figure 8: Stability property of characteristic polynomial for ξ1 against w. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 15 of 45 Figure 9: Stability property of characteristic polynomial for ξ2 against w. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 16 of 45 Figure 10: Periodicity stability region for z against w. Definition 3 ([14]). An interval I = (−w̃, 0) is said to be an interval of stability of the TFETDRKN(5) method if for all w̃ ∈ I, it is |ξ1,2| < 1. The performance of the range of w when ξ is between -1 and 1 is worth studying, then by calculation via Maple, we obtain the interval of stability for ξ1 is (−1.866, 0) and for ξ2 is (−1.588, 0), so the interval of stability for the TFETDRKN(5) method is (−1.588, 0). 5.3.2. For the case of different frequency Definition 4 ([17]). For the TFETDRKN(5) method with the stability matrix N ( z2;w ) and its eigenvalues, ξi (i = 1, 2), the region in complex plane RP = {(w, z) : |ξi| < 1, i = 1, 2} is called the periodicity region of the TFETDRKN(5) method. The periodicity region of the TFETDRKN(5) method with λ ̸= µ is depicted in Fig. 10. For similar reason to the stability regions mentioned in the subsection 5.3.1, we are more focused on the distribution when w and z take relatively smaller values. Figs. 11-14 illustrate the stability property of characteristic polynomial (34) for differ- ent w, z range and contour on different value of ξ. Note that four different w, z ranges of four figures can fit the frequency, λ, µ and step size, h of most experiments, so we can obtain the approximate w, z regions that meet “|ξi| < 1, i = 1, 2”. 6. Algebraic Order Analysis In order to derive a general formula for the higher-order derivatives of the theoretical solution of the problem (1), we propose expressing the theoretical solution, u (t) at t = t0 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 17 of 45 Figure 11: Stability property of characteristic poly- nomial with w, z ∈ [−5, 5] and contour on ξ = −8. Figure 12: Stability property of characteristic poly- nomial with w, z ∈ [0, 3] and contour on ξ = −8. Figure 13: Stability property of characteristic poly- nomial with w, z ∈ [−1.588, 1.588] and contour on ξ = −2. Figure 14: Stability property of characteristic poly- nomial with w, z ∈ [−0.1, 0.1] and contour on ξ = −2. from first to seventh derivatives u(1) = z, u(2) = f, u(3) = g, u(4) = guz + gzf, u (5) = guuz 2 + 2guzfz + guf + gzzf 2 + gzg, u(6) = guuuz 3 + 3 ( guuzfz 2 + guufz + guzzf 2z + guzgz + guzf 2 + gzzgf ) + gug + gzzzf 3 + gz (guz + gzf) , u(7) =12 ( guzzgfz + guuzf 2z ) + 10guzgf + 6 ( guuufz 2 + guuzzf 2z2 + guuzgz 2 + gzzzgf 2 + guzzf 3 ) + 4 ( guugz + guzzzf 3z + guzguz 2 + guzgzfz + gzzgufz + gzzgzf 2 + guuuzfz 3 ) + 3 ( guuf 2 + gzzg 2 ) + gz ( guuz 2 + 2guzfz + guf + gzzf 2 + gzg ) + ′ u (guz + gzf) + guuuuz 4 + gzzzzf 4. (41) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 18 of 45 Then we demonstrate the LTEs for u and u′ generated by the ETDRKN(5) method. The Taylor series expansion is applied over h to the theoretical solution u(tn + h), its derivative u′(tn + h), the numerical solution un+1, and its derivative u′n+1. For the ET- DRKN(5) method, the LTEs at tn+1 for both the theoretical and numerical solutions, along with their respective derivatives, are given as follows LTEno fitted(u) = un+1 − u (tn + h) = ( 1 2400 ( 5 3 − √ 5 )( g2z + gzguz + 4gug )) · w6 +O ( w7 ) , LTEno fitted(u ′) = u′n+1 − u′ (tn + h) = ( 1 1200 (√ 5− 5 3 )( guzgz − 1 2 guu ) z2 + ( √ 5 325731600 ( 271443gzzgzf − 723848guug − 1033931gz 2 ) − 452405gzzgzf + 1925540gz 2 ) z + √ 5 651463200 ( 271443gzzgzf 2 − 1447696guzgf − 2339305gzgzf − 271443gz 2g ) − 1 1440 ( gzzf 2 − 860697 90481 gzf − gzg ) gz ) · w6 +O ( w7 ) . (42) LTEs for u and u′ of the ETDRKN(5) method up to h5 vanish, so the ETDRKN(5) method is algebraic of order 5. 7. Numerical Experiments and Results We solve problems (1) with the theoretical solution generated by the following linear space { 1, t, t2, cos(λt), sin(λt) } (43) We use several typical second-order ODEs with periodic solutions (43) to demonstrate the effectiveness and accuracy of the TFETDRKN(5) method. We conduct the experiments testing and get the numerical results by plotting the tables and figures comparing other selected existing methods in terms of the maximum global truncation error, the num- ber of function evaluations and CPU time in seconds. The numerical methods used for comparison are • TFTDRKN3s5: Existing trigonometrically-fitted two-derivative Runge-Kutta-Nyström method of three stage fifth-order [17]. • STDRKN5(3): Existing efficient two-derivative Runge-Kutta-Nyström Method of three stage fifth-order [18]. • PFAFRKN6-6ER: Existing optimized sixth-order explicit Runge-Kutta-Nyström Method to solve oscillating systems [24]. • TFIRKN5: Existing fifth-order improved Runge-Kutta-Nyström Method using trigonometrically [25]. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 19 of 45 Several notations used are • h: step size • B: endpoint of value t • Time(s): CPU time in seconds • NFE: number of function evaluation • MGTE: maximum global truncation error from numerical experiments • 4.084364(−16): 4.084364× 10−16 The maximum global truncation error (MGTE) is defined by MGTE = max |un − u (tn) |, where un is numerical solution and u (tn) is theoretical solution at endpoint, tn = B. Experiment 1: Consider the homogeneous problem [26]u ′′ (t) = −64u (t) , t ∈ [0, 100], u (0) = −1 4 , u′ (0) = −1 2 , with the theoretical solution: u (t) = − 1 16 sin (8t) + 1 4 cos (8t) . Natural frequency, λ = 8. Experiment 2: Consider the linear second-order problem [27]{ u′′ (t) = −u (t) + 2, t ∈ [0, 100], u (0) = 0, u′ (0) = 1, with the theoretical solution: u (t) = 2 (1− cos (t)) + sin (t) . Fitted frequency, λ = 1. Experiment 3: Consider one typical stiff problem [28]u ′′ (t) = ( −1 2 ( β2 + 1 ) −1 2 ( β2 − 1 ) −1 2 ( β2 − 1 ) −1 2 ( β2 + 1 ) ) u (t) , t ∈ [0, 100], u (0) = (1,−1)T , u′ (0) = (1,−1)T , Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 20 of 45 Figure 15: Experiment 4 with the theoretical solution: u (t) = (cos (t) + sin (t) , − (cos (t) + sin (t)))T . Fitted frequency, λ = 1. We select β = 2, 102 and 104 to test this problem, respectively. Experiment 4 (Fig. 15): Consider one undamped mass-spring system [29] ( m1 0 0 m2 ) u′′ (t) = ( − (k1 + k2) k2 k2 −k2 ) u (t) , t ∈ [0, 100], u (0) = (0, 0)T , u′ (0) = (3,−2)T . (44) We set m1 = 13 3 kg, m2 = 5kg, k1 = 1N ·m−1, k2 = 2N ·m−1, so the theoretical solution is u (t) = (3 sin (t) , −2 sin (t))T . (45) Fitted frequency, λ = 1. Experiment 5 (Fig. 16) : Consider another undamped mass-spring system [29] Figure 16: Experiment 5  ( m1 0 0 m2 ) u′′ (t) = ( −k k k −k ) u (t) , t ∈ [0, 100], u (0) = (0, 0)T , u′ (0) = (10,−1)T . (46) We set m1 = 8kg, m2 = 7kg,k1 = 56N ·m−1, so the theoretical solution is u (t) = ( 73 15 t+ 77 √ 15 225 sin (√ 15t ) , 73 15 t− 88 √ 15 225 sin (√ 15t ))T . (47) Fitted frequency, λ = √ 15. Let t = 0 be the moment when the two cars meet. Then two cars begin to compress the spring till the spring’s elastic potential energy reaches its maximum critical point, the spring stretches and becomes longer, and finally the two cars separate (not considering any car hitting the wall). Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 21 of 45 Figure 17: RK4 method for ex- periment 8 with h = 10−5 and B = 10. Figure 18: TFETDRKN(5) method for experiment 8 with h = 10−3 and B = 10. Figure 19: Numerical simulation for experiment 8 with h = 10−3 and B = 10. Experiment 6: Consider the “two coupled oscillators with different frequencies” sys- tem [30]u ′′ (t) = ( −1 0 0 −2 ) u (t) + ( 2εu1 (t)u2 (t) ε ( u1 2 (t) + 4u2 3 (t) ) ) , t ∈ [0, 5], u (0) = (1, 1)T , u′ (0) = (0, 0)T . (48) The use of a perturbation method yields the first-order fitted frequencies [31] λu1 = 1 · λu2 = √ 2− 3ε√ 2 . (49) We set ε = 10−4. As we are unable to obtain the theoretical solution for experiment 8, we use the numerical solution generated by the classical fourth-order Runge-Kutta method with step size, h = 10−5 as its theoretical solution and compare it with other comparative methods. Figs. 17, 18, 19 display the numerical simulation with B = 10. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 22 of 45 8. Tables and Graphs of Numerical Results Table 5: Experiment 1. h METHOD MGTE NFE Time(s) TFETDRKN(5) 4.084364(-16) 12000 0.519 TFTDRKN3s5 1.802952(-13) 12000 0.539 0.025 STDRKN5 1.562715(-5) 12000 0.587 PFAFRKN6-6ER 5.960929(-13) 24000 0.719 TFIRKN5 1.912294(-4) 32000 0.632 TFETDRKN(5) 1.144546(-17) 15000 0.595 TFTDRKN3s5 1.934776(-14) 15000 0.632 0.020 STDRKN5 5.115801(-6) 15000 0.749 PFAFRKN6-6ER 1.039706(-13) 30000 0.856 TFIRKN5 6.171445(-5) 40000 0.754 TFETDRKN(5) 1.142831(-19) 20000 0.756 TFTDRKN3s5 1.089515(-15) 20000 0.932 0.015 STDRKN5 1.214425(-6) 20000 0.938 PFAFRKN6-6ER 1.140415(-14) 40000 1.145 TFIRKN5 1.450237(-5) 53333 1.014 TFETDRKN(5) 1.737694(-22) 30000 1.250 TFTDRKN3s5 7.568127(-17) 30000 1.260 0.010 STDRKN5 1.599277(-7) 30000 1.309 PFAFRKN6-6ER 5.498041(-16) 60000 1.769 TFIRKN5 1.895651(-6) 80000 1.509 TFETDRKN(5) 2.648241(-27) 60000 2.447 TFTDRKN3s5 1.847588(-20) 60000 2.601 0.005 STDRKN5 4.996194(-9) 60000 2.688 PFAFRKN6-6ER 3.887940(-18) 120000 3.390 TFIRKN5 5.894875(-8) 160000 3.072 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 23 of 45 Table 6: Experiment 2. h METHOD MGTE NFE Time(s) TFETDRKN(5) 1.529179(-30) 12000 0.426 TFTDRKN3s5 1.790330(-22) 12000 0.470 0.025 STDRKN5 5.168409(-10) 12000 0.467 PFAFRKN6-6ER 1.998151(-18) 24000 0.675 TFIRKN5 6.103274(-9) 32000 0.595 TFETDRKN(5) 4.303836(-32) 15000 0.517 TFTDRKN3s5 1.922284(-23) 15000 0.593 0.020 STDRKN5 1.693915(-10) 15000 0.594 PFAFRKN6-6ER 4.654337(-19) 30000 0.794 TFIRKN5 1.999367(-9) 40000 0.798 TFETDRKN(5) 4.313117(-34) 20000 0.768 TFTDRKN3s5 1.082461(-24) 20000 0.771 0.015 STDRKN5 4.016910(-11) 20000 0.856 PFAFRKN6-6ER 7.367519(-20) 40000 1.083 TFIRKN5 4.739915(-10) 53333 1.087 TFETDRKN(5) 6.566061(-37) 30000 1.022 TFTDRKN3s5 1.877093(-26) 30000 1.079 0.010 STDRKN5 5.295473(-12) 30000 1.142 PFAFRKN6-6ER 5.770130(-21) 60000 1.547 TFIRKN5 6.245495(-11) 80000 1.381 TFETDRKN(5) 1.001838(-41) 60000 2.075 TFTDRKN3s5 1.833039(-29) 60000 2.340 0.005 STDRKN5 1.655133(-13) 60000 2.502 PFAFRKN6-6ER 8.089612(-23) 120000 3.285 TFIRKN5 1.951468(-12) 160000 3.022 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 24 of 45 Table 7: Experiment 3. h METHOD MGTE NFE Time(s) TFETDRKN(5) 4.256602(-21) 3000 1.194 TFTDRKN3s5 1.213163(-16) 3000 1.126 0.1 STDRKN5(3) 3.294553(-7) 3000 1.167 PFAFRKN6-6ER 1.998695(-14) 6000 1.296 TFIRKN5 3.907366(-6) 8000 1.202 TFETDRKN(5) 6.479610(-26) 6000 2.293 TFTDRKN3s5 1.184842(-19) 6000 2.369 0.05 STDRKN5(3) 1.029791(-8) 6000 2.433 PFAFRKN6-6ER 1.175856(-16) 12000 2.609 TFIRKN5 1.213642(-7) 16000 2.288 TFETDRKN(5) 9.881954(-31) 12000 4.717 TFTDRKN3s5 1.156993(-22) 12000 4.636 0.025 STDRKN5(3) 3.217851(-10) 12000 4.643 PFAFRKN6-6ER 1.187549(-18) 24000 5.514 TFIRKN5 3.787327(-9) 32000 4.811 TFETDRKN(5) 1.507622(-35) 24000 9.141 TFTDRKN3s5 1.129826(-25) 24000 9.395 0.0125 STDRKN5(3) 1.005527(-11) 24000 8.970 PFAFRKN6-6ER 1.701612(-20) 48000 10.141 TFIRKN5 1.183191(-10) 64000 9.670 TFETDRKN(5) 2.300318(-40) 48000 18.993 TFTDRKN3s5 1.103319(-28) 48000 18.874 0.00625 STDRKN5(3) 3.142220(-13) 48000 18.869 PFAFRKN6-6ER 2.679244(-22) 96000 20.045 TFIRKN5 3.697447(-12) 128000 19.941 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 25 of 45 Table 8: Experiment 4. h METHOD MGTE NFE Time(s) TFETDRKN(5) 8.759277(-21) 6000 0.349 TFTDRKN3s5 2.498800(-16) 6000 0.403 0.1 STDRKN5(3) 7.048727(-7) 6000 0.387 PFAFRKN6-6ER 7.562986(-14) 12000 0.594 TFIRKN5 8.102380(-6) 16000 0.504 TFETDRKN(5) 1.334936(-25) 12000 0.829 TFTDRKN3s5 2.440855(-19) 12000 0.922 0.05 STDRKN5(3) 2.202624(-8) 12000 1.075 PFAFRKN6-6ER 4.602406(-16) 24000 1.146 TFIRKN5 2.514992(-7) 32000 0.951 TFETDRKN(5) 2.036323(-30) 24000 1.611 TFTDRKN3s5 2.384085(-22) 24000 1.700 0.025 STDRKN5(3) 6.882515(-10) 24000 1.937 PFAFRKN6-6ER 4.373250(-18) 48000 2.187 TFIRKN5 7.849286(-9) 64000 2.040 TFETDRKN(5) 3.106811(-35) 48000 3.079 TFTDRKN3s5 2.328292(-25) 48000 3.517 0.0125 STDRKN5(3) 2.150783(-11) 48000 3.398 PFAFRKN6-6ER 5.732428(-20) 96000 4.433 TFIRKN5 2.452453(-10) 128000 4.163 TFETDRKN(5) 4.740479(-40) 96000 6.506 TFTDRKN3s5 2.273722(-28) 96000 6.956 0.00625 STDRKN5(3) 6.721108(-13) 96000 6.572 PFAFRKN6-6ER 8.526926(-22) 192000 9.334 TFIRKN5 7.664112(-12) 256000 8.725 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 26 of 45 Table 9: Experiment 5. h METHOD MGTE NFE Time(s) TFETDRKN(5) 4.622706(-11) 6000 0.366 TFTDRKN3s5 3.732735(-10) 6000 0.378 0.1 STDRKN5(3) 1.207628(-3) 6000 0.449 PFAFRKN6-6ER 1.599813(-9) 12000 0.514 TFIRKN5 1.614405(-2) 16000 0.450 TFETDRKN(5) 6.864049(-16) 12000 0.725 TFTDRKN3s5 3.684791(-13) 12000 0.789 0.05 STDRKN5(3) 3.769888(-5) 12000 0.894 PFAFRKN6-6ER 6.523295(-12) 24000 1.134 TFIRKN5 4.568589(-4) 32000 1.355 TFETDRKN(5) 1.039558(-20) 24000 1.416 TFTDRKN3s5 3.597012(-16) 24000 1.495 0.025 STDRKN5(3) 1.177672(-6) 24000 1.769 PFAFRKN6-6ER 2.986460(-14) 48000 2.146 TFIRKN5 1.388403(-5) 64000 1.829 TFETDRKN(5) 1.582759(-25) 48000 2.941 TFTDRKN3s5 3.512880(-19) 48000 3.156 0.0125 STDRKN5(3) 3.680162(-8) 48000 3.282 PFAFRKN6-6ER 1.852287(-16) 96000 4.082 TFIRKN5 4.308162(-7) 128000 3.997 TFETDRKN(5) 2.414464(-30) 96000 6.112 TFTDRKN3s5 3.430763(-22) 96000 6.285 0.00625 STDRKN5(3) 1.149993(-9) 96000 6.428 PFAFRKN6-6ER 1.795039(-18) 192000 8.336 TFIRKN5 1.344022(-8) 256000 7.670 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 27 of 45 Table 10: Experiment 6 h METHOD MGTE NFE Time(s) TFETDRKN(5) 5.677404(-8) 300 0.015 TFTDRKN3s5 7.829360(-9) 300 0.018 0.1 STDRKN5(3) 8.618116(-8) 300 0.020 PFAFRKN6-6ER 3.727733(-12) 600 0.074 TFIRKN5 1.424941(-1) 800 0.023 TFETDRKN(5) 1.717588(-9) 600 0.037 TFTDRKN3s5 2.441104(-10) 600 0.039 0.05 STDRKN5(3) 2.677502(-9) 600 0.084 PFAFRKN6-6ER 2.800499(-14) 1200 0.105 TFIRKN5 7.075647(-2) 1600 0.047 TFETDRKN(5) 5.289724(-11) 1200 0.088 TFTDRKN3s5 7.621687(-12) 1200 0.074 0.025 STDRKN5(3) 8.348074(-11) 1200 0.121 PFAFRKN6-6ER 6.531133(-16) 2400 0.148 TFIRKN5 3.528249(-2) 3200 0.135 TFETDRKN(5) 1.641778(-12) 2400 0.134 TFTDRKN3s5 2.380660(-13) 2400 0.167 0.0125 STDRKN5(3) 2.605529(-12) 2400 0.208 PFAFRKN6-6ER 4.074754(-17) 4800 0.262 TFIRKN5 1.759057(-2) 6400 0.234 TFETDRKN(5) 5.113543(-14) 4800 0.350 TFTDRKN3s5 7.437993(-15) 4800 0.394 0.00625 STDRKN5(3) 8.136990(-14) 4800 0.376 PFAFRKN6-6ER 2.546304(-18) 9600 0.523 TFIRKN5 8.750423(-3) 12800 0.382 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 28 of 45 The following tables and figures present the numerical results, each demonstrating the performance of five different methods, with the exception of experiment 6.The model of computer for experiments is Lenovo ideapad 330 Intel Core i5-8050U (1.8GHz). Ev- idently, for the MGTE of the first seven experiments, the TFETDRKN(5) method sig- nificantly outperforms other existing methods. This is attributable not only to the in- corporation of trigonometrically-fitted terms, χ̄ and χ̃, but also to the application of the trigonometrically-fitting technique, as outlined in section 4, to a larger set of coefficients (compared to conventional trigonometrically-fitting techniques in existing literature), ren- dering them frequency-dependent due to the inclusion of these terms. Furthermore, our proposed method is specifically tailored for solving special second-order ODEs with peri- odic solutions like (1). Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 29 of 45 Figure 20: Numerical graph for the experiment 1 with B = 100 and h = 0.025 − 0.005i(i = 0, 1, 2, 3, 4). Figure 21: Numerical graph for the experiment 2 with B = 100 and h = 0.025 − 0.005i(i = 0, 1, 2, 3, 4). Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 30 of 45 Figure 22: Numerical graph for the experiment 3 with B = 100 and h = 0.1/2i(i = 0, 1, 2, 3, 4). Figure 23: Numerical graph for the experiment 4 with B = 100 and h = 0.1/2i(i = 0, 1, 2, 3, 4). Figure 24: Numerical graph for the experiment 5 with B = 100 and h = 0.1/2i(i = 0, 1, 2, 3, 4). Figure 25: Numerical graph for the experiment 6 with B = 5 and h = 0.1/2i(i = 0, 1, 2, 3, 4). Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 31 of 45 For the experiments 1 to 5, it is worth noting that the order of increasing values of the MGTE is the TFETDRKN(5), TFTDRKN3s5, PFAFRKN6-6ER, STDRKN, and TFIRKN5 method. An analysis of this phenomenon is provided below. The trigonometrically- fitting technique of the TFETDRKN(5) method is superior to that of the TFTDRKN3s5 for solving first seven special second-order ODEs with periodic solutions like (1), so TFEDTRKN(5) method outperforms TFTDRKN3s5 method. Subsequently, the STDRKN5(3) method is just an efficient TDRKNmethod without the implementation of trigonometrically- fitting technique, so it certainly performs worst among three TDRKNmethods. The reason the PFAFRKN6-6ER method (despite being a sixth-order algebraic method, which would typically result in smaller errors than fifth-order methods) and the TFIRKN5 method (where improvements are generally expected to reduce errors) perform worse than the three above mentioned TDRKN methods is that the TDRKN methods incorporate the second derivative directly into the numerical scheme (3). This integration allows for higher- order accuracy without requiring additional function evaluations, leading to more accurate solutions for the same computational effort compared to the standard RKN methods (see [15]). Unfortunately, for the experiment 6, from Fig. 25, the TFETDRKN(5) method does not show the best error accuracy but ranks in the middle among the five methods, lag- ging behind the PFAFRKN6-6ER and the TFTDRKN3s5 method, and better than the STDRKN and the TFIRKN5 method. The PFAFRKN6-6ER method, which performs poorly in error accuracy in the first seven experiments, performs the best in experiment 6, while the TFTDRKN3s5 method performs slightly better than the TFETDRKN(5) method. The problem in experiment 6 has more secular terms, u1 (t)u2 (t) , u1 2 (t) and u2 3 (t), than the other five experiments. The PFAFRKN6-6ER method implements the phase- and amplification-fitted technique during the construction of their method (see further in [13]), which is specifically designed to handle nonlinearities in oscillatory systems. It more accurately captures the phase and amplitude changes caused by the nonlinear term. In the equations shown in the experiment 6, it is the nonlinear terms that causes the amplitude and frequency changes (which is why the frequency λ in Eq 49 depends on ϵ). While, the TFETDRKN(5) method accurately captures periodic behavior through trigonometric fitting and uses second-order derivatives to improve accuracy. However, it is not specifically optimized for nonlinear terms. The TFETDRKN(5) method is not highly effective for solving second-order systems of ODEs with many nonlinear terms. Consequently, in the simulation of Experiment 6, error accumulation is relatively large compared with the PFAFRKN6 method, and the expected accuracy is not fully achieved. Nevertheless, although the proposed method does not yield the best performance in this case, it still demonstrates significance and value, owing to its lower NTE and reduced computational time (see Table 10) for simulating special second-order systems of ODEs without theoretically periodic solutions. From all numerical tables, we obtain that when B and h are fixed, the number of function evaluations of the three TDRKN methods are the same, but smaller than that of the PFAFRKN6-6ER (also smaller than the TFIRKN5 method) and the TFIRKN5 Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 32 of 45 method, which is consistent with the relatively low computational complexity of the three TDRKNmethods. Also, the computational time of the three TDRKNmethods is generally shorter than that of the PFAFRKN6-6ER and TFIRKN5 methods in all experiments, which the reason of is the three TDRKN methods are all two-derivative Runge-Kutta- Nyström methods, where the complexity of a single function evaluation is lower than that of the PFAFRKN6-6ER method which is of sixth order and the TFIRKN5 method which requires the most function evaluations during process. Overall, the TFETDRKN(5) method outperforms other compared methods when solv- ing special second order (system) ODEs which have theoretical solutions. This is due to its lowest MGTE, which translates to the best accuracy, as well as its relatively better convergence, as well as low number function evaluations and computational time. Hence, the TFETDRKN(5) method outperforms other compared methods about the effectiveness and computational time. 9. Further Study for Real-World Experiment 6 9.1. Theoretical verification of energy preservation for TFETDRKN(5) Background and notation Consider the second-order Hamiltonian system q′′(t) = f(q(t)) = −∇U(q(t)), p := q′. The Hamiltonian (total energy) is H(q, p) = 1 2 pT p+ U(q). We denote by qn, pn the numerical approximations at time tn and by qn+1, pn+1 their one-step updates by the TFETDRKN(5) method with stepsize h. The goal is to prove ∆H := H(qn+1, pn+1)−H(qn, pn) = O(h6), i.e. the local energy error is of order h6 (so energy is preserved up to the method order). The three-stage two-derivative Runge–Kutta–Nyström update has the form (matching the equation in Hamiltonian system) qn+1 = qn + hpn + h2 2 fn + h3 3∑ i=1 d̄i ℓi, pn+1 = pn + hfn + h2 3∑ i=1 d̃i ℓi, (50) where fn := f(qn) = −∇U(qn) and each stage evaluation ℓi is a shorthand for the two- derivative evaluation at the i-th stage, ℓi = g ( tn + cih, Ui, U ′ i ) , Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 33 of 45 and the stages Ui, U ′ i depend on qn, pn, fn, {ĀijandÃij}. For the algebra below it suffices that the method coefficients {ci, Āij , Ãij , d̄i, d̃i} satisfy the ETDRKN(5) order conditions (these will be used to cancel coefficient combinations). Step 1: Expand the kinetic energy term. Compute the kinetic part at n+ 1: 1 2 pTn+1pn+1 = 1 2 ( pn + hfn + h2 ∑ i d̃iℓi )T( pn + hfn + h2 ∑ j d̃jℓj ) = 1 2 pTnpn + h pTnfn + h2 2 fTn fn + h2pTn ∑ i d̃iℓi + h3fTn ∑ i d̃iℓi + h4 2 (∑ i d̃iℓi )T(∑ j d̃jℓj ) . Retain all terms up to order h5 (the last displayed term is O(h4 ·ℓ2); since we assume ℓ = O(1) these terms are of order O(h4) or higher). Step 2: Expand the potential energy term. Taylor-expand U(qn+1) about qn. Using (50) write ∆q := qn+1 − qn = hpn + h2 2 fn + h3 ∑ i d̄iℓi, and expand: U(qn+1) = U(qn) +∇U(qn) T∆q + 1 2 ∆qT∇2U(qn)∆q + 1 6 D3U(qn)[∆q,∆q,∆q] +O(∥∆q∥4). Substitute ∆q; using fn = −∇U(qn) we get U(qn+1) = U(qn) + h∇U(qn) T pn + h2 2 ∇U(qn) T fn + h3∇U(qn) T ∑ i d̄iℓi + h2 2 pTn∇2U(qn)pn + h3pTn∇2U(qn) h 2 fn +O(h3 ·h3). Retain terms up to order h5; higher-order terms will be grouped later into O(h6). Step 3: Assemble ∆H and identify grouped terms. Using the expansions above, ∆H = ( 1 2 pTn+1pn+1 − 1 2 pTnpn ) + ( U(qn+1)− U(qn) ) = [ hpTnfn + h2 2 fTn fn + h2pTn ∑ i d̃iℓi + h3fTn ∑ i d̃iℓi + h4 2 ∥ ∑ i d̃iℓi∥2 ] + [ h∇U(qn) T pn + h2 2 ∇U(qn) T fn + h3∇U(qn) T ∑ i d̄iℓi + h2 2 pTn∇2U(qn)pn + h3(other pn, fn mixed terms) ] − h pTn∇U(qn) +O(h6). Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 34 of 45 Now use fn = −∇U(qn). Several cancellations are immediate: - The O(h) terms cancel: hpTnfn + h∇U(qn) T pn − hpTn∇U(qn) = 0. - Combine O(h2) terms. Replace ∇U(qn) by −fn where convenient: h2 2 fTn fn + h2 2 ∇U(qn) T fn + h2pTn ∑ i d̃iℓi + h2 2 pTn∇2U(qn)pn = h2 2 fTn fn − h2 2 fTn fn + h2pTn ∑ i d̃iℓi + h2 2 pTn∇2U(qn)pn = h2pTn ∑ i d̃iℓi + h2 2 pTn∇2U(qn)pn. The two first terms cancel exactly. The remaining h2 expression involves pn and stage contributions ℓi. These are not yet zero but will be converted (via integration-by-parts like manipulations using the stage definitions) into a symmetric bilinear form in the stage forces f(Q·) with kernel given below. After collecting all terms up to order h5 and doing the routine algebra (expand stages Ui, U ′ i in Taylor series about (qn, pn), replace stage ℓi by their expansions in terms of fn and derivatives), the change in energy can be written in the compact bilinear form ∆H = h2 2 5∑ r,s=0 hrhs Crs(qn, pn) = h2 2 ∫∫ [0,1]2 f(Qτ ) T D(τ, σ; v) f(Qσ) dσdτ + O(h6), (51) where D(τ, σ; v) is the kernel D(τ, σ; v) := Bτ (v)Bσ(v)− ∂τAτσ(v)− ∂σAστ (v), and Aτσ, Bτ are the continuous kernels associated with the method coefficients (discrete collocation/quadrature reproduces these at nodes ci). The representation (51) is stan- dard in continuous-stage proofs: the h2-term collects into a symmetric double integral in the stage forces f(Q·) with kernel D; higher-order contributions (from cubic and quartic expansions) are absorbed into O(h6) once the method order conditions are used. Step 4: Interpretation of the kernel D and sufficient condition for energy preservation. Equation (51) shows that the leading nontrivial contribution to ∆H is the h2-term with kernel D. Therefore a sufficient and (for arbitrary f) necessary condition to annihilate that leading term is Bτ (v)Bσ(v) = ∂τAτσ(v) + ∂σAστ (v) for all τ, σ ∈ [0, 1]. (52) If (52) holds pointwise, the h2-term vanishes. (This is the continuous-stage identity analogous to the discrete condition you quoted; for discrete nodes it becomes BiBj = ∂τAτσ|(τ,σ)=(ci,cj) + · · ·.) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 35 of 45 Step 5: Cancellation of Higher-Order Terms via Order Conditions After the h2-term is removed by (52), the remaining contributions in ∆H are of higher powers of h. A careful expansion shows that the next orders that may survive are h3, h4, h5; each of these is a finite linear combination of the method coefficients multiplying deriva- tives/compositions of f and U evaluated at the stages. Crucially, those linear combinations are exactly the order conditions for the ETDRKN(5) method (the conditions you solved to obtain the coefficients: the third–sixth order conditions for q and second–sixth for p, plus the simplifying assumptions). Because TFETDRKN(5) was constructed so that these order conditions hold, each of the algebraic coefficient combinations multiplying the h3, h4, h5 contributions is zero. Therefore, when the method coefficients satisfy: • the trigonometrically-fitted stage/exponential relations (which algebraically imply (52)), and • the ETDRKN(5) order conditions (which annihilate the h3, h4, h5 coefficient combi- nations), then all contributions to ∆H up to order h5 vanish identically. The first possibly nonzero term is of order h6. Symbolically, ∆H = O(h6). Step 6: Discrete-stage remark (collocation/nodes). In practice we work with the discrete three-stage TFETDRKN(5) tableau: the continuous kernels Aτσ, Bτ are sampled at the nodes τ, σ ∈ {c1, c2, c3} and replaced by the discrete coefficients Āij(w), Ãij(w), d̄i(w), d̃i(w) (frequency dependent). The same chain of algebra applies: substitute the discrete formu- las (for example Ā21(w) = (c2w − sin(c2w))/w 3, Ã21(w) = (1 − cos(c2w))/w 2, etc.) into the discrete analog of the kernel Dij(w): Dij(w) := Bi(w)Bj(w)− ( ∂τAτσ(w) + ∂σAστ (w) )∣∣∣ (τ,σ)=(ci,cj) . Expanding Dij(w) in a Taylor series in w shows Dij(w) = O(w8) (i.e. zero up to w6) be- cause the trigonometrical fitting identities were enforced exactly and the order conditions annihilate the lower-order terms. Hence the discrete local energy error is O(h6) as well. Conclusion. Combining the continuous algebra with the discrete collocation sampling used to obtain TFETDRKN(5) coefficients, we conclude that the TFETDRKN(5) method satisfies ∆H = H(qn+1, pn+1)−H(qn, pn) = O(h6). Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 36 of 45 provided the method coefficients are chosen to satisfy the trigonometrical reproduction identities and the ETDRKN(5) order conditions. This proves that TFETDRKN(5) pre- serves the Hamiltonian up to (and including) order five; the leading energy error term is of order h6. Definition 5 ([32]). A Hamiltonian system with n-degree of freedom is characterized by dv1 dt = ∂H ∂u1 , dv2 dt = ∂H ∂u2 , · · · , dvn dt = ∂H ∂un ; du1 dt = ∂H ∂v1 , du2 dt = ∂H ∂v2 , · · · , dun dt = ∂H ∂vn , (53) where H is the Hamiltonian energy of the system. In real-word applications, ui (i = 1, · · · , n) represent the generalized coordinates and vi (i = 1, · · · , n) are generalized momentum. Definition 6. The Hamiltonian and Lagrangian energy can be expressed as, respectively H (u, v; t) = TE (v; t) + VE (u; t) , L (u, v; t) = TE (v; t)− VE (u; t) , (54) where TE and VE are the kinetic and potential energies, respectively. So it is convenient to characterized a Hamiltonian system with 2-degree of freedom, the Hamiltonian and Lagrangian energy, respectively H (u1, u2, v1, v2; t) = TE (v1, v2; t) + VE (u1, u2; t) , L (u1, u2, v1, v2; t) = TE (v1, v2; t)− VE (u1, u2; t) . (55) Definition 7. If the kinetic energy TE (v1, v2; t) of a Hamiltonian system with 2-degree of freedom (55) has a standard quadratic form such as: TE (v1, v2; t) = 1 2m1 v1 2 + 1 2m2 v2 2, (56) where v1 and v2 are the generalized momentum mentioned above (55), m1 and m2 are the masses associated with each degree of freedom, respectively. Then the concrete expression of momentum of the Hamiltonian system with 2-degree of freedom (55) can be expressed as p1 = m1 du1 dt , p2 = m2 du2 dt , (57) where du1 dt and du2 dt are the corresponding velocities, then p1 and p2 represent the linear momentum corresponding to these velocities. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 37 of 45 Figure 26: 3D theoretical phase for experiment 6 with h = 0.125 and t ∈ [0, 100]. Figure 27: 3D numerical phase for experiment 6 with h = 0.125 and t ∈ [0, 100]. Figure 28: LTE for experiment 6 with h = 0.125 and t ∈ [0, 100]. Figure 29: Comparison phase of two methods for experiment 6 with h = 0.125 and t ∈ [0, 100]. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 38 of 45 In this section, we conduct a detailed investigation of real-world experiment 6 to demonstrate the superiority of the TFETDRKN(5) method in terms of accuracy and en- ergy conservation, achieved through the use of the trigonometrical integration technique. Figs. 26-28 visualize the theoretical and numerical solutions of experiment 6 and the abso- lute value of their differences. Fig. 28 displays how the LTE changes during the integration process from t = 0 to t = 100. Note that LTE doesn’t decrease linearly but with oscillatory behavior, in the meanwhile, the oscillation amplitude is regular. The reasons for these two phenomena are that the displacement change of two cars in experiment 6 behaves in an oscillatory manner and the two ODEs of experiment 6 is linear. Also, Fig. 29 illustrates the superiority of the accuracy of trigonometrical integration technique through comparing the LTEs generated by the ETDRKN(5) and TFETDRKN(5) methods, respectively. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 39 of 45 Figure 30: Hamiltonian energy conversation for the experiment 6 with h = 0.125 and t ∈ [0, 100]. Figure 31: Lagrangian energy conversation for the experiment 6 with h = 0.125 and t ∈ [0, 100]. Figure 32: Momentum conversation for experi- ment the 6 with h = 0.125 and t ∈ [0, 100]. Figure 33: Relative energy exchange for the exper- iment 6 with h = 0.125 and t ∈ [0, 5]. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 40 of 45 Subsequently, Figs. 30-32 visually illustrate the superior performance of the Hamilto- nian energy, Lagrangian energy, and momentum conservation in experiment 6, using the TFETDRKN(5) method with h = 0.125 and t ∈ [0, 100]. It should be noted that these three errors do not exhibit a linear increase but rather display an oscillatory behavior due to the oscillatory manner of the system’s ODEs in experiment 6. As time t progresses from 1 to 100, the already small relative errors increase in an oscillatory manner, further high- lighting their stability. At last, Fig. 33 shows the relative energy exchange process with h = 0.125 and t ∈ [0, 5] for the experiment 6. Throughout the process, the three quan- tities exhibit a regular, periodic behavior, most notably following an oscillatory pattern resembling sine and cosine functions like sin (t) and cos (t), as they evolve over time. 10. Conclusions In this paper, we proposed a new three-stage and fifth-order TFETDRKN(5) method, based initially on the construction of the ETDRKN(5) method. We performed an detailed stability analysis of the TFETDRKN(5) method including stability matrix, dispersion and dissipation error analysis and stability properties and regions. We derived the stability matrix and decided to conduct the stability analysis considering two cases, λ = µ and λ ̸= µ. In this case, the dispersion and the dissipation error were introduced, followed by the presentation of the necessary and sufficient conditions. It was then proven that the TFETDRKN(5) method achieves dispersion of order 2 and dissipation of order 3. For the case, λ = µ, we defined and plotted the stability region, then illustrated the stability property for the TFETDRKN(5) method. Stability interval was introduced and calculated, (−1.588, 0). For the other case, λ ̸= µ, we similarly defined and plotted the pe- riodicity stability region, then analysed the stability property of characteristic polynomial for the TFETDRKN(5) method. After that, we proved that the ETDRKN(5) method is algebraic of order 5. Importantly, to evaluate the numerical effectiveness and accuracy of the TFETDRKN(5) method, we performed numerical tests on the first seven experi- ments with theoretically periodic solutions, and on experiment 6, which lacks a periodic solution. The new proposed method was compared with four selected existing methods by plotting the MGTE against the step size, h. The results of tables and graphs and discussion demonstrated that the TFETDRKN(5) method outperforms other compared methods by producing nearly the least absolute maximum global truncation error (except for experiment 6) , relatively small time of computations and relatively less (the same as the other TDRKN methods) number of function evaluations in experiments 1 to 5. Despite not achieving the lowest numerical MGTE, the new method proved crucial for simulating special second order (system) ODEs like those in experiment 6, which lack theoretical solutions. We further extended our study of experiment 6, and the results show that the TFETDRKN(5) method significantly outperforms the EDTKRN(5) method, following the application of the trigonometrically-fitting technique. By comparing the global truncation errors over time, the superiority of the TFETDRKN(5) becomes evident. Moreover, the proposed method successfully preserves the energy conservation property. For future research, some new two-derivative Runge-Kutta-Type methods can be de- Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 41 of 45 veloped specifically for solving third-order ODEs. Even though some third-order ODEs can be decomposed into first-order or second-order ODEs to be solved using the proposed method in this paper (or other articles), this will increase computational complexity and make the process more cumbersome. Meanwhile, no one has studied the two-derivative Runge-Kutta-Nyström methods for solving third-order ODEs. So in the future, we can focus on developing the two-derivative Runge-Kutta-Nyström methods for solving third- order ODEs with periodic solutions, even exponential solutions based on our study. Acknowledgements This study is supported by the Grant Scheme (Ref. No. GGPM-2023-029) from Universiti Kebangsaan Malaysia. Declaration of competing interest The authors declare no conflict of interest. References [1] R. P. Agarwal, S. R. Grace, and D. O’Regan. Oscillation theory for second order dynamic equations. CRC Press, 2002. [2] M. Condon, A. Deaño, and A. Iserles. On highly oscillatory problems arising in electronic engineering. ESAIM: Mathematical Modelling and Numerical Analysis, 43(4):785–804, 2009. [3] R. Ciii. About the numerical solution of differential equations. Mathematical Annals, 46(2):167–178, 1895. [4] K. Wiii. Contribution to the approximate integration of total differential equations. Mathematical Annals, 45(1):435–453, 1901. [5] E. J. Nyström. About the numerical integration of differential equations. Fennic Society of Sciences, 50(13):1–55, 1925. [6] B. Paternoster. Runge-Kutta-Nyström methods for ODEs with periodic solutions based on trigonometric polynomials. Applied Numerical Mathematics, 28(2-4):401– 412, 1998. [7] P. Albrecht. The extension of the theory of a-methods to rk methods. In Numerical Treatment of Differential Equations, Proceedings 4th Seminar NUMDIFF, volume 4, pages 9–18, 1987. [8] J. M. Franco. Runge-Kutta-Nyström methods adapted to the numerical integration of perturbed oscillators. Computer Physics Communications, 147(3):770–787, 2002. [9] J. Li. Trigonometrically fitted three-derivative runge–kutta methods for solving os- cillatory initial value problems. Applied Mathematics and Computation, 330:103–117, 2018. [10] J. Li, S. Deng, and X. Wang. Extended explicit pseudo two-step rkn methods for oscillatory systems y′′ +my = f(y). Numerical Algorithms, 78:673–700, 2018. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 42 of 45 [11] M. A. Demba, N. Senu, and F. Ismail. Trigonometrically-fitted explicit four-stage fourth-order Runge-Kutta-Nyström method for the solution of initial value problems with oscillatory behavior. Global Journal of Pure and Applied Mathematics, 12(1):67– 80, 2016. [12] M. A. Demba, N. Senu, and F. Ismail. An embedded 4 (3) pair of explicit trigonometrically-fitted Runge-Kutta-Nyström method for solving periodic initial value problems. Appl. Math. Sci, 11:819–838, 2017. [13] M. A. Demba, H. Ramos, P. Kumam, and W. Watthayu. A phase-fitted and amplification-fitted explicit Runge-Kutta-Nyström pair for oscillating systems. Math- ematical and Computational Applications, 26(3):59, 2021. [14] M. A. Demba, N. Senu, H. Ramos, P. Kumam, and W. Watthayu. A phase-and amplification-fitted 5 (4) diagonally implicit Runge-Kutta-Nyström pair for oscilla- tory systems. Bulletin of the Iranian Mathematical Society, 49(3):24, 2023. [15] Z. Chen, Z. Qiu, J. Li, and X. You. Two-derivative Runge-Kutta-Nyström methods for second-order ordinary differential equations. Numerical Algorithms, 70:897–927, 2015. [16] S. N. Jator. Implicit third derivative Runge-Kutta-Nyström method with trigono- metric coefficients. Numerical Algorithms, 70:133–150, 2015. [17] Z. Chen, L. Shi, S. Liu, and X. You. Trigonometrically fitted two-derivative runge- kutta-nyström methods for second-order oscillatory differential equations. Applied Numerical Mathematics, 142:171–189, 2019. [18] T. S. Mohamed, N. Senu, Z. B. Ibrahim, and N. M. A. Nik Long. Efficient two- derivative Runge-Kutta-Nyström methods for solving general second-order ordinary differential equations. Discrete Dynamics in Nature and Society, 2018, 2018. [19] J. O. Ehigie, M. Zou, X. Hou, and X. You. On modified TDRKN methods for second-order systems of differential equations. International Journal of Computer Mathematics, 95(1):159–173, 2018. [20] K. C. Lee, M. A. Alias, N. Senu, and A. Ahmadian. On efficient frequency-dependent parameters of explicit two-derivative improved Runge-Kutta-Nyström method with application to two-body problem. Alexandria Engineering Journal, 72:605–620, 2023. [21] G. Xue and Y. Ye. An efficient algorithm for minimizing a sum of euclidean norms with applications. SIAM Journal on Optimization, 7(4):1017–1036, 1997. [22] H. Ramos and J. Vigo-Aguiar. On the frequency choice in trigonometrically fitted methods. Applied Mathematics Letters, 23(11):1378–1381, 2010. [23] P. J. Van der Houwen and B. P. Sommeijer. Diagonally implicit Runge-Kutta-nyström methods for oscillatory problems. SIAM Journal on Numerical Analysis, 26(2):414– 429, 1989. [24] M. A. Demba, H. Ramos, P. Kumam, and W Watthayu. An optimized sixth-order ex- plicit RKN method to solve oscillating systems. In Proceedings of the XXVI Congreso de Ecuaciones Diferenciales y Aplicaciones. XVI Congreso de Matemática Aplicada. Servicio de Publicaciones de la Universidad de Oviedo, 2021. [25] W. J. Hasan and K. A. Hussain. Fifth order improved Runge-Kutta-Nystrom method using trigonometrically-fitting for solving oscillatory problems. Al-Nahrain Journal Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 43 of 45 of Science, 25(4):63–67, 2022. [26] B. S. Attili, K. Furati, and M. I. Syam. An efficient implicit Runge-Kutta method for second order systems. Applied Mathematics and Computation, 178(2):229–238, 2006. [27] D. O. Awoyemi. A p-stable linear multistep method for solving general third or- der ordinary differential equations. International Journal of Computer Mathematics, 80(8):985–991, 2003. [28] J. M. Franco, I. Gómez, and L. Rández. Four-stage symplectic and p-stable SDIRKN methods with dispersion of high order. Numerical Algorithms, 26:347–363, 2001. [29] J. Lebl. Notes on diffy Qs: differential equations for engineers. Independent, 2014. [30] X. Wu and B. Wang. Geometric integrators for differential equations with highly oscillatory solutions. Springer, 2021. [31] J. Vigo-Aguiar, T. E. Simos, and J. M. Ferrándiz. Controlling the error growth in long–term numerical integration of perturbed oscillations in one or several frequencies. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460(2042):561–567, 2004. [32] S. Lynch. Dynamical systems with applications using MAPLE. Springer, 2010. Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 44 of 45 Appendices Frequency-Dependent Coefficients of Two-Derivative Runge-Kutta-Nyström Method Ā21 (w) = 1 30 + √ 5 75 + ( − 1 1200 − 11 √ 5 30000 ) w2 + ( 13 1260000 + 29 √ 5 6300000 ) w4 + ( − 17 226800000 − 19 √ 5 567000000 ) w6 + ( 89 249480000000 + 199 √ 5 1247400000000 ) w8 +O ( w10 ) , Ā32 (w) = 98209 2714430 − √ 5 75 + ( 136849 108577200 + 476881 √ 5 2714430000 ) w2 + ( 11338513 57003030000 + 4109663 √ 5 57003030000 ) w4+( 79811269 2052109080000 + 1733239283 √ 5 102605454000000 ) w6 + ( 184081799677 22573199880000000 + 409862403109 √ 5 112865999400000000 ) w8 +O ( w10 ) , χ̄2 (w) = 1 + 1 240000 ( 5 + √ 5 )4 w4 − 1 720000000 ( 5 + √ 5 )6 w6 + 1 4032000000000 (5 + √ 5)8w8 +O ( w10 ) , χ̄3 (w) = 1 + ( − 121393 21715440 + 59569 √ 5 108577200 ) w4 + ( − 1019767 2035822500 − 214129 √ 5 1628658000 ) w6+( − 394695169 4560242400000 − 163481333 √ 5 4560242400000 ) w8 +O ( w10 ) , Ã21 (w) = 3 20 + √ 5 20 + ( − 7 1200 − √ 5 400 ) w2 + ( 1 10000 + √ 5 22500 ) w4 + ( − 47 50400000 − √ 5 2400000 ) w6+( 41 7560000000 + 11 √ 5 4536000000 ) w8 +O ( w10 ) , Ã32 (w) = 3 20 − √ 5 20 + ( 1 240 + √ 5 400 ) w2 + ( 11 10000 + 41 √ 5 90000 ) w4 + ( 2281 10080000 + 241 √ 5 2400000 ) w6+( 361961 7560000000 + 97009 √ 5 4536000000 ) w8 +O ( w10 ) , χ̃2 (w) = 1− 1 600 ( 5 + √ 5 )2 w2 + 1 1200000 ( 5 + √ 5 )4 w4 − 1 5040000000 ( 5 + √ 5 )6 w6+ 1 36288000000000 (5 + √ 5)8w8 +O ( w10 ) , χ̃3 (w) = 1 + ( 1 20 + √ 5 60 )2 w2 + ( −29 √ 5− 75 −15000 + 3000 √ 5 ) w4 + ( −683 √ 5− 1560 −1575000 + 315000 √ 5 ) w6+ 1 226800000 −105209 √ 5− 235985 −5 + √ 5 w8 +O ( w10 ) . (58) Z. Sun et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6042 45 of 45 d̄2 (w) = − √ 5 48 + 1 16 + ( √ 5 504000 − 1 201600 ) w4 + ( − √ 5 181440000 − 1 181440000 ) w6+( √ 5 66528000000 − 17 39916800000 ) w8 +O ( w10 ) , d̄3 (w) = √ 5 48 + 1 16 + ( − √ 5 504000 − 1 201600 ) w4 + ( √ 5 181440000 − 1 181440000 ) w6+( − √ 5 66528000000 − 17 39916800000 ) w8 +O ( w10 ) , d̃2 (w) = − √ 5 24 + 5 24 + 1 252000 √ 5w4 + ( − √ 5 90720000 + 1 6048000 ) w6 + ( √ 5 33264000000 + 1 362880000 ) w8 +O ( w10 ) , d̃3 (w) = √ 5 24 + 5 24 − 1 252000 √ 5w4 + ( √ 5 90720000 + 1 6048000 ) w6 + ( − √ 5 33264000000 + 1 362880000 ) w8 +O ( w10 ) . (59)