Novel methodology of a moving cart on a horizontal plane coupled with an inclined damped nonlinear spring

Abstract

This work introduces a comprehensive framework in analysing nonlinear dynamics of a moving cart with an inclined smooth surface on a horizontal plane featuring a damped nonlinear spring. It integrates nonlinear elastic behaviour and damping effects to depict mechanical engineering systems. The study aims to enhance the comprehension of intricate dynamic reactions affected by gravity, nonlinear stiffness, and energy dissipation, which are relevant to vibration isolation and mechanical transport. The issue involves a two-degrees-of-freedom of dynamical model with parametric excitation. The methodology employs He’s frequency formula that implemented via a non-perturbative approach. This method aims to convert ordinary differential equations that are nonlinear into linear ones. To confirm the significant accuracy between nonlinear and linear/nonlinear ODEs, Mathematica Software was adopted. The distinctive methodology is defined by its simplicity, attractiveness, potential, and efficiency, rendering it beneficial for highly nonlinear ordinary differential equations. It can be employed in different classes’ fields of dynamical systems as well as fluid mechanics. The stability analysis illustrates several non-dimensional physical concerns. The system’s erratic behaviour is analysed to understand its dynamics and transition from a periodic to a chaotic state. This organization utilizes bifurcation diagrams and the largest Lyapunov exponent which communally provide a comprehensive depiction of the system’s behaviour across several stages. Moreover, the system’s complexity landscape is visualized through 2D chaos maps based on Kaplan–Yorke dimension, pinpointing the exact parameter shifts where simple order collapses into high-dimensional fractal chaos.

1. Introduction

In engineering systems, joints and connectors serve a variety of purposes, including increasing total cost efficiency, enabling assembly of large arrangements from simple, easily manufactured, and adding flexibility to precise DOF. The use of connections makes it easier to build large buildings in their assigned operational locations, while allowing of disassembly when necessary of effective packing and practical transportation. As seen by interference fits used between blade and disk in bladed disk assemblies, they may contain other components, such as bolts or rivets common in many technical applications, or they may exist without any parts at all. Furthermore, they are utilized in buildings, ¹ bladed disk assemblies,^2,3 aircraft, ⁴ and automobile applications.^5,6 Joints cause localized changes in the system, such as increased damping and changes in the stiffness characteristics of the structure, which can have a significant impact on the dynamics of the interconnected structures. The development of a reliable mathematical model of jointed systems is hampered by the noticeable nonlinear behaviour of contacting surfaces caused by friction, particularly under large-amplitude vibrations. While there are analytical techniques, such as lumped models for bolted flanges, ⁷ it is very difficult, if not impossible, to capture joint dynamics accurately analytically. Accordingly, experimental parametric identification has become a conventional method, either via static testing, e.g., tensile tests, ⁸ or dynamic measures derived from modal and vibrational responses.^9,10 Hammer strikes Ref. ¹¹ or shakers Ref. ¹² can produce dynamic excitation, and accelerometers, ¹³ laser vibrometers, ¹⁴ and microphones ¹⁵ can be used to record reactions. A new method for determining the nonlinear properties of a connection used to couple two structures was proposed. ¹⁶ Using the measured nonlinear responses of system in the frequency domain, the proposed method in identifying nonlinear internal forces calculates the nonlinear internal forces and their corresponding relative displacements at the nonlinear connections. The technique used in the present study is completely different from that used in other studies. As previously mentioned, this study investigates the nonlinear 2DOF of a unique method known as NPA. The NPA diverges from previous models by incorporating and thus addressing certain deficiencies discovered in past studies. This differentiation is not merely incremental but represents a substantial methodological advancement. Furthermore, while Refs. 1–16 provide a valuable basis, their techniques do not cover the whole spectrum of situations and scenarios that NPA aims to address. We hope this explanation highlights the novelty and significance of the proposed methodology, and we are prepared to revise the work to make this distinction more apparent if necessary.

An analysis of more than 200 major problems with nonlinear oscillation computational and analytical techniques was conducted. ¹⁷ The oscillation of a pendulum on an extended, elastic thread was examined. ¹⁸ The relationship between linear/nonlinear stabilities in the propagation of explosive waves in complex sensitive systems was examined via linear stability analysis. ¹⁹ The magnetic forces and swirl currents affect the nonlinear vibrational characteristics of conductive beams was studied. ²⁰ New correlations between electromagnetic interaction forces and ferromagnetic materials were made possible via Maxwell’s equations and Lorentz effects. Analysis was performed on the transient response characteristics with initial geometric mistakes. ²¹ An articulated ring truss antenna’s nonlinear dynamic response under different 1:1 internal resonances during style transition was reported. ²² A strong magnet interacting with two symmetrically positioned loops close to the equilibrium angular position has an impact on the nonlinear dynamics of a magnetic pendulum. ²³ A system consisting of two coupled van der Pol oscillators impacted by an external force supplied to one oscillator was examined for nonlinear dynamics and attraction basin features. ²⁴ Recent single wave findings were obtained via the modified formulation. ²⁵ The complex wave patterns were examined via an efficient auxiliary formulation mapping technique. It was shown that incompressible fluids have less intermittency and dimensional dissipation. ²⁶

Extra extensive difficulties in computational and analytical techniques of nonlinear oscillations were examined. ²⁷ The correlation between linear/nonlinear stabilities in the propagation of explosive waves within intricately sensitive systems was examined via linear stability analysis. ²⁸ The influences of eddy currents and magnetic forces on the nonlinear vibrational characteristics of conductive beams were examined. ²⁹ Innovative correlations between electromagnetic interaction forces and ferromagnetic materials were made possible via Maxwell’s equations and Lorentz effects. The transient response characteristics with initial geometric defects were examined. ³⁰ The nonlinear dynamic behaviour of an articulated ring truss antenna under various 1:1 internal resonances during style swaps was documented. ³¹ A powerful magnet interacting with two symmetrically organized loops close to the equilibrium angular position affects the nonlinear dynamics of a magnetic pendulum. ³² A system consisting of two coupled van der Pol oscillators, affected by an external force supplied to one oscillator, was studied for its nonlinear dynamics and attraction basin behaviour. ³³ Recent solitary wave results were obtained via the updated Ginzburg–Landau formulation. ³⁴ The complex wave forms were examined via a potent auxiliary formulation mapping approach. There were reports of decreased dimensional dissipation and intermittency in incompressible fluids. ³⁵ A crucial component of nonlinear dynamics is the stability analysis of 2DOF oscillators with Duffing–Helmholtz-type nonlinearities that are controlled by coupled nonlinear ODEs with damping and periodic forcing. The study Ref. ³⁶ investigated the stability and nonlinear dynamics of a coupled parametric oscillator using a non-perturbation analytical approach. This method effectively captures intense nonlinear behaviors and detects transitions from periodic motion to chaos by analyzing stability and bifurcations.

The HFF was introduced by Professor Ji-Huan He to analyse nonlinear oscillations. He offered improvements to the initial prototype and summarized the amplitude‒frequency representations of nonlinear oscillators. The choice of an assessment point is crucial to both the original and improved versions, but there is no standardized process for making this choice. ³⁷ Investigators have utilized the weighted residual method to facilitate swift estimation of the amplitude-frequency curves of oscillators displaying discontinuities. ³⁸ The method yields precise findings with minimum processing, and a specific criterion for identifying assessment locations was proposed. ³⁹ An association of researchers has successfully implemented HFF. A study demonstrated a novel trial function, determined that it was a workable solution for a nonlinear oscillator at a given frequency, and confirmed its effectiveness on the DO. ⁴⁰ The HFF, which is based on an ancient Chinese mathematical technique, is still a reliable framework in dealing with nonlinear vibration problems. The method’s effectiveness in real-world applications was confirmed by an additional theoretical evaluation, which also recommended a changed formulation. ⁴¹ Additionally, researchers have found a precise link between the amplitude and period of the controlling ODE by comparing a nonlinear oscillator with its linear equivalent in the relevant domain Ref. ⁴² . An integer parameter that can be appropriately modified of each instance is included in the amplitude-period ODE; several illustrative examples compared this approximation with the conventional HFF. ⁴³ By choosing a suitable weighting function inside the objective function, a formula of ideal formative period of a conservative nonlinear oscillator was developed; in a number of common situations, the estimated period closely resembles the precise value. ⁴⁴ The precision of the improved methodology is assessed via highly nonlinear DO as a benchmark. ⁴⁵ Researchers in the domains of fluid mechanics and dynamical systems use NPA to overcome the drawbacks of conventional perturbation techniques.^46–48 This innovative method offers a profound comprehension of complex dynamical systems and transcends conventional limitations.

The dynamic study of a rolling cart on a horizontal plane coupled with a damped nonlinear spring on its inclined surface constitutes a significant experimental framework in comprehending complex mechanical systems characterized by nonlinearity, energy dissipation, and gravity-driven motion. Researchers can examine the interplay of nonlinear restoring forces, damping processes and slope-induced gravitational components on system stability, oscillatory behaviour, and transient responses through the integration of theoretical modelling and controlled testing. These studies have substantial practical applications in engineering disciplines, encompassing the design of vibration isolation systems, vehicle suspension components, and energy-dissipating structures. Furthermore, the empirical validation of nonlinear dynamic models enhances predictive simulations utilized in robotics, automated transport systems on inclined surfaces, and mechanical control systems where motion transpires under fluctuating loads and angles. Comprehending the interrelated impacts of damping and nonlinear elasticity aids in optimizing energy efficiency and improving safety in mechanical devices exposed to repetitive or impact-driven motion, rendering this system a significant reference of both fundamental research and practical engineering applications. Considering the aforementioned features, the present situation is examined. To clarify the presentation of the manuscript, its remainder is structured as follows: The construction of the issue § 2. The NPA of the coupled ODEs of the examined system is detailed in § 3. § 4 analyses the stability analysis, utilizing Arnold tongues, Floquet multipliers. § 5 is devoted to analyse BDs, LLE, and KYD to thoroughly elucidate its dynamics. Final observations are encapsulated in § 6.

2. Formation of issue

In this paper, we analyse the 2DOF dynamical system involving two interacting subsystems. First, the first subsystem includes the cart having a mass $m_o$ which slides horizontally on the plane. The cart is attached to the nonlinear suspension system with stiffness factors $k_2$ and $k_4$ , damping factor $b_1$ , and reference length $l_2$ . The inclined guide mounted on the cart makes an angle ${\theta }_o$ with the horizontal line. Second, the second subsystem involves the uniform disc having a mass $m$ which rolls along the inclined plane but still maintains contact with the cart. The mass $m$ is restricted by a damper with a coefficient $b_2$ and the non-linear spring with stiffness coefficients $k_1$ and $k_3$ , and reference length $l_1$ . The coefficients $k_1$ and $k_2$ are for the linear stiffness, $k_3$ and $k_4$ are for the nonlinear stiffness. Moreover, a periodic force $F\left(t\right)$ acts on the mass $m$ of the disc at its centre, as depicted in Figure 1.

Figure 1.

Displays the structure of the analyzed model.

The system under consideration is a coupled nonlinear dynamical system of 2DOF, where the generalized displacements of two interacting modes of motion in a mechanical structure are represented by the time-dependent functions $\mathrm{x}\left(\mathrm{t}\right)$ and $\mathrm{y}\left(\mathrm{t}\right)$ . The system’s inertial effects are denoted by $\ddot{x}$ and $\ddot{y}$ , whereas the damping that dissipates energy during motion is represented by $2{\mathrm{c}}_1\dot{x}$ and $2c_2\dot{y}$ . The cross-coupling damping is characterized by values of $d_1$ and $d_2$ , which show that the velocity of one coordinate affects the dynamics of the other. The natural frequencies of the two modes are denoted by ${\omega }_1$ and ${\omega }_2$ , whereas the linear restoring forces are represented by ${\omega }_1^{2}x$ and ${\omega }_2^{2}y$ . The Duffing-type nonlinear stiffness, which simulates the geometric or material nonlinearities that emerge for large oscillations, is described by the cubic nonlinear terms ${\mu }_1{x}^3$ and ${\mu }_2{y}^3$ . The amplitudes of a parametric periodic excitation are denoted by the terms $f_1$ and $f_2$ , which act through the factors $x\hspace{0.17em}{f}_1\hspace{0.17em}\cos \hspace{0.17em}\sigma t$ and $y\hspace{0.17em}{f}_2\hspace{0.17em}\cos \hspace{0.17em}\sigma t$ , where $\sigma$ is the excitation frequency. Physically, this means that the system’s effective stiffness fluctuates over time, which could result in complicated nonlinear behaviour and parametric resonance. In many engineering applications, such as flexible structures, vibration of coupled beams or plates, rotating machinery, and nonlinear vibration absorbers, interactions between various modes of motion and nonlinear stiffness effects are crucial to the system response, such as coupled nonlinear oscillator models, which frequently appear. Consequently, the system under study is depicted in Figure 1.

The coupled governing ODEs of motion for the present system in the dimensionless form may be formulated as follows Ref. ²⁹ :

$\ddot{x}+d_1\ddot{y}+2c_1\dot{x}+{\omega }_1^{2}x+{\mu }_1{x}^3=f_1\hspace{0.17em}\cos \hspace{0.17em}\sigma t,$

(1)

and

$\ddot{y}+d_2\ddot{x}+2c_2\dot{y}+{\omega }_2^{2}y+{\mu }_2{y}^3=f_2\hspace{0.17em}\cos \hspace{0.17em}\sigma t.$

(2)

Quadratic terms can be converted into other odd terms via multiplication, division, integration, or differentiation within the context of perturbation techniques. The categorization of odd functions resulting from damping and other odd terms associated with stiffness effects is the main basis for the application of NPA. Damping or stiffness forces are addressed by including quadratic elements. The equivalent damping is obtained from the first of the two main integrals, whereas the equivalent frequency is obtained from the second. Either the primary integral or the alternate integral can include the integral of the even terms. The thorough efforts of El-Dib,^49,50 who authored multiple articles on the subject, clarified and analysed the specific approach of the NPA. Only two items from El-Dib are chosen.^49,50 Additionally, for convenience, we replace $\sigma$ with $2\sigma$ , and then, the coupled system of nonlinear ODE given by Eqs. (3) and (4). It should be noted that the NPA primarily addresses secular terms, which pertain exclusively to the odd functions. Since the function $f_1\hspace{0.17em}\cos \hspace{0.17em}\sigma t$ can be expressed as a constant function with a time-dependent parameter $x$ . Similarly the function $f_2\hspace{0.17em}\cos \hspace{0.17em}\sigma t$ can be expressed as a constant function with a time-dependent parameter $y$ . This action may be referred to as the “masking technique”. In accordance with El-Dib, ⁴⁹ the utilization of the masking technique in nonlinear dynamics, particularly in scenarios involving periodic forces, constitutes a strategic and distinctive scientific methodology. This method primarily focuses on altering the governing ODE of the system to conceal the explicit influence of the periodic force. This reformulation optimizes the analytical process while efficiently maintaining the fundamental attributes of the system’s reaction to periodic input. The masking technique circumvents the direct complexities of periodic forces while still providing a comprehensive understanding of their impact on the system’s behaviour. This method is especially beneficial in clarifying the complex dynamics of nonlinear systems influenced by periodic forces, hence offering a more efficient analytical framework. The application of the masking technique commences with the initial formulation, as illustrated by Eqs. (1) and (2), in our work advance to effectively, remove the explicit expression of the periodic forces. This crucial transformation involves devising an alternative equation that obscures the influence of the periodic force on the system. The objective is to develop a comparable system that emulates the behaviour of the original system under periodic forcing, while omitting the force from the equation itself. Eq. (9) is derived by a sequence of mathematical transformations or approximations to accomplish this. These steps aim to subtly integrate the impact of the periodic force into the system’s dynamics. Conversely, the inclusion of $f_1\hspace{0.17em}\cos \hspace{0.17em}\sigma t$ and $f_2\hspace{0.17em}\cos \hspace{0.17em}\sigma t$ , as a non-homogeneous component in Eqs. (1) and (2), prevent the excitation force from the stability requirement. Furthermore, El-Dib Ref. ⁵⁰ presented a comprehensive study of the frequency-amplitude formulae for nonlinear oscillators and its developments.

$\ddot{x}+d_1\ddot{y}+2c_1\dot{x}+{\omega }_1^{2}x+{\mu }_1{x}^3-x{f}_1\hspace{0.17em}\cos \hspace{0.17em}2\sigma t=0,$

(3)

and

$\ddot{y}+d_2\ddot{x}+2c_2\dot{y}+{\omega }_2^{2}y+{\mu }_2{y}^3-y{f}_2\hspace{0.17em}\cos \hspace{0.17em}2\sigma t=0,$

(4)

with ICs: $x\left(0\right)=A,y\left(0\right)=B,\dot{x}\left(0\right)=\dot{y}\left(0\right)=0$ , the dimensionless coefficients are defined as:

$\begin{array}{l}{c}_1=\frac{b_1}{2{\omega }_o\left(m+m_o\right)},c_2=\frac{b_2}{2{\omega }_{o}m},d_1=\frac{m\hspace{0.17em}\cos \hspace{0.17em}{\theta }_o}{m+m_o},d_2=\cos \hspace{0.17em}{\theta }_o,{\omega }_1^2=\frac{k_2}{{\omega }_o^2\left(m+m_o\right)},\\ {\omega }_2^2=\frac{k_1}{m{\omega }_o^2},{\mu }_1=\frac{k_4{l}_1^2}{{\omega }_o^2\left(m+m_o\right)},{\mu }_2=\frac{k_3{l}_1^2}{m{\omega }_o^2},f_1=\frac{F}{{\omega }_o^2\left(m+m_o\right)},\text{and}\hspace{0.5em}{f}_2=\frac{F\hspace{0.17em}\cos \hspace{0.17em}{\theta }_o}{m{\omega }_o^2}.\end{array}$

(5)

3. NPA methodology

It is possible to simplify and transform the non-autonomous system into an analogous autonomous form, which facilitates the analysis of its dynamic behaviour and stability. The following formula can be used to approximate these natural frequencies Ref. ⁵¹ :

${\omega }_{1i}={\int }_0^{\pi /2\sigma }\hspace{0.17em}{\cos }^2\hspace{0.17em}\tau \left(f_i\hspace{0.17em}\cos \hspace{0.17em}2\sigma \tau \right)d\tau /{\int }_0^{\pi /2\sigma }\hspace{0.17em}{\cos }^2\hspace{0.17em}\tau d\tau =\frac{f_i\sigma \hspace{0.17em}\sin \left(\pi /\sigma \right)}{\left({\sigma }^2-1\right)\left(\pi +\sigma \hspace{0.17em}\sin \left(\pi /\sigma \right)\right)},i=1,2.$

(6)

By removing the challenges associated with variable coefficients, the coupled system shown in ODEs as given in Eqs. (3) and (4) significantly simplifies the understanding of the dynamic behaviour. Important phenomena, including resonance, balance, and matched nonlinear interactions, can be examined in greater detail by simplifying the underlying ODEs. Automating the process makes it easier to use sophisticated analytical approaches, ensuring precise dynamic behaviour characterization and enabling a deeper understanding of its performance and balance. The equations depict the coupled system as specified in Eqs. (3) and (4) turn into:

$\ddot{x}+d_1\ddot{y}+2c_1\dot{x}+\left({\omega }_1^2-{\omega }_{11}^2\right)x+{\mu }_1{x}^3=0,$

(7)

and

$\ddot{y}+d_2\ddot{x}+2c_2\dot{y}+\left({\omega }_2^2-{\omega }_{12}^2\right)y+{\mu }_2{y}^3=0.$

(8)

A comparison of the two nonlinear ODEs before and after the application of excitation forces is carried out for improved clarity to validate the preceding methods. For a chosen system with the following information:

$\begin{array}{l}A=0.7,B=0.8,c_1=0.05,c_2=0.3,d_{1\hspace{0.17em}}=0.05,d_2=0.08,\sigma =\pi /4,\\ f_1=0.02,f_2=0.015,{\omega }_1=0.8,{\omega }_2=0.9,{\mu }_1=0.13;\hspace{0.5em}\text{and}\hspace{0.5em}{\mu }_2=0.12.\end{array}$

A comparison of the two curves is shown in Figure 2, which reveals that adding the periodic term $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ causes slight changes in the oscillation amplitude and phase during the first time interval. However, the system’s general qualitative features are essentially unaltered. For the chosen parameter values, the red curve (with $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ ) closely resembles the blue curve (without $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ ), suggesting that the parametric forcing has little effect on the long-term dynamics. Figure 3 demonstrates that the evolution of v(t) is largely unaffected by the inclusion of the periodic forcing term cos(2σt), with both solutions remaining almost indistinguishable over the full time interval. Any difference are confined to the early transient phase, where only a very slight phase shift and negligible amplitude variation can be detected before rapidly fading. The reported maximum absolute error of approximately 0.02 further confirms the weak influence of the forcing on this state variable, indicating a higher level of robustness compared to the previously discussed dynamics. Overall, both responses converge to the same damped oscillatory decay, showing that the long-term behaviour is preserved and that the periodic perturbation introduces only a brief and minor transient effect.

Figure 2.

Illustrates an excellent comparison for $u\left(t\right)$ with respect to the nonlinear ODE (3) when $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ is present and its comparable ODE (7) when $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ is not. A maximum absolute error is approximately $0.03$ .

Figure 3.

Displays an excellent comparison for $v\left(t\right)$ with respect to the nonlinear ODE (4) when $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ is present and its comparable ODE (8) when $\cos \hspace{0.17em}2\sigma \hspace{0.17em}t$ is not. A maximum absolute error is approximately $0.02$ .

The excellent correlation between two nonlinear ODEs suggests that the solutions exhibit almost identical behavior throughout the relevant area, even though their mathematical representations may differ. Both ODEs yield nearly identical trajectories with nearly comparable growth, decay, oscillation patterns, or steady-state behavior when solved numerically or analytically under identical ICs. This suggests that one ODE may serve as an accurate approximation or simplified representation of the other while preserving the essential system behavior since both models adequately capture the basic nonlinear dynamics, such as feedback mechanisms, stability features, and responses to disturbances.

The strong correlation between two nonlinear ODEs indicates that the solutions they produce have almost identical behavior within the relevant domain, even though their mathematical forms may differ. Both ODEs yield trajectories that substantially match when solved numerically or analytically with the same beginning conditions, exhibiting almost identical patterns of growth, decay, oscillation, or steady-state behavior. This suggests that one ODE may function as an accurate approximation or simplified representation of the other while maintaining the basic system behavior since both models adequately capture the essential nonlinear dynamics, such as feedback mechanisms, stability characteristics, and responses to disturbances.

A schematic flowchart illustrating the sequential process of the suggested NPA framework is shown in Figure 4. The diagram shows the key steps of the methodology, beginning with the choice of a suitable trial function and continuing with the estimation of effective system characteristics, including frequency and damping, using integral-based formulations. The main dynamical performance of the original nonlinear ODE is then captured by using these parameters to create a consistent linearized representation. HFF is used to enhance frequency estimation, and numerical simulations using Mathematica’s NDSolve are used to validate the analytical approximation that results.^46–48

Figure 4.

Outlines the methodology of the solution, emphasizing the NPA processes.

There are no periodic coefficients in the nonlinear ODEs as agreed in Eqs. (7) and (8). Converting these two ODEs into two linear ones is the goal. Thus, the NPA of a connected system is expressed in accordance with Moatimid et al.^46–48 and El-Dib.^49–51 Creating similar linear ODEs that correlate to the nonlinear ODEs as approved in Eqs (7) and (8) is our main goal. The following is an articulation of the proposed trial solutions:

$u=A\hspace{0.17em}\cos \hspace{0.17em}{\mathrm{\Omega }}_{1}t\hspace{0.17em}and\hspace{0.17em}v=B\hspace{0.17em}\cos \hspace{0.17em}{\mathrm{\Omega }}_{2}t,$

(9)

with initial amplitudes $A$ and $B$ where the parameters ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_2$ address the total frequencies to be determined later.

The previous trial solutions satisfy the following ICs:

$u\left(0\right)=A,\dot{u}\left(0\right)=0,v\left(0\right)=B,\text{and}\hspace{1em}\dot{v}\left(0\right)=0.$

(10)

Notably, NPA has three primary limitations. The following is a possible list of these limitations:

1. Only weakly nonlinear oscillators are addressed by the method.

2. There are no changes made to the ICs.

3. The initial amplitudes need to be less than unity to improve accuracy.

Define $H\left(x\right)$ and $G\left(y\right)$ which represent the odd functions in ODEs as established in Eqs (7) and (8), respectively, as follows:

$H\left(x\right)=\left({\omega }_1^2-{\omega }_{11}^2\right)x+{\mu }_1{x}^3,$

(11)

and

$G\left(y\right)=\left({\omega }_2^2-{\omega }_{12}^2\right)y+{\mu }_2{y}^3,$

(12)

the equivalent frequencies are evaluated using MS as follows Ref. ^45–51:

${\omega }_{eqv1}^2={\int }_0^{2\pi /{\mathrm{\Omega }}_1}uH\left(u\right)dt/{\int }_0^{2\pi /{\mathrm{\Omega }}_1}{u}^{2}dt={\omega }_1^2-{\omega }_{11}^2+\frac{3}{4}{A}^2{\mu }_1,$

(13)

and

${\omega }_{eqv2}^2={\int }_0^{2\pi /{\mathrm{\Omega }}_2}vG\left(v\right)dt/{\int }_0^{2\pi /{\mathrm{\Omega }}_2}{v}^{2}dt={\omega }_2^2-{\omega }_{12}^2+\frac{3}{4}{B}^2{\mu }_2.$

(14)

Therefore, the comparable linear ODEs are written as follows

$\ddot{u}+d_1\ddot{v}+2c_1\dot{u}+{\omega }_{eqv1}^{2}u=0,$

(15)

and

$\ddot{v}+d_2\ddot{u}+2c_2\dot{v}+{\omega }_{eqv2}^{2}v=0.$

(16)

To decouple the system and facilitate resolution, we rearrange linear ODEs as decided in Eqs (15) and (16) as follows:

$\ddot{u}+2c_1\dot{u}+g_1\left(u,\ddot{v}\right)u=0,$

(17)

and

$\ddot{v}+2c_2\dot{v}+g_2\left(v,\ddot{u}\right)v=0,$

(18)

where the introduced nonlinear stiffening functions $g_1\left(u,\ddot{v}\right)$ and $g_2\left(v,\ddot{u}\right)$ are defined as follows:

$g_1\left(u,\ddot{v}\right)={\omega }_{eqv1}^2+\frac{d_1\ddot{v}}{u},$

(19)

and

$g_2\left(v,\ddot{u}\right)={\omega }_{eqv2}^2+\frac{d_2\ddot{u}}{v}.$

(20)

Using the trial (guessing) solutions that are provided in equation (9), the nonlinear stiffening functions $g_1\left(u,\ddot{v}\right)$ and $g_2\left(v,\ddot{u}\right)$ are given by

$g_1\left(u,\ddot{v}\right)={\omega }_{eqv1}^2-\frac{{\mathrm{\Omega }}_2^2{d}_{1}v}{u},$

(21)

and

$g_2\left(v,\ddot{u}\right)={\omega }_{eqv2}^2-\frac{{\mathrm{\Omega }}_1^2{d}_{2}u}{v}\hspace{0.17em}.$

(22)

The proposed square expressions support the calculation of coupled nonlinear interactions by providing an analytical framework in comprehending the system’s dynamic behaviour. Significant dynamic parameters are measured and complex nonlinear interactions are decoupled to produce mean square forms, which provide a scientific approach to addressing the problematic behaviour inherent in such systems. This method is extensively employed for analysing nonlinear vibration phenomena due to its ability to generate sufficiently accurate solutions and enhance the evaluation process, for more details, see El-Dib and Al-Ghamdi. ⁵¹ In other words, the obtained mean-square formulas establish a foundation for analysing nonlinear dynamical systems characterized by complex coupling. Such formulas establish an appropriate foundation for assessing various parameters of nonlinear dynamics and for separating nonlinear couplings. The proposed method constitutes a reliable and efficient approach for examining nonlinear vibrations.

${\overline{u}}^2=\frac{1}{4T_1^2}{\int }_0^{T_1}u{\left(t\right)}^{2}d\tau =\frac{A^2}{8},$

(23)

and

${\overline{v}}^2=\frac{1}{4T_2^2}{\int }_0^{T_2}v{\left(t\right)}^{2}d\tau =\frac{B^2}{8},$

(24)

where $T_i=\frac{\pi }{2{\mathrm{\Omega }}_i},i=1,2$ . In this approach, the shifted functions $g_1\left(u,\ddot{v}\right)$ and $g_2\left(v,\ddot{u}\right)$ are used to establish the new additional frequencies ${\varpi }_1^2$ and ${\varpi }_2^2$ defined as follows:

${\varpi }_1^2=\underset{u\to \overline{u},v\to \overline{v}}{\lim }{g}_1\left(u,\ddot{v}\right)={\omega }_{eqv1}^2-\frac{{\mathrm{\Omega }}_2^2{d}_{1}B}{A},$

(25)

and

${\varpi }_2^2=\underset{u\to \overline{u},v\to \overline{v}}{\lim }{g}_2\left(v,\ddot{u}\right)={\omega }_{eqv2}^2-\frac{{\mathrm{\Omega }}_1^2{d}_{2}A}{B}.$

(26)

The coupled linear ODEs as shown in Eqs (15) and (16) are finally transformed to decoupled linear ODEs that can be expressed as follows:

$\ddot{\mathrm{u}}+2{\mathrm{c}}_1\dot{\mathrm{u}}+{\mathrm{\varpi }}_1^2\mathrm{u}=0,$

(27)

and

$\ddot{v}+2c_2\dot{v}+{\varpi }_2^{2}v=0\hspace{0.17em}.$

(28)

To remove the damping terms $2c_1\dot{u}$ and $2c_2\dot{v}$ from the decoupled linear ODE (27) and (28), one assumes that $u\left(t\right)=e^{-c_{1}t}{R}_1\left(t\right)$ and $u\left(t\right)=e^{-c_{2}t}{R}_2\left(t\right)$ . Then, the corresponding standard normal form reads as:

${\ddot{R}}_1+{\mathrm{\Omega }}_1^2{R}_1=0,{\mathrm{\Omega }}_1^2={\varpi }_1^2-c_1^2,$

(29)

and

${\ddot{R}}_2+{\mathrm{\Omega }}_2^2{R}_2=0,{\mathrm{\Omega }}_2^2={\varpi }_2^2-c_2^2,$

(30)

which satisfy the same ICs as in coupled ODE (1) and (2) and has the following solution:

$R_1\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}{\mathrm{\Omega }}_{1}t,$

(31)

and

$R_2\left(t\right)=B\hspace{0.17em}\cos \hspace{0.17em}{\mathrm{\Omega }}_{2}t\hspace{0.17em}.$

(32)

Using equations (25) and (26) and ODEs as agreed in Eqs (29) and (30), the total frequencies ${\mathrm{\Omega }}_1,\text{and}\hspace{0.5em}{\mathrm{\Omega }}_2$ are given by

${\mathrm{\Omega }}_1^2=\frac{A\left(c_1^2-{\omega }_{eqv1}^2\right)+d_{1}B\left({\omega }_{eqv2}^2-c_2^2\right)}{A\left(d_1{d}_2-1\right)},$

(33)

and

${\mathrm{\Omega }}_2^2=\frac{B\left(c_2^2-{\omega }_{eqv2}^2\right)+d_{2}A\left({\omega }_{eqv1}^2-c_1^2\right)}{B\left(d_1{d}_2-1\right)}\hspace{0.17em}.$

(34)

Comparing the nonlinear ODE as agreed in Eq (3) with the corresponding linear ODE as agreed in Eq. (27) is helpful for validating the NPA. A selected sample of the relevant constant coefficients may be used in the following way to accomplish this goal:

$\begin{array}{l}A=0.7,B=0.8,c_1=0.05,c_2=0.3,d_{1\hspace{0.17em}}=0.05,d_2=0.08,\sigma =\pi /4,\\ f_1=0.02,f_2=0.015,{\omega }_1=0.8,{\omega }_2=0.9,{\mu }_1=0.13;\hspace{0.17em}\text{and}\hspace{0.17em}{\mu }_2=0.12.\end{array}$

The numerical solution for the nonlinear coupled ODE as agreed in Eq. (3) obtained via MS with the command NDSolve is compared with the analogous linear ODE as approved in Eq. (27) in Figure 5. These ODEs clearly indicate that the two plane curves exhibit a high degree of concordance over a long time span. Additionally, the absolute error is approximately $0.03$ , supporting the accuracy and reliability of the proposed analytical approximation. The close agreement between the solutions of the linearized ODE as settled in Eq. (27) and the nonlinear ODE as decided in Eq. (1) demonstrates how well the approach used captures the essential dynamical behaviour of the system.

Figure 5.

Shows an excellent comparison of the linear and nonlinear ODEs provided by equations (3) and (27), respectively. About $0.03$ is the maximum absolute error.

The classical idea of local linearization around an equilibrium configuration can be used to explain the remarkable agreement between the nonlinear model and its linear equivalent. The nonlinear restoring forces can be approximated in the neighbourhood of a stable equilibrium point, resulting in a linear ODE that precisely characterizes the system’s dominant dynamical properties. Consequently, the nonlinear dynamics resemble those of a linear oscillator when the system undergoes slight perturbations about the equilibrium point. NDSolve in the MS is used to numerically solve the nonlinear system to quantitatively validate the suggested method. Consequently, the nonlinear solution of ODE as granted in Eq. (4) and the corresponding linear approximation provided by ODE as approved in Eq. (28) are directly compared in Figure 6. With an absolute error of approximately $0.02$ , the results show strong agreement between the two plane curves across a large time range. This close connection demonstrates how well the linearized oscillator captures the fundamental dynamical behaviour of the system and validates the high correctness and dependability of the suggested analytical representation.

Figure 6.

Displays an excellent comparison concerning the nonlinear and linear ODEs as specified by equations (4) and (28), respectively, with maximum absolute error $\approx 0.02$ .

Additional numerical validations utilizing various sets of ICs and parameter values are shown in Figures 5 and 6 to further evaluate the robustness and reliability of the suggested methodology. These findings demonstrate that even in more challenging dynamical scenarios; the suggested method retains a high degree of accuracy.

It is helpful to compare the nonlinear ODE as shown in Eq. (4) and the corresponding linear ODE as agreed in Eq. (28) as shown Figure 6 to validate the NPA. A chosen sample of the pertinent constant coefficients may be used in the following ways to achieve this goal:

$\begin{array}{l}A=0.7,B=0.8,c_1=0.04,c_2=0.05,d_{1\hspace{0.17em}}=0.05,d_2=0.05,\sigma =\pi /4,\\ f_1=0.02,f_2=0.015,{\omega }_1=0.5,{\omega }_2=1.9,{\mu }_1=0.013;\hspace{0.5em}\text{and}\hspace{0.5em}{\mu }_2=0.01.\end{array}$

The temporal history of the Absolute errors related to the state variables $x\left(t\right)$ and $y\left(t\right)$ is shown in Figure 7. Subfigures 7(a) and 7(b) both show an oscillatory error profile that rapidly decays, indicating that the suggested approach has significant convergence features. The errors show comparatively larger amplitudes (approximately $0.025$ for $x\left(t\right)$ and $0.022$ for $y\left(t\right)$ ) at the start stage $\left(t\approx 0\right)$ , which is predicted because of IC mismatches or transient dynamics. The error magnitudes, however, drastically diminish over time in an exponentially damped oscillatory fashion. This behaviour demonstrates that both transient and steady-state mistakes are effectively suppressed by the approach. Interestingly, within a short time span (approximately $t<40$ ), the error decreases to almost zero values, after which the system maintains a minor error level. The robustness and numerical stability of the proposed method are demonstrated by the lack of persistent oscillations or divergence. Moreover, the approach appears to be devoid of numerical aberrations such as instability or noise amplification on the basis of smooth decay without unusual oscillations. The ssubfigures 7(a) and 7(b) show consistent performance across both state variables, suggesting that the technique is not biased in favor of one system component over another. The consistent convergence demonstrated the effectiveness and reliability of the suggested approach.

Figure 7.

Demonstrates the temporal evolution of the absolute errors for $x\left(t\right)$ and $y\left(t\right)$ in subfigures (a) and (b), respectively.

4. Stability configuration

The stability criterion is plotted below. The stability distributions of the nonlinear system in the $A-{\mathrm{\Omega }}_1^2$ plane for different values of the governing parameters are thus shown in Figure 8. The stable oscillatory regime is represented by the shaded regions, whereas unstable responses are represented by the unshaded parts. These stability maps show how changes in the system parameters affect the range of acceptable oscillation amplitudes by changing the stability boundary.

Figure 8.

Illustrates the stability zone of $x$ corresponding to different values of system parameters.

The impact of the nonlinear parameter ${\mu }_1$ on the stability characteristics is shown in Figure 8(a). When ${\mu }_1$ increases, the stability boundary clearly shifts upward, causing the stable region to contract. This behaviour indicates that the nonlinear damping parameter has a destabilizing effect. The system becomes more prone to instability at lower amplitude levels consequently of the nonlinear damping’s physical modification of the amplitude-dependent dissipation process, which changes the equilibrium between energy dissipation and nonlinear stiffness effects.

Figure 8(b) shows the impact of the parameter $c_1$ . The stable region expands as $c_1$ increases because the stability boundary moves downward. This pattern suggests that the oscillatory motion is stabilized by the parameter $c_1$ . Physically speaking, the presence of $c_1$ alters the system’s effective dynamic properties, enhancing its capacity to sustain constrained oscillations throughout a larger amplitude range.

Figure 8(c) shows the impact of the excitation frequency parameter ${\omega }_1$ . The stability border clearly shifts upward as ${\omega }_1$ increases. Consequently, the permissible stable region gradually shrinks, suggesting that the excitation frequency has a destabilizing effect on the dynamics of the system. Higher excitation frequencies physically change the phase connection between the excitation and the system response by intensifying the dynamic interplay between the nonlinear restorative forces and the external forcing. Consequently, the system’s effective damping capacity is diminished, increasing the vulnerability of the oscillatory motion to instability at lower amplitude levels.

Figure 8(d) shows how the external forcing parameter $f_1$ affects the stability characteristics. As $f_1$ rises, the stability boundary distinctly swings downward, leading to a noticeable expansion of the stable oscillatory region. This behaviour suggests that the system dynamics are stabilized by the external excitation parameter. Physically, the increase in $f_1$ alters the equilibrium between the energy provided by the excitation and the nonlinear restoring forces, improving the effective dynamic response and enabling the system to maintain constrained oscillations over a larger amplitude range.

Overall, the findings shown in Figure 8 illustrate how sensitive the system stability is to changes in the governing parameters. By moving the stability boundary southward, the parameters $c_1$ and $f_1$ are found to increase the admissible stable oscillatory domain and stabilize the dynamics of the system. On the other hand, increasing the parameters ${\mu }_1$ and ${\omega }_1$ causes the stability boundary to move upward, which causes the stable region to contract and ultimately encourages instability. These results demonstrate how important parameter selection is in regulating the system’s dynamic behaviour. These insights are especially helpful for the correct tuning of nonlinear systems in real-world engineering applications, such as vibration suppression, nonlinear structural dynamics, and delayed feedback control mechanisms.

Furthermore, the stability maps of the system response $y$ for various parameter values as functions of the amplitude $B$ are shown in Figure 9. The stable domain is indicated by the shaded area, whereas unstable oscillatory behaviour is represented by the area below the stability bounds. As shown in Figure 9(a), increasing the parameter ${\mu }_2$ narrows the stable region by shifting the stability border toward larger values of ${\mathrm{\Omega }}_2^2$ . This suggests that the range of parameters for which stable oscillations can occur is reduced with greater values of ${\mu }_2$ . Similar trends can be seen in Figure 9(b), where an increase in ${\omega }_2$ lifts the stability boundary and causes the allowable stability domain to contract further. These findings show that the stability characteristics are significantly influenced by both ${\mu }_2$ and ${\omega }_2$ . These findings show that the nonlinear system’s stability properties are significantly influenced by both ${\mu }_2$ and ${\omega }_2$ .

Figure 9.

Demonstrates the stability zone of $y$ corresponding to different parameters of the system.

4.1 Stability tongues

The parametric regions where nonlinear dynamical systems display frequency locking and synchronized periodic responses under external excitation are defined by stability tongues, also known as Arnold tongues. Resonance interactions between intrinsic system frequencies and stimulation give rise to these wedge-shaped formations. Their shape offers important information on energy transfer methods, modal interactions, and nonlinear coupling. Changes in system stability and bifurcation behaviour are reflected in the evolution and deformation of these regions. Consequently, stability tongues are effective tools for forecasting instability limits and synchronization thresholds. They are frequently used to improve performance, stability, and robustness in the analysis and design of nonlinear engineering systems.

The stability structure of the coupled nonlinear dynamical system of 2DOF is shown in Figure 10(a) and 10(b) in two distinct parametric planes, emphasizing the creation of Arnold tongues under changes in the excitation amplitudes $f_1$ and $f_2$ with respect to the excitation frequency $\sigma$ . The blue areas in both subfigures represent a weak system response, in which there is no frequency locking and the dynamics are primarily quasi-periodic. Wedge-shaped areas appear as the excitation amplitude increases, indicating the beginning of synchronization between the intrinsic system frequencies and the stimulation frequency. These areas, which resemble typical Arnold tongues, gradually increase as the forcing intensity increases, indicating improved resonance capture and stability locking.

Figure 10.

Investigates stability tongues (Arnold tongues) of coupled nonlinear dynamical system of 2DOF in the $\left(f_1,\sigma \right)$ and $\left(f_2,\sigma \right)$ parameter spaces as shows in subfigures 10(a) and 10(b), respectively.

In Figure 10(a), (a)well-defined and broad Arnold tongue structure is observed for the system parameters:

$c_1=0.15,c_2=0.14,d_{1\hspace{0.17em}}=0.15,d_2=0.2,f_2=0.01,{\omega }_1=1.67,{\omega }_2=0.4,{\mu }_1=0.03,{\mu }_2=0.15.$

The synchronization region grows asymmetrically as $f_1$ increases, starting at a crucial frequency close to the first mode’s natural frequency. The surrounding yellow region depicts intermediate nonlinear dynamics linked to partial synchronization and the transition to resonance, whereas the red region shows severe frequency locking and large-amplitude periodic responses. This prominent asymmetry and downward bending of the tongue’s lower border are noteworthy features. Strong nonlinear interactions and modal coupling effects are suggested by this behaviour, which causes the resonance condition to change. The parametric excitation applied through the first mode greatly increases the resilience of synchronization and expands the stability domain by widening the tongue with increasing $f_1$ .

In contrast, Figure 10(b) shows a noticeably narrower and more confined Arnold tongue structure for the system parameters:

$c_1=0.15,c_2=0.14,d_{1\hspace{0.17em}}=0.1,d_2=0.45,f_1=0.01,{\omega }_1=0.4,{\omega }_2=1.6,{\mu }_1=0.15,{\mu }_2=0.05.$

As $f_2$ increases, the expansion rate decreases, and the synchronization region becomes more localized. This suggests that, in comparison with the first mode, stimulation via the second mode is less successful in producing strong resonance. The red area is still small, indicating poorer energy transfer efficiency and weaker frequency locking. Furthermore, the tongue is less deformed and more symmetric, suggesting that coupling contributions and nonlinear effects are less significant in this arrangement. The system transitions into synchronized motion more slowly and weakly, as further supported by the decreased spread of the yellow transition region.

The influence of the two stimulation channels clearly differs when Figure 10(a) and 10(b) are directly compared:

• Stronger, wider, and more distorted Arnold tongues are produced by excitation via $f_1$ (Figure 10(a)), suggesting dominant nonlinear coupling and improved resonance processes.

• Weaker and more confined synchronization regions result after excitation via $f_2$ (Figure 10(b)), indicating the secondary involvement of the second mode in controlling the system dynamics.

This asymmetry draws attention to the uneven contribution of the coupled modes and implies that the first mode serves as the main channel for energy injection and synchronization, with the second mode having a more passive or secondary function.

The observed behaviour shows that transitions between quasi-periodic motion, partial synchronization, and strong nonlinear resonance are governed by the interaction between the excitation amplitude and frequency. Arnold tongue expansion and deformation are unmistakable markers of the system’s energy transfer efficiency, modal dominance, and nonlinear coupling strength. Overall, the findings demonstrate that the path to synchronization is strongly influenced by the mode through which the parametric excitation is supplied and offers a thorough global understanding of the system’s bifurcation structure.

4.2 Stability multipliers via floquet theory

A basic framework for examining the stability of linear ODEs with periodic coefficients is provided by Floquet theory. By calculating Floquet multipliers, which establish whether perturbations increase or decrease over time, the behaviour of the system can be described. Without the need for lengthy simulations, this method allows for the accurate identification of stable and unstable periodic solutions. When examining parametric resonance, bifurcations, and time-periodic dynamical systems, the Floquet theory is especially effective. This approach provides an in-depth understanding of how system factors affect transition mechanisms and stability bounds. Consequently, it is frequently used in physics and engineering to forecast, regulate, and maximize the stability of oscillatory and periodically driven systems.

Furthermore, the present coupled nonlinear dynamical system of 2DOF can be linearized via Eqs. (13) and (14), leading to a system of linear ODEs with time-periodic coefficients as follows:

$\ddot{x}+d_1\ddot{y}+2c_1\dot{x}+\left({\varpi }_{eqv1}^2-f_1\hspace{0.17em}\cos \hspace{0.17em}2\sigma \hspace{0.17em}t\right)x=0,$

(35)

and

$\ddot{y}+d_2\ddot{x}+2c_2\dot{y}+\left({\varpi }_{eqv2}^2-f_2\hspace{0.17em}\cos \hspace{0.17em}2\sigma \hspace{0.17em}t\right)y=0.$

(36)

Floquet stability analysis is performed via linearization of the nonlinear ODEs about the periodic solution. As such, the Floquet multipliers can be derived using the monodromy matrix and used in assessing the stability properties of the system. In this regard, the system is considered stable provided that all the Floquet multipliers remain inside the unit circle of the complex plane; if at least one of them lies beyond it, the system will be deemed unstable.

A generalized Mathieu-type system is comparable to such systems. Consequently, a rigorous framework for analyzing the stability of the ensuing periodic solutions is provided by Floquet theory. By adding the state variables, the governing ODEs are rewritten in first-order form:

${\lambda }_1=x,{\lambda }_2=\dot{x},{\lambda }_3=y,{\lambda }_4=\dot{y}.$

(37)

Thus, the system can be expressed in compact form as follows:

$\frac{d\lambda }{dt}=A\left(t\right)\lambda ,$

(38)

where $\lambda =\left({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4\right)$ and $A\left(t\right)$ is a time-periodic matrix with period $T=2\pi /\sigma$ and is given by,

$A\left(t\right)=\left(\begin{array}{llll}0& 1& 0& 0\\ a_{21}\left(t\right)& a_{22}\left(t\right)& a_{23}\left(t\right)& a_{24}\left(t\right)\\ 0& 0& 0& 1\\ a_{41}\left(t\right)& a_{42}\left(t\right)& a_{43}\left(t\right)& a_{44}\left(t\right)\end{array}\right),$

(39)

where the matrix elements are given by:

$a_{21}\left(t\right)={\mathrm{\Gamma }}^{-1}\left(-{\omega }_{eqv1}^2+f_1\hspace{0.17em}\cos \hspace{0.17em}2\sigma t\right),a_{22}=-2c_1{\mathrm{\Gamma }}^{-1},$

(40)

$a_{23}\left(t\right)=d_1{\mathrm{\Gamma }}^{-1}\left({\omega }_{eqv2}^2-f_2\hspace{0.17em}\cos \hspace{0.17em}2\sigma t\right),a_{24}=2c_2{d}_1{\mathrm{\Gamma }}^{-1},$

(41)

$a_{41}\left(t\right)=d_2{\mathrm{\Gamma }}^{-1}\left({\omega }_{eqv1}^2-f_1\hspace{0.17em}\cos \hspace{0.17em}2\sigma t\right),a_{42}=2c_1{d}_2{\mathrm{\Gamma }}^{-1},$

(42)

and

$a_{43}\left(t\right)={\mathrm{\Gamma }}^{-1}\left(-{\omega }_{eqv2}^2+f_2\hspace{0.17em}\cos \hspace{0.17em}2\sigma t\right),a_{44}=-2c_2{\mathrm{\Gamma }}^{-1},$

(43)

with $\mathrm{\Gamma }=1-d_1{d}_2$ .

The monodromy matrix is the state transition matrix calculated for one excitation period and is given by:

$M=\mathrm{\Phi }\left(T\right),$

(44)

where the fundamental solution matrix $\mathrm{\Phi }\left(T\right)$ is defined as:

$\dot{\mathrm{\Phi }}\left(t\right)=A\left(t\right)\mathrm{\Phi }\left(t\right),\mathrm{\Phi }\left(0\right)=I_4,$

(45)

where $I_4$ is the square identity matrix of order $4\times 4$ and the periodic time $T=2\pi /\sigma$ .

Owing to the nonlinear and coupled nature of the system, monodromy matrix cannot be derived in analytical form and is therefore computed numerically by integrating the variational system over one excitation period. The computed Floquet multipliers are as follows:

${\rho }_{1,2}=0.521686 \pm 0.52577\mathrm{i},{\rho }_{3,4}=-0.691398\pm 0.214363\mathrm{i},$

(46)

that satisfies

$\left|{\rho }_i\right|<1;\hspace{0.5em}i=1,2,3,4.$

(47)

This confirms that the periodic response is asymptotically stable, as any infinitesimal perturbation diminishes exponentially with time.

The stability of the periodic solutions is quantitatively assessed via Floquet theory. The linearized variational equations are integrated over one excitation period ( $T=2\pi /\sigma$ ) to create the analytic continuation behaviour matrix. The Floquet multipliers, which describe how perturbations change over subsequent periods, are defined by the eigenvalues of this matrix. Instability occurs when at least one Floquet multiplier in the complex plane intersects the unit circle. Floquet multipliers obtained through calculation have been observed to be conjugate pairs of complex numbers with real and imaginary parts less than unity, lying strictly within the unit circle in the complex plane, as illustrated in Figure 11. Thus, the periodic solution is stable. From a physical standpoint, this suggests that minor perturbations, such as brief excitations or changes in the surroundings, gradually diminish with time, restricting the associated oscillatory motion. Moreover, heavily damped oscillatory disturbances are reflected by complex-conjugate Floquet multipliers with very modest magnitudes, demonstrating the efficiency of the system damping and energy dissipation mechanisms under the examined excitation conditions.

Figure 11.

Depicts the Floquet on the complex plane.

5. Bifurcation analysis

In this Section, we investigated the analysis of the BDs and LLE,^52–54 which is essential in understanding the chaotic behaviour of dynamical systems, as these methods offer valuable insights into their intricate dynamics. BDs illustrate how a system transitions between states, such as periodic, periodic doubling, quasi-periodic, or chaotic, by mapping changes in behaviour as parameters vary, helping identify stability boundaries.

In this study, we took the excitation amplitudes $f_1$ and $f_2$ as bifurcation parameters to acquire the needed BDs and LLE. The upcoming Figs. provided illustrate the transition of a dynamic system from stable, predictable behavior to a state of total chaos as the control parameter $f_1$ increases. In Figure 12, the bifurcation diagrams for $x\left(t\right)$ and $y\left(t\right)$ act as a visual map of this journey. Initially, for values of $f_1$ between 0 and approximately 2, the system stays on a very narrow, singular path, indicating a stable periodic state. However, once the excitation amplitude $f_1$ crosses the threshold near 2.2, these lines “explode” into a dense cloud of points. This sudden expansion represents a BD, where the system no longer follows a single predictable loop but instead begins to oscillate across a wide range of values in an unpredictable, chaotic manner.

Figure 12.

Shows BDs of (a) $x\left(t\right)$ and (b) $y\left(t\right)$ vs. $f_1$ .

Figure 13 provides the mathematical validation for this transition by examining the LLE and classifying the motion types. In the first plot, the LLE starts below zero, which confirms that the system is naturally stable and resistant to small disturbances. As $f_1$ approaches 2, the LLE trends upward and eventually crosses the zero line into positive territory. A positive LLE is the “smoking gun” for chaos, indicating that even the tiniest change in starting conditions will cause the system’s trajectory to diverge exponentially. The accompanying classification chart simplifies this by color-coding the behavior: green represents the initial periodic stability, blue marks the transitionary “quasi-periodic” phase, and red highlights the dominant chaotic regime, though small “windows of order” still occasionally reappear even at high $f_1$ values.

Figure 13.

Shows a simulation of (a) LLE and its motion type classification in (b).

Moreover, a two-dimensional chaos map based on the KYD provides a global view of how the system behaves when two parameters vary simultaneously, rather than along a single path as in classical bifurcation diagrams. It allows us to clearly identify regions of periodic, quasi-periodic, and chaotic motion, while also quantifying the intensity of chaos through the attractor’s fractal dimension. In general, when the KYD is close to 1, the system exhibits stable periodic motion; values around 2 indicate quasi-periodic behaviour on a torus; while values greater than 2 reflect chaotic dynamics with increasing complexity as the dimension grows. This makes it a powerful tool in understanding the interplay between excitation and nonlinearity, and for selecting parameter ranges that yield stable or desirable system performance. So, Figure 14 expands this analysis into a two-dimensional chaos map using the KYD to show how the system evolves when both $f_1$ and the nonlinear parameter ${\mu }_1$ are adjusted. This heat map visualizes the “complexity” of the chaos. The deep blue regions on the left correspond to a low KYD, representing simple, non-chaotic behavior. As you move to the right and $f_1$ increases, the map shifts into vibrant teals, yellows, and reds. These warmer colors indicate a higher fractal dimension, meaning the system’s chaotic attractor has become more structurally complex and “fills” more of the available phase space. While $f_1$ is clearly the primary driver of the chaos, the vertical axis ( ${\mu }_1$ ) subtly shifts the intensity and exact boundaries of these chaotic zones, creating a complex landscape of unpredictability.

Figure 14.

Shows 2D Chaos map using KYD.

On the other hand, the new set of Figures. 15–17 reveal a system that is significantly more “stubborn” and resistant to chaos than the one we looked at previously. In Figure 15, the BDs show that for a large portion of the range, specifically from $f_2=0.0$ all the way to about 6, the system manages to stay relatively contained. While there is a slight thickening of the lines after $f_2=2$ , suggesting some complexity, the dramatic “ink-spill” effect of total chaos does not really take over until much later. Even when the system breaks down into chaos around the $f_2=6$ to 8 mark, it appears much more fragmented and less “dense” than before, frequently trying to pull itself back into a stable rhythm.

Figure 15.

Shows BDs of (a) $x\left(t\right)$ and (b) $y\left(t\right)$ vs. $f_2$ .

Figure 16.

Shows simulation of (a) LLE and its motion type classification in (b).

Figure 17.

Displays 2D Chaos map using KYD.

The data in Figure 16 confirms that this system is more stable. The LLE remains mostly in negative values, only briefly rising above zero during sharp spikes. This is clearly shown in the motion classification chart, which is mostly green—a sign of periodic (stable) behavior. Unlike the initial case, which was filled with red chaotic points, this system shows only a few “mood swings” into quasi-periodic (blue) and chaotic (red) states. Likewise, the 2D Chaos Map in Figure 17 appears overall “cooler”. Large areas are deep blue, with the high-complexity chaotic regions forming narrow, vertical “flames” instead of expansive heatwaves, as seen in the previous map.

When we compare the two cases, the difference in “volatility” is striking. The first system (controlled by $f_1$ ) was highly sensitive; once it hit the threshold of 2, it essentially surrendered to chaos and stayed there. It was a high-energy, highly unpredictable environment. In contrast, this second system (controlled by $f_2$ ) acts more like a dampened or more robust mechanism. It has a much higher “breaking point” and exhibits what we call intermittent chaos; it does not stay chaotic for long before settling back into a predictable periodic state. If the first case was a storm that never ended, this second case is more like a generally calm day with a few sudden, short-lived gusts of wind.

Knowing that the nonlinear governing equations were numerically solved using the adaptive fourth–fifth order Runge–Kutta method (ode45) in MATLAB. The excitation period was defined as $T=\frac{2\pi }{s}$ . The total simulation time for the bifurcation and Lyapunov analyses was taken as $1200T$ , while approximately 70% of the computed response was discarded to remove transient dynamics. For the two-dimensional chaos maps, the transient response was first eliminated over $150T$ . In the first case, the bifurcation parameter $f_1$ was varied within $0\le f_1\le 10$ using resolutions of $\mathbf{\Delta }{f}_1=0.03$ for the one-parameter bifurcation diagrams and $\mathbf{\Delta }{f}_1=0.1$ for the two-dimensional chaos maps, while ${\mu }_1$ was varied within $0.05\le {\mu }_1\le 0.25$ using $\mathbf{\Delta }{\mu }_1=0.02$ . Similarly, in the second case, the same numerical procedure was repeated using $f_2$ and ${\mu }_2$ as the varying control parameters.

6. Conclusions

This study aimed to examine the issue of a cart on a horizontal plane coupled with a damped nonlinear spring on an inclined surface. The study investigated the interplay between linear/nonlinear dynamics within a system, concentrating on a cart, inclined surface, and sliding disc, and offering insights into stability, control, and energy dynamics in both theoretical and practical situations. This dynamical system comprises 2DOFs. The methodology implemented through the NPA. The aim of NPA was to transform the governing ODEs into linear ones. MS was employed to validate the significant accuracy between the nonlinear and linear/nonlinear ODEs. The NPA was completely distinct from all traditional perturbations. It was necessary to analyse the stability criterion in both the theoretical and numerical simulations. The NPA was distinguished by its simplicity, attractiveness, potential, and effectiveness, rendering it beneficial for weakly nonlinear ODEs spanning several categories within dynamic systems as well as fluid mechanics. The temporal variation in the results is depicted to emphasize the influence of several impacted parameters. The chaotic behaviour of the system was investigated via BDs and LLE. The analysis revealed that the system transitions from periodic motion to a chaotic state. Additionally, the complexity of the system is illustrated using two-dimensional chaos maps derived from the KYD. These maps clearly identify the specific parameter changes at which the system transitions from regular, predictable behaviour to intricate, high-dimensional chaotic dynamics.

As future work, we aim to apply the NPA to the following works:

The (2+1)-dimensional nonlinear damped Klein–Gordon Fock equation was analysed via the classical Lie symmetry approach. ³¹ A comprehensive discrimination framework of such systems facilitates qualitative and quantitative evaluations of a given observable model. The method’s efficacy was substantiated by precise findings obtained from the model under the requisite constraint conditions. ³² The comprehensive version of the Zakharov system, delineating dispersive and ion acoustic wave propagation in plasma was analysed. ³³ The dynamics of nonlinear waves in a modified Zakharov–Kuznetsov equation and their interactions with discrete electric lattice structures were investigated. ³⁴