Analytical modelling of a nonlinear 2DOF Mathieu

Abstract

The nonlinear Mathieu-Duffing oscillator is a basic model of parametrically driven dynamical systems with nonlinear coefficients, and it serves as an important model for the study of such phenomena as bifurcation, chaotic transitions and resonance instabilities in engineering structures and physical systems subjected to periodic external excitations. In this work, a non-perturbative analytical approach is proposed to analyse a coupled parametric nonlinear oscillatory system and to determine its dynamical features with much more computing efficiency. The proposed methodology is based on He’s frequency formulation, which systematically linearises weakly nonlinear ordinary differential oscillators, avoiding the fundamental limitations of classical perturbation methods and being effective for large amplitude nonlinear oscillatory regimes. A prominent contribution of this work is the extension of the non-perturbative framework to the consideration of linked oscillator systems, which is a unique application domain for the method. The analytical conclusions are validated by symbolic computation in Mathematica. It is found that there is strong conformity with the original governing equations for different parametric configurations. The stability study is conducted under different operating situations, and it always shows the accuracy, analytical tractability, and numerical accuracy of the procedure. The global dynamical behavior of the system is also characterized by bifurcation diagrams, which identify critical parameter-dependent transitions, and the largest Lyapunov exponent, which quantitatively distinguishes between periodic and chaotic oscillatory regimes and defines the threshold conditions for the onset of chaos. The results agree with the conclusion that the non-perturbative technique is a powerful, and efficient for the analysis of parametrically stimulated nonlinear dynamical systems.

Keywords

nonlinear differential equations Mathieu Duffing oscillators parametric excitation He’s frequency formula non-perturbative approach bifurcation diagrams largest Lyapunov exponent

1. Introduction

Dynamical systems whose governing ODEs contain time–periodic coefficients are generally classified as parametrically excited systems. Since the eighteenth century, such systems are extensively studied to understand a wide range of wave phenomena occurring in both fluid and solid media. In contemporary research, parametric excitation has become increasingly important in several areas of science and engineering, including mechanical resonators, optical systems, and microwave devices.¹ One of the simplest and most representative mathematical models used to describe these systems is the Mathieu oscillator. The DO serves as a classical model of nonlinear ODEs and is widely used to describe nonlinear dynamic behavior in many physical systems, particularly micro-mechanical resonators. When parametric excitation was incorporated into the DO, the resulting system is often referred to as the Mathieu–DO.² Both linear/nonlinear ODEs of the Mathieu type attract significant attention due to their wide range of theoretical and practical applications. Several of their important properties can be understood within the general framework of ODEs with periodic coefficients. Over the years, numerous perturbation techniques have been developed and applied to obtain approximate analytical solutions for these equations.^3–6 Among these techniques, the MTSM is commonly employed, as it allows the analysis of weakly nonlinear ODEs by introducing a small parameter into the formulation.⁷ Finding reliable analytical solutions of nonlinear problems has long attracted considerable attention, mainly because nonlinear ordinary differential ODEs are inherently more complex than linear ones. In many cases, perturbation techniques are traditionally applied to study ODEs subjected to periodic excitation. However, when nonlinear effects are present, the analysis becomes considerably more challenging, particularly when the nonlinear terms interact with periodic components of the system. The NPA deviates from all earlier models by integrating and thus rectifying specific shortcomings identified in prior research. This divergence is not only incremental but signifies a significant methodological progress. Moreover, whereas references^1–7 offer a significant foundation, their methodologies do not encompass the entirety of situations and scenarios that NPA seeks to address. We trust this explanation underscores the innovation and importance of the proposed methodology, and we are willing to amend the work to clarify this distinction if required.

The 2DOF Mathieu–DO is significant as it provides more precise representation of complex nonlinear dynamical systems defined by linked oscillations, parametric excitation, and notable nonlinear stiffness effects. Unlike 1DOF models, the 2DOF system can encompass substantial phenomena such as internal resonance, mode coupling, energy transfer, bifurcation, and chaotic motion, commonly seen in engineering structures, MEMS devices, aerospace systems, and mechanical vibrations. The study is essential in understanding stability characteristics and developing suitable vibration control and nonlinear dynamic analysis methods. The fundamentals and applications of nonlinear systems in science and engineering were analyzed by exploring the dynamic behavior of a connected cubic-quintic damping fractal system.⁸ The analytical and numerical approaches to analyze the dynamic properties of coupled dynamical systems, incorporating two periodic excitations and distributed time delays, were investigated.⁹ A 2DOF asymmetric dynamic oscillator with fractional damping was investigated through the analysis of the Jeffcott rotor system with horizontal support, focusing on the dynamic behavior of its resonance and attraction domain structures.¹⁰ A novel methodology in investigating the nonlinear dynamic behaviors of a belt-driven, two-coupled friction oscillator was analyzed.¹¹ Initially, Buckingham’s theorem was employed to generate several dimensionless and physically significant parameters. Their associations with the steady-state reactions of the system were elucidated by numerical simulations. A novel mapping technique called the sticking phase plot was created to document the transition process of stick–slip motion. A 2DOF piecewise nonlinear system model incorporating a fractional differential term was developed.¹² The harmonic balance approach was employed to analyze the 2DOF piecewise nonlinear system with a fractional differential term, resulting in the achievement of the system’s amplitude–frequency curve.

Despite the wide range of practical applications of nonlinear phenomena, deriving an exact analytical solution, or even a close approximation, remains a challenging task of theoretical physicists. Nonlinear ODEs are inherently more complex than the linear ones, which has encouraged the development of various analytical alternatives. Consequently, the asymptotic analysis of nonlinear ODEs has attracted considerable attention from Mathematicians. Techniques such as the small parameter method and the averaging approach were applied to study weakly nonlinear equations.¹³ The HPM was used to obtain accurate asymptotic approximations, particularly in problems involving low-intensity acoustics. In addition, MTSM was widely employed to determine the behavior of oscillatory systems. However, the reliance on a limited number of assumptions or parameters may sometimes lead to inconsistent results.⁷ In most perturbation or asymptotic techniques, identifying an appropriate small parameter is an essential step toward constructing a realistic mathematical representation of the governing equations. The increasing popularity of HPM-based approaches has improved the prediction of solutions for many nonlinear ODEs.¹⁴ For example, HPM was applied to obtain analytical approximations for magnetic spherical pendulum systems¹⁵; meanwhile Moatimid et al.¹⁶ used the method to address various problems in fluid mechanics and dynamical systems. Nevertheless, despite these developments, applying HPM to non-conservative oscillators still presents certain challenges. Further progress was therefore made to enhance the analytical capability of HPM in nonlinear vibration theory, especially for systems of coupled damped nonlinear oscillators. Within this context, HFF was introduced as a simple and efficient technique in studying conservative nonlinear oscillators.¹⁷ This formulation was originally proposed by the Chinese Mathematician Prof. J-H He. The DO was examined in relation to periodic vibration behavior under generalized ICs.¹⁸ To further extend the HFF, oscillators were classified into two limiting cases.¹⁹ The cubic–quintic DO was used to demonstrate the accuracy and computational simplicity of this approach, and a direct technique in analysing this system was proposed.²⁰ This method offers an efficient and reasonably accurate way to estimate the frequency of nonlinear conservative oscillators. Several studies have also presented demonstrations and improvements of the optimized HFF of nonlinear oscillators.^21,22 Moreover, a direct frequency expression of fractal systems was derived based on the HFF, providing a useful tool in investigating fractal vibration phenomena through straightforward calculations and reliable results.

The examination of 2DOF of MDO yields significant experimental and practical understanding of nonlinear dynamics of interconnected mechanical, electrical, and structural systems influenced by periodic and parametric excitations. This model experimentally summarizes the intricate interactions of nonlinearity, coupling, and parametric resonance, allowing researchers to replicate and scrutinize phenomena such as mode coupling, internal resonance, bifurcation, chaos, and amplitude modulation in laboratory environments utilizing coupled beams, mass-spring systems, or MEMS resonators. These experimental validations are essential in enhancing analytical and computational predictions of dynamic behavior under multi-frequency excitations. Insights from this oscillator are crucial of the design and vibration management of aerospace structures, vehicle suspensions, precision instruments, and energy harvesting systems, where the optimization of nonlinear stiffness and damping is essential to prevent instability or improve performance. Moreover, comprehending the parametric domains of stability and chaotic transitions facilitates the formulation of resilient control techniques and adaptive materials that can preserve structural integrity under fluctuating loads. The 2DOF of MDO functions as a theoretical framework and a practical instrument in enhancing nonlinear system analysis and engineering design in several fields. A study closely related to MDO was recently carried out by Barakat et al. 2²³ in the case of the cubic DO. In their work, the stability of the trivial solution was examined through perturbation analysis of the nonlinear system. The analysis relied mainly on the MTSM, while different resonance mechanisms were considered, including parametric resonance and internal resonance under parametric excitation. The authors also investigated the stability criteria that are associated with combined parametric resonances. Additionally, the behavior of the system was discussed in the non-resonant region, where the characteristics of the resulting limit cycles were analyzed. More recently, the NPA has attracted attention in the analysis of dynamical systems. Several studies have reported the use of NPA in stability investigations across different areas of dynamical system research.^24–30

The authors of^31–34 concentrate on innovative advances in nonlinear fractional dynamical systems and mathematical physics, especially in the fields of variational principles, sensitivity analysis, Hamiltonian structures, and soliton dynamics. In addition to examining nonlinear fractional equations, these studies offer analytical methods for obtaining precise and approximative wave solutions with applications in the engineering and physical sciences. Additionally, the works enhance the theoretical foundation of fractal and fractional calculus and promote its applications in contemporary nonlinear research by helping to comprehend the stability and dynamical behavior of complicated fractional systems. Furthermore, analytical and numerical studies on nonlinear interactions and internal resonance in mechanical systems such as n-layer X-shape vibration isolators have advanced the understanding and prediction of nonlinear dynamic performance for engineering applications greatly.³⁵

This paper incorporates the NPA of coupled nonlinear ODEs, specifically the Mathieu oscillators, which extends to this research³⁰ and represents a novel approach in the NPA methodology. Consequently, two coupled linear ODEs are derived in contrast to the nonlinear ones. To elucidate the presentation of the essay, the remainder of the material is structured as follows: § 2 provides a concise overview of the coupled system. The numerical and tabular representations of the theoretical outcomes are delineated in § 3. § 4 presents the BDs, LE and PMs. The final observations are encapsulated in § 5.

2. Problem construction

This Section is devoted to providing a theoretical overview of the model by describing its core assumptions, describing the principal variables and parameters, and outlining the mathematical framework that dictates its behavior. This formulation arranges the groundwork of the ensuing study and offers the essential context in comprehending the model’s characteristics, constraints, and relevance to the issue.

2.1. Governing ODEs

The following ODEs represent a set of coupled, damped, nonlinear oscillators influenced by parametric (time-periodic) stimulation. Such systems frequently occur in mechanical vibrations, MEMS resonators, nonlinear optics, and structural dynamics. We will elucidate the physical significance of each term individually, followed by a comprehensive physical interpretation. The nonlinear coupled system under consideration is given as:

x ″ + ω_{1}^{2} x + μ_{1} x' + α_{4} x'^{3} + α_{1} x^{3} + α_{2} x^{5} + α_{3} x^{7} + α_{5} x \cos (2 ω t) + α_{6} y \cos (2 ω t) = 0,

(1)

and

y ″ + ω_{2}^{2} y + μ_{2} y' + β_{4} y'^{3} + β_{1} y^{3} + β_{2} y^{5} + β_{3} y^{7} + β_{5} y \cos (2 ω t) + β_{6} x \cos (2 ω t) = 0 .

(2)

The physical significance of this equation has been elucidated as follows:

This system visually represents two interconnected vibrating entities that exchange energy while oscillating due to internal restorative forces, damping effects, and periodic external stimulation. The damping terms denote energy loss due to friction or resistance, whereas the restoring forces function to return each body to its equilibrium position. The nonlinear effects indicate that the system’s stiffness fluctuates with displacement, signifying that the response is not proportional to the applied motion, especially at high amplitudes. The interrelationship between the two entities indicates that vibrations in one directly affect the motion of the other, leading to energy transfer and synchronized or complex oscillations. Periodic excitation refers to an external time-dependent energy source, such as rotating machinery, electromagnetic forces, or structural vibrations, that can induce resonance and intricate dynamic behavior. This system can physically replicate consistent mechanical structures, nonlinear vibration dampers, aeronautical frameworks, MEMS systems, or other engineering systems where nonlinear interactions and time-dependent excitations are critical.

Furthermore, equations (1) and (2) contain parametric excitation terms of the form $\cos (2 ω t)$ , which make the system non-autonomous and introduce a Mathieu-type structure. This structure is known to generate instability regions and complicates the application of perturbation. To eliminate the explicit time depend and avoid the Mathieu formulation, the parametric terms must be replaced by equivalent constant coefficients that represent their effective contribution over one oscillation cycle.

1. General structure

The variables $x (t)$ and $y (t)$ represent generalized coordinates (displacements or amplitudes).

2. (a) Inertia $x^{″}$ and $y^{″}$ are standard Newtonian inertia, where mass is normalized to unity.

(b) The linear stiffness (linear restoring forces) are represented by: $ω_{1} x$ and $ω_{2} y$ (represent springs or elastic restoring forces), where $ω_{1}$ and $ω_{2}$ are natural frequencies.

(c) The linear viscous damping is representing by: $μ_{1} x^{'}$ and $μ_{2} y^{'}$ . They are energy dissipation due to friction, are resistance, or material damping.

3. The nonlinear stiffness terms (enharmonic effects) are:

(a) $α_{1} x^{3} + α_{2} x^{5} + α_{3} x^{7}$ and $β_{1} y^{3} + β_{2} y^{5} + β_{3} y^{7}$ . These signify forceful nonlinear restorative forces, originating from: significant distortions, nonlinear material properties, geometric nonlinearities.

Interpretation in a physical context: Cubic word indicates Duffing-type nonlinearity. Quintic and septic words represent higher-order corrections. May induce frequencies contingent upon amplitude jump phenomena, and multi-stability

4. Parametric excitation (time-periodic coefficients):

(a) Direct parametric terms are: $α_{5} x \cos (2 ω t)$ , and $β_{5} y \cos (2 ω t)$ are time-periodic modulation of stiffness. The frequency $2 ω$ signifies parametric resonance. These terms are analogous to: Periodically varying spring stiffness or vertical oscillation of pendulum’s pivot.

(b) Cross-parametric coupling are: $β_{6} x \cos (2 ω t)$ , and $α_{6} y \cos (2 ω t)$

External excitation of one mode affects the other. Simultaneous external field acting on both coordinates.

5. Comprehensive physical interpretation

This system simulates: Two weakly damped, highly nonlinear oscillators exhibiting nonlinear coupling and parametric excitation near double a reference frequency. Common bodily manifestations encompass: Coupled nonlinear beams or plates. They are nonlinear vibration dampers. Mechanically systems with parametric excitation

6. Anticipated Key Phenomena

This setup can demonstrate: Internal resonance.

Finally, these ODEs characterize a 2DOF nonlinear vibrating system with damping, pronounced polynomial stiffness nonlinearities, nonlinear modal coupling, and parametric excitation, encapsulating intricate energy exchange and resonance phenomena observed in sophisticated mechanical and micro-scale systems.

To achieve this, we adopt the introduced transformation; the averaging operator is defined by³⁰:

K = \int_{0}^{π / 2 ω} ψ (t) χ^{2} (t) d t / \int_{0}^{π / 2 ω} χ^{2} (t) d t,

(3)

The periodic stiffness is defined as function $ψ (t)$ . This function represents time-dependent or parametrically varying part of the system’s stiffness. It typically contains terms such as: $ψ (t) = \cos (2 ω t)$ . The weighting function is given as $χ (t)$ . This is a chosen comparison (or test) function used to project periodic stiffness onto an equivalent constant value. In the original formulation, it is assumed as: $χ (t) = \cos t$ . It ensures orthogonality and allows accurate averaging over one oscillation cycle.

After substituting the average parametric coefficients obtained from equation (3) into the original governing system, the parametric excitation terms $\cos (2 ω t)$ are replaced by their corresponding constant effective contributions. Consequently, the system can be rewritten in the following autonomous form:

x ″ + (ω_{1}^{2} + ϖ_{11}^{2}) x + μ_{1} x' + α_{4} x'^{3} + α_{1} x^{3} + α_{2} x^{5} + α_{3} x^{7} + ϖ_{12}^{2} y = 0,

(4)

and

y ″ + (ω_{2}^{2} + ϖ_{21}^{2}) y + μ_{2} y' + β_{4} y'^{3} + β_{1} y^{3} + β_{2} y^{5} + β_{3} y^{7} + ϖ_{22}^{2} x = 0,

(5)

where

ϖ_{1 j}^{2} = \int_{0}^{π / 2 ω} α_{5, 6} \cos (2 ω t) \cos^{2} (t) d t / \int_{0}^{π / 2 ω} \cos (t) d t = \frac{ω \sin (π / ω) α_{5, 6}}{(ω^{2} - 1) (ω \sin (π / ω) + π)}, j = 1, 2

(6)

and

ϖ_{2 j}^{2} = \int_{0}^{π / 2 ω} β_{5, 6} \cos (2 ω t) \cos^{2} (t) d t / \int_{0}^{π / 2 ω} \cos (t) d t = \frac{ω \sin (π / ω) β_{5, 6}}{(ω^{2} - 1) (p \sin (π / ω) + π)} j = 1, 2 .

(7)

Applying the adopted analytical technique to the original nonlinear coupled ODEs with periodic parametric excitation, the system can be reduced to the following form. A graphical comparison is presented in Figure 1, to demonstrate the equivalence between the original and the transformed systems. The following figures are plotted in accordance with the subsequent data sample:

Figure 1.

Shows graphical validation of the transformation performance the agreement between the original systems (equations (1) and (2) red curve) and the reduced system (equations (4) and (5) blue curve).

α_{1} = 0.1, α_{2} = 2.1, α_{3} = 1.1, α_{4} = 0.01, α_{5} = 0.1, α_{6} = 0.1, β_{1} = 1, β_{2} = 0.5, β_{3} = 0.5, β_{4} = 0.3, β_{5} = 0.1, β_{6} = 0.2, ω_{1} = 1.87, ω_{2} = 1.41, μ_{1} = 0.2, μ_{2} = 0.3, ω = 5, A = 0.1, B = 0.2 .

Following the same methodology that was previously adopted,³⁰ the next step is to convert the newly obtained autonomous system into an equivalent linear ODEs. Consequently, the governing ODEs are reduced to the following autonomous linear ones:

x ″ + σ_{1} x' + λ_{1} x + ϖ_{12}^{2} y = 0,

(8)

and

y ″ + σ_{2} y' + λ_{2} y + ϖ_{22}^{2} x = 0 .

(9)

As earlier stated, the objective is to convert the coupled nonlinear ODEs as presented in equations (1) and (2) into an equivalent set of linear ODEs. To achieve this, two assumed (trial) forms of the solution are introduced of the nonlinear ODEs as specified in equations (1) and (2), as follows:

\tilde{x} (t) = A \cos Ω_{1} t, and \tilde{y} (t) = B \cos Ω_{2} t,

(10)

where

Ω_{1}

and

Ω_{2}

represent the resulting frequencies that will be evaluated later.

Equation (10) is considered as a trial solution, which represents an assumed form of the system response used as an initial approximation in the analysis. This form may later be modified if necessary. The expressions in equation (10) are obtained using the following standard ICs:

\tilde{x} (0) = A, \tilde{\dot{x}} (0) = 0, \tilde{y} (0) = B, and \tilde{\dot{y}} (0) = 0 .

(11)

The next step is to derive the corresponding linear ODEs. Following the approach adopted in the previous studies, the equivalent damping term can be expressed as follows³⁰:

σ_{1} = \int_{0}^{2 π / Ω_{1}} \tilde{\dot{x}} (μ_{1} \dot{x} + α_{4} {\dot{x}}^{3}) d t / \int_{0}^{2 π / Ω_{1}} {\tilde{\dot{x}}}^{2} d t = \frac{3}{4} α_{4} A^{2} Ω_{1}^{2} + μ_{1},

(12)

and

σ_{2} = \int_{0}^{2 π / Ω_{2}} \tilde{\dot{y}} (μ_{2} \dot{y} + β_{4} {\dot{y}}^{3}) d t / \int_{0}^{2 π / Ω_{2}} {\tilde{\dot{y}}}^{2} d t = \frac{3}{4} β_{4} B^{2} Ω_{2}^{2} + μ_{2} .

(13)

Conversely, the associated equivalent frequencies can be obtained as follows:

λ_{1}^{2} = \int_{0}^{2 π / Ω_{1}} \tilde{x} [(ω_{1} + ϖ_{11}) x + α_{1} x^{3} + α_{2} x^{5} + α_{3} x^{7}] d t / \int_{0}^{2 π / Ω_{1}} {\tilde{x}}^{2} d t = ϖ_{11}^{2} + ω_{1}^{2} + \frac{3}{4} A^{2} α_{1} + \frac{5}{8} A^{4} α_{2} + \frac{35}{64} A^{6} α_{3},

(14)

and

λ_{2}^{2} = \int_{0}^{2 π / Ω_{2}} \tilde{y} [(ω_{2}^{2} + ϖ_{21}^{2}) y + β_{1} y^{3} + β_{2} y^{5} + β_{3} y^{7}] d t / \int_{0}^{2 π / Ω_{2}} {\tilde{y}}^{2} d t = ϖ_{21}^{2} + ω_{2}^{2} + \frac{3}{4} B^{2} β_{1} + \frac{5}{8} B^{4} β_{2} + \frac{35}{64} B^{6} β_{3} .

(15)

To facilitate a clearer analysis of the system, the objective is to transform the coupled ODEs as agreed in equations (8) and (9) into simplified linear and decoupled forms. This procedure aims to represent the system in a more manageable structure by separating the two variables $x$ and $y$ . Accordingly, equations (8) and (9) can be reformulated by expressing one variable in terms of the other, which allows a more straightforward examination of the dynamical behavior of the system.

x ″ + σ_{1} x' + f_{1} (x, y) x = 0,

(16)

and

y ″ + σ_{2} y' + f_{2} (x, y) y = 0,

(17)

where the nonlinear stiffening functions

f_{1} (x, y)

and

f_{2} (x, y)

are defined as:

f_{1} (x, y) = λ_{1}^{2} + ϖ_{12}^{2} (y / x),

(18)

and

f_{2} (x, y) = λ_{2}^{2} + ϖ_{22}^{2} (x / y) .

(19)

The obtained mean–square expressions provide a practical way to address the complex behavior of such systems by allowing the evaluation of key dynamic parameters and simplifying the nonlinear interactions between the variables. This approach is widely used in the analysis of nonlinear vibration problems, as it simplifies the calculations while still producing reasonably accurate results.

{\bar{x}}^{2} = \frac{1}{4 T_{1}} \int_{0}^{T_{1}} x_{0}^{2} d t = \frac{1}{8} A^{2},

(20)

and

{\bar{y}}^{2} = \frac{1}{4 T_{1}} \int_{0}^{T_{1}} y_{0}^{2} d t = \frac{1}{8} B^{2}

(21)

In this approach, the stiffened functions

f_{1} (x, y)

and

f_{2} (x, y)

will be used to establish the additional national frequencies

ω_{e q v 1}^{2}

and

ω_{e q v 2}^{2}

are defined as follows:

ω_{e q v 1}^{2} = \lim_{x \to \bar{x}, y \to \bar{y}} f_{1} (x, y) = λ_{1}^{2} + ϖ_{12}^{2} (B / A),

(22)

and

ω_{e q v 2}^{2} = \lim_{x \to \bar{x}, y \to \bar{y}} f_{2} (x, y) = λ_{2}^{2} + ϖ_{22}^{2} (A / B) .

(23)

Accordingly, the system can be represented by the following decoupled equations:

x ″ + σ_{1} x' + ω_{e q v 1}^{2} x = 0,

(24)

and

y ″ + σ_{2} y^{'} + ω_{e q v 2}^{2} y = 0 .

(25)

The obtained ODEs are now uncoupled, which simplifies the study of the system dynamics. This separation allows stability and resonance characteristics to be examined more easily. Accordingly, the total frequencies $Ω_{1}$ and $Ω_{2}$ are determined from the linear damped harmonic ODEs as follows:

Ω_{1}^{2} = ω_{e q v 1}^{2} - \frac{1}{4} σ_{1}^{2},

(26)

and

Ω_{2}^{2} = ω_{e q v 2}^{2} - \frac{1}{4} σ_{2}^{2} .

(27)

Employing equations (22) and (23) into equations (26) and (27) then using equations (12) to (15), one gets

Ω_{1}^{2} = ω_{1}^{2} + \frac{3}{4} A^{2} α_{1} + \frac{5}{8} A^{4} α_{2} + \frac{35}{64} A^{6} α_{3} + ω \sin (π / ω) α_{5} / (ω^{2} - 1) (ω \sin (π / ω) + π) + B ω \sin (π / ω) α_{6} / A (ω^{2} - 1) (ω \sin (π / ω) + π) - \frac{1}{4} {(\frac{3}{4} α_{4} A^{2} Ω_{1}^{2} + μ_{1})}^{2},

(28)

and

Ω_{2}^{2} = ω_{2}^{2} + \frac{3}{4} B^{2} β_{1} + \frac{5}{8} B^{4} β_{2} + \frac{35}{64} B^{6} β_{3} + ω \sin (π / ω) β_{5} / (ω^{2} - 1) (ω \sin (π / ω) + π) + A ω \sin (π / ω) β_{6} / B (ω^{2} - 1) (ω \sin (π / ω) + π) - \frac{1}{4} {(\frac{3}{4} β_{4} B^{2} Ω_{2}^{2} + μ_{2})}^{2} .

(29)

We can write:

Ω_{1}^{4} + c_{1} Ω_{1}^{2} + c_{2} = 0,

(30)

and

Ω_{2}^{4} + c_{3} Ω_{2}^{2} + c_{4} = 0 .

(31)

As Ω_{1}^{2} = - \frac{c_{1}}{2} + \frac{\sqrt{c_{1}^{2} - 4 c_{2}}}{2},

(32)

and

Ω_{2}^{2} = - \frac{c_{3}}{2} + \frac{\sqrt{c_{3}^{2} - 4 c_{4}}}{2} .

(33)

where values of

c_{i}, i = 1, 2, 3, 4

are moved to the Appendix.

Accordingly, the solutions of equations (24) and (25) have the form

x (t) = A e^{- \frac{1}{2} μ_{1} t} (\cos (Ω_{1} t) + \frac{μ_{1}}{2 Ω_{1}} \sin (Ω_{1} t)),

(34)

and

y (t) = B e^{- \frac{1}{2} μ_{2} t} (\cos (Ω_{2} t) + \frac{μ_{2}}{2 Ω_{2}} \sin (Ω_{2} t)) .

(35)

3. Results and discussion

In this Section, the dynamic behavior of the system is examined through graphical pictures. The analysis focuses on time histories, phase-plane trajectories, and stability characteristics to provide a clear understanding of how the system response progresses under different parameter settings. These graphical results are used to illustrate the influence of the governing parameters on the overall dynamic behavior, particularly in terms of transient response, steady-state motion, and stability limits. By comparing the responses obtained from the original system, the analytical formulation, and the corresponding linear model, the reliability of the adopted approach can be clearly assessed. The following figures present detailed insights into the system’s behavior and highlight the agreement between different solution methods, as well as the conditions under which stable motion is achieved.

In Figure 2 the responses of $x (t)$ and $y (t)$ obtained from the original system equations (1) and (2), the analytical approach equations (24), and (25), and the equivalent exact linear ODEs as given in equations (34) and (35), are compared in this figure. At the initial stage, the system shows oscillations with relatively large amplitudes due to the ICs.

Figure 2.

(a), and (b): Shows time histories of the responses $x (t)$ and $y (t)$ obtained from the original system (red curve), the analytical solution (blue curve), and the exact linear differential solution (green curve).

As time progresses, the oscillations decrease gradually, indicating stable behavior governed by damping. The analytical results (blue curves) follow the original system response (red curves) closely over the entire time interval, showing no noticeable difference in either amplitude or frequency. The exact linear solution (green curves) shows a similar decay pattern, the three solutions become almost identical. This confirms that the analytical and linear models provide an accurate representation of the time-history behavior of the original system. For more convenience, consider the following data:

α_{1} = 0.1, α_{2} = 2.1, α_{3} = 1.1, α_{4} = 0.01, α_{5} = 0.1, α_{6} = 0.1, β_{1} = 1, β_{2} = 0.5, β_{3} = 0.5, β_{4} = 0.3, β_{5} = 0.1, β_{6} = 0.2, ω_{1} = 1.87, ω_{2} = 1.41, μ_{1} = 0.2, μ_{2} = 0.3, ω = 5, A = 0.1, B = 0.2 .

Figure 3(a) and (b) show the absolute errors between the numerical and the approximate analytical solution for $x$ and $y$ correspondingly. The inaccuracy for the displacement $x$ is quite modest in the whole examined time range with a maximum value about 2.92348× $10^{- 4}$ . The error presents an oscillating behavior with a slow decline with time until it becomes zero. This shows a very good agreement between the numerical and the approximate solution and validates the capacity of the approximation method to accurately represent the system dynamics. The error for displacement $y$ shows a similar fading oscillatory behavior, but with a considerably larger magnitude compared to $x$ . In the transient phase the error initially grows up to a maximum value of about $3.64461 \times 1 0^{- 3}$ and then is decreasing gradually, tending towards zero as the duration increases. This behavior indicates that the approximate solution can be used to characterise the transient and long-term response of the system. In conclusion, the small scale and decreasing trend of the errors in both $x$ and $y$ , validate the approximate solution method and show great accuracy in anticipating the nonlinear response of the coupled system.

Figure 3.

(a), and (b): Absolute error between the numerical and approximate analytical solutions for $x$ and $y$ versus time $t$ .

To gain further insight into the system behavior, the influence of selected parameters on the dynamic response is examined. Figure 4 and Figure 5 show how the system response changes when the parameter $μ_{1}$ and $μ_{2}$ are varied. Different values of $μ_{1}$ and $μ_{2}$ leads to different response patterns.

Figure 4.

Shows time histories and phase portraits for different values of $μ_{1}$ : (a and b) $μ_{1} = 0.5$ , (c and d) $μ_{1} = 0.1$ , and (e and f) $μ_{1} = 0.03$ .

Figure 5.

Displays time histories and phase portraits for different values of $μ_{2}$ : (a and b) $μ_{2} = 0.5$ , (c and d) $μ_{2} = 0.1$ , and (e and f) $μ_{2} = 0.03$ .

From time histories, the oscillation amplitude decreases with time in all cases. When $μ_{1}$ or $μ_{2}$ are small, the oscillations last longer. For larger values of $μ_{1}$ or $μ_{2}$ , the motion disappears more quickly and the response settles earlier. This effect also appears in the phase-plane plots. Moreover, for higher damping, the trajectories shrink faster toward the center. At lower damping levels, the spirals remain wider and take more time to decay. These plots clearly show the influences of $μ_{1}$ and $μ_{2}$ on the system motion.

Figures 6(a), (c) and 6(e) show the time domain responses of the displacement variables $x (t)$ and $u (t)$ for different values of equivalent natural frequency $ω_{e q 1} = 1.0, 2.0, and 4.0$ , correspondingly. As $ω_{e q 1}$ increases, the oscillation density becomes much higher within the same period, indicating that the system experiences more cycles of vibration per unit time. Even though the oscillation frequency increases, the response amplitudes clearly exhibit an exponential decay due to the damping effect, where the vibration amplitudes gradually decrease until they completely vanish with time. This confirms that the dissipative forces continuously remove energy from the system and drive the motion toward equilibrium. The accompanying phase portraits in Figures 6(b), (d), and (f) depict the relation between displacement $u (t)$ and velocity $u^{'} (t)$ .The trajectories form inward spiral patterns converging toward the equilibrium point $(0, 0)$ , which is a characteristic feature of an asymptotically stable damped system. Furthermore, increasing $ω_{e q 1}$ makes the spirals denser and more compact, which agrees with the higher rotational speed of the system trajectory in the phase plane. However, all trajectories eventually converge to the origin with time, demonstrating the strong stabilizing effect of damping and confirming the asymptotic stability of the system.

Figure 6.

Time-domain responses and corresponding phase portraits of the variables $x (t)$ and $u (t)$ for different values of the equivalent natural frequency $ω_{e q 1}$ : (a and b) $ω_{e q 1} = 1.0$ , (c and d) $ω_{e q 1} = 2.0$ , and (e and f) $ω_{e q 1} = 4.0$ .

Figures 7(a), (c), and (e) show the time histories of $y (t)$ and $v (t)$ for different values of the equivalent natural frequency $ω_{e q 2} = 1.0, 2.0, and 4.0$ , respectively. The present responses indicate a faster decay rate compared to the first system, where the oscillations practically vanish at approximately $t \approx 25$ . At lower frequency values, the oscillations appear relatively well separated and smooth, whereas increasing $ω_{e q 2}$ leads to tightly packed high-frequency oscillations accompanied by rapid decay. Such behavior suggests that higher equivalent frequencies enhance the dynamic interaction and increase the rate of energy dissipation within the system. The corresponding phase-plane plots shown in Figures 7(b), (d), and (f) illustrate inward spiral trajectories that start from relatively larger initial amplitudes and gradually converge toward the equilibrium point at the origin. It is also observed that the velocity magnitude $v^{'} (t)$ increases significantly with increasing $ω_{e q 2}$ . For $ω_{e q 2} = 4.0$ , the velocity range becomes considerably larger compared to the lower-frequency cases, indicating that higher frequencies generate faster instantaneous motion of the vibrating body. Despite the increase in both velocity magnitude and oscillation density, all trajectories eventually collapse into the equilibrium point, confirming the global asymptotic stability of the considered nonlinear system.

Figure 7.

Time-domain responses and corresponding phase portraits of the variables $y (t)$ and $v (t)$ for different values of the equivalent natural frequency $ω_{e q 2}$ : (a and b) $ω_{e q 2} = 1.0$ , (c and d) $ω_{e q 2} = 2.0$ , and (e and f) $ω_{e q 2} = 4.0$ .

Figure 8 illustrates the stability boundaries in $(ω_{m}, Ω_{1}^{2})$ and $(ω_{m}, Ω_{2}^{2})$ planes under systematic variations of damping and nonlinear parameters. The shaded regions denote stable steady-state responses, whereas the lower un-shaded areas correspond to instability zone. In Figure 8(a) and (b), the analysis is carried out for the parameter set:

α_{1} = α_{2} = α_{3} = 0.1, α_{4} = 4.01, α_{5} = 7.1, α_{6} = 0.1, β_{1} = 0.1, β_{2} = 0.5, β_{3} = 0.5, β_{4} = 0.3, β_{5} = 0.1, β_{6} = 0.2, ω_{1} = 1.87, ω_{2} = 2, A = 0.1, B = 0.2 .

Within this configuration, the increasing of the primary damping

μ_{1}

and the secondary damping

μ_{2}

from lower to higher values results in a steady upward shift of the stability boundary throughout the considered excitation-frequency range. Although the overall contour of the transition curve remains similar, the critical amplitude associated with instability onset increases as the damping coefficients grow. This behavior indicates that enhanced dissipation suppresses the amplification mechanisms induced by parametric excitation and limits the expansion of instability regions near critical frequency ratios. A distinct effect appears in Figure 8(c), where the nonlinear stiffness parameter is varied as

α_{1} = - 15, - 5, 5, and 15

; meanwhile the remaining parameters are fixed at

α_{2} = α_{3} = 0.1, α_{4} = 4.01, α_{5} = 7.1, α_{6} = 0.1, β_{1} = 0.1, β_{2} = 0.5, β_{3} = 0.5, β_{4} = 0.3, β_{5} = 4.1, β_{6} = 6.2, ω_{1} = 1.87, ω_{2} = 2, μ_{1} = 0.1, μ_{2} = 0.3, A = 0.1, B = 0.2 .

In this case, the variation of

α_{1}

does not merely produce a vertical displacement of the boundary; instead, it modifies its curvature and alters the position of the maximum stability peak. Positive values of

α_{1}

(hardening behavior) slightly elevate the peak and shift it toward higher frequencies, whereas negative values (softening behavior) reduce the peak and move it toward lower frequencies. These observations confirm that cubic nonlinear stiffness influences the effective frequency–amplitude relationship through amplitude-dependent stiffness modulation, which is characteristic of Duffing-type systems.

Figure 8.

Shows stability curves versus $ω$ for different parameter values: (a) $μ_{1}$ , (b) $μ_{2}$ , (c) $α_{1}$ , (d) $β_{1}$ , (e) $α_{4}$ ,and (f) $β_{4}$ .

Figure 8(d) examines the influence of the nonlinear coupling parameter $β_{1} = - 15, - 5, 5, and 15$ . An increase in $β_{1}$ leads to a clear expansion of the stable domain, accompanied by an upward movement of the transition boundary and a widening of the frequency interval that supports bounded oscillations. This trend suggests that stronger nonlinear modal interaction promotes energy exchange between coupled modes and mitigates instability growth. The impact of the higher-order stiffness parameter is presented in Figure 8(e), where $α_{4} = - 10, - 1, 5, and 10$ while.

α_{1} = α_{2} = α_{3} = 0.1, α_{5} = 7.1, α_{6} = 0.1, β_{1} = 0.1, β_{2} = 0.5, β_{3} = 0.5, β_{4} = 0.3, β_{5} = 4.1, β_{6} = 6.2, ω_{1} = 1.87, ω_{2} = 2, μ_{1} = 0.3, μ_{2} = 1.5, A = 0.1, B = 0.2 .

Variations in $α_{4}$ primarily produce a moderate elevation of the maximum stability level without significantly altering the global shape of the boundary, indicating that higher-order nonlinear stiffness contributes additional resistance at larger oscillation amplitudes.

Finally, Figure 8(f) considers the nonlinear control parameter $β_{4} = - 10, - 1, 5, an d 10$ . Among all investigated coefficients, $β_{4}$ exhibits the strongest stabilizing effect. Increasing its magnitude generates noticeable separation between the stability curves and significantly reduces the extent of the unstable region within the operational frequency band. The consistent elevation of the transition boundary demonstrates that $β_{4}$ plays a dominant role in shaping the overall stability topology of the nonlinear excited system. Overall, the results indicate that the stability characteristics arise from the combined action of damping-induced dissipation and nonlinear stiffness–coupling interactions, where damping primarily suppresses parametric amplification, while nonlinear coefficients reshape resonance behavior and determine the extent of the stable operating domain.

4. Chaotic response

In this Section, we show that a dynamical model’s chaotic behavior^36,37 must be examined to comprehend how sensitive the system is to changes in the parameters and ICs. The BD shows how variations in a control parameter led to qualitative changes in the dynamics of the system. To do this, we build the BD for the nonlinear equations (1) and (2) to distinguish between the various motion types and to show the parameter ranges where the system behaves in a periodic or chaotic manner.

Solving ODEs alone in studying nonlinear dynamical systems is insufficient; it is essential to study the nature of the system and how it responds to the changes of the parameters. Here, we employ a combination of some analytical tools as BDs, which illustrate changes in system behavior, and LEs, which act as a measuring of the disorder of the system, in which a positive value indicates chaotic behavior. Additionally, we use PMs and phase portraits to realize the nature of stability and disorder. Furthermore, we investigate the effect of different frequencies $ω_{1}$ and $ω_{2}$ : meanwhile monitoring changes in the system’s fundamental variables $x$ and $y$ . The BDs and the corresponding LEs versus $α_{5}$ are showed in Figure 9 and Figure 10 for the variables $x$ and $y$ respectively. We clearly observe a high degree of agreement in the dynamic paths. It shows that the system starts with stable periodic behavior at low values of the excitation coefficient. As the value increases, the system undergoes periods of period-doubling until it reaches a peak of disorder and chaos, as confirmed by the accompanying Lyabonov curves below each diagram. Analysis of phase paths in Figure 11 and Figure 12, when the value of $α_{5} = 0.05 .$ Figure 11 show that the system exhibits high discipline, with the PM showing limited points indicating a stable periodic orbit. However, when moving to $α_{5} = 4.4$ as in Figure 12, we observe an explosion of paths and complex interference, where the Poincaré points transform into a chaotic cloud, indicating that the system has entered the “strange attractor” region. The second case study (at $ω_{1} = 2$ ), the Figure 13, Figure 14, Figure 15, and Figure 16, further illustrate this view. In the BDs as given in Figure 13 and Figure 14, we observe how the system reconfigures itself in response to frequency. The phase portrait and PMs in Figure 15 and Figure 16 show profound path shifts, where the complexity of the variables increases proportionally with the excitation strength, reflecting the system’s extreme sensitivity to ICs. The Third case study the effect of increasing frequency (at $ω_{2} = 3$ ), in Figure 17, Figure 18, Figure 19, and Figure 20. It must note that the system moves to a more dramatic stage. When the frequency increased to $ω_{2} = 3$ . In Figure 17 and Figure 18, we observe that the regions of stability become narrower, and the system tends to slide into disorder more rapidly than in previous cases. On the other hand, the phase portrait and PMs, as shown in Figure 19 and Figure 20, provide that the paths become denser and more widespread, suggesting that the higher frequency has expanded the range of disorders and made it difficult to predict its future behavior without rigorous statistical analysis. Finally, the system does not have a fixed behavior, but rather a “dynamic change” that follows each increase in arousal.

Figure 9.

Shows the BDs of $x$ and the corresponding LE versus $α_{5}$ .

Figure 10.

Shows the BDs of $y$ and the corresponding LE versus $α_{5}$ .

Figure 11.

Shows the phase portrait and Pioncaré maps for the variables $x$ and $y$ for at $α_{5} = 0.05$ .

Figure 12.

Shows the phase portrait and Pioncaré maps for the variables $x$ and $y$ for at $α_{5} = 4.4$ .

Figure 13.

The bifurcation diagrams of $x$ and the regards the LE versus $α_{5}$ at $ω_{1} = 2$ .

Figure 14.

The bifurcation diagrams of $y$ and the regards the LE versus $α_{5}$ at $ω_{1} = 2$ .

Figure 15.

The phase portrait and Pioncaré maps for the variables $x$ and $y$ for at $α_{5} = 0.05$ .

Figure 16.

The phase portrait and Pioncaré maps for the variables $x$ and $y$ for at $α_{5} = 4.4$ .

Figure 17.

The bifurcation diagrams of $x$ and the regards the LE versus $α_{5}$ at $ω_{2} = 3$ .

Figure 18.

The bifurcation diagrams of $y$ and the regards the LE versus $α_{5}$ at $ω_{2} = 3$ .

Figure 19.

The phase portrait and Pioncaré maps for the variables $x$ and $y$ for at $α_{5} = 0.05$ .

Figure 20.

Shows the phase portrait and Pioncaré maps for the variables $x$ and $y$ for at $α_{5} = 4.4$ .

5. Conclusions

The dynamics of a coupled nonlinear parametric oscillatory system were examined in the current work using an NPA to evaluate the efficacy of the suggested methodology. By converting the original nonlinear ordinary differential equations into comparable linear forms, the main idea of the chosen method avoids the drawbacks frequently connected with conventional perturbation approaches. The current formulation is especially useful for systems displaying strong nonlinearities and large-amplitude oscillations because it does not rely on the existence of a small perturbation parameter, unlike traditional asymptotic procedures, except for the use of appropriate series expansions.

Additionally, the established framework makes it possible to evaluate the frequency-amplitude relationship efficiently, which is essential for building successive approximations of the system response under parametric excitation. Additionally, with improved mathematical simplicity, numerical efficiency, and practical application, the suggested methodology offers a fresh analytical viewpoint for researching coupled nonlinear oscillators.

Stability maps in many parameter planes were used to examine the system’s stability properties under different parameter settings. The obtained results showed that raising the damping coefficients significantly decreases the instability caused by parametric excitation by expanding the stable operating zones and shifting the stability boundaries upward. Furthermore, because the effective stiffness is amplitude-dependent, it was discovered that the nonlinear stiffness coefficients, especially the cubic Duffing-type nonlinearity, considerably altered the location and form of stability areas. Additionally, through nonlinear energy exchange between the interacting oscillatory modes, the nonlinear coupling parameters helped to expand the stable region.

Phase portraits, PMs, BDs, and LLEs were used to further investigate the system’s global dynamics. As the excitation parameters increase, the system response shifts from regular periodic oscillations to period-doubling bifurcations and ultimately into chaotic motion, according to these nonlinear dynamical tools. The long-term stability and the beginning of chaos were better understood because of the LLE analysis’s successful differentiation between periodic and chaotic oscillatory regimes. All things considered, the current findings offer insightful knowledge of the stability and nonlinear dynamical behavior of parametrically excited coupled oscillatory systems and could be a helpful foundation for examining comparable nonlinear models found in contemporary scientific and engineering applications.

Footnotes

Acknowledgments

The researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University for its financial support (QU-APC-2026).

ORCID iD

T. S. Amer

Authors contributions

Ahmad Almutlg: Investigation, Conceptualization, Validation, and Reviewing.

Galal M. Moatimid: Investigation, Methodology, Data curation, Conceptualization, Validation, Reviewing.

T. S. Amer: Resources, Methodology, Conceptualization, Validation, Formal analysis, Visualization, Reviewing and Editing.

Y. Y. Ellabban: Methodology, Resources, Data curation, Validation, Writing- Original draft preparation.

Funding

The researchers would like to thank the Deanship of Graduate Studies and Scientific Research at Qassim University for its financial support (QU-APC-2026).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data Availability Statement

The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.*

Appendix

The coefficients appearing in equations (30) and (31) may be formulated as follows:

c_{1} = \frac{24 α_{4} A^{2} μ_{1} + 64}{9 α_{4}^{2} A^{4}},

c_{2} = - \frac{64}{9 α_{4}^{2} A^{4}} (ω_{1}^{2} - \frac{1}{4} μ_{1}^{2} + \frac{3}{4} A^{2} α_{1} + \frac{5}{8} A^{4} α_{2} + \frac{35}{64} A^{6} α_{3} + \frac{p \sin (π / p)}{(p^{2} - 1) (p \sin (π / p) + π)} (α_{5} + \frac{B}{A} α_{6})),

c_{3} = \frac{24 β_{4} B^{2} μ_{2} + 64}{9 β_{4}^{2} B^{4}},

and

c_{4} = - \frac{64}{9 β_{4}^{2} B^{4}} (ω_{2}^{2} - \frac{1}{4} μ_{2}^{2} + \frac{3}{4} B^{2} β_{1} + \frac{5}{8} B^{4} β_{2} + \frac{35}{64} B^{6} β_{3} + \frac{p \sin (π / p)}{(p^{2} - 1) (p \sin (π / p) + π)} (β_{5} + \frac{A}{B} β_{6})) .

References

Aumentado

. Superconducting Parametric Amplifiers: The State of the Art in Josephson Parametric Amplifiers. IEEE Microwave Magazine 2020; 21(8): 45–59. https://doi.org/10.1109/mmm.2020.2993476

Rhoads

Shaw

Turner

. Nonlinear dynamics and its applications in micro-and nanoresonators. Journal of Dynamic Systems, Measurement, and Control 2010; 132(3): 1–14. https://doi.org/10.1115/1.4001333

Kovacic

Rand

Sah

. Mathieu’s equation and its generalizations: Overview of stability charts and their features. Applied Mechanics Reviews 2018; 70(2): 020802. https://doi.org/10.1115/1.4039144

McLachlan

. Theory and Applications of Mathieu Functions. Clarendon Press, 1947.

Achala

. Mathematical analysis and applications of Mathieu’s equation revisited. International Journal of Mathematics and its Applications 2021; 9(2): 49–54.

Ramani

Keith

Rand

. Perturbation solution for secondary bifurcation in the quadratically-damped Mathieu equation. International Journal of Non-Linear Mechanics 2004; 39(3): 491–502. https://doi.org/10.1016/s0020-7462(02)00218-4

Nayfeh

Zavodney

. The response of two-degree-of-freedom systems with quadratic non-linearities to a combination parametric resonance. Journal of Sound and Vibration 1986; 107(2): 329–350. https://doi.org/10.1016/0022-460x(86)90242-7

Xiu

J-H

El-Dib

, et al. Fractal Duffing oscillators with two degrees of freedom and cubic–quintic nonlinear stiffness. Fractal and Fractional 2026; 10(4): 265. https://doi.org/10.3390/fractalfract10040265

Guo

Liu

. Bifurcation and chaos analysis of coupled Duffing oscillators with two-periodic excitations and distributed time delays. International Journal of Bifurcation and Chaos 2026; 36(03): 2650029. https://doi.org/10.1142/s021812742650029x

10.

Xie

, et al. Resonance and attraction domain analysis of asymmetric Duffing systems with fractional damping in two degrees of freedom. Chaos, Solitons & Fractals 2024; 187: 115440. https://doi.org/10.1016/j.chaos.2024.115440

11.

Wang

Long

Yue

, et al. Bifurcation analysis of stick–slip vibration in a 2-DOF nonlinear dynamical system with dry friction. Communications in Nonlinear Science and Numerical Simulation 2022; 111: 106475. https://doi.org/10.1016/j.cnsns.2022.106475

12.

W-L

Wang

. Study on the dynamics of two-degree-of-freedom fractional differential piecewise nonlinear systems under harmonic excitation. International Journal of Non-Linear Mechanics 2022; 149(1): 104302. https://doi.org/10.1016/j.ijnonlinmec.2022.104302

13.

Wang

Liu

. Dynamics analysis of an impulsive stochastic model for spruce budworm growth. Applied Mathematics and Computations 2021; 19: 336–359.

14.

. Homotopy perturbation technique. Computer Methods in Applied Mechanics and Engineering 1999; 178: 257–262. https://doi.org/10.1016/s0045-7825(99)00018-3

15.

Moatimid

Amer

. Analytical approximate solutions of a magnetic spherical pendulum: Stability analysis. Journal of Vibration Engineering & Technologies 2023; 11: 2155–2165. https://doi.org/10.1007/s42417-022-00693-8

16.

Moatimid

Amer

. Analytical solution for the motion of a pendulum with rolling wheel: Stability Analysis. Scientific Reports 2022; 12(1): 12628. https://doi.org/10.1038/s41598-022-15121-w

17.

Zhang

J-G

Song

Q-R

Zhang

J-Q

, et al. Application of He’s frequency formula to nonlinear oscillators with generalized initial conditions. Facta Universitatis Series: Mechanical Engineering 2023; 21(4): 701–712. https://doi.org/10.22190/fume230909047z

18.

. A short remark on He’s frequency formulation. Journal of Low Frequency Noise, Vibration and Active Control 2022; 41(4): 1380–1385. https://doi.org/10.1177/14613484221107291

19.

. Simplified Hamiltonian-based frequency-amplitude formulation for nonlinear vibration systems. Facta Universitatis Series Mechanical Engineering 2022; 20(2): 445–455. https://doi.org/10.22190/fume220420023m

20.

J-H

. The simplest approach to nonlinear oscillators. Results in Physics 2019; 15: 102546. https://doi.org/10.1016/j.rinp.2019.102546

21.

C-H

Liu

. A modified frequency-amplitude formulation for fractal vibration systems. Fractals 2022; 30(03): 2250046. https://doi.org/10.1142/s0218348x22500463

22.

Niu

J-Y

Feng

G-Q

Gepreel

. A simple frequency formulation for fractal-fractional non-linear oscillators: A promising tool and its future challenge. Frontiers in Physics 2023; 11: 1158121. https://doi.org/10.3389/fphy.2023.1158121

23.

Barakat

Weig

Hagedorn

. Non-trivial solutions and their stability in a two-degree-of-freedom Mathieu–Duffing system. Nonlinear Dynamics 2023; 111: 22119–22136. https://doi.org/10.1007/s11071-023-08659-5

24.

El-Dib

. Insightful and Comprehensive Formularization of Frequency-Amplitude Formula for Strong or Singular Nonlinear Oscillators. Journal of Low Frequency Noise, Vibration and Active Control 2022; 42(1): 89–109. https://doi.org/10.1177/14613484221118177

25.

El-Dib

. The Frequency Estimation for Non-Conservative Nonlinear Oscillation. ZAMM: Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 2021; 101(12): e202100187. https://doi.org/10.1002/zamm.202100187

26.

El-Dib

. The damping Helmholtz-Rayleigh-Duffing oscillator with the non-perturbative approach. Mathematics and Computers in Simulation 2022; 194: 552–562. https://doi.org/10.1016/j.matcom.2021.12.014

27.

El-Dib

. The periodic property of Gaylord’s oscillator with a non-perturbative method. Archive of Applied Mechanics 2022; 92(11): 3067–3075. https://doi.org/10.1007/s00419-022-02269-0

28.

El-Dib

. Properties of complex damping Helmholtz-Duffing oscillator arising in mechanics, Journal of Low Frequency Noise. Vibration and Active Control 2023; 42(2): 589–597. https://doi.org/10.1177/14613484221138560

29.

El-Dib

. The simplest approach to solving the cubic nonlinear Jerk oscillator with the non-perturbative method. Mathematical Methods in the Applied Sciences 2022; 45: 1–19. https://doi.org/10.1002/mma.8099

30.

El-Dib

Al-Ghamdi

. Renormalization method for variable stiffness in fractal 2DOF nonlinear parametric oscillators. Journal of Vibration Engineering & Technologies 2025; 13: 550. https://doi.org/10.1007/s42417-025-02108-w

31.

Liu

Wang

. Variational principle, Hamiltonian, sensitivity analysis and new soliton solutions for the nonlinear fractional Schrödinger–Hirota equation. Fractals 2026; 34(03): 2650002. https://doi.org/10.1142/s0218348x26500027

32.

Wang

. Dynamical behaviors and diverse soliton solutions for the fractional Yajima–Oikawa equation. Fractals 2026; 34(05): 2650020. https://doi.org/10.1142/s0218348x26500209

33.

Wang

. Bifurcation, chaotic, variational principle and new soliton solutions for the nonlinear fractional Drinfeld–Sokolov system. Fractals 2026; 2650062. https://doi.org/10.1142/S0218348X26500623

34.

Wang

. Investigation of soliton structures and dynamical behaviors for fractional Kaup–Newell system. Fractals 2026; 2650081. https://doi.org/10.1142/S0218348X26500817

35.

Saeed

EllabbanY

Moatimid

, et al. Nonlinear interactions of an n-layer X-shape low-frequency vibration isolator equipped with a nonlinear vibration absorber at 1:1 internal resonance: analytical and numerical investigations. Physica Scripta 2024; 99: 105207. https://doi.org/10.1088/1402-4896/ad6f53

36.

Amer

Moatimid

Zakria

, et al. Dynamical analysis of a four-degree-of-freedom vibratory structure: Bifurcation, stability, and resonance exploration, Journal of Low Frequency Noise. Journal of Low Frequency Noise, Vibration and Active Control 2025; 44(3): 1555–1575. https://doi.org/10.1177/14613484251328786

37.

Alanazy

Elbaz

Amer

, et al. Analyzing stability and bifurcation of a controlled dynamical system: A comparative study of negative velocity feedback and positive position feedback controllers. European Journal of Mechanics /A Solids 2026; 118: 106053. https://doi.org/10.1016/j.euromechsol.2026.106053

Analytical modelling of a nonlinear 2DOF Mathieu–Duffing oscillator

Abstract

Keywords

1. Introduction

2. Problem construction

2.1. Governing ODEs

3. Results and discussion

4. Chaotic response

5. Conclusions

Footnotes

Acknowledgments

ORCID iD

Authors contributions

Funding

Declaration of conflicting interests

Data Availability Statement

Appendix

References