Examining the behavior of a two-connected-spring pendulum according to non-resonance and resonance conditions

Abstract

This paper focuses on the planar forced oscillations of a particle linked to its support point through a system comprising two nonlinear springs arranged in series with two viscous dampers. The motion of this point is restricted to being on an elliptic route. Three degrees-of-freedom (DOF) characterize the system’s motion, explained by two differential equations (DEs) and another algebraic equation derived using Lagrange’s equations. The traditional method of multiple-time scales (MMTS) is applied in the context of the time realm. The dynamics of forced and damped oscillations is explored in this work, highlighting non-resonant conditions along with two distinct external resonances. Furthermore, stationary periodic states influenced by external resonances are examined, and the stability of each case is evaluated. The asymptotic solutions of the controlling equations are solved up to the third approximation. The classification of resonance cases is based on the second and third order of the approximate solutions. Therefore, the solvability criteria and the modulation equations (ME) are obtained. The Routh–Hurwitz conditions (RHCs) serve as the basis for assessing the stability and instability criteria. A discussion and graphical representations of the time histories, frequency response curves (FRC), and stability regions are provided to demonstrate the beneficial effects of different physical parameter inputs on the system’s performance. The obtained outcomes can improve the performance of spring-based mechanical systems like vibrational isolation devices, shock absorbers, and suspension mechanisms. Additionally, they have the potential to enhance stability and control in robotic limbs and aircraft landing gear systems, ensuring better stress absorption and stability.

Introduction

The study of nonlinear dynamics in spring pendulum systems has garnered significant interest due to its rich mathematical structure and wide-ranging applications in engineering, physics, and beyond. Nonlinear systems, characterized by their sensitivity to initial conditions and complex behavior, often require sophisticated analytical techniques to understand and predict their dynamics. A notable example of such a system is the spring pendulum, which exhibits complex motions such as chaos and resonance under various excitations. The MMTS is a perturbation approach for analyzing systems with oscillatory dynamics across varying time scales. It involves introducing separate time variables for the fast and slow dynamics of the system. By expanding the solution as a series in terms of a small parameter, the method helps to systematically solve for the behavior at each time scale, allowing for a more accurate approximation of solutions, especially in nonlinear systems or near resonances. This approach is beneficial for handling complex interactions between fast and slow oscillations.

Early investigations into the chaotic responses of harmonic excited spring pendulums (SPs) provided foundational insights into the behavior of these systems. Lee and Park ¹ developed a second-order approximation method to analyze chaotic responses in SP systems, highlighting the intricate interplay between harmonic excitation and system nonlinearity. Building on this, in Ref. 2, the authors explored vibration reduction in a 3DOF nonlinear SP, demonstrating techniques to mitigate undesirable oscillations through strategic system design. Moreover, it was focused on managing the examined model at the primary resonance, with its stability using both the phase plane graphs and FRC. The chaotic response of an SP whose pivot point moves along a circular path at a consistent angular speed was examined, emphasizing the diverse dynamic behavior that arises under different motion constraints. ³ In Ref. 4, the authors delved into asymptotic analysis near resonances in kinematically excited dynamical systems, providing a deeper understanding of the resonance phenomena. ⁴ For the scenario of a stationary pivot point of an SP, their work in Ref. 5 extended this analysis to both stationary and transient resonant responses, offering comprehensive insights into the SP’s dynamic stability when the motion was affected by radial and transverse forces.

Analytical approaches to investigating spring pendulum dynamics have evolved significantly, as the authors in Ref. 6 presented two distinct analytical methods for studying the system, showcasing the versatility and depth of analytical techniques in addressing complex dynamical problems. Perturbation methods, as detailed in Ref. 7, remain crucial tools in the analysis of nonlinear systems, allowing for the systematic approximation of solutions in the presence of small parameters.

Several studies have continued to expand our understanding of the dynamics of SP, for example, Refs. 8–11. In Ref. 8, the dynamics of rigid bodies in relative equilibrium positions, providing a broader context for pendulum-like systems, was explored in Ref. 8. In Ref. 9, the movement of a 3DOF dynamical system near resonances, contributing to the understanding of multi-DOF systems’ behavior in resonant conditions, was analyzed. In Ref. 10, the authors likewise examined a 2DOF dampened SP, where its support point is confined to an elliptic path, furthering the knowledge of complex oscillatory systems. In Ref. 11, a restriction was considered for the motion of a simple pendulum’s pivot in an elliptical trajectory. An absorber was considered to connect with the pendulum’s rigid arm in its longitudinal direction. The movement was studied considering the effects of two harmonic excitation forces and a single viscous moment. The MMST was used to construct the analytic solutions of these works up to various degrees of approximations, in which the used time has been split into several time scales; the first one is known as a fast time scale, while the others are defined by the slow ones. Secular terms were avoided from the solutions, to be valid through the examined time interval, to achieve the solvability criteria.

According to Ref. 12, the authors achieved the decoupling of fast and slow scales by incorporating local averages and substituting fast scale contributions with time-localized periodic problems. An a posteriori error estimator was devised using the method of dual-weighted residual, facilitating the decomposition of errors into averaging error, slow scale error, and fast scale error. The accuracy of the error estimator was validated by the authors, who also showcased its application in the adaptive control of a numerical MMST. However, in Ref. 13, the authors introduced a method utilizing MMST, which combines the homotopy perturbation method^14,15 and the classic perturbation method. This approach has proven to be highly effective for solving various nonlinear equations, especially in the context of forced nonlinear oscillators. The MMST was applied in Ref. 16 to systematically simplify thermo-mechanical structural models where temperature distributions change slowly. The study in Ref. 17 conducted a bifurcation analysis on a nonlinear railway vehicle equipped with dual bogies to examine the impact of bogie interaction on the vehicle’s hunting behavior. Solutions near the hunting speed were derived using an asymptotic expansion of a small perturbation parameter within the framework of the MMST.

Further works on this topic are outlined in Refs. 18–24. The case of a damped SP rather than an unstretched one was explored in Ref. 18. Mechanical energies are calculated and solved using a numerical method to analyze the regions of stability and instability concerning RHC. Additionally, constrained motions are taken into account, with the spring’s suspension point tracing a Lissajous pattern. The study in Ref. 19 employed the MMST to analyze a two-scale command shaping feed-forward control strategy designed to minimize residual oscillations in both conventional and unconventional Duffing systems. The stationary and nonstationary oscillatory behaviors of a parametric pendulum by applying the concept of limiting phase trajectories were investigated in Ref. 20, in which a simplified model that could predict highly modulated patterns beyond the usual range of initial conditions has been suggested. The solvability criteria and a multi-scale dynamic analysis of one-dimensional structures with nonuniform boundaries, facilitating the examination of strings, cables, beams, and arches, were established in Ref. 21. In Ref. 22, perturbation techniques combined with the MMST were used to derive and analyze a simplified third-order model of microscale oscillators with thermos-optical feedback, and the authors also provided the system’s bifurcation diagram. Ref. 23 focused on examining the planar forced oscillations of a particle connected to a support via two nonlinear springs and viscous dampers. The study utilized the MMST in the time domain. However, a modified version of this approach with three scales has been developed with new algorithms specifically tailored to solve problems involving both algebraic and differential equations. However, in Ref. 24, the authors examined how vibrating rods with nonlinear array springs respond harmonically. The analysis was multi-mode, with response functions plotted across broad frequency ranges to observe various resonances.

In Refs. 25–29, the authors used the amplitude-frequency formulation for nonlinear oscillators to uncover the key mechanism behind pseudo-periodic motion. They demonstrated that quadratic nonlinear forces play a role in the pull-down phenomenon during each cycle of periodic motion. When these forces exceed a certain threshold, pull-down instability arises. The authors examined the periodic vibration characteristics of the Duffing oscillator under generalized initial conditions. ²⁶ In Ref. 27, the authors conducted a numerical analysis of fractional nonlinear vibrations in a restrained cantilever beam with an intermediate lumped mass. This study focused on the complex vibrational dynamics that arise in such systems due to the nonlinearity introduced by the fractional-order derivatives, contributing to the field of low-frequency noise and active control.

Further advancing the exploration of nonlinear oscillators, an analytical approach is developed in Ref. 28 based on the global residue harmonic balance method to investigate the nonlinear differential equations of the circular sector oscillator, a system characterized by complex geometrical and physical constraints.

The vibration frequencies of tapered beams, used in many structural applications, have also been extensively studied. The spreading residue harmonic balance method is employed in Ref. 29 to analyze the vibrational behavior of tapered beams, offering a detailed examination of how geometry influences vibrational frequencies.

This study explores the movement of a system with two nonlinear springs arranged in series, exposed to external forces in both longitudinal and transverse directions. The system operates under the assumption that the suspension point of one spring traces an elliptical path with a stationary angular velocity. By employing the second type of Lagrange’s equations, the system’s equations of motion (EOM) are derived and solved asymptotically using the MMST up to the approximation of third-order. The determination of solvability conditions involves removing secular terms, which leads to equations for modulation in various resonance scenarios. The RHC is utilized to analyze the stability of fixed points at steady-state solutions, investigating regions of stability and instability. Plots of time histories, resonance, and stability boundaries are used to depict the system’s behavior.

The obtained results can be applied to enhance the design and efficiency of mechanical systems such as conveyor belt systems, hydraulic dampers, and seismic isolators that rely on spring dynamics. Such advancements contribute to smoother, more precise operations, minimizing wear and tear while prolonging the service life of components. Moreover, they can improve the precision and reliability of automated manufacturing equipment and prosthetic devices that use spring mechanisms to mimic natural joint movements, offering more lifelike motion and increased functionality in dynamic environments. Additionally, optimizing automotive suspension systems leads to smoother rides, increased stability on uneven roads, and better handling, all of which enhance both safety and passenger comfort.

Beyond these uses, the results can be applied to aerospace engineering, where the precise regulation of spring-damped systems is vital for landing gear performance, vibration suppression, and space mission equipment. In robotics, enhanced spring dynamics can enable more adaptive and responsive systems, benefiting tasks that require fine motor control, such as robotic surgery or dexterous handling in hazardous environments. Finally, these innovations can also enhance energy harvesting systems by increasing the efficiency of devices that transform mechanical vibrations into electrical power, facilitating the development of sustainable, self-powered technologies for remote or difficult-to-reach areas.

Description of the vibrating system

This section presents a full description of the structure of the examined dynamical system and the derivation of its EOM is the dimensionless form. This system consists of two tiny sliding particles on a bar, with masses $m_0$ and $m$ , in which the motion is limited to be in the vertical plane, as depicted in Figure 1. These particles are connected by a spring, and the pin point of the bar $O_1$ is restricted to move in an elliptic trajectory with the constancy of its angular velocity $\Omega$ . The major and minor axes of this trajectory have semi-lengths of $a$ and $b$ to axes $OX$ and $OY$ . It must be noted that on the auxiliary circle $a$ , the point $Q$ corresponds to $O_1$ .

Figure 1.

Non-linear dynamical system.

One of these masses is connected to a point $O_1$ via the second spring, where the springs’ lengths at $t=0$ are $l_{01}$ and $l_{02}$ , respectively. The system’s motion is subjected to two applied harmonic forces $F_s\left(t\right)=F_s^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{s}t\right);(s=1,2)$ on the mass $m$ , along the radial and transversal directions of the spring, and an anticlockwise harmonic moment $M\left(t\right)=M_0^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{0}t\right)$ at $O_1$ . Here ${\Omega }_s$ , ${\Omega }_0$ , and $F_s^{*}$ , $M_0^{*}$ , are the frequencies and amplitudes of $F_s\left(t\right)$ and $M\left(t\right)$ . It is presumed that the nonlinear effects on the total elasticity force are negligible, and the bar, the springs, and the dampers are considered to be massless. In this system, the coefficients of the two viscous dampers are $C_1$ and $C_2$ .

It is noted that 3DOF are available for the system, where the balls’ positions $X_s\left(t\right)$ and the rotatory angle $\phi \left(t\right)$ are defined by the system’s generalized coordinates. For every spring, the correlation between axial extension and the elastic force $R_s$ can be expressed as follows:

$R_s=k_s\left({\mathrm{\Delta }}_s+{\Lambda }_s{\mathrm{\Delta }}_s^3\right),s=1,2$

(1)

where $k_s$ and ${\Lambda }_s$ are the elastic coefficients of the springs and ${\mathrm{\Delta }}_s$ is its total elongation.

The following is the system’s kinetic $T$ and potential $V$ energies:

$\begin{array}{l}T=\frac{1}{2}{m}_0[{\Omega }^2\left(a^2\hspace{0.17em}{\mathit{cos}}^2\hspace{0.17em}\Omega t+b^2\hspace{0.17em}{\mathit{sin}}^2\hspace{0.17em}\Omega t\right)+{\left(L_{01}+X_1\right)}^2{\phi }^{\prime 2}+{X^{\prime }}_1^2+2\Omega \left(L_{01}+X_1\right){\phi }^{\prime }\times \left(a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \right)+2\Omega X_1^{\prime }\left(a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi -b\hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \right)+\frac{1}{2}m[{\Omega }^2\left(a^2\hspace{0.17em}{\mathit{cos}}^2\hspace{0.17em}\Omega t+b^2\hspace{0.17em}{\mathit{sin}}^2\hspace{0.17em}\Omega t\right)+{\left(L_{01}+L_{02}+X_2\right)}^2{\phi }^{\prime 2}+{X^{\prime }}_2^2+2\Omega \left(L_{01}+L_{02}+X_1\right)\times {\phi }^{\prime }\left(a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi \right)+2\Omega X_2^{\prime }\left(a\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{sin}\hspace{0.17em}\phi -b\hspace{0.17em}\mathit{sin}\hspace{0.17em}\Omega t\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \right),\\ V=\frac{1}{2}{k}_1{X}_1^2\left(1+\frac{1}{2}{\Lambda }_1{X}_1^2\right)-m_{0}g\left[\left(L_{01}+X_1\right)\mathit{cos}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t\right]+\frac{1}{2}{k}_2{\left(X_2-X_1\right)}^2\times \left[1+\frac{1}{2}{\Lambda }_2{\left(X_2-X_1\right)}^2\right]-mg\left[\left(L_{01}+L_{02}+X_2\right)\mathit{cos}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{cos}\hspace{0.17em}\Omega t\right].\end{array}$

(2)

Here, the derivatives regarding time $t$ are represented by dashes.

Referring to this equation, we can establish the system’s Lagrangian $L=T-V$ , leading to the subsequent derivation of the mentioned equations:

$\frac{d}{dt}\left(\frac{\partial L}{\partial {\rho }^{\prime }}\right)-\frac{\partial L}{\partial \rho }=Q_{\rho };\rho \left(=X_1,X_2,\phi \right),$

(3)

where $Q_{\rho }$ represent the generalized forces, which consist of the harmonic forces and the damping ones. These forces have the following forms:

$\begin{array}{l}{Q}_{X_1}=0,\\ Q_{X_2}=F_1^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{1}t\right)-C_1{X}_2^{\prime },\\ Q_{\phi }=M_0^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{0}t\right)+\left(L_{01}+L_{02}+X_2\right)F_2^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{2}t\right)-C_2{\phi }^{\prime }.\end{array}$

(4)

Furthermore, the symbols $\left(X_1,X_2,\phi \right)$ and $\left({X_1}^{\prime },{X_2}^{\prime },{\phi }^{\prime }\right)$ denote, respectively, the generalized coordinates and velocities.

Substituting (2) and (4) into (3), we obtain the EOM:

$m_0\left[X_1^{″}-{\Omega }^2\left(a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{cos}\hspace{0.17em}\phi \right)-\left(L_{10}+X_1\right){\phi }^{\prime 2}-g\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \right]+k_1\left(X_1+{\Lambda }_1{X}_1^3\right)-k_2\left[\left(X_2-X_1\right)+{\Lambda }_2{\left(X_2-X_1\right)}^3\right]=0,$

(5)

$m\left\{X_2^{″}-{\Omega }^2\left[a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{cos}\hspace{0.17em}\phi \right]-\left(L_{10}+L_{20}+X_2\right){\phi }^{\prime 2}-g\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \right\}+C_1{X}_2^{\prime }+k_2\left[\left(X_2-X_1\right)+{\Lambda }_2{\left(X_2-X_1\right)}^3\right]=F_1^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{1}t\right),$

(6)

$m_0\{{\left(L_{10}+X_1\right)}^2{\phi }^{″}+2\left(L_{10}+X_1\right)X_1^{\prime }{\phi }^{\prime }+{\Omega }^2\left(L_{10}+X_1\right)[b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi -a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\times \mathit{cos}\hspace{0.17em}\phi ]+g\left(L_{10}+X_1\right)\mathit{sin}\hspace{0.17em}\phi \}+m\{{\left(L_{10}+L_{20}+X_2\right)}^2{\phi }^{″}+2\left(L_{10}+L_{20}+X_2\right)X_2^{\prime }{\phi }^{\prime }+{\Omega }^2\left(L_{10}+L_{20}+X_2\right)\left[b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi -a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\mathit{cos}\hspace{0.17em}\phi \right]+g\left(L_{10}+L_{20}+X_2\right)\mathit{sin}\hspace{0.17em}\phi \}{C}_2{\phi }^{\prime }=\left(L_{10}+L_{20}+X_2\right)F_2^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{2}t\right)+M_0^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{0}t\right).$

(7)

According to equations (5)–(7), the governing EOM for the spring pendulum consists of connected two nonlinear springs in series, can be derived. Figure 2 illustrates such a system, which represents the paper’s topic. The symbol $C$ designates the spring’s connecting point, which can be thought of as the particle’s massless counterpart with a mass that tends toward zero. As a result, the point maintains its mobility along the bar. Formally, the only difference between the drawn systems in Figures 1 and 2 is the inclusion of the ball with mass $m_0$ . Accordingly, if we assume that $m_0=0$ in equations (5)–(7), we obtain the following:

$k_1\left(X_1+{\Lambda }_1{X}_1^3\right)-k_2\left[X_2-X_1+{\Lambda }_2{\left(X_2-X_1\right)}^3\right]=0,$

(8)

$m\left\{X_2^{″}-{\Omega }^2\left[a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi +b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{cos}\hspace{0.17em}\phi \right]-\left(L_{10}+L_{20}+X_2\right){\phi }^{\prime 2}-g\hspace{0.17em}\mathit{cos}\hspace{0.17em}\phi \right\}+C_1{X}_2^{\prime }+k_2\left[\left(X_2-X_1\right)+{\Lambda }_2{\left(X_2-X_1\right)}^3\right]=F_1^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{1}t\right),$

(9)

$m_0\{{\left(L_{10}+X_1\right)}^2{\phi }^{″}+2\left(L_{10}+X_1\right)X_1^{\prime }{\phi }^{\prime }+{\Omega }^2\left(L_{10}+X_1\right)[b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi -a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\mathit{cos}\hspace{0.17em}\phi ]+g\left(L_{10}+X_1\right)\mathit{sin}\hspace{0.17em}\phi \}+m\{{\left(L_{10}+L_{20}+X_2\right)}^2{\phi }^{″}+2\left(L_{10}+L_{20}+X_2\right)X_2^{\prime }{\phi }^{\prime }+{\Omega }^2\left(L_{10}+L_{20}+X_2\right)[b\hspace{0.17em}\mathit{cos}\left(\Omega t\right)\mathit{sin}\hspace{0.17em}\phi -a\hspace{0.17em}\mathit{sin}\left(\Omega t\right)\mathit{cos}\hspace{0.17em}\phi ]+g\left(L_{10}+L_{20}+X_2\right)\mathit{sin}\hspace{0.17em}\phi \}+C_2{\phi }^{\prime }=\left(L_{10}+L_{20}+X_2\right)\times F_2^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{2}t\right)+M_0^{*}\hspace{0.17em}\mathit{sin}\left({\Omega }_{0}t\right).$

(10)

Figure 2.

Three DOF dynamical system of a serial pendulum’s configuration.

The equilibrium circumstance for both springs’ forces is described by the algebraic equation (8). The differential-algebraic equation (8) through (10) regulates a class of dynamical systems, which represents the corresponding model of these equations. For computational reasons, eliminating the $X_1$ coordinate strictly is not a good idea; therefore, both coordinates are required to clearly represent the system’s state. If equation (8) represents constraints, the coordinate $X_1$ can be considered complementary rather than independent. Consequently, the mechanism depicted in Figure 2 has 3DOF.

The dimensionless version of the controlling framework for EOM is structured as follows:

$\begin{array}{l}{x}_1+{\alpha }_1{x}_1\left(x_1^2+3x_{1e}^2+3x_1{x}_{1e}\right)-\lambda \left(x_2-x_1\right)\{1+{\alpha }_2[{\left(x_2-x_1\right)}^2+3\left(x_2-x_1\right)\left(x_{2e}-x_{1e}\right)+3{\left(x_{2e}-x_{1e}\right)}^2]\}=0,\\ {\ddot{x}}_2-p^2\left[\frac{a}{L}\mathit{sin}\left(p\tau \right)\mathit{sin}\hspace{0.17em}\phi +\frac{b}{L}\mathit{cos}\left(p\tau \right)\mathit{cos}\hspace{0.17em}\phi \right]-\left(1+x_2\right){\dot{\phi }}^2+w^2\left(1-\mathit{cos}\hspace{0.17em}\phi \right)+c_1{\dot{x}}_2+\left(1+\lambda \right)\times \left(x_2-x_1\right)\left\{1+{\alpha }_2\left[{\left(x_2-x_1\right)}^2+3\left(x_2-x_1\right)\left(x_{2e}-x_{1e}\right)+3{\left(x_{2e}-x_{1e}\right)}^2\right]\right\}=f_1^{*}\hspace{0.17em}\mathit{sin}\left(p_1\tau \right),\\ {\left(1+x_2\right)}^2\ddot{\phi }+2\left(1+x_2\right){\dot{x}}_2\dot{\phi }-p^2\left(1+x_2\right)\left[\frac{a}{L}\mathit{sin}\left(p\tau \right)\mathit{cos}\hspace{0.17em}\phi -\frac{b}{L}\mathit{cos}\left(p\tau \right)\mathit{sin}\hspace{0.17em}\phi \right]+w^2\left(1+x_2\right)\mathit{sin}\hspace{0.17em}\phi +c_2\dot{\phi }=f_0^{*}\hspace{0.17em}\mathit{sin}\left(p_0\tau \right)+\left(1+x_2\right)f_2^{*}\hspace{0.17em}\mathit{sin}\left(p_2\tau \right),x_{1e}+{\alpha }_1{x}_{1e}^3-\lambda \left[x_{2e}-x_{1e}+{\alpha }_2{\left(x_{2e}-x_{1e}\right)}^3\right]=0,\end{array}$

(11)

where $x_{1e}$ and $x_{2e}$ are dimensionless static extensions

$\begin{array}{c}\left(1+\lambda \right)\left[x_{2e}-x_{1e}+{\alpha }_2{\left(x_{2e}-x_{1e}\right)}^3\right]=w^2,\\ \tau =\omega t,x_1=\frac{X_1-X_{1e}}{L},x_{1e}=\frac{X_{1e}}{L},x_2=\frac{X_2-X_{2e}}{L},x_{2e}=\frac{X_{2e}}{L},k_e=\frac{k_1{k}_2}{k_1+k_2},\\ \lambda =\frac{k_2}{k_1},\omega =\sqrt{\frac{k_e}{m}},{\alpha }_1={\Lambda }_1{L}^2,{\alpha }_2={\Lambda }_2{L}^2,w^2=\frac{g}{L{\omega }^2},c_1=\frac{C_1}{m\omega },c_2=\frac{C_2}{m{L}^2\omega },\\ f_0^{*}=\frac{M_0^{*}}{m{L}^2{\omega }^2},f_1^{*}=\frac{F_1^{*}}{mL{\omega }^2},f_2^{*}=\frac{F_2^{*}}{mL{\omega }^2},p=\frac{\Omega }{\omega },p_0=\frac{{\Omega }_0}{\omega },p_1=\frac{{\Omega }_1}{\omega },p_2=\frac{{\Omega }_2}{\omega };\end{array}$

The dimensionless parameters are indicated, and the differentiation with respect to $\tau$ is represented by the dots. In this context, the initial conditions for both the generalized coordinates and their first derivatives are provided below

$\begin{array}{l}{x}_1\left(0\right)=x_{10},x_2\left(0\right)=x_{20},\phi \left(0\right)={\phi }_0,\\ {\dot{x}}_1\left(0\right)=0,{\dot{x}}_2\left(0\right)=0,\dot{\phi }\left(0\right)=0.\end{array}$

Perturbation’s methodology

In this section, the governing system of the EOM outlined in equation (11) for the examined problem is solved analytically in the form of an approximate solution using the MMTS in the time domain. To achieve this assignment, let us approximate the trigonometric functions in $\phi$ till the third order by selecting only the first two terms of their Taylor series around zero as below

$\mathit{sin}\hspace{0.17em}\phi \approx \phi -\frac{{\phi }^3}{6},\mathit{cos}\hspace{0.17em}\phi \approx 1-\frac{{\phi }^2}{2}.$

(12)

Assuming that the system under consideration is only weakly nonlinear, we limit the analysis to tiny values of certain parameters that describe the pendulum and its loading. Assumptions on smallness are expressed as follows:

$\begin{array}{l}{c}_s={\epsilon }^2{\stackrel{\sim }{c}}_s,{\alpha }_s={\epsilon }^2{\stackrel{\sim }{\alpha }}_s,f_s^{*}={\epsilon }^3{\stackrel{\sim }{f}}_s^{*}\left(s=1,2\right),\\ f_0^{*}={\epsilon }^3{\stackrel{\sim }{f}}_0^{*},a={\epsilon }^2\stackrel{\sim }{a},b={\epsilon }^2\stackrel{\sim }{b},\end{array}$

(13)

where $0<\epsilon \ll 1$ is a small parameter. If $\epsilon$ tends to zero, then the coefficients ${\stackrel{\sim }{c}}_s,{\stackrel{\sim }{\alpha }}_s,{\stackrel{\sim }{f}}_s^{*},{\stackrel{\sim }{f}}_0^{*},\stackrel{\sim }{a},$ and $\stackrel{\sim }{b}$ are considered finite.

Consequently, the vibration amplitudes $x_1,x_2,$ and $\phi$ are thought to be on the order of $\epsilon$ . Thus, let’s consider the new variables ${\stackrel{\sim }{x}}_1,{\stackrel{\sim }{x}}_2,$ and $\stackrel{\sim }{\phi }$ as follows:

$x_1=\epsilon {\stackrel{\sim }{x}}_1,x_2=\epsilon {\stackrel{\sim }{x}}_2,\phi =\epsilon \stackrel{\sim }{\phi }.$

(14)

According to the procedure of MMST, we can write the asymptotic expansion of these functions as follows:

$\begin{array}{l}{\stackrel{\sim }{x}}_1\left(\tau ,\epsilon \right)={\displaystyle \sum }_{j=1}^3{\epsilon }^{j-1}{\stackrel{\sim }{x}}_{1j}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),\\ {\stackrel{\sim }{x}}_2\left(\tau ,\epsilon \right)={\displaystyle \sum }_{j=1}^3{\epsilon }^{j-1}{\stackrel{\sim }{x}}_{2j}\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),\\ \stackrel{\sim }{\phi }\left(\tau ,\epsilon \right)={\displaystyle \sum }_{j=1}^3{\epsilon }^{j-1}{\stackrel{\sim }{\phi }}_j\left({\tau }_0,{\tau }_1,{\tau }_2\right)+O\left({\epsilon }^3\right),\end{array}$

(15)

where ${\tau }_n={\epsilon }^n\tau \left(n=0,1,2\right)$ are three variables, that are connected to the time $\tau$ . The upper orders have been left out because of $\epsilon$ 's modest size.

Operators having partial derivatives of the time scales ${\tau }_n$ are used in place of the differential operators related to time $\tau$ . According to the below chain rule, one writes

$\begin{array}{l}\frac{d}{d\tau }=D_0+\epsilon D_1+{\epsilon }^2{D}_2,\\ \frac{d^2}{d{\tau }^2}=D_0^2+2\epsilon D_0{D}_1+{\epsilon }^2\left(D_1^2+2D_0{D}_2\right),\end{array}$

(16)

where $D_n=\left(\frac{\partial }{\partial {\tau }_n}\right);n=0,1,2.$

Owing to the presence of algebraic and differential equations, one can deal with the MMST procedure to provide the desired solutions.

Non-resonant oscillations

Substituting from (12)–(16) into (11) to generate equations involving a small parameter raised to various powers. It should be possible to satisfy all three equations for any value of $\epsilon$ . Upon organizing the equations’ terms by the magnitude of $\epsilon$ , eliminating terms of $O\left({\epsilon }^4\right)$ order or greater, the method of undetermined coefficients enables compliance with the condition. This process yields a collection of nine equations involving unspecified functions ${\stackrel{\sim }{x}}_{1j},{\stackrel{\sim }{x}}_{2j},$ and ${\stackrel{\sim }{\phi }}_j\left(j=1,2,3\right)$ , in which three distinct sets of partial DEs comprise the entire ones. The equations in the first group are those that follow, accompanied by $\epsilon$ .

$\begin{array}{l}\lambda {\stackrel{\sim }{x}}_{21}-\left(1+\lambda \right){\stackrel{\sim }{x}}_{11}=0,\\ D_0^2{\stackrel{\sim }{x}}_{21}+\left(1+\lambda \right)\left({\stackrel{\sim }{x}}_{21}-{\stackrel{\sim }{x}}_{11}\right)=0,\\ D_0^2{\stackrel{\sim }{\phi }}_1+w^2{\stackrel{\sim }{\phi }}_1=0.\end{array}$

(17)

These equations are designated as the first-order approximation equations, while the terms containing ${\epsilon }^2$ lead to the derivation of equations representing the second-order approximation

$\begin{array}{c}\lambda {\stackrel{\sim }{x}}_{22}-\left(1+\lambda \right){\stackrel{\sim }{x}}_{12}=0,\\ D_0^2{\stackrel{\sim }{x}}_{22}+\left(1+\lambda \right)\left({\stackrel{\sim }{x}}_{22}-{\stackrel{\sim }{x}}_{12}\right)=\frac{p^2\stackrel{\sim }{b}}{L}\mathit{cos}\left(p{\tau }_0\right)-\frac{1}{2}{w}^2{\stackrel{\sim }{\phi }}_1^2+{\left(D_0{\stackrel{\sim }{\phi }}_1\right)}^2-2D_0{D}_1{\stackrel{\sim }{x}}_{21},\\ D_0^2{\stackrel{\sim }{\phi }}_2+w^2{\stackrel{\sim }{\phi }}_2=\frac{p^2\stackrel{\sim }{a}}{L}\mathit{sin}\left(p{\tau }_0\right)-w^2{\stackrel{\sim }{\phi }}_1{\stackrel{\sim }{x}}_{21}-2\left(D_0{\stackrel{\sim }{\phi }}_1{D}_0{\stackrel{\sim }{x}}_{21}+D_0{D}_1{\stackrel{\sim }{\phi }}_1+{\stackrel{\sim }{x}}_{21}{D}_0^2{\stackrel{\sim }{\phi }}_1\right).\end{array}$

(18)

The coefficients associated with ${\epsilon }^3$ constitute the equations representing the third-order approximation

$\begin{array}{l}\lambda {\stackrel{\sim }{x}}_{23}-\left(1+\lambda \right){\stackrel{\sim }{x}}_{13}=3x_{1e}^2{\stackrel{\sim }{\alpha }}_1{\stackrel{\sim }{x}}_{11}+3\lambda {\stackrel{\sim }{\alpha }}_2{\left(x_{1e}-x_{2e}\right)}^2\left({\stackrel{\sim }{x}}_{11}-{\stackrel{\sim }{x}}_{21}\right),\\ D_0^2{\stackrel{\sim }{x}}_{23}+\left(1+\lambda \right)\left({\stackrel{\sim }{x}}_{23}-{\stackrel{\sim }{x}}_{13}\right)={\stackrel{\sim }{f}}_1^{*}\hspace{0.17em}\mathit{sin}\left(p_1{\tau }_0\right)+\frac{p^2\stackrel{\sim }{a}}{L}\mathit{sin}\left(p{\tau }_0\right){\stackrel{\sim }{\phi }}_1-w^2{\stackrel{\sim }{\phi }}_1{\stackrel{\sim }{\phi }}_2+3{\stackrel{\sim }{\alpha }}_2\left(1+\lambda \right){\left(x_{1e}-x_{2e}\right)}^2\left({\stackrel{\sim }{x}}_{11}-{\stackrel{\sim }{x}}_{21}\right)-D_1^2{\stackrel{\sim }{x}}_{21}+2\left[D_0{\stackrel{\sim }{\phi }}_1(D_1{\stackrel{\sim }{\phi }}_1+D_0{\stackrel{\sim }{\phi }}_2)-D_0{D}_2{\stackrel{\sim }{x}}_{21}-D_0{D}_1{\stackrel{\sim }{x}}_{22}\right]+{\stackrel{\sim }{x}}_{21}{\left(D_0{\stackrel{\sim }{\phi }}_1\right)}^2-{\stackrel{\sim }{c}}_1{D}_0{\stackrel{\sim }{x}}_{21}.\\ D_0^2{\stackrel{\sim }{\phi }}_3+w^2{\stackrel{\sim }{\phi }}_3={\stackrel{\sim }{f}}_0^{*}\hspace{0.17em}\mathit{sin}\left(p_0{\tau }_0\right)+{\stackrel{\sim }{f}}_2^{*}\hspace{0.17em}\mathit{sin}\left(p_2{\tau }_0\right)-\frac{p^2}{L}\left[\stackrel{\sim }{b}\mathit{cos}\left(p{\tau }_0\right){\stackrel{\sim }{\phi }}_1-\stackrel{\sim }{a}\mathit{sin}\left(p{\tau }_0\right){\stackrel{\sim }{x}}_{21}\right]+w^2\left(\frac{{\stackrel{\sim }{\phi }}_1^3}{6}-{\stackrel{\sim }{\phi }}_2{\stackrel{\sim }{x}}_{21}-{\stackrel{\sim }{\phi }}_1{\stackrel{\sim }{x}}_{22}\right)-{\stackrel{\sim }{c}}_2{D}_0{\stackrel{\sim }{\phi }}_1-D_1^2{\stackrel{\sim }{\phi }}_1-2D_0{\stackrel{\sim }{\phi }}_1\left(D_1{\stackrel{\sim }{x}}_{21}+D_0{\stackrel{\sim }{x}}_{22}\right)-2D_0{\stackrel{\sim }{x}}_{21}\left(D_1{\stackrel{\sim }{\phi }}_1+{\stackrel{\sim }{x}}_{21}{D}_0{\stackrel{\sim }{\phi }}_1+D_0{\stackrel{\sim }{\phi }}_2\right)-2(D_0{D}_2{\stackrel{\sim }{\phi }}_1+D_0{D}_1{\stackrel{\sim }{\phi }}_2)-2{\stackrel{\sim }{x}}_{21}\times \left(D_0^2{\stackrel{\sim }{\phi }}_2+2D_0{D}_1{\stackrel{\sim }{\phi }}_1\right)-D_0^2{\stackrel{\sim }{\phi }}_1\left({\stackrel{\sim }{x}}_{21}^2+2{\stackrel{\sim }{x}}_{22}\right)\end{array}$

(19)

In contrast with systems regulated purely by DEs, every stage of the asymptotic approximation is expressed via a single algebraic equation. The way in which these equations are approximately solved is intricately linked to the linearity of the differential operator that arises in other equations. According to this requirement, ${\stackrel{\sim }{x}}_{1j}$ and ${\stackrel{\sim }{x}}_{2j}\left(j=1,2,3\right)$ are linearly related. Assumption (13) justifies this rough linearization.

One should note that equations (17)–(19) can be solved sequentially. In the process of estimating the solution, we invariably begin with solving the algebraic equation. Thus, the general solutions for the system of equation (17) of the first group take the following form:

$\begin{array}{l}{\stackrel{\sim }{x}}_{11}=\frac{\lambda }{1+\lambda }\left(y_1{e}^{i{\tau }_0}+{\overline{y}}_1{e}^{-i{\tau }_0}\right),\\ {\stackrel{\sim }{x}}_{21}=y_1{e}^{i{\tau }_0}+{\overline{y}}_1{e}^{-i{\tau }_0},\\ {\stackrel{\sim }{\phi }}_1={\overline{y}}_2{e}^{iw{\tau }_0}+{\overline{y}}_2{e}^{-iw{\tau }_0}.\end{array}$

(20)

In this context, $i$ represents the imaginary unit, and $y_1$ together with its complex conjugate ${\overline{y}}_1$ are complex-valued functions that are not yet determined, dependent on the slower time scales ${\tau }_1$ and ${\tau }_2$ . In view of the previous solutions (20) and the system of equation (18), secular terms have been obtained. The elimination of these terms required that

$D_1{y}_s=0;\left(s=1,2\right).$

(21)

As a consequence, the solutions of second order are given by

$\begin{array}{l}{\stackrel{\sim }{x}}_{12}=\frac{\lambda }{1+\lambda }\left[\frac{p^2\stackrel{\sim }{b}}{2L\left(1-p^2\right)}{e}^{ip{\tau }_0}-\frac{3w^2{y}_2^2}{2\left(1-4w^2\right)}{e}^{2iw{\tau }_0}+w^2{y}_2{\overline{y}}_2\right]+c.c.,\\ {\stackrel{\sim }{x}}_{22}=\left[\frac{p^2\stackrel{\sim }{b}}{2L\left(1-p^2\right)}{e}^{ip{\tau }_0}-\frac{3w^2{y}_2^2}{2\left(1-4w^2\right)}{e}^{2iw{\tau }_0}+w^2{y}_2{\overline{y}}_2\right]+c.c.,\\ {\stackrel{\sim }{\phi }}_2=\frac{\left(2-w\right)w{\overline{y}}_1{y}_2}{1-2w}{e}^{i\left(w-1\right){\tau }_0}-\frac{\left(2+w\right)w{y}_1{y}_2}{1+2w}{e}^{i\left(1+w\right){\tau }_0}-\frac{i{p}^2\stackrel{\sim }{a}}{2L\left(w^2-p^2\right)}{e}^{ip{\tau }_0}+c.c.\end{array}$

(22)

In this instance, $c.c.$ signifies the complex conjugate of the earlier terms.

Substituting (20) and (22) into the third group (19) and removing terms that cause secular ones, the solvability conditions for the third-order approximation are obtained as follows:

$\begin{array}{l}2i{D}_2{y}_1+\left\{i{\stackrel{\sim }{c}}_1+\frac{6\left(1-w^2\right)w^2{y}_2{\overline{y}}_2}{1-4w^2}+\frac{3}{1+\lambda }\left[\lambda {{\stackrel{\sim }{\alpha }}_1{x}_{1e}^2+\left(x_{2e}-x_{1e}\right)}^2{\stackrel{\sim }{\alpha }}_2\right]\right\}{y}_1=0,\\ 2iw{D}_2{y}_2+\left[iw{\stackrel{\sim }{c}}_2+\frac{6\left(1-w^2\right)w^2}{1-4w^2}{y}_1{\overline{y}}_1-\frac{\left(1+7w^2-8w^4\right)w^2}{2\left(1-4w^2\right)}{y}_2{\overline{y}}_2\right]y_2=0.\end{array}$

(23)

Leveraging the preceding context, the solutions of third order can be written as follows:

$\begin{array}{l}{\stackrel{\sim }{x}}_{13}=\frac{\lambda }{1+\lambda }\{-\frac{i{\stackrel{\sim }{f}}_1^{*}}{2\left(1-p_1^2\right)}{e}^{i{p}_1{\tau }_0}-\frac{3\left(1+w\right)w}{4\left(1+2w\right)}{y}_1{y}_2^2{e}^{i\left(1+2w\right){\tau }_0}+\frac{3\left(1-w\right)w}{4\left(1-2w\right)}{y}_1{\overline{y}}_2^2{e}^{i\left(1-2w\right){\tau }_0}-\frac{i\left(p+2w\right)p^3\stackrel{\sim }{a}{y}_2}{2L\left(p^2-w^2\right)\left[1-{\left(p+w\right)}^2\right]}{e}^{i\left(p+w\right){\tau }_0}+\frac{i\left(p-2w\right)p^3\stackrel{\sim }{a}{y}_2}{2L\left(p^2-w^2\right)\left[1-{\left(w-p\right)}^2\right]}{e}^{i\left(w-p\right){\tau }_0}\}-\frac{3\lambda }{{\left(1+\lambda \right)}^2}\left[x_{1e}^2\left({\stackrel{\sim }{\alpha }}_1-{\stackrel{\sim }{\alpha }}_2\right)+x_{2e}{\stackrel{\sim }{\alpha }}_2\left(2x_{1e}-x_{2e}\right)\right]y_1{e}^{i{\tau }_0}+c.c.,\\ {\stackrel{\sim }{x}}_{23}=-\frac{i{\stackrel{\sim }{f}}_1^{*}}{2\left(1-p_1^2\right)}{e}^{i{p}_1{\tau }_0}-\frac{3\left(1+w\right)w}{4\left(1+2w\right)}{y}_1{y}_2^2{e}^{i\left(1+2w\right){\tau }_0}+\frac{3\left(1-w\right)w}{4\left(1-2w\right)}{y}_1{\overline{y}}_2^2{e}^{i\left(1-2w\right){\tau }_0}-\frac{i\left(p+2w\right)p^3\stackrel{\sim }{a}{y}_2}{2L\left(p^2-w^2\right)\left[1-{\left(p+w\right)}^2\right]}{e}^{i\left(p+w\right){\tau }_0}+\frac{i\left(p-2w\right)p^3\stackrel{\sim }{a}{y}_2}{2L\left(p^2-w^2\right)\left[1-{\left(w-p\right)}^2\right]}{e}^{i\left(w-p\right){\tau }_0}+c.c,\\ {\stackrel{\sim }{\phi }}_3=-\frac{i{\stackrel{\sim }{f}}_0^{*}}{2\left(w^2-{p2}_0\right)}{e}^{i{p}_0{\tau }_0}-\frac{i{\stackrel{\sim }{f}}_2^{*}}{2\left(w^2-{p2}_2\right)}{e}^{i{p}_2{\tau }_0}+\frac{i\left(2+p\right)p^3\stackrel{\sim }{a}{y}_1}{2L\left(p^2-w^2\right)\left[w^2-{\left(1+p\right)}^2\right]}{e}^{i\left(1+p\right){\tau }_0}+\frac{i\left(2-p\right)p^3\stackrel{\sim }{a}{y}_1}{2L\left(p^2-w^2\right)\left[w^2-{\left(1-p\right)}^2\right]}{e}^{i\left(1-p\right){\tau }_0}+\frac{\left(1-p-w\right)\left(1+p+w\right)p\stackrel{\sim }{b}{y}_2}{2L\left(1-p^2\right)\left(p+2w\right)}{e}^{i\left(p+w\right){\tau }_0}+\frac{\left(1-{\left(p-w\right)}^2\right)p\stackrel{\sim }{b}{y}_2}{2L\left(1-p^2\right)\left(p-2w\right)}{e}^{i\left(w-p\right){\tau }_0}+\frac{\left(2+w\right)\left(3+w\right)}{4+8w}{y}_1^2{y}_2{e}^{i\left(2+w\right){\tau }_0}-\frac{\left(1-49w^2\right)}{48\left(1-4w^2\right)}{y}_2^3{e}^{3iw{\tau }_0}-\frac{\left(3-w\right)\left(2-w\right)w}{4\left(1-2w\right)}{y}_1^2{\overline{y}}_2{e}^{i\left(2-w\right){\tau }_0}+c.c.\end{array}$

(24)

By leveraging the assumptions outlined in hypothesis (14) and the solutions provided in series (15), along with the attained solutions in equations (20), (22), and (24), we can promptly derive the required approximate expressions of $x_1,x_2,$ and $\phi$ up to the third approximation.

Taking into account the removal conditions (21) and (23) targeting secular terms, the functions $y_s\left(s=1,2\right)$ can be estimated. From condition (21), we can deduce that the functions $y_s$ do not have any reliance on the variable ${\tau }_1$ . Bearing this circumstance in mind, we can display these functions in the next exponential form

$y_s\left({\tau }_2\right)=\frac{1}{2}{\stackrel{\sim }{y}}_s\left({\tau }_2\right)e^{i{\theta }_s\left({\tau }_2\right)};\left(y_s=\epsilon {\stackrel{\sim }{y}}_s\right).$

(25)

Both functions ${\stackrel{\sim }{y}}_s\left({\tau }_2\right)$ and ${\theta }_s\left({\tau }_2\right)$ assigned to $s$ be real-valued.

Utilizing the analysis presented above, the first-order derivatives of the functions $d_s$ can be summarized as follows:

$\frac{d{y}_s}{d\tau }={\epsilon }^2{D}_2{y}_s.$

(26)

Thus, the substitution into equations (25) and (26) enables the conversion of the partial DEs (21) and (23) into ordinary DEs. By distinguishing the real components from the imaginary ones, one obtains the below system of four first-order ordinary DEs concerning $y_s$ and ${\theta }_s$ .

$\begin{array}{l}{\dot{y}}_1=-\frac{1}{2}{c}_1{y}_1,\\ y_1{\dot{\theta }}_1=\left[\frac{3\left(1-w^2\right)w^2}{4\left(1-4w^2\right)}{y}_2^2+\frac{3}{2\left(1+\lambda \right)}\left(\lambda {\alpha }_1{x}_{1e}^2+{\left(x_{2e}-x_{1e}\right)}^2{\alpha }_2\right)\right]y_1,\\ {\dot{y}}_2=-\frac{1}{2}{c}_2{y}_2,\\ y_2{\dot{\theta }}_2=-\left[\frac{\left(1+7w^2-8w^4\right)w}{16\left(1-4w^2\right)}{y}_2^2-\frac{3\left(1-w^2\right)w}{4\left(1-4w^2\right)}{y}_1^2\right]y_2.\end{array}$

(27)

This system can be solved numerically for both $y_s$ and ${\theta }_s$ using the algorithms of Runge–Kutta in light of Mathematica software 12.3. Therefore, the next data are used

$\begin{array}{l}g=9.8m/s^2,a=0.7m,b=0.4m,l_0=0.9m,m=1kg,f_0^{*}=0.0004,\\ f_1^{*}=0.0001,f_2^{*}=0.008,k_1=30,{\Lambda }_1=0.1,{\Lambda }_2=0.1,\Omega =0.1,{\Omega }_0=0.01,\\ {\Omega }_1=0.0001,{\Omega }_2=0.1.\end{array}$

Distinct plots visually portray the solutions of the aforementioned system for the chosen values of the physical parameters, as illustrated in Figures 3–8. These figures are evaluated when the values of both damping parameters are as follows:

$c_1\left(=0.01,0.015,0.02\right),$

Figure 3.

(a)–(d) shows the time variation of $y_s\left(s=1,2\right)$ and the phases ${\theta }_s$ at $c_1\left(=0.01,0.015,0.02\right),$ $c_2=0.01$ , and $\omega =4.242$ , while (e) and (f) present their projections in the planes ${\theta }_s{y}_s$ .

Figure 4.

(a)–(d) shows the time variation of $y_s\left(s=1,2\right)$ and the phases ${\theta }_s$ at $c_2\left(=0.01,0.02,0.03\right),$ $c_1=0.01$ , and $\omega =4.242$ , while (e) and (f) present their projections in the planes ${\theta }_s{y}_s$ .

Figure 5.

(a)–(d) presents the behavior of the functions $y_s\left(s=1,2\right)$ and the phases ${\theta }_s$ at $\omega \left(=4.472,4.242,3.872\right),$ and $c_1=c_2=0.01$ , while (e) and (f) illustrate their projections in the planes ${\theta }_s{y}_s$ .

Figure 6.

Solutions’ time histories when ${\mathit{c}}_{\mathbf{1}}\left(=\mathbf{0.01},\mathbf{0.015},\mathbf{0.02}\right),$ ${\mathit{c}}_{\mathbf{2}}=\mathbf{0.01}$ , and $\mathit{\omega }=\mathbf{4.242}$ : (a) and (b) at $\tau \in \left[0,200\right]$ ; (c)–(e) at $\tau \in \left[0,2500\right]$ .

Figure 7.

Solutions’ time histories when ${\mathit{c}}_{\mathbf{2}}\left(=\mathbf{0.01},\mathbf{0.02},\mathbf{0.03}\right),$ ${\mathit{c}}_{\mathbf{1}}=\mathbf{0.01}$ , and $\mathit{\omega }=\mathbf{4.242}$ : (a) and (b) at $\tau \in \left[0,200\right]$ ; (c)–(e) at $\tau \in \left[0,2500\right]$ .

Figure 8.

Solutions’ time histories when $\mathit{\omega }\left(=\mathbf{4.472},\mathbf{4.242},\mathbf{3.872}\right)$ and ${\mathit{c}}_{\mathbf{1}}={\mathit{c}}_{\mathbf{2}}=\mathbf{0.01}$ : (a) and (b) at $\tau \in \left[0,200\right]$ ; (c)–(e) at $\tau \in \left[0,2500\right]$ .

$c_2\left(=0.01,0.02,0.03\right),$ and the frequency $\omega \left(=4.472,4.242,3.872\right)$ is included. It is apparent that the plotted curves of the amplitude $y_1$ and the phase plane ${\theta }_2$ are positively impacted by the variation of $c_1$ values, and they have decay forms, as explored in Figure 3(a) and (d), respectively. The underlying reason is that their equations are based on the values of $c_1$ . On the other hand, the temporal curves in the planes $\tau y_2$ and $\tau {\theta }_1$ are not impacted by the change in the values of $c_1$ , as indicated in Figure 3(b) and (c), because the first-order equations of $y_2$ and ${\theta }_1$ don’t depend on the parameter $c_1$ . The projections of the included curves in parts (a), (b), (c), and (d) of Figure 3 on the phase planes ${\theta }_1{y}_1$ and ${\theta }_2{y}_2$ are graphed in Figure 3(e) and (f).

Figure 4 examines the influence of the change in amplitudes and phase angles $y_s$ and ${\theta }_s$ , respectively, when the values of $c_2$ have different values, as well as the study of their phase planes. It is obvious that $y_1$ and ${\theta }_1$ don’t impact with the change in $c_2$ values, while $y_2$ and ${\theta }_2$ vary to have decay forms over time. Furthermore, when the dimensionless time parameter is excluded, the projections of $y_s$ and ${\theta }_s$ in the plane ${\theta }_s{y}_s$ are graphed.

The plotted curves in Figure 5 demonstrate that they are drawn for various values of the frequency $\omega \left(=4.472,4.242,3.872\right)$ . Decay behaviors are noted in Figure 5(a), (b), and (d) when time goes on. Furthermore, the amplitudes $y_s$ don’t impact by the change in $\omega$ values, while ${\theta }_s$ are influenced by the same values, as seen in Figure 5(a)–(d). The curves in Figure 5(e) and (f) show the curves of the phases via amplitudes in the planes ${\theta }_s{y}_s$ .

Now, we will display the identified solutions $x_1,x_2,$ and $\phi$ when $c_1,c_2,$ and $\omega$ have various values, as in Figures 6–8. Their portions (a) and (b) show the solutions of $x_1$ and $x_2$ in short-range, that is, $\tau \in \left[0,200\right]$ , while portions (c) and (d) presented these solutions for long ranges, that is, $\tau \in \left[0,2500\right]$ . These solutions have decay forms through the time interval $\left[0,200\right]$ and then will be periodic till the end of the long-range interval. The analysis of the drawn curves in these figures expresses that the solutions $x_1$ and $x_2$ are influenced by the change in $c_1$ and $\omega$ , as indicated in Figures 6 and 8. On the other hand, the achieved solutions $x_1,x_2,$ and $\phi$ don’t impact with the change in $c_2$ values, as drawn in Figure 7. This is attributed to the mathematical formulations used for the first-order equations in system (26). Whereas, the periodic solution $\phi$ is influenced by the change in $\omega$ only, as explored in Figure 8(c), where the amplitudes of the presented waves of this solution increase with the increase in $\omega$ values.

It must be mentioned that periodic solutions imply that the system’s behavior repeats over time, making future states predictable. This predictability is crucial for planning and control in various applications. These solutions can indicate stability in a dynamic system, meaning small disturbances or changes in initial conditions do not lead to unpredictable behavior. This stability is critical in control systems, ensuring that the system remains in a desirable state. Moreover, it can be represented as optimal operating points for dynamic systems, allowing for efficient and effective performance. For example, in power systems, periodic solutions can help optimize load balancing and energy distribution. These solutions can serve as a benchmark for understanding more complex, non-periodic behaviors, providing a stepping stone for advanced research and development.

Classification of resonance cases

In this part, we will examine how to categorize the different cases of resonance using the solutions mentioned in (22) and (24), which are applicable as long as their denominators are non-zero. Accordingly, these cases are divided into primary external resonance at $p_1=1,p_2=w,p_0=w,$ $p=1,p=w,p=\pm 2w$ ; internal resonance at $w=\pm \left(\frac{1}{2}\right),w=-\left(\frac{1}{4}\right)$ ; and combined resonance at ${\left(p\pm w\right)}^2=1,$ ${\left(1\pm p\right)}^2=w^2$ .

Our focus will be on two primary external resonances. Thus, we include the occurrence of both cases simultaneously.

$p_1\approx 1,p_0\approx w.$

(28)

Relation (28) reveals the closeness of $p_1$ to $1,$ and $p_0$ to $w$ , respectively. To accomplish this, one can introduce the dimensionless values defined by the detuning parameters ${\sigma }_s\left(s=1,2\right)$ (indicating the gap between the oscillations and the precise resonance) as follows:

$\begin{array}{l}{p}_1=1+{\sigma }_1,\\ p_0=w+{\sigma }_2.\end{array}$

(29)

Under the assumption that the parameters ${\sigma }_s$ are of order $O\left({\epsilon }^2\right)$ , we can write

${\sigma }_1={\epsilon }^2{\stackrel{\sim }{\sigma }}_1,{\sigma }_2={\epsilon }^2{\stackrel{\sim }{\sigma }}_2,$

(30)

where the coefficients ${\stackrel{\sim }{\sigma }}_s$ are finite as $\epsilon$ tends to zero.

Substituting (29) and (30) into (11) and using (13)–(16), and then equalling the coefficients of the $\epsilon$ , we obtain three groups of equations according to the various powers of $\epsilon$ . Two of them are the same of the systems of equations (17) and (18), while the latter group has the below form

$\begin{array}{l}\lambda {\stackrel{\sim }{x}}_{23}-\left(1+\lambda \right){\stackrel{\sim }{x}}_{13}=3x_{1e}^2{\stackrel{\sim }{\alpha }}_1{\stackrel{\sim }{x}}_{11}+3\lambda {\stackrel{\sim }{\alpha }}_2{\left(x_{1e}-x_{2e}\right)}^2\left({\stackrel{\sim }{x}}_{11}-{\stackrel{\sim }{x}}_{21}\right),\\ D_0^2{\stackrel{\sim }{x}}_{23}+\left(1+\lambda \right)\left({\stackrel{\sim }{x}}_{23}-{\stackrel{\sim }{x}}_{13}\right)={\stackrel{\sim }{f}}_1^{*}\hspace{0.17em}\mathit{sin}\left[\left(1+{\sigma }_1\right){\tau }_0\right]+\frac{p^2\stackrel{\sim }{a}}{L}\mathit{sin}\left(p{\tau }_0\right){\stackrel{\sim }{\phi }}_1-w^2{\stackrel{\sim }{\phi }}_1{\stackrel{\sim }{\phi }}_2+3{\stackrel{\sim }{\alpha }}_2\left(1+\lambda \right){\left(x_{1e}-x_{2e}\right)}^2\left({\stackrel{\sim }{x}}_{11}-{\stackrel{\sim }{x}}_{21}\right)-D_1^2{\stackrel{\sim }{x}}_{21}+2\left[D_0{\stackrel{\sim }{\phi }}_1(D_1{\stackrel{\sim }{\phi }}_1+D_0{\stackrel{\sim }{\phi }}_2)-D_0{D}_2{\stackrel{\sim }{x}}_{21}-D_0{D}_1{\stackrel{\sim }{x}}_{22}\right]+{\stackrel{\sim }{x}}_{21}{\left(D_0{\stackrel{\sim }{\phi }}_1\right)}^2-{\stackrel{\sim }{c}}_1{D}_0{\stackrel{\sim }{x}}_{21},\\ D_0^2{\stackrel{\sim }{\phi }}_3+w^2{\stackrel{\sim }{\phi }}_3={\stackrel{\sim }{f}}_0^{*}\hspace{0.17em}\mathit{sin}\left[\left(w+{\sigma }_2\right){\tau }_0\right]+{\stackrel{\sim }{f}}_2^{*}\hspace{0.17em}\mathit{sin}\left(p_2{\tau }_0\right)-\frac{p^2}{L}\left[\stackrel{\sim }{b}\mathit{cos}\left(p{\tau }_0\right){\stackrel{\sim }{\phi }}_1-\stackrel{\sim }{a}\mathit{sin}\left(p{\tau }_0\right){\stackrel{\sim }{x}}_{21}\right]+w^2\left(\frac{{\stackrel{\sim }{\phi }}_1^3}{6}-{\stackrel{\sim }{\phi }}_2{\stackrel{\sim }{x}}_{21}-{\stackrel{\sim }{\phi }}_1{\stackrel{\sim }{x}}_{22}\right)-{\stackrel{\sim }{c}}_2{D}_0{\stackrel{\sim }{\phi }}_1-D_1^2{\stackrel{\sim }{\phi }}_1-2D_0{\stackrel{\sim }{\phi }}_1\left(D_1{\stackrel{\sim }{x}}_{21}+D_0{\stackrel{\sim }{x}}_{22}\right)-2D_0{\stackrel{\sim }{x}}_{21}\left(D_1{\stackrel{\sim }{\phi }}_1+{\stackrel{\sim }{x}}_{21}{D}_0{\stackrel{\sim }{\phi }}_1+D_0{\stackrel{\sim }{\phi }}_2\right)-2(D_0{D}_2{\stackrel{\sim }{\phi }}_1+D_0{D}_1{\stackrel{\sim }{\phi }}_2]-2{\stackrel{\sim }{x}}_{21}\times \left(D_0^2{\stackrel{\sim }{\phi }}_2+2D_0{D}_1{\stackrel{\sim }{\phi }}_1\right)-D_0^2{\stackrel{\sim }{\phi }}_1\left({\stackrel{\sim }{x}}_{21}^2+2{\stackrel{\sim }{x}}_{22}\right).\end{array}$

(31)

Based on the solutions (20) and (22), one can obtain the solutions of the third group of approximation (31) after removing the yielded secular terms in the form

$\begin{array}{l}{\stackrel{\sim }{\mathit{x}}}_{\mathbf{13}}=\frac{\mathit{\lambda }}{\mathbf{1}+\mathit{\lambda }}\left\{\frac{\mathbf{3}\left(\mathbf{1}-\mathit{w}\right)\mathit{w}}{\mathbf{4}\left(\mathbf{1}-\mathbf{2}\mathit{w}\right)}{\mathit{y}}_{\mathbf{1}}{\overline{\mathit{y}}}_{\mathbf{2}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathbf{1}-\mathbf{2}\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}-\frac{\mathbf{3}\left(\mathbf{1}+\mathit{w}\right)\mathit{w}}{\mathbf{4}\left(\mathbf{1}+\mathbf{2}\mathit{w}\right)}{\mathit{y}}_{\mathbf{1}}{\mathit{y}}_{\mathbf{2}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathbf{1}+\mathbf{2}\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}-\frac{\mathit{i}\left(\mathit{p}+\mathbf{2}\mathit{w}\right){\mathit{p}}^{\mathbf{3}}\stackrel{\sim }{\mathit{a}}{\mathit{y}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathit{p}+\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\mathit{L}\left({\mathit{p}}^{\mathbf{2}}-{\mathit{w}}^{\mathbf{2}}\right)\left[\mathbf{1}-{\left(\mathit{p}+\mathit{w}\right)}^{\mathbf{2}}\right]}+\frac{\mathit{i}\left(\mathit{p}-\mathbf{2}\mathit{w}\right){\mathit{p}}^{\mathbf{3}}\stackrel{\sim }{\mathit{a}}{\mathit{y}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathit{w}-\mathit{p}\right){\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\mathit{L}\left({\mathit{p}}^{\mathbf{2}}-{\mathit{w}}^{\mathbf{2}}\right)\left[\mathbf{1}-{\left(\mathit{w}-\mathit{p}\right)}^{\mathbf{2}}\right]}\right\}-\frac{\mathbf{3}\mathit{\lambda }}{{\left(\mathbf{1}+\mathit{\lambda }\right)}^{\mathbf{2}}}\left[{\mathit{x}}_{\mathbf{1}\mathit{e}}^{\mathbf{2}}\left({\stackrel{\sim }{\mathit{\alpha }}}_{\mathbf{1}}-{\stackrel{\sim }{\mathit{\alpha }}}_{\mathbf{2}}\right)+{\mathit{x}}_{\mathbf{2}\mathit{e}}{\stackrel{\sim }{\mathit{\alpha }}}_{\mathbf{2}}\left(\mathbf{2}{\mathit{x}}_{\mathbf{1}\mathit{e}}-{\mathit{x}}_{\mathbf{2}\mathit{e}}\right)\right]{\mathit{y}}_{\mathbf{1}}{\mathit{e}}^{\mathit{i}{\mathit{\tau }}_{\mathbf{0}}}+\mathit{c}.\mathit{c}.,\\ {\stackrel{\sim }{\mathit{x}}}_{\mathbf{23}}=\frac{\mathbf{3}\left(\mathbf{1}-\mathit{w}\right)\mathit{w}}{\mathbf{4}\left(\mathbf{1}-\mathbf{2}\mathit{w}\right)}{\mathit{y}}_{\mathbf{1}}{\overline{\mathit{y}}}_{\mathbf{2}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathbf{1}-\mathbf{2}\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}-\frac{\mathbf{3}\left(\mathbf{1}+\mathit{w}\right)\mathit{w}}{\mathbf{4}\left(\mathbf{1}+\mathbf{2}\mathit{w}\right)}{\mathit{y}}_{\mathbf{1}}{\mathit{y}}_{\mathbf{2}}^{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathbf{1}+\mathbf{2}\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}-\frac{\mathit{i}\left(\mathit{p}+\mathbf{2}\mathit{w}\right){\mathit{p}}^{\mathbf{3}}\stackrel{\sim }{\mathit{a}}{\mathit{y}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathit{p}+\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\mathit{L}\left({\mathit{p}}^{\mathbf{2}}-{\mathit{w}}^{\mathbf{2}}\right)\left[\mathbf{1}-{\left(\mathit{p}+\mathit{w}\right)}^{\mathbf{2}}\right]}+\frac{\mathit{i}\left(\mathit{p}-\mathbf{2}\mathit{w}\right){\mathit{p}}^{\mathbf{3}}\stackrel{\sim }{\mathit{a}}{\mathit{y}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathit{w}-\mathit{p}\right){\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\mathit{L}\left({\mathit{p}}^{\mathbf{2}}-{\mathit{w}}^{\mathbf{2}}\right)\left[\mathbf{1}-{\left(\mathit{w}-\mathit{p}\right)}^{\mathbf{2}}\right]}+\mathit{c}.\mathit{c},\\ {\stackrel{\sim }{\mathit{\phi }}}_{\mathbf{3}}=\frac{\mathit{i}\left(\mathbf{2}+\mathit{p}\right){\mathit{p}}^{\mathbf{3}}\stackrel{\sim }{\mathit{a}}{\mathit{y}}_{\mathbf{1}}{\mathit{e}}^{\mathit{i}\left(\mathbf{1}+\mathit{p}\right){\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\mathit{L}\left({\mathit{p}}^{\mathbf{2}}-{\mathit{w}}^{\mathbf{2}}\right)\left[{\mathit{w}}^{\mathbf{2}}-{\left(\mathbf{1}+\mathit{p}\right)}^{\mathbf{2}}\right]}-\frac{\mathit{i}{\stackrel{\sim }{\mathit{f}}}_{\mathbf{2}}^{*}{\mathit{e}}^{\mathit{i}{\mathit{p}}_{\mathbf{2}}{\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\left({\mathit{w}}^{\mathbf{2}}-{\mathit{p}\mathbf{2}}_{\mathbf{2}}\right)}+\frac{\mathit{i}\left(\mathbf{2}-\mathit{p}\right){\mathit{p}}^{\mathbf{3}}\stackrel{\sim }{\mathit{a}}{\mathit{y}}_{\mathbf{1}}{\mathit{e}}^{\mathit{i}\left(\mathbf{1}-\mathit{p}\right){\mathit{\tau }}_{\mathbf{0}}}}{\mathbf{2}\mathit{L}\left({\mathit{p}}^{\mathbf{2}}-{\mathit{w}}^{\mathbf{2}}\right)\left[{\mathit{w}}^{\mathbf{2}}-{\left(\mathbf{1}-\mathit{p}\right)}^{\mathbf{2}}\right]}+\frac{\left(\mathbf{1}-\mathit{p}-\mathit{w}\right)\left(\mathbf{1}+\mathit{p}+\mathit{w}\right)\mathit{p}\stackrel{\sim }{\mathit{b}}{\mathit{y}}_{\mathbf{2}}}{\mathbf{2}\mathit{L}\left(\mathbf{1}-{\mathit{p}}^{\mathbf{2}}\right)\left(\mathit{p}+\mathbf{2}\mathit{w}\right)}{\mathit{e}}^{\mathit{i}\left(\mathit{p}+\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}+\frac{\left[\mathbf{1}-{\left(\mathit{p}-\mathit{w}\right)}^{\mathbf{2}}\right]\mathit{p}\stackrel{\sim }{\mathit{b}}{\mathit{y}}_{\mathbf{2}}}{\mathbf{2}\mathit{L}\left(\mathbf{1}-{\mathit{p}}^{\mathbf{2}}\right)\left(\mathit{p}-\mathbf{2}\mathit{w}\right)}{\mathit{e}}^{\mathit{i}\left(\mathit{w}-\mathit{p}\right){\mathit{\tau }}_{\mathbf{0}}}+\frac{\left(\mathbf{2}+\mathit{w}\right)\left(\mathbf{3}+\mathit{w}\right)}{\mathbf{4}+\mathbf{8}\mathit{w}}\times {\mathit{y}}_{\mathbf{1}}^{\mathbf{2}}{\mathit{y}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathbf{2}+\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}-\frac{\left(\mathbf{1}-\mathbf{49}{\mathit{w}}^{\mathbf{2}}\right)}{\mathbf{48}\left(\mathbf{1}-\mathbf{4}{\mathit{w}}^{\mathbf{2}}\right)}{\mathit{y}}_{\mathbf{2}}^{\mathbf{3}}{\mathit{e}}^{\mathbf{3}\mathit{i}\mathit{w}{\mathit{\tau }}_{\mathbf{0}}}-\frac{\left(\mathbf{3}-\mathit{w}\right)\left(\mathbf{2}-\mathit{w}\right)\mathit{w}}{\mathbf{4}\left(\mathbf{1}-\mathbf{2}\mathit{w}\right)}{\mathit{y}}_{\mathbf{1}}^{\mathbf{2}}{\overline{\mathit{y}}}_{\mathbf{2}}{\mathit{e}}^{\mathit{i}\left(\mathbf{2}-\mathit{w}\right){\mathit{\tau }}_{\mathbf{0}}}+\mathit{c}.\mathit{c}.\end{array}$

(32)

The eliminating conditions of secular terms are given below

$\begin{array}{l}{D}_1{y}_1=0,\\ D_1{y}_2=0,\\ 2i{D}_2{y}_1=-\frac{i}{2}{\stackrel{\sim }{f}}_1^{*}{e}^{i{\epsilon }^2{\stackrel{\sim }{\sigma }}_1{\tau }_0}-\left[i{\stackrel{\sim }{c}}_1+\frac{6\left(1-w^2\right)w^2{y}_2{\overline{y}}_2}{1-4w^2}+\frac{3}{1+\lambda }\left(\lambda {\stackrel{\sim }{\alpha }}_1{x}_{1e}^2+{\left(x_{2e}-x_{1e}\right)}^2{\stackrel{\sim }{\alpha }}_2\right)\right]y_1,\\ 2iw{D}_2{y}_2=-\frac{i}{2}{\stackrel{\sim }{f}}_0^{*}{e}^{i{\epsilon }^2{\stackrel{\sim }{\sigma }}_2{\tau }_0}-\left[iw{\stackrel{\sim }{c}}_2+\frac{6\left(1-w^2\right)w^2}{1-4w^2}{y}_1{\overline{y}}_1-\frac{\left(1+7w^2-8w^4\right)w^2}{2\left(1-4w^2\right)}{y}_2{\overline{y}}_2\right]y_2=0.\end{array}$

(33)

An in-depth analysis of the aforementioned solvability conditions indicates the presence of a system of four nonlinear partial DEs related to the functions $y_s\left(s=1,2\right)$ that vary with the slow scales ${\tau }_2$ . After that, we may present these functions’ polar forms and their first derivatives as in (25) and (26), respectively.

Hence, by applying equation (25) and the subsequent modified phases, we can convert partial DEs (33) into ordinary ones

$\begin{array}{l}{\psi }_1\left({\tau }_2\right)={\stackrel{\sim }{\sigma }}_1{\tau }_2-{\theta }_1\left({\tau }_2\right),\\ {\psi }_2\left({\tau }_2\right)={\stackrel{\sim }{\sigma }}_2{\tau }_2-{\theta }_2\left({\tau }_2\right),\end{array}$

(34)

when the imaginary and real parts are separated, the resulting system consists of four first-order ODEs involving $y_s$ and ${\psi }_s$ . These equations are known by the ME of the examined model when two resonance scenarios are taken into account

$\begin{array}{l}{\dot{y}}_1=-\frac{1}{2}\left(f_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\psi }_1+c_1{y}_1\right),\\ y_1{\dot{\psi }}_1=\frac{1}{2}{f}_1^{*}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\psi }_1+\left\{{\sigma }_1-\left[\frac{3\left(1-w^2\right)w^2}{4\left(1-4w^2\right)}{y}_2^2+\frac{3}{2\left(1+\lambda \right)}\left(\lambda {\alpha }_1{x}_{1e}^2+{\left(x_{2e}-x_{1e}\right)}^2{\alpha }_2\right)\right]\right\}{y}_1,\\ w{\dot{y}}_2=-\frac{1}{2}\left(f_0^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\psi }_2+w{c}_2{y}_2\right),\\ w{y}_2{\dot{\psi }}_2=\frac{1}{2}{f}_0^{*}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\psi }_2+w\left\{{\sigma }_2+\left[\frac{\left(1+7w^2-8w^4\right)w}{16\left(1-4w^2\right)}{y}_2^2-\frac{3\left(1-w^2\right)w}{4\left(1-4w^2\right)}{y}_1^2\right]\right\}{y}_2.\end{array}$

(35)

According to the system of equation (35), one can solve it numerically using the Mathematica 12.3 to obtain the solutions $y_s$ and ${\psi }_s$ when the following data are considered:

$\begin{array}{l}g=9.8m/s^2,a=0.7m,b=0.4m,l_0=0.9m,m=1kg,f_0^{*}=0.0004,\\ f_1^{*}=0.0001,f_2^{*}=0.008,k_1=30,{\Lambda }_1=0.1,{\Lambda }_2=0.1,\Omega =0.1,{\Omega }_0=0.01,\\ {\Omega }_1=0.0001,{\Omega }_2=0.1,{\sigma }_1=0.01,{\sigma }_2=0.01.\end{array}$

These solutions are graphed beside the phase plane diagrams for the projections of these solutions in the planes $y_s{\psi }_s$ , as in Figures 9–11. These graphs are calculated when the $c_1,c_2,$ and $\omega$ have the values $\left(0.01,0.015,0.02\right),$ $\left(0.01,0.02,0.03\right),$ and $\left(4.472,4.242,3.872\right)$ , respectively. It is noted that the temporal history of the amplitude $y_1$ has been impacted by the variation of $c_1$ values, while $y_2,{\psi }_1,$ and ${\psi }_2$ don’t have any change with the same variation, as seen in Figure 9(a)-(d). It is observed that the wave sequence ranges between decreasing waves, which are represented by straight lines, as in Figure 9(a), and between ascending waves for a certain period of time, then they begin to decrease, and then they are stable until the end of the time period, as in Figure 9(b) and (d). Moreover, we also observe that one of the waves tends to decrease sharply until it reaches near the horizontal axis, which represents time, and the wave then stabilizes, as in Figure 9(c). When the time $\tau$ is excluded from the amplitudes $y_s$ and the related phases ${\psi }_s$ , one obtains a projection of them in the planes $y_s{\psi }_s$ , as seen in Figure 9(e) and (f).

Figure 9.

(a)–(d) shows time change of the functions $y_s\left(\tau \right)\left(s=1,2\right)$ and ${\psi }_s\left(\tau \right)$ while (e) and (f) show the curves in the planes ${\psi }_s{y}_s$ at $c_1\left(=0.01,0.015,0.02\right)$ when $c_2=0.01$ , and $\omega =4.242$ .

Figure 10.

(a)–(d) reveals the behavior of $y_s\left(\tau \right)$ and ${\psi }_s\left(\tau \right)$ , while (e) and (f) reveal the curves in the planes ${\psi }_s{y}_s$ when $c_2\left(=0.01,0.02,0.03\right)$ at $c_1=0.01$ and $\omega =4.242$ .

Figure 11.

(a)–(d) portrays the change of the time-dependent functions $y_s\left(\tau \right)$ and ${\psi }_s\left(\tau \right)$ , while (e) and (f) portray the curves in the planes ${\psi }_s{y}_s$ when $\omega \left(=4.472,4.242,3.872\right)$ at $c_1=c_2=0.01$ .

With an insight into the curves drawn in Figure 10, we find that these curves, which represent the amplitude $y_2$ and the phase ${\psi }_2$ , were affected by the change in the values of $c_2$ , as in Figure 10(b) and (d), respectively. Therefore, their phase plane curves undergo a similar change as in Figure 10(f). On the other hand, we find that the curves representing $y_1$ and ${\psi }_1$ did not have any significant effect on changing the values of $c_2$ , as in Figure 10(a) and (c). Consequently, the corresponding phase plane plots for them do not change with the variation of $c_2$ . The behavior of the drawn curves in portions of Figure 10 is considered very close to the behavior of the plotted ones in Figure 9. This can be traced back to the mathematical representation found in the system of equation (35).

By examining the curves in various parts of Figure 11, we observe that the previously mentioned frequency values significantly impact both $y_2$ and ${\psi }_1$ , as shown in Figure 11(b) and (c). Conversely, the effect on $y_1$ is minimal or almost imperceptible, as indicated in Figure 11(a). Additionally, there is no observable effect of changing frequency values on the behavior of the wave representing ${\psi }_2$ , as noticed in Figure 11(d). Consequently, the projections of these curves in the ${\psi }_s{y}_s$ planes are influenced by the changing frequency values, as illustrated in Figure 11(e) and (f). The behavior of the drawn waves is consistent with the behavior observed in the previous Figures 9 and 10.

Taking a closer look at the curves in Figure 12, one can observe that the solutions $x_1,x_2,$ and $\phi$ exhibit a wavelike nature, with all their behaviors consisting of wave packets. The shape of these packets varies depending on the specific solution they represent. These waves are calculated when $c_2=0.01$ and $\omega =4.242$ at aforementioned various values of $c_1$ . Furthermore, changing the values of the damping coefficient $c_1$ does not affect the drawn waves, which is expected since these solutions do not explicitly depend on $c_1$ .

Figure 12.

Solutions’ time histories $x_1,x_2$ , and $\phi$ in parts (a), (b), and (c), respectively, when ${\mathit{c}}_{\mathbf{1}}=\left(\mathbf{0.01},\mathbf{0.015},\mathbf{0.02}\right),$ ${\mathit{c}}_{\mathbf{2}}=\mathbf{0.01}$ , and $\mathit{\omega }=\mathbf{4.242}$ .

In the same context, parts of Figure 13 are plotted at ${\mathit{c}}_{\mathbf{1}}=\mathbf{0.01}$ and $\mathit{\omega }=\mathbf{4.242}$ while ${\mathit{c}}_{\mathbf{2}}$ varying as previously mentioned in Figure 10. A noticeable change is observed in the waves representing the solution $\mathit{\phi }$ as a result of changing the values of ${\mathit{c}}_{\mathbf{2}}$ , as indicated in Fig. (13(c). Moreover, the amplitudes of the waves decrease as ${\mathit{c}}_{\mathbf{2}}$ values increase, while the wavelength remains constant. There is a slight change in the behavior of the waves representing the solutions $x_1$ and $x_2$ when the damping parameter $c_2$ varies, as presented in Figure 13(a) and (b), respectively.

Figure 13.

Time behavior of the solutions $x_1,x_2,$ and $\phi$ in parts (a), (b), and (c), respectively, when ${\mathit{c}}_{\mathbf{2}}=\left(\mathbf{0.01},\mathbf{0.02},\mathbf{0.03}\right),$ ${\mathit{c}}_{\mathbf{1}}=\mathbf{0.01},$ and $\mathit{\omega }=\mathbf{4.242}$ .

The advantageous effect of varying frequency values $\omega$ on the behavior of the solutions $x_1,x_2,$ and $\phi$ is illustrated in the curves drawn in Figure 14. These curves are graphed considering the values $\omega =\left(4.472,4.242,3.872\right),$ $c_1=0.01$ and $c_2=0.01$ . It is evident that the wave amplitudes rise as the values of $\omega$ increase, while the number of waves and the corresponding lengths remain unchanged.

Figure 14.

Behavior of the time dependent functions $x_1,x_2,$ and $\phi$ in parts (a), (b), and (c), respectively, at $\mathit{\omega }=\left(\mathbf{4.472},\mathbf{4.242},\mathbf{3.872}\right),$ ${\mathit{c}}_{\mathbf{1}}=\mathbf{0.01},$ and ${\mathit{c}}_{\mathbf{2}}=\mathbf{0.01}$ .

Steady-state solutions

The key objective of this section is to analyze the oscillations of the examined system at the steady-state scenario when two resonance cases are examined simultaneously, that is, $p_1\approx 1,p_0\approx w$ . Accordingly, we assume the left-hand sides of the equations in the system (35) to be zero ³⁰ in order to find the fixed points for the amplitudes $y_s\left(s=1,2\right)$ and the phases ${\psi }_s$ . Consequently, setting ${\dot{y}}_s={\dot{\psi }}_s=0$ , we establish the subsequent algebraic system, which includes four equations involving the functions $y_s$ and ${\psi }_s$ .

$\begin{array}{l}{f}_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\psi }_1+c_1{y}_1=0,\\ f_1^{*}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\psi }_1+2\left\{{\sigma }_1-\left[\frac{3\left(1-w^2\right)w^2}{4\left(1-4w^2\right)}{y}_2^2+\frac{3}{2\left(1+\lambda \right)}\left(\lambda {\alpha }_1{x}_{1e}^2+{\left(x_{2e}-x_{1e}\right)}^2{\alpha }_2\right)\right]\right\}{y}_1=0,\\ f_0^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\psi }_2+w{c}_2{y}_2=0,\\ f_0^{*}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\psi }_2-2w\left\{{\sigma }_2-\left[\frac{\left(1+7w^2-8w^4\right)w^2}{16\left(1-4w^2\right)}{y}_2^2-\frac{3\left(1-w^2\right)w^2}{4\left(1-4w^2\right)}{y}_1^2\right]\right\}{y}_2=0.\end{array}$

(36)

Thus, we exclude the modified phases ${\psi }_s$ from the previous system, resulting in the below two non-linear algebraic equations regarding $h_s$ and ${\sigma }_s$

$\begin{array}{l}{f_1^{*}}^2=y_1^2\left\{c_1^2+{\left[2{\sigma }_1-\frac{3\left(1-w^2\right)w^2}{2\left(1-4w^2\right)}{y}_2^2-\frac{3}{\left(1+\lambda \right)}\left({\lambda \alpha }_1{x}_{1e}^2+{\left(x_{2e}-x_{1e}\right)}^2{\alpha }_2\right)\right]}^2\right\},\\ {f_0^{*}}^2=w^2{y}_2^2\left\{c_2^2+{\left[2{\sigma }_2+\frac{\left(1+7w^2-8w^4\right)w}{8\left(1-4w^2\right)}{y}_2^2-\frac{3\left(1-w^2\right)w}{2\left(1-4w^2\right)}{y}_1^2\right]}^2\right\}.\end{array}$

(37)

The actual solutions of system (37), with specified values of $f_0^{*}$ and $f_1^{*}$ , as well as ${\sigma }_s$ , represent the steady oscillation amplitudes. By varying the values of ${\sigma }_s$ , one can identify the resonance response curves. However, the assessment of stability plays a crucial role in studying vibrations at steady-state. To investigate this further, the system’s behavior near its fixed points is examined. Consequently, the substitutions detailed below are applied in equation (35) for this specific purpose

$\begin{array}{l}{y}_1=h_{10}+h_{11},{\psi }_1={\theta }_{10}+{\theta }_{11},\\ y_2=h_{20}+h_{21},{\psi }_2={\theta }_{20}+{\theta }_{21}.\end{array}$

(38)

Here, $h_{s0}$ and ${\theta }_{s0}$ represent the steady-state solutions, while $h_{s1}$ and ${\theta }_{s1}$ correspond to relatively small disturbances relative to $h_{s0}$ and ${\theta }_{s0}$ . Consequently, due to linearization and the fixed points of (35), one obtains

$\begin{array}{l}\frac{d{h}_{11}}{d\tau }=\frac{1}{2}\left(f_1^{*}{\theta }_{11}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}-c_1{h}_{11}\right),h_{10}\frac{d{\theta }_{11}}{d\tau }=\frac{f_1^{*}}{2}{\theta }_{11}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}+\left\{{\sigma }_1-[\frac{3\left(1-w^2\right)w^2}{4\left(1-4w^2\right)}\left(h_{20}^2{h}_{11}+2h_{10}{h}_{20}{h}_{21}\right)+\frac{3}{2\left(1+\lambda \right)}\left(\lambda {\alpha }_1{x}_{1e}^2+{\left(x_{2e}-x_{1e}\right)}^2{\alpha }_2\right)]\right\}{h}_{11},\\ w\frac{d{h}_{21}}{d\tau }=\frac{1}{2}\left(f_0^{*}{\theta }_{21}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{20}-w{c}_2{h}_{21}\right),\\ w{h}_{20}\frac{d{\theta }_{21}}{d\tau }=\frac{f_0^{*}}{2}{\theta }_{21}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}+w\left\{{\sigma }_2{h}_{21}+\frac{3\left(1+7w^2-8w^4\right)w}{16\left(1-4w^2\right)}{h}_{20}^2{h}_{21}-\frac{3\left(1-w^2\right)w}{4\left(1-4w^2\right)}\left(h_{10}^2{h}_{21}+2h_{10}{h}_{20}{h}_{11}\right)\right\}.\end{array}$

(39)

Given that $h_{s1}$ and ${\theta }_{s1}\left(s=1,2\right)$ are perturbed functions of the amplitudes and phases in the linear system mentioned above, one expresses their solutions using the linear function $q_k{e}^{\xi \tau }\left(k=1,2,3,4\right)$ in exponential form, where $q_s$ and $\xi$ are constants and the perturbation’s eigenvalues. To ensure the asymptotic stability of the steady-state solutions, the roots’ real components of the characteristic equation in (39) are negative ³¹

${\xi }^4+{\Gamma }_1{\xi }^3+{\Gamma }_2{\xi }^2+{\Gamma }_3\xi +{\Gamma }_4=0,$

(40)

where ${\Gamma }_k\left(k=1,2,3,4\right)$ are given in terms of $h_{s0},{\theta }_{s0},$ and $q_s$ , as follows:

$\begin{array}{l}{\Gamma }_1=\frac{1}{2}\left(c_1+c_2-\frac{f_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}}{h_{10}}-\frac{f_0^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}}{h_{20}w}\right),\\ {\Gamma }_2=\frac{1}{32h_{10}{h}_{20}(\lambda +1)w(4w^2-1)}\left\{{4h}_{20}\left[2w\left(1-4w^2\right)\left(\lambda +1\right)\mathit{cos}\hspace{0.17em}{\theta }_{10}\left(c_1+c_2\right)f_1^{*}+2w{c}_1{c}_2{h}_{10}\left(4w^2-1\right)\times \left(\lambda +1\right)+f_1^{*}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}\left(3h_{20}^2\left(\lambda +1\right)\left(w^2-1\right)w^3+2w\left(4w^2-1\right)\left[3{\alpha }_1\lambda x_{1e}^2-2\left(\lambda +1\right){\sigma }_1+3{\alpha }_2{\left(x_{1e}-x_{2e}\right)}^2\right]\right)\right]+f_0^{*}\left(\lambda +1\right)\left(8\left(4w^2-1\right)\mathit{cos}\hspace{0.17em}{\theta }_{20}\left[f_1^{*}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}-\left(c_1+c_2\right)h_{10}\right]+h_{10}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{20}\left[3w\left(1-w^2\right)\times \left(h_{20}^2\left(8w^2+1\right)-4h_{10}^2\right)+16{\sigma }_2\left(1-4w^2\right)\right]\right)\right\},\\ {\Gamma }_3=\frac{1}{64h_{10}{h}_{20}(\lambda +1)w(4w^2-1)}\left\{\{\left(\lambda +1\right)c_1(8c_2{f}_1^{*}{h}_{20}w\left(1-4w^2\right)\mathit{cos}\hspace{0.17em}{\theta }_{10}+f_0^{*}[8f_1^{*}\left(4w^2-1\right)\times \mathit{cos}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}+h_{10}(\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{20}-8c_2\left(4w^2-1\right)\mathit{cos}\hspace{0.17em}{\theta }_{20}+12h_{10}^{2}w\left(w^2-1\right)-3h_{20}^{2}w\left(8w^4-7w^2-1\right)+16{\sigma }_2\left(1-4w^2\right))])+f_1^{*}[4\hspace{0.17em}{\mathrm{c}}_2(2f_0^{*}(\lambda +1)\left(4w^2-1\right)\mathit{cos}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}+h_{20}w\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}(3h_{20}^2\left(\lambda +1\right)w^2\times \left(w^2-1\right)+2\left(4w^2-1\right)\left(3{\alpha }_2{x}_{2e}^2-6{\alpha }_2{x}_{1e}{x}_{2e}+3x_{1e}^2\left({\alpha }_1\lambda +{\alpha }_2\right)-2\left(\lambda +1\right){\sigma }_1\right)))+f_0^{*}(12h_{10}^2\left(\lambda +1\right)\times w\left(1-w^2\right)\mathit{sin}\hspace{0.17em}{\theta }_{20}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}+3\left(\lambda +1\right)w\left(w^2-1\right)\left[\left(8w^2+1\right)\mathit{sin}\hspace{0.17em}{\theta }_{20}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}-4w\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}\right]h_{20}^2-8\left(4w^2-1\right)(3{\alpha }_2{x}_{2e}^2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}-6{\alpha }_2{x}_{1e}{x}_{2e}\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}+3x_{1e}^2\left({\alpha }_1\lambda +{\alpha }_2\right)\mathit{sin}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}-2\left(\lambda +1\right)\left({\sigma }_1\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{10}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{20}+{\sigma }_2\hspace{0.17em}\mathit{sin}\hspace{0.17em}{\theta }_{20}\hspace{0.17em}\mathit{cos}\hspace{0.17em}{\theta }_{10}\right)))]\right\},\\ {\Gamma }_4=\frac{1}{256w{\left(1-4w^2\right)}^2\left(1+\lambda \right)h_{10}{h}_{20}}\left\{f_0^{*}{f}_1^{*}\left(2\left(4w^2-1\right)\left(1+\lambda \right)\mathit{cos}\hspace{0.17em}{\theta }_{10}{c}_1\left[8\left(4w^2-1\right)\mathit{cos}\hspace{0.17em}{\theta }_{20}{c}_2+{\mathit{sin}\hspace{0.17em}\theta }_{20}\left(3w\left(w^2-1\right)\left[\left(1+8w^2\right)h_{20}^2-4h_{10}^2\right]+16\left(4w^2-1\right){\sigma }_2\right)\right]-{\mathit{sin}\hspace{0.17em}\theta }_{10}\left(8\left(4w^2-1\right){\mathit{cos}\hspace{0.17em}\theta }_{20}{c}_2\times \left(3w^2\left(w^2-1\right)\left(1+\lambda \right)h_{20}^2+2\left(4w^2-1\right)\left({3\lambda x_{1e}^2{\alpha }_1+3\left(x_{1e}-x_{2e}\right)}^2{\alpha }_2-2\left(1+\lambda \right){\sigma }_1\right)\right)+2\left(4w^2-1\right)\times (3\lambda x_{1e}^2{\alpha }_1+3{\left(x_{1e}-x_{2e}\right)}^2{\alpha }_2-2\left(1+\lambda \right){\sigma }_1\left(9w^2\left(w^2-1\right)\left(1+\lambda \right)h_{20}^2-2\left(4w^2-1\right)\left[{3\lambda x_{1e}^2{\alpha }_1+3{\left(x_{1e}-x_{2e}\right)}^2\alpha }_2-2\left(1+\lambda \right){\sigma }_1\right]\right)+\left(3w^2\left(w^2-1\right)\left(1+\lambda \right)h_{20}^2+2\left(4w^2-1\right)\right)\left[3\lambda x_{1e}^2{\alpha }_1+3{\left(x_{1e}-x_{2e}\right)}^2{\alpha }_2-2\left(1+\lambda \right){\sigma }_1\right]\left[3w\left(8w^4-7w^2-1\right)h_{20}^2+16\left(4w^2-1\right){\sigma }_2\right]\right)\right)\right\}.\end{array}$

For certain steady-state solutions, the conditions required and sufficient for achieving stability at the steady-state solution can be outlined according to RHC ³² as follows:

$\begin{array}{l}{\Gamma }_1>0,{\Gamma }_1{\Gamma }_2-{\Gamma }_3>0,\\ {\Gamma }_4>0,{\Gamma }_3\left({\Gamma }_1{\Gamma }_2-{\Gamma }_3\right)-{\Gamma }_4{\Gamma }_1^2>0.\end{array}$

(41)

The most direct feedback from applying the RHC is determining whether a system is stable or unstable. The stability of the system is assured if all elements in the first column of the Routh array have the same sign and none are zero, while the system is unstable if any element in the first column of the Routh array changes sign. By applying these criteria to systems with varying parameters, one can analyze how changes in system parameters affect stability. This sensitivity analysis helps in designing systems that remain stable under different operating conditions. The benefit from the stability RHC is valuable in control system design. It helps in tuning the controller parameters to achieve the desired stability margins and ensure robust performance. By systematically applying the RHC, researchers can make informed decisions about the design, analysis, and tuning of dynamic systems to ensure desired stability and performance characteristics.^33,34

The stability analysis

In this section, we apply nonlinear stability analysis for the RHC to investigate the dynamic behavior of the studied oscillating system. ³⁵ The system’s motion is examined by considering the force $F\left(t\right)$ acting on the spring’s direction and the moment $M\left(t\right)$ at the pivot point $O$ . Stability criteria are implemented along with simulations of the system’s equations. It is found that the parameters of detuning ${\sigma }_s$ , frequency $\omega$ , and coefficients of damping $c_s$ are crucial factors influencing stability. To determine the stability ranges of system (37), various system settings were tested through a predetermined process.

Figures 15–17 show the FRC of $y_s$ against ${\sigma }_s$ concerning the system’s fixed points. They are plots depicting the time-dependent changes in amplitudes $y_s$ across various parametric regions when the following data are considered:

$\begin{array}{l}g=9.8m/s^2,a=0.7m,b=0.4m,l_0=0.9m,m=1kg,f_0^{*}=0.0004,\\ f_1^{*}=0.0001,f_2^{*}=0.008,k_1=30,{\Lambda }_1=0.1,{\Lambda }_2=0.1,\Omega =0.1,{\Omega }_0=0.01,\\ {\Omega }_1=0.0001,{\Omega }_2=0.1,{\sigma }_2=0.01.\end{array}$

Figure 15.

FRC of the planes at ${\mathit{c}}_{\mathbf{1}}\left(=\mathbf{0.01},\mathbf{0.015},\mathbf{0.02}\right)$ when ${\mathit{\sigma }}_{\mathbf{2}}=\mathbf{0.01},{\mathit{c}}_{\mathbf{2}}=\mathbf{0.01},$ and $\mathit{\omega }=\mathbf{4.242}$ : (a) in the plane ${\sigma }_1{y}_1$ , (b) in the plane ${\sigma }_1{y}_2$ .

Figure 16.

FRC of the planes at $c_2\left(=0.01,0.011,0.012\right)$ when ${\sigma }_2=0.01,c_1=0.01,$ and $\omega =4.242$ : (a) in the plane ${\sigma }_1{y}_1$ , (b) in the plane ${\sigma }_1{y}_2$ .

Figure 17.

FRC of the planes at $\omega \left(=4.472,4.242,3.872\right)$ when ${\sigma }_2=0.01,c_1=0.01,$ and $c_2=0.01$ : (a) in the plane ${\sigma }_1{y}_1$ , (b) in the plane ${\sigma }_1{y}_2$ .

The presented stability and instability areas in Figures 15–17 are plotted when the $c_1,c_2,$ and $\omega$ have the aforementioned various values. The range of ${\sigma }_1$ is chosen to be during the interval $\left[-0.1,0.1\right]$ . Based on the change in $c_1$ , that is, the plotted curves in Figure 15, one observes the stability zones at $c_1$ equals $0.01,0.15,$ and $0.2$ are being in the ranges ${\sigma }_1>0.005,{\sigma }_1>0.001,$ and ${\sigma }_1>-0.003$ , respectively. Whereas the instability regions at the same considered values of $c_1$ satisfy the inequalities ${\sigma }_1\le 0.005,{\sigma }_1\le 0.001,$ and ${\sigma }_1\le -0.003$ . It is noted that each curve in the parts of this figure has one critical fixed point and one peak point.

The inspection of the graphed FRC in portions of Figure 16 shows that they are calculated at ${\sigma }_2=0.01,c_1=0.01$ , and $\omega =4.242$ when $c_2\left(=0.01,0.011,0.012\right)$ . It is seen that one drawn curve is calculated with the variation of $c_2$ values, as explored in Figure 16(a), which means that the values of this parameter don’t affect on the variation of the FRC in the plane ${\sigma }_1{y}_1$ , while the same values of $c_2$ have an excellent influence on the FRC in the plane ${\sigma }_1{y}_2$ , as shown in Figure 16(b). The explanation lies in the first and second equations in the system (37), where the first one is dependent on $c_1$ , while the second depends on $c_2$ . Stability zones and instability ones are matched with the inequalities ${\sigma }_1>0.003$ and ${\sigma }_1\le 0.003$ , respectively. As previously mentioned, two fixed points are observed: one peak and one critical.

The plotted FRC in Figure 17 is determined when at ${\sigma }_2=0.01,c_1=0.01,$ and $c_2=0.01$ at $\omega \left(=4.472,4.242,3.872\right).$ Analyzing the curves in this figure shows that these curves are impacted by the change $\omega$ values, in which the stable area and unstable one are found in the ranges ${\sigma }_1>0.005$ and ${\sigma }_1\le 0.005$ , respectively.

The drawn response curves in Figure 18 are calculated when at $\omega =4.242,c_1=0.01,$ and when ${\sigma }_2$ has the values $\left(-0.01,0,0.01\right)$ . An examination of these curves explores that these curves are not impacted by the change ${\sigma }_2$ values, in which the stable areas are found when ${\sigma }_2=0.01$ at the ranges ${\sigma }_1\ge 0.005$ , while the unstable regions are found in the ranges ${\sigma }_1<0.005$ . On the other hand, all the represented areas when ${\sigma }_2=0$ and ${\sigma }_2=-0.01$ are unstable during the plotted domain of ${\sigma }_1$ . The reason goes back to the fact that at least one condition of the RHC is not satisfied.

Figure 18.

FRC of the planes ${\sigma }_1{y}_s$ when $c_1=0.01,c_2=0.01$ , and $\omega =4.242$ : (a) and (b) at ${\sigma }_2=0.01$ ; (c) and (d) at ${\sigma }_2=0$ ; and (e) and (f) at ${\sigma }_2=-0.01$ .

The graphed curves in Figure 19 show the FRC in the plane ${\sigma }_2{y}_s\left(s=1,2\right)$ when ${\mathit{c}}_{\mathbf{1}}\left(=0.01,0.015,0.02\right)$ have various values. These curves are calculated at ${\mathit{\sigma }}_{\mathbf{1}}=0.01,{\mathit{c}}_{\mathbf{2}}=0.01,$ and $\mathit{\omega }=4.242$ . It is noted that the domains of stable regions are smaller than those of unstable regions. The stable’s regions are found at $0.003\le {\sigma }_2\le 0.016$ , while the unstable ones are presented at $-0.1\le {\sigma }_1<0.003$ and $0.016<{\sigma }_1\le 0.1$

Figure 19.

FRC of the planes at ${\mathit{c}}_{\mathbf{1}}\left(=\mathbf{0.01},\mathbf{0.015},\mathbf{0.02}\right)$ when ${\mathit{\sigma }}_{\mathbf{1}}=\mathbf{0.01},{\mathit{c}}_{\mathbf{2}}=\mathbf{0.01},$ and $\mathit{\omega }=\mathbf{4.242}$ : (a) in the plane ${\sigma }_2{y}_1$ , (b) in the plane ${\sigma }_2{y}_2$ .

Conclusion

Exploring the planar forced oscillations, a particle attached to its support point through a series of two nonlinear springs and two viscous dampers has been examined. This point follows an elliptic route due to the presence of an external harmonic force along the spring’s elongation and a viscous harmonic moment applied at the support point. The governing two differential equations have been formulated using Lagrange’s equations besides an algebraic equation. These equations have been solved for the non-resonant and resonant scenarios up to the third order of approximation utilizing the MMTS, in light of three distinct time scales. For the resonant scenarios, two cases of external resonance have been examined simultaneously. Consequently, the solvability criteria and the ME have been obtained. The RHCs have been employed to analyze the stability and instability areas for the drawn FRC. Detailed discussions and visual depictions of the motion temporal histories and stability regions are provided. These visuals highlight the advantageous effects of varying physical parameters on the system’s operational performance. The achieved outcomes are considered generalizations of those which were gained in Ref. 23 when the springs’ suspension point was stationary. These outcomes can be used to improve the design and performance of mechanical systems such as shock absorbers, suspension systems, and vibration isolators that rely on spring dynamics. Moreover, they enhance the control and stability of robotic arms and joints that use spring mechanisms to mimic muscle movements. They also optimize the behavior of landing gear systems in aircraft, ensuring better shock absorption and stability during landing.

Future work

Future work could branch out in the following several important directions to enhance both the theoretical understanding and practical applications of the studied system.

• The current study focuses can be extended to be in three dimensions, which provide a more comprehensive understanding of the system’s dynamics. Studying the role of spatial motion in the performance of nonlinear springs and dampers might reveal further resonance phenomena and intricate dynamic patterns.

• Further work could investigate the system’s behavior through higher-order asymptotic approximations, providing deeper insights into subtle nonlinear effects, especially in highly nonlinear regimes or under stronger excitations where more intricate interactions may occur.

• An obvious extension would be to explore internal resonances, where interactions between multiple natural frequencies occur, and potentially reveal novel dynamic behaviors such as mode-to-mode energy transfer or intricate bifurcation patterns.

• Exploring the effects of various damping models, such as nonlinear damping or fractional damping,^36,37 could be a focus of future research. The He’s frequency and the semi-inverse method^38–41 can be used to establish future results. This would improve the fidelity of real-world damping approximations and enable more accurate predictions for practical scenarios.

• Examining the system’s behavior in response to parametric excitations, where parameters such as stiffness or damping coefficients are periodically adjusted, would enable a more thorough investigation of the system’s stability limits and dynamic behavior in practical contexts like engineering applications that deal with time-varying external conditions.

• To confirm the theoretical results, experimental studies could be conducted, or the findings could be applied to real-world applications, such as vibration absorbers, mechanical linkages, or suspension systems.

• Future work could include large-scale numerical simulations that explore a wider range of parameter values, beyond the scope of the analytical approximations used here. Exploring the system’s capability to demonstrate complex nonlinear phenomena, like bifurcations, period-doubling, or chaotic dynamics, could be a focus of future studies, particularly under conditions of varying external forcing amplitudes or frequencies.

Focusing on these aspects in future work can substantially improve both the theoretical comprehension and practical implementation of the forced oscillatory system discussed in this paper.