EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 5662 ISSN 1307-5543 – ejpam.com Published by New York Business Global Insights in Inspecting Two-Degree-of-Freedom of Mathieu-Cubic-Quintic Duffing Oscillator Asma Alanazy1, Galal M. Moatimid2, and Mona A.A. Mohamed2,∗ 1 Department of Mathematics, Northern Border University, Arar, Kingdom of Saudi Arabia 2 Department of Mathematics, Faculty of Education, Ain Shams University, Roxy, Cairo, Egypt Abstract. Mathieu–cubic-quintic Duffing oscillator exemplifies a fundamental model of a para- metric stimulated system characterized by nonlinearities. This issue is essential in understanding intricate, nonlinear dynamical behavior, including bifurcations, chaos, and resonance phenomena, pertinent to engineering and physics for systems influenced by periodic external incentives. A uni- versal two-degree-of-freedom of Mathieu–cubic-quintic Duffing oscillator is examined. This study utilizes an adaptation of the non-perturbative approach, to encompass the coupled system and inspect a parametric nonlinear oscillatory system and evaluate its efficacy. The non-perturbative approach is based mainly on the He’s frequency formula. The main aim of this technique is to convert any nonlinear ordinary differential equation into a linear one. The approximate solutions are derived independently of conventional perturbation methods, excluding the series expansion. Consequently, the objective of this study is to depart from traditional perturbation methods and derive estimated solutions for tiny amplitude parametric components without restrictions. Further- more, this technique is extended to ascertain optimal responses for the nonlinear large amplitude fluctuations. The capacity to swiftly assess the frequency-amplitude correlation is essential for deriving successive estimates of the solutions to parametric nonlinear fluctuations. In addition, the discussion of the coupled oscillators is considered as a new direction in the usage of the current approach. The Mathematica Software is employed to verify the agreement parametric equation, which demonstrates significant concordance with the original equation. The stability behavior is examined in multiple instances. The existing methodology is characterized by clear principles, ap- propriateness, user-friendliness, and exceptional numerical precision. The current method reduces mathematical complexity, making it advantageous in addressing nonlinear parametric problems. It likewise discussed the influence of multiple coefficients on solution behavior, ensuring stability through time-dependent plots and PolarPlots, and detecting periodic behavior with parameter variations. 2020 Mathematics Subject Classifications: 34C15, 34C25, 34D20, 37D45, 70K30, 70H06 Key Words and Phrases: Nonlinear Mathieu oscillator, Cubic-qunitic Duffing oscillator, Cou- pled nonlinear system, Parametric excitation, He’s frequency formula, Non-perturbative approach ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.5662 Email addresses: Asma.Jazza2@gmail.com (A. Alanazy), gal moa@edu.asu.edu.eg (G. M. Moatimid), monaali@edu.asu.edu.eg (M. A. A. Mohamed) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 2 of 36 1. Introduction Dynamical systems with time-periodic coefficients are commonly considered as para- metric stimulated systems. They have been studied since the eighteenth century to elu- cidate various wave processes in fluids and solids. Currently, parametric excitation is fundamental in numerous domains of research and engineering, such as mechanical res- onators [1], optics [2], microwave systems [3], and atomic physics [4]. The prototypical and most fundamental mathematical representation of these parametrically energized sys- tems is found in the Mathieu equation. Consequently, considerable focus was directed towards these systems in examining the stability of their trivial solutions. In a nonlinear system, parametric resonance may result in steady-state non-trivial solutions and limit cycles, where the amplitudes are determined by the system’s nonlinear characteristics. The Duffing oscillator (DO) is a conspicuous example of nonlinear ordinary differential equations (ODEs), particularly in simulating nonlinear systems [5], including micro and nanomechanical resonators. Particular attention is devoted to the parametrically stim- ulated DO, sometimes referred to as the Mathieu–DO [6, 7]. The parametric ODEs of the Mathieu type, whether nonlinear or linear, have attracted significant attention due to their numerous potential applications. Numerous characteristics of the Mathieu ODE can be inferred from the overarching theory of ODEs with periodic coefficients. Mathieu’s ODEs are prevalent in various domains of engineering, physics, and applied mathematics. This renowned equation necessitates careful examination and comprehensive investigation across several contexts due to its extensive applications in physics and chemistry. The theory of ODEs with periodic coefficients, sometimes referred to as flux theory, produces specific elements of the Mathieu ODE. Many properties of the Mathieu ODE have been established through the general theory of periodic coefficients in ODEs. Scientists have long employed various perturbation methods to investigate the analytical solutions related to these equations [8-11]. The perturbation methodology, according to the many-scales approach, is predicated on the concept of a minor parameter [12]. The pursuit of a correct analytical solution to nonlinear problems was a study focus since nonlinear ODEs proved more challenging to solve than their linear counterparts. Accordingly, the perturbation approach is the exclusive method in analyzing any linear ODE subjected to periodic forces. Moreover, the existence of a nonlinear state in ODEs presents many difficulties, especially when periodic components are incorporated into the nonlinear state. One of the primary objectives of this work is to employ an alternative method to attain an estimated solution without utilizing any existing perturbation techniques. The perturbation approaches are challenging for obtaining an approximate solution. This analysis focuses on Mathieu’s nonlinear oscillator, with/without damping forces. Disregarding the ramifications of non- linearity, it can be elucidated using perturbation systems, especially due to the presence of periodic activities. The introduction of a unique method in achieving a quasi-perfect response is regarded as a substantial achievement that advantages both technological and physical fields. Inspired by the aforementioned theoretical discoveries and the importance of para- metric excitation techniques in micro-and nano-applications, the two-degree-of-freedom A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 3 of 36 (TDOF) of The Mathieu–cubic-quintic Duffing oscillator (MDO) is examined in a general configuration, particularly where the broadband destabilization impact is demonstrated. Parametric stimulated nonlinear TDOF systems were previously examined in various con- figurations. Besides cubic nonlinearities, certain systems featuring either quadratic or coupling nonlinear components were analyzed analytically. The TDOF’s response to a combination of parametric resonance, influenced by quadratic nonlinearity and two-to-one internal resonances was examined [13]. The normal forms approach was employed to an- alyze the nonlinear response of a TDOF system exhibiting recurrent natural frequencies and cubic nonlinearity under principal parametric excitation [14]. The linear component of the system exhibits a non-semi-simple one-to-one resonance. The findings were utilized for the interest analysis of a simply supported panel in supersonic airflow. The multiple time- scales method [7, 12] is employed to derive a first-order uniform expansion, resulting in four first-order nonlinear ODEs that govern the modulation of the amplitudes and phases of the two modes. The steady-state responses and their stability are calculated for specific values of the system parameters. The impacts of altering the internal resonance, adjusting the parametric resonance, the phase and amplitude of the second mode parametric excita- tion, and the beginning conditions were examined. Certain studies concentrated on TDOF systems for particular applications utilizing numerical or experimental techniques, incor- porating various nonlinear terms according to the individual application. A novel efficient numerical method for the stability analysis of linear systems with periodic parameters was proposed [15]. The methodology was founded on the premise that the state vector and the periodic matrix of the system can be expressed using Chebyshev polynomials across the major period. The method integrates Floquet theory to produce the transition matrix at the conclusion of one period and establish stability requirements through an eigenvalue analysis procedure. Transient vibrations in mechanical systems represent a prevalent issue in engineering. Numerous theoretical investigations have demonstrated both analytic and computational analysis that a vibrating system can be stabilized or its vibrations dimin- ished when subjected to excitation around a certain parametric combination resonance frequency. In this procedure, transitory vibrations are efficiently mitigated using para- metric excitation. The fundamental stage in utilizing this technology is its experimental application in mechanical systems. This review discusses recent experiments conducted on a basic chain mass system, a continuous cantilever, and a flexible rotor system [16, 17]. Despite the nonlinearity’s extensive practical uses, researchers and theoretical physi- cists have found it difficult to obtain an exact solution or one that closely approximates the correct solution. They typically appear to possess options, even if their drive stems from their objectives. As a consequence, asymptotic analyses of several nonlinear ODEs were the focus of numerous mathematicians. The smallest factor method and the averag- ing technique were employed to illustrate weak nonlinear equations [18]. The Homotopy Perturbation Method (HPM) is employed to achieve accurate asymptotic computations in scenarios involving low-intensity sounds. Solutions of oscillation systems are also ob- tained through the application of multiple time-scales method. Regrettably, utilizing this minimal quantity in both methods yielded inconsistent results [12]. In any asymptotic or perturbation method, identifying the small parameter necessary for a more practical A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 4 of 36 and realistic expression of the fundamental equations is a crucial beginning step. There was a heightened appeal of numerous HPM-based methodologies for forecasting various nonlinear ODEs, hence approximating their solutions more closely [19]. If the beginning assumption did not align with the employed approach for the inquiry, the procedure would diverge and fail to yield the intended results. These strategies depended on the initial esti- mation of the solution. The HPM is employed in the analysis of analytical approximations for magnetic spherical pendulums [20]. Moatimid et al. [21-24] employed the HPM in an- alyzing different problems in Fluid Mechanics as well as Dynamical Systems. Despite the evident advancements, the HPM demonstrated difficulties in use with non-conservative oscillators. The enhanced analytical capacity of the HPM for nonlinear vibration theory is assured. It concentrated on interconnected damping nonlinear oscillators across sev- eral categories. The examples within each category may function as heuristic elucidations or as prototypes for further applications. He’s frequency formula (HFF) is a straight- forward and efficient method for constructing a conservation nonlinear oscillator in the context of nonlinear oscillator challenges [25]. The establishment of HFF is credited to the Chinese mathematician Prof. He. A study of the DO in periodic vibrational behav- ior was conducted under generalized starting conditions [26]. The HFF is determined to be mathematically simple, biologically insightful, and practically beneficial by numerical evaluation. Engineers may swiftly and precisely assess nonlinear vibration systems uti- lizing the HFF through the distinctive methodology outlined in the research. To modify the HFF creatively, it is proposed to categorize the oscillators into two distinct extremes [27]. There is substantial concordance when comparing the approximate and exact fre- quencies across various amplitudes. The Hamiltonian function derived from the HFF has garnered significant attention due to its simplification of the computational analysis of a complex nonlinear vibration system. The cubic-quintic differential operator was utilized to illustrate the precision and simplicity of the calculation. A direct approach for resolving the cubic-quintic differential operator is introduced [28]. The method provided a highly effective and relatively precise means to compute the frequency of a nonlinear conservative oscillator. An exhibition and demonstration of the optimized HFF for nonlinear oscillators are provided [29]. A direct frequency formula in fractal systems is presented and derived from the HFF [30]. The direct calculation and dependable results amalgamated to produce a valuable instrument for a comprehensive investigation of fractal vibration phenomena. The non-Newtonian fluids are essential in various areas, including technology and indus- try. Consequently, the study of these fluids is highly intriguing. The nonlinear stability analysis concentrated on several non-Newtonian fluids. An illustration was created for the creation of some dynamical systems with a flat disturbed interface. The main aim of this theoretical analysis is to utilize the non-perturbative approach (NPA) to augment our recent findings [30-44]. An overview of the actual implementation and physical applications of the MDO can be summarized in the following points: A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 5 of 36 Pragmatic Execution: 1. The MDO is frequently applied in systems that demonstrate both parametric and nonlinear characteristics. This necessitates the development of an oscillator that is: i. They are influenced by periodic or parametric forcing, resulting in Mathieu-like behavior. ii. The restoring force incorporates nonlinear components, reflecting real-world sys- tems where higher-order effects gain prominence. 2. The restoring force of the oscillator is not proportional to displacement and in- cludes nonlinear components. This leads to phenomena such as bi-stability or multi- stability when the system can exhibit several stable and unstable states. 3. To emulate the parametric characteristic, the system is generally exposed to peri- odic or time-dependent forcing, frequently in contexts such as mechanical vibrations or electrical circuits. This forcing interacts with the nonlinear terms to generate intricate dynamical phenomena such as chaotic motion, bifurcations, and resonance. 4. Due to the complexity of the MDO, numerical simulations are frequently necessary to evaluate its behavior. This entails integrating the equations of motion, frequently employing techniques such as the Runge-Kutta method, to analyze the system’s response to various beginning circumstances, forcing frequencies, and amplitudes. Physical Applications: 1. This model describes the motion of beams or plates subjected to periodic external loads. In these instances, the system demonstrates nonlinear vibrations resulting from both parametric excitation and inherent material nonlinearities. 2. The oscillator is applicable to specific electromechanical devices, including microelec- tromechanical systems and nanomechanical systems, where nonlinearities develop from geometric limitations or material characteristics. For instance, it can delineate the movement of a micro-beam subjected to electrostatic actuation. 3. In nonlinear optics, the MDO can simulate the dynamics of light propagation across nonlinear media, whereas the parametric driving may be from medium modulations or external influences such as fluctuating laser intensities. 4. In structural dynamics, it can describe structures or bridges subjected to time- varying stresses, such as wind or earthquakes, when nonlinear behavior significantly influences the system’s reaction to external excitations. 5. Certain biological systems, such as brain oscillations or rhythmic physiological pro- cesses, can be represented by this type of oscillator when they display nonlinear, time-dependent characteristics. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 6 of 36 Key Phenomena: 1. The interplay of parametric forcing and nonlinearities can lead the system to display complex dynamic behavior, including bifurcations and chaos. 2. The parametric driving can induce resonance at particular frequencies, resulting in an amplified system response, as well as sub-harmonic oscillations, when the apparatus fluctuates at a fraction of the controlling frequency. Recently, Barakat et al. [45] studied a similar problem to MDO in the case of cubic DO. In their work, they examined each of the stability of the trivial solution-perturbation anal- ysis of the nonlinear system, via the multiple-time scales-principle parametric resonance- internal resonance under parametric excitation and the stability analysis-combination parametric resonances. Finally, the non-resonant limit cycles are discussed. The NPA of a coupled nonlinear oscillator like the Mathieu oscillator is included in this paper, which builds on our earlier work [46] and is considered as a new direction in the procedure of the NPA. Therefore, two coupled linear ODEs are obtained as compared with the nonlinear ones. To crystalize the presentation of the article, the rest of the paper is organized as follows: A brief description of the NPA of the coupled system is presented in § 2. The NPA of the coupled system is obtainable in § 3. In this Section, an excellent between the nonlinear in addition to the linear ODEs is announced. The damped TDOF of the MDO configuration is accessible in § 4. The un-damped TDOF of the MDO configuration is offered in § 5. The concluding remarks are summarized in § 6. 2. A Brief Description of NPA in the Coupled System This section is devoted to developing the previous procedure of the NPA for a single nonlinear ODE as previously given [30-44 and 46] considering the coupled nonlinear ODEs. Simply, this coupled system may be written as follows: ẍ+ f1(x, ẋ, ẍ, y, ẏ, ÿ) + f2(x, ẋ, ẍ, y, ẏ, ÿ) = f3(x, ẋ, ẍ, y, ẏ, ÿ), (1) and ÿ + h1(y, ẏ, ÿ, x, ẋ, ẍ) + h2(y, ẏ, ÿ, x, ẋ, ẍ) = h3(y, ẏ, ÿ, x, ẋ, ẍ), (2) where f1 as well as h1represent odd functions that relate to the damping forces, f2 and h2 refer to some other odd functions. Away from the damping forces and the restoring forces, the functions f3 and h3 signify even functions. As well from the perturbation approaches [7, 12]; the odd functions produce secular terms. In contrast, the even functions do not yield any secular terms. For straightforwardness, the given procedure is not greatly different from the case of single linear ODE. As previously revealed, the aim is to transform the coupled nonlinear Eqs. (1) and (2) to another equivalent linear equation in the form of the simple harmonic motion. For this purpose, consider two guessing (trial) solutions of the given nonlinear ODEs that are given in Eqs. (1) and (2) as follows: A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 7 of 36 ũ(t) = A cosΩ1t, and ṽ(t) = B cosΩ2t, (3) where Ω1 and Ω2 characterize the total frequencies that are determined later. The solutions given in Eq. (3) obey the following standard initial conditions (ICs): ũ(0) = A, ˙̃u(0) = 0, ṽ(0) = B, and ˙̃v(0) = 0. (4) Now, the goal here is to obtain the corresponding linear ODEs. For this aim, following the aforementioned works [30-44] for the single nonlinear ODE. The equivalent damping term is given as follows: σeqv1 = ∫ 2π/Ω1 0 ˙̃u f1(ũ, ˙̃u, ¨̃u, ṽ, ˙̃v, ¨̃v)dt/ ∫ 2π/Ω1 0 ˙̃u2dt = g1(Ω1, Ω2). (5) It would be mentioned that the function g1 includes all the constants that appear in the function f1 in addition to the unknown frequencies. On the other hand, the equivalent frequency may be evaluated as follows: ϖ2 eqv1 = ∫ 2π/Ω1 0 ũ f2(ũ, ˙̃u, ¨̃u, ṽ, ˙̃v, ¨̃v)dt/ ∫ 2π/Ω1 0 ũ2dt = g2(Ω1, Ω2). (6) The function g2 incorporates all the parameters that perform the function f2 and uniden- tified frequencies that exist within it. Concerning the quadratic function f3, actually, the corresponding term represents a non- homogeneous part. As previously shown [30-44 and 46], this non-homogeneous part can be determined by the following direct substitutions: Replace: u → A 2 √ 2 , u̇ → AΩ1 2 √ 2 , ü → AΩ2 1 2 √ 2 , v → B 2 √ 2 , v̇ → BΩ2 2 √ 2 , and v̈ → BΩ2 2 2 √ 2 . The non-homogeneous part of the corresponding linear ODE may be estimated as follows: λ1 = f3 ( A 2 √ 2 , AΩ1 2 √ 2 , AΩ2 1 2 √ 2 , B 2 √ 2 , BΩ2 2 √ 2 , BΩ2 2 2 √ 2 ) . (7) Finally, the corresponding linear ODE that corresponds to Eq. (1) may be formulated as follows: ü+ σeqv1u̇+ϖ2 eqv1u = λ1. (8) In accordance with the standard normal form, the middle term in Eq. (8) may be removed via the transformation: u(t) = ũ(t)Exp(−σeqv1t/2), (9) where F1(t) is an arbitrary function. The new function F1(t) will be determined via the linear ODE as given in Eq. (8): ¨̃u+Ω2 1ũ = λ1Exp(σeqv1/2), (10) where Ω2 1 = ϖ2 eqv1 − 1 4 σ2 eqv1, (11) A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 8 of 36 which signifies the square of the total frequency of the first linear ODE. Additionally, the stability criteria are evaluated through the following: Ω2 1 > 0, and σeqv1 > 0. (12) Using similar arguments that are previously given, the second comparable linear ODE is given as: v̈ + σeqv2v̇ +ϖ2 eqv2v = λ2. (13) As previously shown the parameters σeqv2 and ϖ2 eqv2 are determined as shown in Eqs. (5) and (6). Remember that the integration is evaluated over the interval 0 → 2π/Ω2. Once more, the middle term in Eq. (13) could be ignored via the following normal form: v(t) = ṽ(t)Exp(−σeqv2t/2), (14) The unknown function F2(t) is determined as: ¨̃v +Ω2 2ṽ = λ2Exp(σeqv2t/2), (15) where Ω2 2 = ϖ2 eqv2 − 1 4 σ2 eqv2, (16) which denotes the square of the total frequency of the second linear ODE. Furthermore, the stability requirements are assessed using the following: Ω2 2 > 0, and σeqv2 > 0. (17) To this end, it should be mentioned that the total frequencies Ω1 and Ω2are evaluated via the Mathematica software (MS) during the command FindRoot for the coupled equations (11) and (16). 3. Significance of the Problem The TDOF of MDO is an extension of the classical DO model that describes a more complex and nonlinear dynamical system of TDOF. This oscillator typically arises in the study of nonlinear vibrations, resonance phenomena, and stability analysis in various engineering and physical systems. In what follows, some features of the given problem will be raised. 1. Mathieu term supplements parametric excitation to the ODE, which leads a system to be non-autonomous. Parametric excitation indicates that some system parameters differ with time, frequently in sinusoidal waves, giving a more complex dynamic performance. This term is a reason for phenomena such as parametric resonance. 2. The cubic-DO, which presents a restoring force that is proportional to the cube of the displacement, makes the system’s behavior depend strongly on amplitude. The Duffing term can result in hardening or softening stiffness features. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 9 of 36 3. The quantic-DO signifies a higher-order restoring force. This enhances additional complexity, permitting more applicable dynamical behavior, such as bifurcations, multiple stable equilibria, and chaotic motion under definite conditions. 4. Mathieu’s term presents the probability of parametric resonance, where definite mo- tivating frequencies lead to enormous amplitude oscillations. 5. The mixture of cubic and quintic nonlinearities leads to enhancement dynamic be- haviors, such as bifurcations, chaotic motion, and limit cycles. 6. The TDOF systems may be combined in several ways, which increases the complexity of the system’s structure. Some Application 1. Mechanical systems (suspension systems, pendulums with elastic supports). 2. Structural engineering (vibration of beams, plates, and shells). 3. Control systems and robotics. 4. Nonlinear optics (wave propagation in nonlinear media). 4. Damped TDOF of MDO Configuration The current study investigates a general version of the damping TDOF of MDO, which deliberates the synchronous parametric excitation in a TDOF nonlinear system with a filled excitation matrix. The synchronous parametric excitation denotes the state where each degree of freedom includes the same or harmonically related frequencies that excite all degrees of freedom simultaneously. This indicates that the parametric excitation terms influencing both degrees of freedom in a TDOF of MDO are driven by the same or synchronized periodic function, usually a sinusoidal function of time. This opposes the asynchronous parametric excitation which indicates that the excitation acting on the two degrees of freedom occurs at different frequencies or phases as discussed by [38]. What follows, are some applications of parametric excitation: 1. Light transmission in nonlinear media can be affected by synchronous parametric stimulation in optical systems. 2. Synchronous parametric excitation can be utilized in structural engineering vibration control systems. 3. Systems that maximize energy reduction by using parametric resonance, in which the excitation indicates enormous oscillations. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 10 of 36 In order to accomplish that, the nonlinear TDOF of the MDO system is suggested as follows [45]: ẍ+ ω2 1x+ µ1ẋ+ γ1x 3 + λ1x 5 + α1ẋ 3 + η11x cosΩpt+ η12y cosΩpt = 0, (18) and ÿ + ω2 2y + µ2ẏ + γ2y 3 + λ2y 5 + α2ẏ 3 + η21x cosΩpt+ η22y cosΩpt = 0, (19) where ω1, and ω2 are the natural frequencies of the two oscillators, µ1, and µ2 are the linear damping coefficients, γ1, and γ2 are the nonlinear cubic coefficients, λ1, and λ2 are the nonlinear quintic coefficients, α1, and α2 are the nonlinear damping coefficients, η11, η21, η12, and η22 are the joint excited amplitudes of the two oscillators, and Ωp is the excited frequency. The suggested ICs of Eq. (18) and (19) are regularly assumed as [30-44 and 46]: x(0) = A, y(0) = B, ẋ(0) = 0, and ẏ(0) = 0. (20) Eqs. (18) and (19) could be re-written as: ẍ+ F1(x, ẋ, ẍ, y, ẏ, ÿ) + F2(x, ẋ, ẍ, y, ẏ, ÿ) = 0, (21) and ÿ +H1(y, ẏ, ÿ, x, ẋ, ẍ) +H2(y, ẏ, ÿ, x, ẋ, ẍ) = 0, (22) where F1(x, ẋ, ẍ, y, ẏ, ÿ)and H1(x, ẋ, ẍ, y, ẏ, ÿ) are the odd secular terms in the two oscillators, while F2(x, ẋ, ẍ, y, ẏ, ÿ) andH2(x, ẋ, ẍ, y, ẏ, ÿ) are the damping secular terms in them. These terms are defined as: F1(x, ẋ, ẍ, y, ẏ, ÿ) = ω2 1x+ γ1x 3 + λ1x 5 + η11x cosΩpt+ η12y cosΩpt, (23) F2(x, ẋ, ẍ, y, ẏ, ÿ) = µ1ẋ+ α1ẋ 3, (24) H1(x, ẋ, ẍ, y, ẏ, ÿ) = ω2 2y + γ2y 3 + λ2y 5 + η21x cosΩpt+ η22y cosΩpt, (25) and H2(x, ẋ, ẍ, y, ẏ, ÿ) = µ2ẏ + α2ẏ 3. (26) The trial solutions of the system as given in Eqs. (21) and (22) are assumed, as per the NPA method, as: ũ = A cosΩ1t, ṽ = B cosΩ2t, (27) where A and B denote the amplitudes of the equivalent solutions, whereas Ω1 and Ω2 are the total frequencies to be estimated later. The functions u and v satisfy the equivalent system:{ ü+ ξ1u̇+ϖ2 1u = 0 v̈ + ξ2v̇ +ϖ2 2v = 0 . (28) Through the inspection of this system, one can expect that the solutions’ behavior for their equations has damping manners throughout the explored time interval. The factors A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 11 of 36 ω1, and ω2 and ξ1, and ξ2 are known by equivalent frequencies and equivalent damping coefficients, which are defined as follows: (i) The Equivalent Frequencies To determine the equivalent frequency, the HFF technique can be utilized with the odd terms F1(x, ẋ, ẍ, y, ẏ, ÿ) and H1(x, ẋ, ẍ, y, ẏ, ÿ) as: ϖ2 1 = ∫ 2π/Ω1 0 ũF1(ũ, ˙̃u, ¨̃u, ṽ, ˙̃v, ¨̃v) dt/ ∫ 2π/Ω1 0 ũ2dt = a1 + ω2 1 + a2 sin ( 2πΩp Ω1 ) + a3 sin ( 2π(Ω2 − Ωp) Ω1 ) + a4 sin ( 2π(Ω2 − Ωp) Ω1 ) , (29) ϖ2 2 = ∫ 2π/Ω1 0 ũH1(ũ, ˙̃u, ¨̃u, ṽ, ˙̃v, ¨̃v) dt/ ∫ 2π/Ω1 0 ũ2dt = a5 + ω2 2 + a6 sin ( 2πΩp Ω2 ) + a7 sin ( 2π(Ω1 − Ωp) Ω2 ) + a8 sin ( 2π(Ω1 − Ωp) Ω2 ) , (30) where the coefficients a1 → a8 are defined through the Appendix. (ii) The Equivalent Damping: The calculation of the equivalent damping coefficients utilizing the HFF and by employing the odd terms F2(x, ẋ, ẍ, y, ẏ, ÿ) and H2(x, ẋ, ẍ, y, ẏ, ÿ) through the integrations: ξ1 = ∫ 2π/Ω1 0 ũF2(ũ, ˙̃u, ¨̃u, ṽ, ˙̃v, ¨̃v) dt/ ∫ 2π/Ω1 0 ˙̃u2dt = µ1 + 3 4 A2α1Ω 2 1, (31) ξ2 = ∫ 2π/Ω1 0 ṽH2(ũ, ˙̃u, ¨̃u, ṽ, ˙̃v, ¨̃v) dt/ ∫ 2π/Ω1 0 ˙̃v2dt = µ2 + 3 4 B2α2Ω 2 2. (32) In light of this, it should be noted that the corresponding frequencies ϖ1, and ϖ2 and the damping coefficients ξ1, and ξ2 are calculated as functions of the total frequencies Ω1, and Ω2 using the MS; hence, these forms are relegated to the Appendix because of their complex formulas. Assume the following standard normal forms: u(t) = ũ(t)Exp(−ξ1 t/2). (33) v(t) = ṽ(t)Exp(−ξ2 t/2). (34) Therefore, substituting by Eqs. (33) and (34) into Eqs. (28), one gets: ¨̃u(t) + Ω2 1ũ(t) = 0, (35) and ¨̃v(t) + Ω2 2ṽ(t) = 0, (36) A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 12 of 36 where Ω2 1 = ϖ2 1 − ξ21/4, (37) Ω2 2 = ϖ2 2 − ξ22/4. (38) The goal of the current method is to comprehend an auxiliary type of the inferred non- linear ODEs. The following Figs and Tables demonstrate high convergence between these two responses, confirming the methodology’s accuracy and insight. It is important to note that the nonlinear ODEs as given in Eqs. (18) and (19) are replaced by alternating linear ODEs. As given in Eq. (28) utilizing the NPA. The (NDSolve) command is used in conjunction with the MS to compare the numerical solutions (NSs) of these two systems in view of the above values of the system coefficients. The NSs of the coupled system, as provided in Eqs. (18), (19), and the NPA solution, as provided by Eqs. (28) are compared for convenience in Figs. 1 and 2 and Tables 1 and 2. The nonlinear stability analysis was the focus of the first model. This new method is thought to be comprehensive, striking, consistent, easy to use, and effective in setting itself apart from the more recent methods. It may be used to examine a number of dynamical system categories, such as hydrodynamic stability. As can be seen, there is a significant convergence between the two results in these figures and Tables. As estimated by the MS, the assessed absolute error equals 0.0064 for x distribution and 0.0021 for y distribution in view of the following data: µ1 = 0.2, µ2 = 0.5, α1 = 0.1, α2 = 0.3, ω1 = 2, ω2 = 1.5, γ1 = 0.3, γ2 = 0.3, η11 = 0.1, η12 = 0.2, η21 = 0.1, η22 = 0.2, Ωp = 0.5, λ1 = 0.2, λ12 = 0.3, A = 0.05 and B = 0.05. Figure 1: Demonstrates the matching between the solutions of the nonlinear dynamical system (18) and (19) and those of the NPA system (28) for the oscillator. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 13 of 36 Table 1: Validates the convergence of the original and NPA solutions for the x oscillator. t Real oscillator x Approximate oscillator Absolute error 0 0.05 0.05 0. 5 -0.025963 -0.0228042 0.00315881 10 0.00792326 0.00153441 0.00638885 15 0.00310068 0.00702038 0.00391969 20 -0.00489663 -0.00667517 0.00177854 25 0.00412394 0.00323986 0.000884078 30 -0.00238374 -0.000375066 0.00200867 35 0.00098347 -0.000860581 0.00184405 40 -0.0000360658 0.000886748 0.000922813 45 -0.000237634 -0.000456704 0.00021907 50 0.000278283 0.0000727985 0.000205485 Figure 2: Demonstrates the matching between the solutions of the nonlinear dynamical system (18) and (19) and those of the NPA system (28) for the oscillator. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 14 of 36 Table 2: Validates the convergence of the original and NPA solutions for the y oscillator. t Real oscillator y Approximate oscillator Absolute error 0 0.05 0.05 0. 5 0.00631459 0.0058608 0.000453797 10 -0.00238827 -0.00323727 0.000849002 15 -0.000359652 -0.000941673 0.000582021 20 0.000675203 0.000129051 0.000546152 25 0.000208416 0.0000952715 0.000113144 30 -0.00019148 3.12676 ∗ 10−6 0.000194607 35 -0.000129185 −7.32034 ∗ 10−6 0.000121864 40 0.0000458095 −1.30373 ∗ 10−6 0.0000471132 45 0.0000670557 4.10079 ∗ 10−7 0.0000666456 50 0.0000201414 1.65056 ∗ 10−7 0.0000199764 The solutions of the system of Eqs. (35) and (36) using the NPA are graphically repre- sented in Figs. 3-16, when the above-considered values of the parameters are used. Firstly, Figs. 3-8 are plotted to examine the influences of the relevant factors ω1, ω2, γ1, γ2, η11, η21, η12, η22, λ1, λ2, and Ωp on the time-dependent variation and the polar plots of u(t) and v(t). These parameters take the abovementioned values that vary from one figure to an- other with the variation of each studied parameter. Examining the waves of time history, it is apparent that these waves exhibit damping periodic shapes across the whole duration. The oscillation numbers for both solutions vary based on the factors that are in effect, as seen by the graphs. The effects of the natural frequencies ω1 and ω2 are shown through Figs. 3. It is seen from this figure that the amplitudes of the waves for the two solutions continue damping, while the number of vibrations increases significantly with the growth of ω1 and ω2. Furthermore, the wavelengths are shrinking with the rise of these two frequencies. An acceleration in the waves is observed also due to their growth, this acceleration is greater for u(t) than v(t) which is due to the bigger of ω1 than ω2. The presence of quintic and cubic nonlinearities in the Mathieu system, especially at large oscillation amplitudes, alters naturally the actual frequencies. This is referred to as amplitude-dependent frequency or nonlinear frequency shift. The NPA is used in the current study to realize these new frequencies, as seen above. In a fully linear system, the natural frequency remains constant; however, in a nonlinear system, such as the DO, the nonlinear restoring forces cause the effective natural frequency to vary with the oscillation amplitude. Additionally, it is seen from Figs. 4 that the amplitudes of the waves for the two solutions are damped, whereas the number of vibrations is approximately constant with the growth of the cubic factors γ1 and γ2. The wavelengths are found unvarying with the impact of γ1 and γ2, meanwhile, a notable delay in the waves is obtained in Figs. 4 with the rise of these cubic coefficients. Systems with cubic nonlinearities do not have a set natural frequency of oscillation, in contrast to linear systems. Instead, the amplitude starts to depend on the inherent frequency. Intricate behaviors and improved system dynamics are A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 15 of 36 provided by cubic nonlinearities, which are not possible in simply linear systems. Figs. 5 are plotted to examine the influence of the coefficients of the periodic terms of x in the nonlinear system of Mathieu Eqs. (18) and (19), η11 and η21, on both u(t) and v(t), respectively. From these figures, one sees that the amplitudes of the waves are decreasing, meanwhile the wavelengths remain slightly stable for u(t) and v(t) with the rise of η11 and η21. Moreover, Figs. 6 are schemed to illustrate the effects of the coefficients of the periodic terms of y in the nonlinear system of Mathieu Eqs. (18) and (19), η12and η22, on both u(t) and v(t), individually. These figures indicate similar impacts as Figs. 4. Physically, parametric resonance is made possible by the rise of the excited amplitudes. This can lead to nonlinear phenomena such as high amplitude oscillations, amplitude- dependent resonance, and perhaps instability. Higher stimulation amplitudes extend the instability zones and increase the system’s susceptibility to parametric resonance and enormous oscillations. Figs. 7 are plotted to inspect the influence of the quinitic coefficients λ1 and λ2 in the coupled system (18) and (19). The amplitudes of the waves as well as the number of vibra- tions for the two solutions declined with the growth of the quinitic factors λ1 and λ2. The wavelengths are found unvarying with the impact of λ1and λ2, meanwhile, an insignificant delay in the waves is obtained in Figs. 7 with the rise of these quintic coefficients. The dynamic behavior of a TDOF of MDO is further complicated by quintic nonlinearities, particularly when paired with parametric excitation and other nonlinear effects like cubic terms. Higher-order stiffness effects brought forth by the quintic term affect the stability, bifurcation structures, and oscillations of the system. At higher oscillation amplitudes, the quintic term changes the behavior of the system and may be combined with cubic nonlinearities and parametric stimulation to produce complex dynamics. Similar to cubic nonlinearities, quintic nonlinearities cause the system’s natural frequencies to change as the amplitude grows. At larger amplitudes, the frequency-amplitude dependency becomes more pronounced due to the higher-order frequency shift introduced by the quintic term. At this end, Figs. 8 show the impact of the excited frequency Ωp, of the two waves x and y which characterizes the two coupled oscillators on the distributions u(t) and v(t). It is seen also from these figures that wavelengths and the number of waves remain unvarying, while significant progress in the waves is obtained in Figs. 8 with the rise of Ωp. The excited frequency in a TDOF of MDO is crucial in determining how the system reacts to parametric stimulation. Bifurcations, stability, oscillation amplitude, the beginning of parametric resonance, and the shift to complex and chaotic dynamics are all influenced by the excited frequency. The excited frequency Ωp has a significant impact on the system’s response amplitude. The large oscillations can result from even tiny parametric stimulation around resonant frequencies. The magnitude of the oscillatory response of the system is determined by the connection between Ωp and its inherent frequencies. Based on the aforementioned findings, one might conclude that these outcomes can aid in the development and regulation of systems that rely on periodic solutions. For example, engineers may utilize periodic behavior in mechanical systems to improve system perfor- mance or design systems that operate effectively at particular frequencies. The interaction between parametric excitation, excitation frequency, and nonlinearities results in compli- A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 16 of 36 cated behavior of the system. The system may exhibit bifurcations, nonlinear frequency changes, enhanced oscillations, or even chaotic motion close to resonances. Especially with large oscillations, quintic terms improve the amplitude dependency of these effects. Figure 3: The effect of the natural frequencies ω1 and ω2 on the time-dependent variation of u(t) and v(t) . Figure 4: The effect of the cubic coefficients γ1 and γ2 on the time-dependent variation of u(t) and v(t) . Figure 5: The effect of the excited amplitudes η11 and η21 on the time-dependent deviation of u(t) and v(t) . A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 17 of 36 Figure 6: The effect of the excited amplitudes η12 and η22 on the time-dependent deviation of u(t) and v(t) . Figure 7: The effect of the quinitic coefficients λ1 and λ2 on the time-dependent distribution of u(t) and v(t). . Figure 8: The effect of the excited frequency Ωp on the time-dependent distribution of u(t)andv(t) . A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 18 of 36 Figs. 9-14 indicate the PolarPlots of the suggested solutions as given in Eq. (27) through the time interval [0, 50π] with the distinction of the parameters ω1, ω2, µ1, µ2, α1, α2, λ1, and λ2 to illustrate the polarized behavior of the functions u(t) and v(t) as defined in Eq. (27) in terms of the total frequencies Ω1 and Ω2, respectively. According to the same above data except that A = B = 0.5, such data vary from one figure to another according to the studied parameter. Figs. 9 and 10 are plotted due to the changed values of the natural frequencies ω1 and ω2. It is exposed from these figures that, with the rising of ω1 and ω2, the curvilinear circles regulate and collect in a more uniform style about the core. These outcomes are approved by the stimulating effects of the parameters ω1 and ω2 that appear in the Figs. 3. Moreover, a symmetric distribution encircles the centers of the rounded, connected curves. This dispersion may be used to determine the stable form in which these curves operate. The impact of the modifying elements determines how much the intersecting plots circulate. Figs. 11 and 12 are plotted discussing the variations of the damping coefficients µ1 andµ2, to illustrate the PolarPlot of the functions u(t) and v(t) as signified in Eq. (27) utilizing the different parameters and the total frequencies Ω1and Ω2. It is detected that the impacts of µ1 and µ2 oppose the influence of ω1 and ω2 , where the loops are dispersed more about the center, but the behavior is still uniformly stable. This model has a symmetric distribution around the centers of the circular, linked curves. These curves could function in a steady state, as indicated by this distribution. The influence of the influencing factors causes the circulation of the connected plots to diminish. Furthermore, the PolarPlots of the recommended solutions u(t) and v(t) conferring to Eq. (28) are considered for altered values of the nonlinear damping coefficients α1 and α2 as detected in Figs. 13 and 14. These illustrations comprehend the PolarPlots of the equal solutions with a small variation in the nonlinear damping coefficients. One can notice that the impacts of α1 and α2 are similar to those of the linear ones µ1 and µ2, approximately they give the same distributed circular curves. These created curves are exposed to spin in repeated pathways, resulting in remarkable spin structures that are symmetrical around their center, which validates the dominance of the stable mode in the current model with the applicable importance of this conclusion. Due to the similarity between the cubic and quintic coefficients impacts, as seen before, it is adequate to study the influence of one of them. Consequently, Figs. 15 and 16 are graphed to illustrate the impacts of the quintic coefficients λ1 and λ2 , as new terms added in the current study contrary to the previous [38]. It is noticed that the number of loops is growing with the rise of λ1 and λ2, which indicates that the current system in general verifies the stability manner. All these curves have a variety of closed or semi-closed elliptical trajectories that distribute steadily around the center point giving these amazing shapes. Because of the similarity of data used to create some figures, it is found that some parts of them are comparable. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 19 of 36 Figure 9: displays the PolarPlot of u(t) for different measures of the x-oscillator natural frequency ω1. Figure 10: displays the PolarPlot of v(t) for different measures of the y-oscillator natural frequency ω2. Figure 11: displays the PolarPlot of u(t) for different measures of the x-oscillator damping coefficient µ1. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 20 of 36 Figure 12: displays the PolarPlot of v(t) for different measures of the y-oscillator damping coefficient µ2. Figure 13: displays the PolarPlot of u(t) for different measures of the x-oscillator cubic coefficient α1. Figure 14: displays the PolarPlot of v(t) for different measures of the y-oscillator cubic coefficient α2. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 21 of 36 Figure 15: displays the PolarPlot of u(t) for different measures of the x-oscillator quinitic coefficient λ1. Figure 16: displays the PolarPlot of v(t) for different measures of the y-oscillator quinitic coefficient λ2. 5. The Un-Damped TDOF of MDO Inspecting the conventional situation of the coupled system as given in Eqs. (18) and (19), in the case of nonexistence of damping terms (µ1 = µ2 = α1 = α2 = 0). The absence of damping in the un-damped situation of a TDOF of MDO prevents energy from dissipating over time. This lack of damping considerably changes the dynamics of the system, resulting in more noticeable nonlinear effects due to the existence of the cubic and quintic terms, magnified reactions close to resonance, and persistent oscillations. In this case, if the system is close to resonance conditions, continuous large-amplitude oscillations will result from oscillations in both degrees of freedom. Furthermore, strong coupling and possible internal resonance effects result from the free flow of energy between modes. For these important aspects, we are going to illustrate the TDOF of MDO without the damping parameters. The coupled harmonic MDO (18) and (19) in this case is identified as follows: ẍ+ ω2 1x+ γ1x 3 + λ1x 5 + η11x cosΩpt+ η12y cosΩpt = 0, (39) and ÿ + ω2 2y + γ2y 3 + λ2y 5 + η21x cosΩpt+ η22y cosΩpt = 0, (40) A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 22 of 36 As usual, the suggested ICs of Eq. (39) and (40) are supposed as: x(0) = C, y(0) = D, ẋ(0) = 0, ẏ(0) = 0 = α1 = 0 (41) Eqs. (39) and (40) could be reconstructed as: ẍ+ F1(x, ẋ, ẍ, y, ẏ, ÿ) = 0, (42) and ÿ +H1(y, ẏ, ÿ, x, ẋ, ẍ) = 0, (43) where F1(x, ẋ, ẍ, y, ẏ, ÿ) and H1(x, ẋ, ẍ, y, ẏ, ÿ) are the odd secular terms in the two oscillators, which are defined previously in Eqs. (23) and (25). The trial solutions of the system of Eqs. (21) and (22) are assumed, as per the NPA method, as: w̃ = C cosΩ3t, g̃ = D cosΩ4t, (44) where C and D represent the amplitudes of the equivalent solutions in this case, whereas Ω3 and Ω4 are the total frequencies to be valued later. Moreover, the functions w and g satisfy the corresponding system: { ẅ +ϖ2 3w = 0 g̈ +ϖ2 4g = 0 . (45) One may anticipate that the solutions’ behavior for the system equations (39) and (40) will exhibit an un-damping form across the investigated period of time based on this system’s examination. The elements ϖ3 and ϖ4 are identified by equivalent frequency coefficients, which are defined as follows: ϖ2 3 = ∫ 2π/Ω3 0 w̃F1(w̃, ˙̃w, ¨̃w, g̃, ˙̃g, ¨̃g) dt/ ∫ 2π/Ω3 0 w̃2dt = b1 + ω2 1 + b2 sin ( 2πΩp Ω3 ) + b3 sin ( 2π(Ω4 − Ωp) Ω3 ) + b4 sin ( 2π(Ω4 − Ωp) Ω3 ) , (46) ϖ2 4 = ∫ 2π/Ω4 0 g̃H1(w̃, ˙̃w, ¨̃w, g̃, ˙̃g, ¨̃g) dt/ ∫ 2π/Ω4 0 g̃2dt = b5 + ω2 2 + b6 sin ( 2πΩp Ω4 ) + b7 sin ( 2π(Ω3 − Ωp) Ω4 ) + b8 sin ( 2π(Ω3 − Ωp) Ω4 ) , (47) where the coefficients b1 → b8 are defined through the Appendix. It should be noted that the equivalent frequencies ϖ3, and ϖ4 are expected to be equal to the total frequencies Ω3, and Ω4 here due to the absence of the damping coefficients. In this case, there is no need to utilize the standard normalized technique and the two equations (45) will be scrutinized directly. The objective of the current case is to follow the above procedures to study other linear ODEs instead of the non-linear ODEs described by Eqs. (39) and (40). This equivalent A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 23 of 36 system as defined by (45) is investigating utilizing the forms (46) and (47) and the following Figs and Tables validate the extreme matching between the two solutions, approving the methodology’s correctness and awareness. Given the aforementioned values of the system coefficients, the numerical solutions of these two systems are compared using the (NDSolve) command in combination with the MS. Figs. 17 and 18 and Tables 3 and 4 are designed to achieve this convergence. As can be noticed through these Figs. and Tables, there is a noteworthy conjunction between the numerical solutions of the original and equivalent systems. As estimated by the MS, the assessed absolute error equals 0.0108 for x distribution and 0.0135 for y distribution in view of the following data: ω1 = 2, ω2 = 1.5, γ1 = 0.3, γ2 = 0.3, η11 = 0.1, η12 = 0.2, η21 = 0.1, η22 = 0.2,Ωp = 0.5, λ1 = 0.2, λ12 = 0.3, A = 0.05 and B = 0.05. Figure 17: Demonstrates the matching between the solutions of the nonlinear dynamical system as given in Eqs. (39) and (40) and those of the NPA system (45) for the un-damped x oscillator. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 24 of 36 Table 3: Validates the convergence of the original and NPA solutions for the un-damped x oscillator. t Real oscillator Approximate oscillator Absolute error 0 0.05 0.05 0. 5 -0.0436541 -0.0425569 0.00109724 10 0.0188486 0.0224435 0.00359486 15 0.00600942 0.00435191 0.00165751 20 -0.0347072 -0.0298516 0.0048556 25 0.0479366 0.0464638 0.00147285 30 -0.0492413 -0.0492425 1.22066 ∗ 10−6 35 0.030144 0.0373605 0.00721643 40 -0.00746833 -0.0143553 0.00688697 45 -0.0230729 -0.0129238 0.010149 50 0.0419213 0.0363552 0.00556605 Figure 18: Demonstrates the matching between the solutions of the nonlinear dynamical system as given in Eqs. (39) and (40) and those of the NPA system (45) for the un-damped y oscillator. Table 4: Validates the convergence of the original and NPA solutions for the un-damped y oscillator. t Real oscillator Approximate oscillator Absolute error 0 0.05 0.05 0. 5 0.0171072 0.0171218 0.0000145838 10 -0.0418844 -0.0382738 0.00361066 15 -0.0495195 -0.0433344 0.00618507 20 0.00346055 0.00859528 0.00513473 25 0.0481834 0.0492211 0.00103769 30 0.0296473 0.0251148 0.00453244 35 -0.0316531 -0.0320206 0.000367506 40 -0.0541948 -0.0470448 0.00714994 45 -0.0105252 -0.00019903 0.0103262 50 0.0428419 0.0469085 0.00406667 A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 25 of 36 The solutions of the un-damped oscillators’ Eqs. (39) and (40) utilizing the NPA are graphically represented in Figs. 19-30, when the previously measured data of the pa- rameters are considered. Initially, Figs. 19-24 are schemed to inspect the impacts of the pertinent factors ω1, ω2, γ1, γ2, η11, η21, η12, η22, λ1, λ2, and Ωp, the same parameters dis- cussed in the damped case, on the time deviation and the PolarPlots of w(t) and g(t). These parameters use the previously stated values, which change depending on how each parameter under study changes from one picture to another. It is evident from analyzing the time history waves that they display damping periodic forms over the whole duration. The graphs show that the oscillation numbers for both solutions differ according to the factors at play. Fig. 19 illustrates the impact of the natural frequencies ω1 and ω2 on the time distribution of w(t) and g(t), respectively. These figures illustrate how the wave amplitudes for the two solutions remain constant but the number of vibrations rises sharply as ω1 and ω2 grows. Additionally, when these two frequencies increase, the wavelengths get smaller. Because of their growth, waves are also seen to accelerate. The presence of quintic and cubic nonlin- earities in the Mathieu system, especially at large oscillation amplitudes, alters naturally occurring frequencies. This is referred to as amplitude-dependent frequency or nonlinear frequency shift. The NPA is used in the current study to realize these new frequencies, as seen above. In a nonlinear system, such as the DO, nonlinear restoring forces cause the effective natural frequency to change with the oscillation amplitude, whereas the natural frequency remains constant in a fully linear system. Furthermore, Figs. 20 show that the number of vibrations and amplitudes of the waves for the two solutions are roughly constant as the cubic factors γ1 and γ2 are increased. The influence of γ1 and γ2 is shown to have little effect on the wavelengths, but the growth of these cubic coefficients in Figs. 4 results in a noticeable delay in the waves. Unlike linear systems, systems with cubic nonlinearities do not have a fixed natural frequency of oscillation. Rather, the intrinsic frequency begins to determine the amplitude. The complicated behavior and improved system dynamics are provided by cubic nonlinearities, which are not possible in simply linear systems. Figs. 21 and 22 are plotted to examine the influence of the coefficients of the excited am- plitudes of x and y, η11, η21,η12, and η22, in the nonlinear system of non-damped Mathieu Eqs. (39) and (40), on both w(t) and g(t), respectively. From these figures, one can see that the amplitudes of the waves are stay stable, meanwhile, the wavelengths are shrink- ing for w(t) and g(t) with the rise of η11, η21, η12, and η22. The oscillations are swinging between delays and enhance with the varying of these parameters. Physically, parametric resonance is made possible by the growth in stimulated amplitudes. Nonlinear phenom- ena including amplitude-dependent resonance, large amplitude oscillations, and perhaps instability can result from this. The instability zones are expanded by higher stimula- tion amplitudes, which also make the system more vulnerable to massive oscillations and parametric resonance. To examine the impact of the quintic coefficients λ1 and λ2 on the linked system as given in Eqs. (39) and (40), Fig. 23 is shown. Both the wave amplitudes and the number of vibrations for the two solutions stay constant as these quinitic factors increase. The A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 26 of 36 influence of and is shown to have little effect on the wavelengths, but the growth of these quintic coefficients in Figs. 23 results in a negligible delay in the waves. Quintic nonlinear- ities complicate the dynamic behavior of a TDOF of MDO, especially when combined with additional nonlinear effects such as cubic terms and parametric excitation. The quintic terms introduce higher-order stiffness effects that impact the system’s oscillations, bifur- cation structures, and stability. The quintic term can interact with cubic nonlinearities and parametric excitation to produce complex dynamics, and they change the behavior of the system at greater oscillation amplitudes. Similar to cubic nonlinearities, quintic non- linearities cause the system’s natural frequencies to alter as the amplitude rises. At larger amplitudes, the frequency amplitude dependency is more apparent due to the higher-order frequency shift produced by the quintic components. Moreover, Figs. 24 illustrate how the distributions are affected by the excited frequency Ωp of the two un-damped waves x and y and that characterizes the two linked oscillators by the functions w(t) and g(t). These figures likewise show that the number of waves and wavelengths do not change, however, Fig. 24 shows a notable improvement in the waves with the increase of Ωp. A TDOF of MDO’s excited frequency plays a critical role in shaping the system’s response to parametric stimulation. The excited frequency affects bifurcations, oscillation amplitude, stability, the onset of parametric resonance, and the transition to complex and chaotic dynamics. The system’s response amplitude is significantly influenced by the stimulated frequency Ωp. Around resonant frequencies, even little parametric stimulation can produce large oscillations. The relationship between and the system’s intrinsic frequencies determines the size of its oscillatory response. Based on the facts presented above and as stated in the damped instance, it is possible to realize the conclusion that these results can help with the construction and management of systems that depend on periodic solutions. To improve system performance or design systems that operate well at particular frequencies, engineers could, for instance, leverage periodic behavior in mechanical systems. Nonlinearities, parametric excitation, and exci- tation frequency interact to produce the system’s complicated behavior. The system may exhibit bifurcations, nonlinear frequency changes, heightened oscillations, or even chaotic motion close to resonances. The quintic term improves the amplitude dependency of these effects, especially when there are large oscillations. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 27 of 36 Figure 19: The effect of the natural frequencies ω1 and ω2 on the time-dependent variation of w(t) and g(t). Figure 20: The effect of the cubic coefficients γ1 and γ2 on the time-dependent variation of w(t) and g(t). Figure 21: The effect of the excited amplitudes η11 and η21 on the time-dependent deviation of w(t) and g(t). A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 28 of 36 Figure 22: The effect of the excited amplitudes η12 and η22 on the time-dependent deviation of w(t) and g(t). Figure 23: The effect of the quinitic coefficients λ1 and λ2 on the time-dependent distribution of w(t) and g(t). Figure 24: The effect of the excited frequency Ωp on the time-dependent distribution of w(t) and g(t). Once more, Figs. 25-30 designate the PolarPlots of the proposed solutions as given in Eq. (45) over the time interval[0, 50π] with the difference of the factors ω1, ω2, α1, α2, λ1, and λ2 to demonstrate the polarization performance of the functions w(t) and g(t) as given in Eq. (45) in terms of the total frequencies Ω3 and Ω4, respectively. Considering the identical data but taking C = D = 0.5, which differ from one figure to another rendering to the deliberated stricture. According to the altered values of the natural frequencies ω1 and ω2, Figs. 25 and 26 are displayed. These figures reveal that as ω1 and ω2 rise, the curvilinear circles control and A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 29 of 36 gather around the center more consistently. These results are supported by the parameters’ stimulating effects that are shown in Fig. 19. In this idea, a symmetric distribution encircles the centers of the rounded, connected curves. This dispersion may be used to determine the stable form in which these curves operate. The impact of the modifying elements determines how much the intersecting plots circulate. Additionally, the PolarPlots of the recommended solutions w(t) and g(t) conferring to Eqs. (45) are considered for altered values of the coefficients of the excited amplitude η11 and η21 as detected in Figs. 27 and 28. These illustrations comprehend the PolarPlots of the equivalent solutions w(t) and g(t) with a small variation in the excited amplitudes coefficients of the x oscillator, where the effects of the other two excited amplitudes of the yoscillator η12 and η22 are excluded here due to their similarity with the effects of η11 and η21. One can notice that the impacts of the excited amplitudes are similar for both w(t) andg(t), almost they give the same dispersed circular curves. When these generated curves are subjected to swirl in recurring tracks, amazing gyration constructions are constructed, which are regular about their center. It is sufficient to examine the impact of one of the cubic and quintic factors because of their similarities, as previously shown. In order to show the effects of the quinitic coefficients λ1 and λ2, a new term included in the current research in contrast to the previous one [38], Figs. 29 and 30 are graphed. The number of loops is observed to increase as and rises, indicating that the existing system generally verifies the stability method. All of these curves have a variety of closed or semi-closed elliptical trajectories that, when they progressively spread out around the center point, produce these amazing shapes. Because the data used to create the figures is similar, some of its components are deemed comparable. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 30 of 36 Figure 25: displays the PolarPlot of w(t) for different measures of the x-oscillator natural frequency ω1. Figure 26: displays the PolarPlot of g(t) for different measures of the y-oscillator natural frequency ω2. Figure 27: displays the PolarPlot of w(t) for different measures of the x-oscillator excited amplitude coefficient η11. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 31 of 36 Figure 28: displays the PolarPlot of g(t) for different measures of the excited amplitude coefficient η21. Figure 29: displays the PolarPlot of w(t) for different measures of the x-oscillator quinitic coefficient λ1. Figure 30: displays the PolarPlot of g(t) for different measures of the y-oscillator quinitic coefficient λ2. 6. Concluding Remarks The MDO served as a foundational model in a parametrically excited system that distinguished by the nonlinearities. This matter is crucial for understanding complex, nonlinear dynamical behavior, including bifurcations, chaos, and resonance phenomena, A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 32 of 36 relevant to engineering and physics for systems affected by periodic external stimuli. A universal TDOF model of dynamic optimization was analyzed. This work employed an application of the NPA to analyze a coupled system of MDO and assess the efficacy of a parametric nonlinear oscillatory system in the damped and non-damped cases. The NPA primarily relies on the HFF. Therefore, the current study aims to move away from con- ventional perturbation techniques and provide approximated solutions for small amplitude parametric components without limitations. Moreover, the technique is expanded to de- termine optimal solutions for nonlinear large amplitude fluctuations. The MS was utilized to validate the resulting parametric equation, which exhibits substantial agreement with the original equation. The current methodology was defined by explicit principles, suit- ability, ease of use, and remarkable numerical accuracy. The existing approach diminishes mathematical complexity, rendering it beneficial for tackling nonlinear parametric issues. The main key outcomes of the current may be summarized as follows: 1. The NPA was adapted to find equivalent linear ODEs. An excellent matching be- tween the original nonlinear system and the corresponding linear one was found and validated through some graphs and error tables. 2. The stability of the various states was guaranteed by time-dependent plots and PolarPlots of these equivalent solutions with the variation of the multiple pertinent coefficients. 3. Regular reduction conduct in the time histories of these solutions indicates stable modes. 4. Additionally, the PolarPlots curves displayed continued uniform spirals meeting at a central point. 5. Periodic behaviors were detected in all cases with the variation of several parameters. The current study will be enhanced in the subsequent work to include the stability con- ditions profile, the resonance, the bifurcation, and the chaotic behavior of coupled and multiple systems. Acknowledgements The authors extend their appreciation to the Deanship of Scientific Research at North- ern Border University, Arar, KSA for funding this research work through the project number ”NBU-FFR-2024-912-04” Author Contributions Asma Alanazy: Resources, Methodology, Formal Analysis, Validation, Visualization and Reviewing. Galal M. Moatimid: Conceptualization, Resources, Methodology, Formal analysis, Val- idation, Writing-Original draft preparation, Visualization and Reviewing. A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 33 of 36 Mona A. A. Mohamed: Investigation, Methodology, Data duration, Examination, Organization, Validation, Reviewing and Editing. All authors have read and agreed to the published version of the manuscript. Data Availability Statement All data generated or analyzed during this study are included in this article. Conflicts of Interest: The authors have disclosed no competing interests. References [1] Rhoads J.F., Guo C., and Fedder G.K., Parametrically excited micro- and nanosys- tems. In: Resonant MEMS, Principles, Modeling, Implementation, and Applications, pp. 73–95. Wiley-VCH Verlag (2015). [2] Cerullo G., De Silvestri S., Ultrafast optical parametric amplifiers, Review of Scientific Instruments, 74(1), 1–18 (2003). [3] Aumentado J., Superconducting Parametric Amplifiers: The State of the Art in Josephson Parametric Amplifiers, IEEE Microwave Magazine, 21(8), 45–59 (2020). [4] Gemelke N., Sarajlic E., Bidel Y., Hong S., and Chu S., Parametric amplification of matter waves in periodically translated optical lattices, Physical Review Letters, 95(17), 170404 (2005). [5] Rhoads J.F., Shaw S.W., and Turner K.L., Nonlinear dynamics and its applications in micro-and nanoresonators, Journal of Dynamic Systems, Measurement, and Control, 132(3), 1–14 (2010). [6] Kovacic I., Rand R., Sah S.M., Mathieu’s equation and its generalizations: Overview of stability charts and their features, Applied Mechanics Reviews, 70(2), 020802 (2018). [7] Nayfeh A.H., and Mook D.T., Nonlinear Oscillations, Wiley & Sons, New York (1979). [8] McLachlan N.W., Theory and Applications of Mathieu Functions, Clarendon Press, Oxford, UK, (1947). [9] Achala L.N., Mathematical analysis and applications of Mathieu’s equation revisited, International Journal of Mathematics And its Applications, 9(2), 49-54 (2021). [10] Kovacic I., Rand R., and Sah S.M., Mathieu’s equation and its generalizations: Overview of stability charts and their features, Applied Mechanics Reviews, 70, 020802 (22 Pages) (2018). [11] Ramani D.V., Keith WL, and Rand R.H.. Perturbation solution for secondary bi- furcation in the quadratically-damped Mathieu equation, International Journal of Non-Linear Mechanics, 39 (3), 491-502 (2004). [12] Nayfeh A.H., Introduction to perturbation techniques, New York, Wiley & Sons (1981). [13] Nayfeh A.H., and Zavodney L.D., The response of two-degree-of-freedom systems with quadratic non-linearities to a combination parametric resonance, Journal of Sound and Vibration, 107(2), 329–350 (1986). [14] Nayfeh A.H., Chin C., and Mook D.T., Parametrically excited nonlinear two-degree- A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 34 of 36 of-freedom systems with repeated natural frequencies, Shock and Vibration, 2(1), 43–57 (1995). [15] Sinha S.C., and Wu D.H., An efficient computational scheme for the analysis of periodic systems, Journal of Sound and Vibration, 151(1), 91–117 (1991). [16] Dohnal F., Experimental studies on damping by parametric excitation using electro- magnets, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 226(8), 2015–2027 (2012). [17] Zaghari B., Dynamic analysis of a nonlinear parametrically excited system using electromagnets, Ph.D. Thesis, University of Southampton (2016). [18] Ji W.M., Wang H., and Liu M., Dynamics analysis of an impulsive stochastic model for spruce budworm growth, Applied Mathematics and Computations, 19, 336–359 (2021). [19] He J.H, Homotopy perturbation technique, Computer Methods in Applied Mechanics and Engineering, 178, 257–262 (1999). [20] Moatimid G.M., and Amer T.S. Analytical approximate solutions of a magnetic spher- ical pendulum: Stability analysis, Journal of Vibration Engineering & Technologies, 11, 2155–2165 (2023). [21] Ghaleb A.F., Abou-Dina M.S., Moatimid G.M., and Zekry M.H., Analytic approxi- mate solutions of the cubic-quintic Duffing Van der Pol equation with two-external periodic forces terms: Stability analysis, Mathematics and Computers in Simulation, 180, 129-151 (2021). [22] Moatimid G.M., El-Dib Y.O., and Zekry M.H., Stability analysis using multiple scales homotopy approach of coupled cylindrical interfaces under the influence of periodic electrostatic fields, Chinese Journal of physics, 56(5), 2507-2522 (2018). [23] Moatimid G.M., Sliding bead on a smooth vertical rotated parabola: Stability con- figuration, Kuwait Journal of Science, 47(2), 6-21 (2020). [24] Moatimid G.M., and Amer T.S., Analytical solution for the motion of a pendulum with rolling wheel: Stability Analysis, Scientific Reports, 12(1), 12628 (2022). [25] Zhang J-G. , Song Q-R. , Zhang J-Q., and Wang F., Application of He’s frequency formula to nonlinear oscillators with generalized initial conditions, Facta Universitatis Series: Mechanical Engineering, 21(4), 701-712 (2023). [26] Ma H., A short remark on He’s frequency formulation, Journal of Low Frequency Noise, Vibration and Active Control, 41(4), 1380–1385 (2022). [27] Ma H., Simplified Hamiltonian-based frequency-amplitude formulation for nonlinear vibration systems, Facta Universitatis Series Mechanical Engineering, 20(2), 445-455 (2022). [28] He J-H., The simplest approach to nonlinear oscillators, Results in Physics, 15, 102546 (2019). [29] He C-H., and Liu C., A modified frequency-amplitude formulation for fractal vibration systems, Fractals, 30 (03), 2250046 (2022). [30] Niu J-Y., Feng G-Q., and Gepreel K.A., A simple frequency formulation for frac- tal–fractional non-linear oscillators: A promising tool and its future challenge, Fron- tiers in Physics, 11, 1158121 (2023). A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 35 of 36 [31] Moatimid G.M., and Amer T.S., Dynamical system of a time-delayed -Van der Pole oscillator: A non-perturbative approach, Scientific Reports, 13, 11942 (2023). [32] Moatimid G.M., Amer T.S., and Ellabban Y.Y., A novel methodology for a time- delayed controller to prevent nonlinear system oscillations, Journal of Low Frequency Noise, Vibration and Active Control, 43(1), 525-542 (2024). [33] Moatimid G.M., Amer T.S., and Galal A.A., Studying highly nonlinear oscillators using the non-perturbative methodology, Scientific Reports, 13, 20288 (2023). [34] Moatimid G.M., Mohamed M. A.A., and Elagamy Kh., Nonlinear Kelvin-Helmholtz instability of a horizontal interface separating two electrified Walters’ B liquids: A new approach, Chinese Journal of Physics, 85, 629-648 (2023). [35] Moatimid G.M., and Sayed A., Nonlinear EHD stability of a cylindrical interface separating two Rivlin-Ericksen fluids: A Novel analysis, Chinese Journal of Physics, 87, 379–397 (2024). [36] Moatimid G.M., El-Sayed A.T., and Salman H.F., Different controllers for suppressing oscillations of a hybrid oscillator via non-perturbative analysis, Scientific Reports, 14, 307 (2024). [37] Moatimid G.M., and M. Mohamed Y.M., A novel methodology in analyzing nonlinear stability of two electrified viscoelastic liquids, Chinese Journal of Physics, 89, 679-706 (2024). [38] Moatimid G.M., and M. Mohamed Y.M., Nonlinear electro-rheological instability of two moving cylindrical fluids: An innovative approach, Physics of Fluids, 36, 024110 (2024). [39] Moatimid G.M., Mostafa D.M., and Zekry M.H., A new methodology in evaluating nonlinear Eelectrohydrodynamic azimuthal stability between two dusty viscous fluids, Chinese Journal of Physics, 90, 134–154 (2024). [40] Moatimid G.M., and Mostafa D.M., Nonlinear stability of two superimposed electri- fied dusty fluids of type Rivlin-Ericksen: Non-perturbative approach, Partial Differ- ential Equations in Applied Mathematics, 10, 100745 (2024). [41] Moatimid G.M., Mohamed M. A.A., and Elagamy Kh., Insightful Inspection of the Nonlinear Instability of an Azimuthal Disturbance Separating Two Rotating Magnetic Liquid Columns, The European Physical Journal Plus, 139, 590 (2024). [42] Ismail G.M., Moatimid G.M., Yamani M.I., Periodic solutions of strongly nonlinear oscillators using He’s frequency formulation, European Journal of Pure and Applied Mathematics, 17(3), 2154-2171 (2024). [43] Alluhydan Kh., Moatimid G.M., and Amer T.S., The non-perturbative approach in examining the motion of a simple pendulum associated with a rolling wheel with a time-delay, European Journal of Pure and Applied Mathematics, 17(4), 3185-3208 (2024). [44] Alluhydan Kh., Moatimid G.M., Amer T.S., and A.A. Galal, A novel inspection of a time-delayed rolling of a rigid rod, European Journal of Pure and Applied Mathe- matics, 17(4), 2878-2895 (2024). [45] Barakat A.A., Weig E.M., and Hagedorn P., Non-trivial solutions and their stabil- ity in a two-degree-of-freedom Mathieu–Duffing system, Nonlinear Dynamics, 111, A. Alanazy, G. M. Moatimid, M. A. A. Mohamed / Eur. J. Pure Appl. Math, 18 (2) (2025), 5662 36 of 36 22119–22136 (2023). [46] Moatimid G.M., Mohamed M.A.A., and Elagamy Kh., An innovative ap- proach in inspecting a damped Mathieu cubic–quintic Duffing oscillator, Jour- nal of Vibration Engineering & Technologies, Published online: 23 July (2024). https://doi.org/10.1007/s42417-024-01506-w Appendix The coefficients defined in Eqs. (29) and (28) as evaluated by MS are given as: a1 = 3 4 A2γ1 + 5 8 A4λ1, a2 = η11Ω1(−2Ω2 1 +Ω2 p) π(−4Ω2 1Ωp +Ω3 p) , a3 = Bη12Ω1 2πA ( (−Ω2 +Ωp) Ω2 1 − Ω2 2 + 2Ω2Ωp − Ω2 p ) , a4 = Bη12Ω1 2πA ( (Ω2 +Ωp) Ω2 2 − Ω2 1 + 2Ω2Ωp +Ω2 p ) , a5 = 3 4 B2γ2 + 5 8 B4λ′ 2a6 = η22Ω2(−2Ω2 2 +Ω2 p) π(−4Ω2 2Ωp +Ω3 p) , a7 = Aη21Ω2 2πB ( (Ω1−Ωp) Ω2 1−Ω2 2−2Ω1Ωp+Ω2 p ) , and a8 = Aη21Ω2 2πB ( (Ω1+Ωp) Ω2 1−Ω2 2+2Ω1Ωp+Ω2 p ) . The coefficients defined in Eqs. (46) and (47) as evaluated by MS are given as: b1 = 3 4 C2γ1 + 5 8 C4λ1, b2 = η11Ω3(−2Ω2 3 +Ω2 p) π(−4Ω2 3Ωp +Ω3 p) , b3 = Dη12Ω3 2πC ( (−Ω4 +Ωp) Ω2 3 − Ω2 4 + 2Ω4Ωp − Ω2 p ) , b4 = Dη12Ω3 2πC ( (Ω4 +Ωp) Ω2 4 − Ω2 3 + 2Ω4Ωp +Ω2 p ) , b5 = 3 4 D2γ2 + 5 8 D4λ′ 2b6 = η22Ω4(−2Ω2 4 +Ω2 p) π(−4Ω2 4Ωp +Ω3 p) , b7 = Cη21Ω4 2πD ( (Ω3−Ωp) Ω2 3−Ω2 4−2Ω3Ωp+Ω2 p ) , and b8 = Cη21Ω4 2πD ( (Ω3+Ωp) Ω2 3−Ω2 4+2Ω3Ωp+Ω2 p ) .