Modeling and Dynamic Analysis for a Motion of Mounted-Based Axisymmetric Rigid Body under Self-Excited Vibrations in an Attractive Newtonian Field

Abstract

In this article, the mathematical modeling and dynamic analysis for a motion of a rigid body mounted on a vibrating base affected by a strongly nonlinear damping and an attractive Newtonian field are investigated. The governing equations are derived, and specific conditions for asymptotic stability of the motion are stated. Results of analytical analysis are used to identify some distinct types of motion, including the cases that are apparently regular or chaotic. A reduction technique via averaging is used to obtain the amplitude equation and the response diagrams to identify regions where the jump phenomenon may occur. Moreover, a simple analytical relationship between the parameters of the system, describing whether the jump phenomenon will be possible, is obtained. Homoclinic bifurcation diagrams are used to argue apparently that the chaotic behavior is slowly approaching a strange attractor at specific values of system parameters. Results of numerical solutions using the fourth-order Runge-Kutta method are closely coincided with the analytical ones to verify all types of motion.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In engineering applications, the vibrational motion is typically produced internally from sources of noise such as motors, bearings, and other moving parts in machines. Therefore, it is necessary to take into account the most significant destabilizing source that can seriously degrade the operation and, in some cases, it leads to chaotic behavior attached by a catastrophic failure of mechanical system devices. Generally, the vibration is described to characterize the machines dynamic force in terms of a single dominant frequency called the operating frequency. This particular modeling, indeed, contains a generic assumption; namely, the excitation force is ideal, which means that the dynamics of excitations are not affected by the dynamics of responses. This assumption becomes more untenable as the operating frequency increases further into the audio or noise ranges, cf. [1].

Clearly, the sources of vibrations output in machines are considered due to the presence of several moving parts such as reciprocating pistons and rotating shafts. For instance, reciprocating pumps and washing machines both generate sinusoidal vibration type due to the reciprocating primary motion or the primary rotational motion, respectively. Indeed, this phenomenon generates dynamic forces which might be sinusoidal, random, or transient. Therefore, adequately, these machines can be modeled by the motion of a rigid body mounted on a vibrating base with the presence of dampers, cf. [2].

In rigid body dynamics, the recent studies have focused on the damped motion of self-excited rigid bodies. This type of motion is described as a body free to rotate about a fixed point in a resistant medium acted upon a time-dependent torque vector. The external torques are often arising from internal reactions which do not appreciably change the distribution of their mass, cf. [3–6]. The interesting engineering examples of this motion are gyroscopes and dynamics of rotating missiles. However, a modern case study of a self-excited motion of a rigid body can arise from the attitude motion of a rigid spacecraft subject to thruster torques, cf. [7–9].

Due to the nature of the vibrational motion, only a few mechanical systems are dominated by motion in one dimension and thus single degree of freedom (SDOF) models can be used. But, in general, most of the real mounted-based systems are always modeled as a piece of equipment vibrating about its center of mass in two or more of degrees of freedom (DOFs). So, as consequences in the presence of all six modes’ modeling, they are coupled resulting in six coupled equations to describe such motion, cf. [9, 10].

Self-excited vibrational systems might be characterized by a principal feature which is resided in the damping behavior due to the dissipative forces. Thus, even an infinitesimal disturbance causes sustained oscillations in the system; however, when the perturbed displacement becomes sufficiently large, the damping becomes positive and suppresses further increase in its amplitude.

In this work, the following specific form of damping function is used, which is considered a generalized form of van der Pol’s damping.$$\begin{array}{}\text{(1)}& hx,y=\epsilon \lambda +\mu x^2+\nu y^{2}y,y=\dot{x}.\end{array}$$where the parameters $\lambda$, $\mu$, and $\nu$ are real constants; $0\le \epsilon \ll 1$.

Using the dimensionless small parameter $\epsilon$, the concept of homoclinic bifurcation and the route to chaos are necessary in the study due to the severe existence in complex systems, in particular, with strong damping forms. Homoclinic bifurcation takes place due to a parameter when the change of the parameter causes a homoclinic path to appear at a certain value and then it disappears again, cf. [11–13].

The generation of chaotic set can be arisen by such parameter to bifurcate at an intermediate value into a solution which circuits the origin point (for example, say twice). Further increasing shows further bifurcation has taken place which results in paths circuiting the origin many times. If the process continues with further period-doubling circuits, then the parameter sequence converges to a certain value beyond which all these periods are presented consisting of a chaotic attracting set. Then, this configuration is known as a strange attractor which is developed parametrically through processes of the period-doubling, cf. [14, 15]. Consequently, at each bifurcation, there exists an unstable periodic orbit that will be always dwelled in, so that the final attractor has to include unbounded sets of unstable limit cycles. This means that the periodic solution becomes closer to the destabilizing effect of the homoclinic paths which appear as the parameter increases. Then, now, the chaos is generated by the effect of that parameter, cf. [16–18].

In common, two main methods are used to measure the existence of chaotic behavior in nonlinear (nonautonomous) systems: Lyapunov exponents and Poincare’s maps. The method of Lyapunov exponents measures the exponential rates of divergence (convergence) of nearby orbits of an attractor in the state space. Therefore, it is the most precisely powerful tool for the identification of motion in dynamical systems. In [19], an algorithm is developed based on tracing the evolution of an initial sphere of small perturbation to a nominal trajectory. In [20, 21], a method to compute Lyapunov exponents from an experimental time series based on the QR methods is presented.

In this work, we study the dynamics of a rigid body mounted on a vibrating base with a generalized nutational damping, using simulation, numerical continuation, and analytical stability theory. We begin by stating the equations of motion, the stability conditions in terms of problem parameters, and then theorems to predict the existence of periodic solutions. Analytical methods and numerical simulation results are used to identify the distinct types of motion, including persistent oscillations that appear to be a stable, periodic, or chaotic in their nature. The theoretical results will be applied to the tackled application beside the study of bifurcation and routes to chaos of the governing system. Melnikov’s method is used to transpire that the ranges of chaotic behavior are affected by the change of some certain parameters. The verification and compatibility of theoretical results with the numerical ones through bifurcation and chaos diagrams as well as phase plane trajectories are taken into consideration in the study.

The rest of the article is summarized as follows: Section 2 presents the modeling of a damped motion of mounted-based rigid body as an SDOF model in attractive Newtonian field. Section 3 discusses the study of stability analysis in both autonomous and nonautonomous cases of motion. Section 4 discusses the existence of periodic solutions for some specific cases of motion. Section 5 presents a reduction technique via averaging to capture the jump phenomenon. Sections 6 and 7 discuss the study of homoclinic bifurcation and routes to chaos. In the last section, the conclusion is given.

2. The Modeling and Governing Equations

Consider that a rigid body of mass $m$ is mounted on a vibrating base of a sinusoidal excitation frequency $\omega$ as shown in Figure 1 with two sets of coordinates: one is fixed in space and the other is moving with the body. The two coordinates are coincident when the body is in equilibrium under the action of gravity force alone. The motion of the body is described by giving it a displacement relative to the inertial fixed axes. Let us choose the axes $OX$, $OY$, and $OZ$ to represent the fixed frame in space with unit vectors $\widehat{\mathsf{I}}$, $\widehat{\mathsf{J}}$, and $\widehat{\mathsf{K}}$, respectively, and the axes $ox$, $oy$, and $oz$ to represent the principal axes constructed on the body at the fixed point $o$ with unit vectors $\widehat{i}$, $\widehat{j}$, and $\widehat{k}$, respectively.

[figure omitted; refer to PDF]

Hence, the translational displacement of the center of mass (C.G.) with respect to the fixed axes $X,Y,Z$ reads$$\begin{array}{}\text{(2)}& r_c=x_c\widehat{I}+y_c\widehat{J}+z_c\widehat{K}+dt,dt=d_{1}t\widehat{I}+d_{2}t\widehat{J}+d_{3}t\widehat{K}.\end{array}$$here, $dt$ is the excitation displacement vector. Then,$$\begin{array}{}\text{(3)}& r_c=d_{1}t+r_c\sin \varphi \sin \theta \widehat{I}+d_{2}t-r_c\sin \varphi \cos \theta \widehat{J}+d_{3}t+r_c\cos \theta \widehat{K}.\end{array}$$

Then, we obtain the velocity of the center of mass:$$\begin{array}{}\text{(4)}& v_c=\frac{d{r}_c}{dt}.\end{array}$$

We assume that the principal axes rotate in space with the same angular velocity of the body $\varpi$ defined by (in terms of Euler’s angles)$$\begin{array}{}\text{(5)}& \varpi t=\dot{\psi }\sin \varphi \sin \theta +\stackrel{.}{\theta }\cos \varphi \widehat{i}+\dot{\psi }\cos \varphi \sin \theta -\stackrel{.}{\theta }\sin \varphi \widehat{j}+\dot{\varphi }+\dot{\psi }\cos \theta \widehat{k},\end{array}$$where $\theta t$, $\varphi t$, and $\psi t$ are the nutation, spin, and precession angles, respectively, and $\cdot$ refers to the derivative with respect to time.

Let $\widehat{K}t$ be the Poisson vector with the following direction cosines with respect to the moving axes:$$\begin{array}{}\text{(6)}& \widehat{K}t={\gamma }_{1}t\widehat{i}+{\gamma }_{2}t\widehat{j}+{\gamma }_{3}t\widehat{k}.\end{array}$$

The relations of the direction cosines with Euler’s angles read$$\begin{array}{}\text{(7)}& {\gamma }_1=\sin \theta \sin \varphi ,{\gamma }_2=\sin \theta \cos \varphi ,{\gamma }_3=\cos \theta .\end{array}$$

The potential energy of a Newtonian force field is considered by the Newtonian attraction of a point $P$, of mass $\mathbb{M}$, on the body mass $m$. If $P$ is located on the negative part of the axis $OZ$, then we have$$\begin{array}{}\text{(8)}& \widehat{R}=d_{3}t-R\widehat{K},d_o\ll Rg=\mathbb{G}\frac{\mathbb{M}}{R^2}.\end{array}$$here, $\mathbb{G}$ is the universal gravitation constant and $R$ is the distance between point $P$ and the fixed point $O$. In general, the fact that $R$ is larger than the dimensions of the body is taken into account, cf. [6]. Hence, the potential energy reads$$\begin{array}{}\text{(9)}& U_N=U_0+U_1+U_2+U_3+\dots \end{array}$$where$$\begin{array}{c}\text{(10)}& U_0=mg{d}_{3}t-R,\text{(11)}& U_1=mg{r}_c\cdot \widehat{K},\text{(12)}& U_2=-\frac{3g}{2R}A{\gamma }_1^2+B{\gamma }_2^2+C{\gamma }_3^2,\text{(13)}& U_3=-\frac{3g}{2R^2}{J}_{xxx}+J_{xyy}+J_{xzz}{\gamma }_1+J_{xxy}+J_{yyy}+J_{yzz}{\gamma }_2+J_{xxz}+J_{yyz}+J_{zzz}{\gamma }_3+\frac{5g}{2R^2}{J}_{xxx}{\gamma }_1^3+J_{yyy}{\gamma }_2^3+J_{zzz}{\gamma }_3^3+3J_{xxy}{\gamma }_1^2{\gamma }_2+3J_{xxz}{\gamma }_1^2{\gamma }_3\\ +3J_{yyx}{\gamma }_2^2{\gamma }_1+3J_{yyz}{\gamma }_2^2{\gamma }_3+3J_{xzz}{\gamma }_1{\gamma }_3^2+3J_{yzz}{\gamma }_2{\gamma }_3^2+6J_{xyz}{\gamma }_1{\gamma }_2{\gamma }_3.\end{array}$$here, $A$, $B$, and $C$ are the principal moments of inertia along the moving axes $ox$, $oy$, and $oz$, respectively, and $J_{\dots }$ is the third moment of inertia of the body. The total potential energy $U$ of the mass with respect to the inertial frame $OXYZ$ is$$\begin{array}{}\text{(14)}& U=U_N.\end{array}$$

Consider the following approximation:$$\begin{array}{}\text{(15)}& U\approx U_0+U_1+U_2,\end{array}$$

and then$$\begin{array}{}\text{(16)}& U=-mgR+mg{x}_c\sin \theta \sin \varphi +y_c\sin \theta \cos \varphi +z_c\cos \theta -\frac{3g}{2R}A{\sin }^2\theta {\sin }^2\varphi +B{\sin }^2\theta {\cos }^2\varphi +C{\cos }^2\theta +mg{d}_{3}t.\end{array}$$

The kinetic energy reads$$\begin{array}{}\text{(17)}& K.E.=\frac{1}{2}m{v}_c^2+\frac{1}{2}\varpi I_c\varpi .\end{array}$$

Then, the system admits the following Lagrangian $\mathbb{L}$ on the Lie group SO (3):$$\begin{array}{}\text{(18)}& \mathbb{L}=\frac{1}{2}m{v}_c^2+\frac{1}{2}\varpi {\mathsf{I}}_c\varpi -U,\end{array}$$where $I_c=I_1,I_2,I_3$ refers to the principal moments of inertia about the equatorial and polar directions through the C.G. of the rigid body.

For simplicity, the mass of rigid body satisfies the symmetric conditions and the excitation displacement vector is only on the fixed vertical direction as follows:$$\begin{array}{}\text{(19)}& I_1=I_2,x_c=y_c=0,d_{1}t=d_{2}t=0,d_{3}t=dt=d_o\sin \omega t,\text{(20)}& A=B=I,I_3=C,I=I_1+m{z}_c^2.\end{array}$$

Then,$$\begin{array}{}\text{(21)}& U=-mgR+mg{z}_c\cos \theta -\frac{3g}{2R}2I{\sin }^2\theta +I_3{\cos }^2\theta +Mgdt,\end{array}$$where $dt=d_0\sin \omega t$ is the time-varying amplitude of the vertical support motion along the vertical fixed axis $OZ$ with constant amplitude $d_o$ and forced frequency $\omega$.

Using Routhian’s function $\Re$,$$\begin{array}{}\text{(22)}& \Re =\mathbb{L}-p_{\psi }\dot{\psi }{p}_{\psi },p_{\varphi },\theta -p_{\varphi }\dot{\varphi }{p}_{\psi },p_{\varphi },\theta ,\end{array}$$where $p_{\varphi }$ and $p_{\psi }$ are the corresponding generalized momenta to the cyclic coordinates $\varphi$ and $\psi$, respectively, cf. [8, 9, 22, 23]:$$\begin{array}{c}\text{(23)}& p_{\varphi }=I_3\dot{\varphi }+\dot{\psi }\cos \theta =\text{const}.,p_{\psi }=I\dot{\psi }{\sin }^2\theta +p_{\varphi }\cos \theta =\text{const.},\end{array}$$

and then the equation of motion reads$$\begin{array}{}\text{(24)}& \frac{d}{dt}\frac{\partial \Re }{\partial \stackrel{.}{\theta }}-\frac{\partial \Re }{\partial \theta }=-\mathfrak{D}\theta ,\stackrel{.}{\theta }+{\tilde{F}}_1\cos \omega t+{\tilde{F}}_2\sin \omega t.\end{array}$$

The right terms in equation (24) represent the damping term $\mathfrak{D}$ and the harmonic forcing terms of amplitudes ${\tilde{F}}_1$ and ${\tilde{F}}_2>0$ with circular frequency ($\omega$) applied to the body on such nonlinear spring.

The damping term $\mathfrak{D}$ is assumed to have the following generalized van der Pol’s damping form:$$\begin{array}{}\text{(25)}& \mathfrak{D}\theta ,\stackrel{.}{\theta }=\epsilon \tilde{\lambda }+\tilde{\mu }{\theta }^2+\tilde{\nu }{\stackrel{.}{\theta }}^2\stackrel{.}{\theta },\tilde{\lambda },\tilde{\mu },\tilde{\nu }\in ℝ,0\le \epsilon <1.\end{array}$$

Using equation (24) and if $p_{\psi }=p_{\varphi }=p_o$, the governing equation becomes$$\begin{array}{}\text{(26)}& \ddot{\theta }+\frac{1}{I}+\frac{p_o^2{1-\cos \theta }^2}{I^2{\sin }^3\theta }-\frac{mg{z}_c}{I}\sin \theta -\frac{3g}{2RI}4I-2I_3\sin \theta \cos \theta -\frac{m{z}_c}{I}\sin \theta \ddot{d}t=\frac{{\tilde{F}}_1}{I}\cos \omega t+\frac{{\tilde{F}}_2}{I}\sin \omega t.\end{array}$$

Equation (26) can be generalized as follows:$$\begin{array}{}\text{(27)}& \ddot{\theta }+h\theta ,\stackrel{.}{\theta }+f\theta ,t=et,\end{array}$$where$$\begin{array}{}\text{(28)}& f\theta ,t=\frac{p_o^2}{I^2}\frac{{1-\cos \theta }^2}{{\sin }^3\theta }-\frac{mg{z}_c}{I}\sin \theta -\frac{3g}{2RI}2I-I_3\sin 2\theta +\frac{1}{I}{\omega }^{2}m{z}_c{d}_o\sin \omega t\sin \theta ,\text{(29)}& h\theta ,\stackrel{.}{\theta }=\epsilon \lambda +\mu {\theta }^2+\nu {\stackrel{.}{\theta }}^2\stackrel{.}{\theta },et=F_1\cos \omega t+F_2\sin \omega t.\end{array}$$

Consider the following notations:$$\begin{array}{}\text{(30)}& \lambda =\frac{\tilde{\lambda }}{I},\mu =\frac{\tilde{\mu }}{I},\nu =\frac{\tilde{\nu }}{I},F_1=\frac{{\tilde{F}}_1}{I},F_2=\frac{{\tilde{F}}_2}{I},F=\sqrt{F_1^2+F_2^2}.\end{array}$$

For simplicity, let us assume that$$\begin{array}{}\text{(31)}& \alpha =\frac{p_o^2}{I^2},\beta =\frac{mg{z}_c}{I},\gamma =\frac{1}{I}{\omega }^{2}m{z}_c{d}_o,\kappa =\frac{3g}{2RI}2I-I_3.\end{array}$$

Then, equation (27) reads$$\begin{array}{}\text{(32)}& \ddot{\theta }+\epsilon \lambda +\mu {\theta }^2+\nu {\stackrel{.}{\theta }}^2\stackrel{.}{\theta }+\alpha \frac{{1-\cos \theta }^2}{{\sin }^3\theta }-\beta \sin \theta -\kappa \sin 2\theta +\gamma \sin \omega t\sin \theta =et.\end{array}$$

If $\theta =x$ and $\stackrel{.}{\theta }=y$, then equation (32) is converted to the following system:$$\begin{array}{}\text{(33a)}& \dot{x}=y,\text{(33b)}& \dot{y}=-\epsilon \lambda +\mu x^2+\nu y^{2}y-\alpha \frac{{1-\cos x}^2}{{\sin }^{3}x}+\beta \sin x+\kappa \sin 2x-\gamma \sin \omega t\sin x+et.\end{array}$$

Use the following approximations:$$\begin{array}{c}\text{(34)}& \sin \theta =\theta -\frac{{\theta }^3}{6},\frac{{1-\cos \theta }^2}{{\sin }^3\theta }=\frac{\theta }{4}+\frac{{\theta }^3}{12},{\omega }_n^2=\frac{\alpha }{4}-\beta -2\kappa ,\Omega =\frac{\alpha }{12}+\frac{\beta }{6}+\frac{4\kappa }{3},\end{array}$$

and then equation (32) reads$$\begin{array}{}\text{(35)}& \ddot{\theta }+\epsilon \lambda +\mu {\theta }^2+\nu {\stackrel{.}{\theta }}^2\stackrel{.}{\theta }+{\omega }_n^2+\gamma \sin \omega t\theta +\Omega -\frac{\gamma }{6}\sin \omega t{\theta }^3=et,\end{array}$$

and the approximate system of equation (35) reads$$\begin{array}{}\text{(36a)}& \dot{x}=y,\text{(36b)}& \dot{y}=-\epsilon \lambda +\mu x^2+\nu y^{2}y-{\omega }_n^2+\gamma \sin \omega tx-\Omega -\frac{\gamma }{6}\sin \omega t{x}^3+et.\end{array}$$

If $\theta \ll 1$, then equation (35) has the following linear form under a proposed fractional damping:$$\begin{array}{}\text{(37)}& \ddot{\theta }+\epsilon {\lambda }^L{D}^q\theta +{\omega }_n^2+\gamma \sin \omega t\theta =et,0<q<2,\end{array}$$where ${{\displaystyle {}^{L}D}}^q\theta$ is Liouville’s fractional derivative defined as$$\begin{array}{}\text{(38)}& {}^L{{\displaystyle D}}^q\theta t=\frac{1}{\Gamma n-q}{\displaystyle {\int }_{t_o}^t\frac{{\theta }^n\xi }{{t-\xi }^{-n+q+1}}}d\xi ,t_o⟶-\infty ,\end{array}$$where $\mathrm{\Gamma }.$ is Euler’s Gamma function and $n-1<q<n$, $n=\mathrm{1,2,3},\dots$. Liouville’s fractional derivative possesses the following fractional derivative properties:$$\begin{array}{}\text{(39a)}& {{\displaystyle {}^{L}D}}^q\cos \omega t={\omega }^q\cos \omega t+\frac{\pi q}{2},\text{(39b)}& {{\displaystyle {}^{L}D}}^q\sin \omega t={\omega }^q\sin \omega t+\frac{\pi q}{2},\text{(39c)}& {{\displaystyle {}^{L}D}}^q\exp \omega t={\omega }^q\exp \omega t.\end{array}$$

3. Stability Analysis

The importance of stability analysis in dynamical systems lies in understanding how changes in design parameters might cause bounded or unbounded responses. Two methods are proposed to study the stability for equation (32): linear analysis to obtain the transition curves of the system taking into account the fractional derivative of damping and the nonlinear analysis to obtain the general conditions of the stability on the whole, cf. [24–30].

3.1. Linear Analysis

The linear approximation of the nonlinear system by linearizing it at the fixed points is mostly useful technique to predict the geometrical nature of the equilibria and the prediction of stability domains separated by boundary curves. In general, equation (32) has the following equilibria for all parameter values:

(i) $x^{\ast },y^{\ast }=\mathrm{0,0}$, which describes a motion in which one of the body principal axes coincides with the local vertical axis

(ii) $x^{\ast },y^{\ast }=\pm \pi ,0$ which describes an inverted motion of the rigid body

Thus, the linearized homogenous form of equation (37) at a general equilibrium point ${\theta }^{\ast }$  = 0 $\pi$ or -$\pi$ for the variable $\xi t=\theta t-{\theta }^{\ast }$ with the time scale $2\tau =\omega t$; then it reads$$\begin{array}{}\text{(40)}& {\xi }^{″}+\frac{2^{q-2}\epsilon \lambda }{{\omega }^{q-2}}{{\displaystyle {}^{L}D}}^q\xi +\frac{4}{{\omega }^2}{\omega }_n^2+\gamma \sin 2\tau \xi =0.\end{array}$$

Then, consider that$$\begin{array}{}\text{(41)}& \omega =\iota {\omega }_n,{\omega }^{q-2}={\iota }^{q-2}{\omega }_n^{q-2}.\end{array}$$

Hence,$$\begin{array}{}\text{(42)}& {\xi }^{″}+\frac{2^{q-2}\epsilon \lambda }{{\iota }^{q-2}{\omega }_n^{q-2}}{{\displaystyle {}^{L}D}}^q\xi +\frac{4}{{\omega }^2}{\omega }_n^2+\gamma \sin 2\tau \xi =0,\text{(43)}& {\xi }^{″}+\frac{2^{q-2}\epsilon \lambda }{{\iota }^{q-2}{\omega }_n^{q-2}}{{\displaystyle {}^{L}D}}^q\xi +\frac{4}{{\iota }^2}+\frac{4{\omega }_n^2\gamma }{{\iota }^2}\sin 2\tau \xi =0.\end{array}$$

Then,$$\begin{array}{c}\text{(44)}& {\xi }^{″}+\zeta {{\displaystyle {}^{L}D}}^q\xi +\delta +\epsilon \sin 2\tau \xi =0.\text{where\hspace{0.17em}}\zeta =\frac{2^{q-2}\epsilon \lambda }{{\iota }^{q-2}{\omega }_n^{q-2}},\delta =\frac{4}{{\iota }^2},\epsilon =\frac{4{\omega }_n^2\gamma }{{\iota }^2},{}^{\prime }=\frac{d}{d\tau }.\end{array}$$

Due to the fact that the fractional derivative is a nonlocal operator, the choice of Liouville’s definition allows us to obtain the periodic solution based on the results provided in [31, 32]. Since the fractional derivative typically depends on the values between the lower limit $t_o$ and the point of interest (t), in this case, we are interested in looking for periodic solutions of equation (44) in $\delta ,ϵ$ domain.

However, we look for the regions of stability using the harmonic balance method for the homogenous linear simplified form of equation (44) under the action of Liouville’s fractional derivative. From Floquet theory, it is clear that this equation admits solutions of period $\pi$ or $2\pi$, cf. [33, 34].

Typically, for given values of the parameters $\delta$ and $ϵ$, all solutions of this equation should be either stable, unstable, or periodic. The $\delta ,ϵ$ plane (Ince’s diagram) is normally divided into stable and unstable regions separated by the transition curves to represent the periodic solution curves. The aim here is to present these transition curves in order to be able to assess the influence of the fractional derivative term on the system.

Now, we seek the stability domains of equation (44) separated by the transition curves using the harmonic balance method. This method was firstly presented to specify the periodic solution curves by employing operational matrices associated with a Fourier basis; and, as a result, the conditions for the existence of nontrivial periodic solution were obtained by setting Hill’s determinant to zero. So, we consider a formal solution in the following Fourier expansion:$$\begin{array}{}\text{(45)}& \xi \tau \approx A_o+{\displaystyle \sum _{\ell =1}^{\infty }{A}_{\ell }\cos \ell \tau +B_{\ell }\sin \ell \tau }.\end{array}$$

Then, inserting equation (45) into equation (44) and collecting terms, we obtain the following recurrence relations:$$\begin{array}{}\text{(46)}& \ell =0,-\frac{\epsilon }{2}{B}_2=\delta A_o,\text{(47)}& \ell =1,1-\zeta \cos \frac{\pi }{2}q{A}_1-\frac{\epsilon }{2}+\zeta \sin \frac{\pi }{2}q{B}_1-\frac{\epsilon }{2}{B}_3=\delta A_1,\text{(48)}& -\frac{\mathrm{ϵ}}{2}+\zeta \sin \frac{\pi }{2}q{A}_1+1-\zeta \cos \frac{\pi }{2}q{B}_1+\frac{\mathrm{ϵ}}{2}{A}_3=\delta B_1,\text{(49)}& \ell =2,2^2-2^q\zeta \cos \frac{\pi q}{2}{A}_2-2^q\zeta \sin \frac{\pi q}{2}{B}_2-\frac{\epsilon }{2}{B}_4=\delta A_2,\text{(50)}& -\epsilon A_o+2^q\zeta \sin \frac{\pi q}{2}{A}_2+2^2-2^q\zeta \cos \frac{\pi q}{2}{B}_2+\frac{\epsilon }{2}{A}_4=\delta B_2,\text{(51)}& \ell \ge 3,\frac{\epsilon }{2}{B}_{\ell -2}+{\ell }^2-\zeta {\ell }^q\cos \frac{\pi }{2}q{A}_{\ell }-\zeta {\ell }^q\sin \frac{\pi }{2}q{B}_{\ell }-\frac{\epsilon }{2}{B}_{\ell +2}=\delta A_{\ell },\text{(52)}& -\frac{\epsilon }{2}{A}_{\ell -2}+\zeta {\ell }^q\sin \frac{\pi }{2}q{A}_{\ell }+{\ell }^2-\zeta {\ell }^q\cos \frac{\pi }{2}q{B}_{\ell }+\frac{\epsilon }{2}{A}_{\ell +2}=\delta B_{\ell },\end{array}$$

As shown in Figure 2, the transition curves of the fractional damped homogenous linear part of equation (32) are obtained, as well as subsequently the stability regions. It is vividly shown that a slight change in the fractional order $q$ or the damping coefficient $\zeta$ affects the shape and location of the transition curves. Clearly, this effect can be characterized by the location of the lowest point on the transition curve, which specifies the minimum value of forcing amplitude $ϵ$, necessary to produce unstable domains in the linearized form. It is clear that the location of the minimum point on the transition curve due to the fractional damping moves this lowest point in a horizontal direction, but due to the damping coefficient it moves in the vertical direction until the unstable region disappears. These results are completely coincided with the fractional parametric excitation of Mathieu’s equation in [35] and Ince’s equation in [36].

[figures omitted; refer to PDF]

3.2. Nonlinear Analysis

In this part, we consider the study of stability on the whole for the governing equation (equation (32)) for both autonomous and nonautonomous forms.

Theorem 1.

Consider the autonomous form of equation (32).$$\begin{array}{}\text{(53)}& \ddot{\theta }+h\theta ,\dot{\theta }+f\theta ,0=0,,\theta \in 0,\pi .\end{array}$$

Then, the sufficient condition for a stable motion is$$\begin{array}{}\text{(54)}& \frac{\alpha }{4}>\beta +2\kappa .\end{array}$$

Proof.

The autonomous form of equation (32), when $\gamma =0$, is equivalent to$$\begin{array}{}\text{(55)}& \dot{x}=y,\dot{y}=-\epsilon \lambda +\mu x^2+\nu y^{2}y-\alpha \frac{{1-\cos x}^2}{{\sin }^{3}x}+\beta \sin x+\kappa \sin 2x.\end{array}$$

We assume that$$\begin{array}{}\text{(56)}& Hx,y=-\epsilon \lambda +\mu x^2+\nu y^{2}y-\alpha \frac{{1-\cos x}^2}{{\sin }^{3}x}+\beta \sin x+\kappa \sin 2x.\end{array}$$

By taking Lyapunov function in the form$$\begin{array}{}\text{(57)}& \mathbb{V}x,y=\frac{1}{2}{y}^2-{\displaystyle {\int }_0^x}Hx,0dx,\end{array}$$

we obtain$$\begin{array}{}\text{(58)}& \stackrel{.}{\mathbb{V}}=yHx,y-Hx,0.\end{array}$$

Hence, we have the following conditions:$$\begin{array}{}\text{(59a)}& i.Hx,0x<0,\text{for}x\ne 0,\text{(59b)}& ii.Hx,y-Hx,0y<0\text{for}y\ne 0,\text{(59c)}& iii.-{\displaystyle {\int }_0^x}Hu,0du⟶\infty \text{as}x⟶\infty .\end{array}$$

Apply the conditions on $Hx,y$; then, the conclusion holds.

Theorem 2.

Consider the following form of equation (32):$$\begin{array}{}\text{(60)}& \ddot{\theta }+\tilde{h}\theta ,\dot{\theta }\dot{\theta }+f\theta ,0=et,,\theta \in 0,\pi .\end{array}$$

Then, the sufficient conditions for a stable motion are$$\begin{array}{}\text{(61)}& \lambda ,\mu ,\nu >0,\frac{\alpha }{4}>\beta +2\kappa .\end{array}$$

Proof.

The form of equation (60) is equivalent to the system$$\begin{array}{}\text{(62)}& \dot{x}=y,\dot{y}=-\tilde{h}x,yy-fx,0+et,\end{array}$$where,$$\begin{array}{}\text{(63)}& \tilde{h}x,y=\lambda +\mu x^2+\nu y^2,fx,0=\alpha \frac{{1-\cos x}^2}{{\sin }^{3}x}-\beta \sin x-\kappa \sin 2x.\end{array}$$

By taking Lyapunov function in the form$$\begin{array}{}\text{(64)}& \mathbb{V}t,x,y=e^{-2Et}\frac{1}{2}{y}^2+{\displaystyle {\int }_0^x}fu,0du+1,\end{array}$$

we obtain$$\begin{array}{}\text{(65)}& \stackrel{.}{\mathbb{V}}<-e^{-2E\infty }\tilde{h}{y}^2,E\infty ={\displaystyle {\int }_0^{\infty }}eudu.\end{array}$$

So, every solution $xt$ exists in the future and $xt⟶0$ and $yt⟶0$ as $t⟶\infty$ if the following conditions are satisfied:$$\begin{array}{c}\text{(66a)}& i.\tilde{h}x,y\text{is\hspace{0.17em}continuous\hspace{0.17em}and\hspace{0.17em}nonnegative\hspace{0.17em}on,}\text{(66b)}& ℝ\times ℝ\text{and\hspace{0.17em}is\hspace{0.17em}bounded\hspace{0.17em}when}{x}^2+y^2\text{is\hspace{0.17em}bounded,}ii.fx,0\text{is\hspace{0.17em}continuous\hspace{0.17em}in}ℝ,xfx,0>0,\text{\hspace{0.17em}for\hspace{0.17em}}x\ne 0,\\ \text{(66c)}& \text{and}{\displaystyle {\int }_0^{x}fu,0}du⟶\infty asx⟶0,iii.et\text{is\hspace{0.17em}continuous\hspace{0.17em}on}0,\infty ,\text{and\hspace{0.17em}}Et={\displaystyle {\int }_0^{t}eu}du<\infty ,\end{array}$$

Apply the conditions; then, the conclusion holds.

Theorem 3.

The governing equation, equation (27),$$\begin{array}{}\text{(67)}& \ddot{\theta }+h\theta ,\dot{\theta }+f\theta ,t=et,,\theta \in 0,\pi ,\end{array}$$has an asymptotically stable solution if$$\begin{array}{}\text{(68)}& \lambda ,\mu ,\nu >0,\frac{\alpha }{4}>\beta +2\kappa .\end{array}$$

Proof.

By considering that$$\begin{array}{}\text{(69a)}& hx,y=\epsilon \lambda +\mu x^2+\nu y^{2}y,\text{(69b)}& fx,t=\frac{1}{2}\alpha \tan \frac{x}{2}{\sec }^2\frac{x}{2}-\beta \sin x-\kappa \sin 2x+\gamma \sin \omega t\sin x.\end{array}$$

Take that Lyapunov function $\mathbb{V}$ is$$\begin{array}{}\text{(70)}& \mathbb{V}x,y,t=\exp \sqrt{y^2+f^{\ast }x,t+k}-{\displaystyle {\int }_0^t|e}u|du-kt,\end{array}$$in the domain $k\in {ℝ}^{+}$, $0\le t\le T$, $|x|<\infty$, $|y|<\infty$, and $x^2+y^2\ge k^2$. Hence, we have$$\begin{array}{}\text{(71)}& \stackrel{.}{\mathbb{V}}\le 0.\end{array}$$

Then, it is easy to obtain the general conditions:

(i) $hx,y$ and $fx,t$ are locally Lipschitzian in $x$ and $y$, respectively, and $et$ is continuous on $ℝ$

(ii) $fx,t$ and $et$ are periodic in $t$ of period $T=2\pi /\omega$

(iii) $f^{\ast }x,t={\displaystyle {\int }_0^{x}f}u,tdu>-k$ for all $t$ and $x$ and $f_t^{\ast }/\sqrt{f^{\ast }x,t+k}$ is bounded

Apply the conditions; then, the conclusion holds.

4. Existence of Periodic Solutions

Let us assume that equation (32) has a periodic solution with period $T=2\pi /\omega$. Also, it is considered that, at the equilibrium point ${\theta }^{\ast }=0$, the following assumption holds:$$\begin{array}{}\text{(72)}& h\mathrm{0,0}+f0,t=0.\end{array}$$

Hence, the only solution of equation (72) is $\theta ={\theta }^{\ast }$. Equation (32) is rewritten as follows ($\theta =x$, $\stackrel{.}{\theta }=y$):$$\begin{array}{}\text{(73)}& \ddot{x}+fx,t=et-hx,y,\end{array}$$where the RHS represents the absorption $hx,y$ and the supply $et$ of the system.

Thus, the following theorem gives for equation (73) the general conditions for the existence of periodic solutions.

Theorem 4.

Equation (73) has at least one periodic solution of period $2\pi /\omega$ if

(i) $hx,y$ and $fx,t$ are locally Lipschitzian in $y$ and $x$, respectively, and $et$ is continuous on ${ℝ}^{+}$

(ii) $fx,t$ and $et$ are periodic in $t$ of period $T=2\pi /\omega$

(iii) $f^{\ast }x,t={\displaystyle {\int }_0^{x}f}u,tdu>-k$ for all $t$ and $f_t^{\ast }/\sqrt{f^{\ast }x,t+k}$ is bounded

(iv) $hx,yy\ge 0$

(v) There exist $s_1,s_2$, $s_1<s_2$, $h{s}_1,0+f{s}_1,t-et\ge 0$, and $h{s}_2,0+f{s}_2,t-et\le 0$

(vi) $fx,t+et\le k_o$

Proof.

For the domain $t\in 0,2\pi /\omega$, let $|x|<\infty$, $|y|<\infty$, and $x^2+y^2\ge k^2$ and $k\in {ℝ}^{+}$. Consider that there exists a continuous function such that$$\begin{array}{}\text{(74)}& |hx,y|\le h_{o}y.\end{array}$$

Assume that the following identities hold:$$\begin{array}{}\text{(75)}& {{\displaystyle \int }}^{-\infty }\frac{u}{h_{o}u+k_o}du⟶\infty ,{{\displaystyle \int }}^{\infty }\frac{u}{h_{o}u+k_o}du⟶\infty ,\end{array}$$

and then let the two functions ${\mathbb{V}}_1$ and ${\mathbb{V}}_2$ be as follows:$$\begin{array}{c}\text{(76)}& {\mathbb{V}}_1=e^{x+{\displaystyle {\int }_k^y\frac{u}{h_{o}u+k_o}}du},{\mathbb{V}}_2=e^{x+{\displaystyle {\int }_{-k}^y\frac{u}{h_{o}u+k_o}}du}.\end{array}$$

Hence, it can be concluded that there exist in the interior of two domains$$\begin{array}{}\text{(77)}& {\stackrel{.}{\mathbb{V}}}_1\ge 0,{\stackrel{.}{\mathbb{V}}}_2\le 0.\end{array}$$

Then, the solutions $xt$ and $yt$ exist and are bounded for all $t\ge 0$. Then, the conclusion holds.

Theorem 5.

In the autonomous case of the approximate governing equation $\gamma =0$, there exists a limit cycle under the following conditions:$$\begin{array}{}\text{(78)}& {\omega }_n=1,\Omega =0,\lambda <0.\end{array}$$

Proof.

By using the polar coordinates$$\begin{array}{}\text{(79)}& {\rho }_o^2=x^2+y^2,\tan \sigma =\frac{y}{x},\end{array}$$

one gets the following:$$\begin{array}{}\text{(80)}& {\dot{\rho }}_o=\frac{x\dot{x}+y\dot{y}}{{\rho }_o},\dot{\sigma }=\frac{x\dot{y}-\dot{x}y}{{\rho }_o^2}.\end{array}$$

After substitutions and using the conditions ${\omega }_n=1$ and $\mathrm{\Omega }=0$, we obtain$$\begin{array}{}\text{(81a)}& {\dot{\rho }}_o=-\epsilon {\rho }_o{\sin }^2\sigma \lambda +{\rho }_o^2\mu {\cos }^2\sigma +\nu {\sin }^2\sigma ,\text{(81b)}& \dot{\sigma }=-1+\frac{\epsilon }{2}\sin 2\sigma \lambda +r^2\mu {\cos }^2\sigma +\nu {\sin }^2\sigma ,\end{array}$$

and one particular solution is$$\begin{array}{}\text{(82)}& {\rho }_o^2=-\frac{\lambda }{\mu {\cos }^2\sigma +\nu {\sin }^2\sigma },\dot{\sigma }=-1.\end{array}$$

Corresponding to the limit cycle,$$\begin{array}{c}\text{(83)}& x=\sqrt{-\frac{\lambda }{\mu {\cos }^2\sigma +\nu {\sin }^2\sigma }}\cos \sigma ,y=\sqrt{-\frac{\lambda }{\mu {\cos }^2\sigma +\nu {\sin }^2\sigma }}\sin \sigma .\end{array}$$

It is noticed that$$\begin{array}{}\text{(84)}& {\dot{\rho }}_o>0\text{if}0<{\rho }_o<\frac{\lambda }{\sqrt{{\mu }^2+{\nu }^2}},\text{and}{\dot{\rho }}_o<0\text{if}{\rho }_o>\frac{\lambda }{\sqrt{{\mu }^2+{\nu }^2}}.\end{array}$$

This means that all paths approach the shown radius of the limit cycle ${\rho }_o$ especially when $\mu =\nu$. Then, the conclusion holds.

5. Averaging Analysis

We assume that the solution of the approximate form of governing equation, equation (35), is still given by$$\begin{array}{}\text{(85)}& xt=at\cos \omega t+bt\sin \omega t,\end{array}$$

With time-varying coefficients $at$ and $bt$, in other words, we consider equation (85) as a transformation from $xt$ to $at$ and $bt$. Hence, by imposing a third condition or a third equation independent of equation (35) and equation (85), assuming that $xt$ and $\dot{x}=yt$ have the same form as the linear case,$$\begin{array}{}\text{(86)}& yt=\omega b\cos \omega t-a\omega \sin \omega t,\end{array}$$

and the following acceleration term is obtained:$$\begin{array}{}\text{(87)}& \dot{y}t=-{\omega }^{2}a+2\omega \dot{b}\cos \omega t-b{\omega }^2+2\omega \dot{a}\sin \omega t.\end{array}$$

Consider the following condition:$$\begin{array}{}\text{(88)}& \dot{a}\cos \omega t+\dot{b}\sin \omega t=0.\end{array}$$

Then, we have to transform the system equation (85) to equation (88) with equation (35) into the two unknowns $at$ and $bt$ as follows:$$\begin{array}{}\text{(89)}& \dot{a}=-\frac{1}{2\omega }{f}_o\sin \omega t,\dot{b}=\frac{1}{2\omega }{f}_o\cos \omega t,\end{array}$$where$$\begin{array}{c}\text{(90)}& f_{o}a,b,t=\frac{3}{48}\gamma b{a}^2+b^2-\frac{1}{2}\gamma b+\cos \omega t{\omega }^2-{\omega }_n^{2}a-\epsilon \lambda \omega b-\frac{1}{2}\epsilon \mu \omega b{a}^2+b^2-\frac{3}{4}\epsilon \nu b{\omega }^3{a}^2+b^2-\frac{3}{4}a\Omega a^2+b^2+F_1\\ +\sin \omega t{\omega }^2-{\omega }_n^{2}b+\epsilon \lambda \omega a+\frac{1}{4}\epsilon \mu \omega a{a}^2+b^2+\frac{3}{4}\epsilon \nu {\omega }^{3}a{a}^2+b^2-\frac{3}{4}\Omega b{a}^2+b^2+F_2\\ +\cos 2\omega t\frac{1}{2}\Gamma b-\frac{3}{48}\Gamma b{a}^2+b^2+\frac{1}{48}\Gamma b3a^2-b^2\\ +\sin 2\omega t-\frac{1}{2}\Gamma b+\frac{3}{48}\Gamma a{a}^2+b^2-\frac{1}{48}\Gamma a{a}^2-3b^2\\ +\text{higher\hspace{0.17em}harmonics}.\end{array}$$

Consequently, by considering that the major parts of $a$ and $b$ are slowly varying functions of time, they change very little during the time period $T=2\pi /\omega$ and, in addition to the first approximation, they can be considered constant in the interval $0,T$. Hence, we average both sides of the system equation (89) for $a$ and $b$ over the interval $0,T$ and obtain$$\begin{array}{}\text{(91)}& \dot{a}=-\frac{1}{4\pi }{\displaystyle {\int }_0^T{f}_o}\sin \omega tdt,\dot{b}=\frac{1}{4\pi }{\displaystyle {\int }_0^T{f}_o}\cos \omega tdt.\end{array}$$

Since $a$ and $b$ can be considered constant in the interval $0,T$, $a$ and $b$ (and so $\dot{a}$ and $\dot{b}$) can be taken outside the integral signs in equation (91). Typically, the set of functions $1,\cos n\omega t,\sin n\omega t$ forms an orthogonal set over the interval $0,T$. Then, we have, for all $m,n=\mathrm{0,1,2},\dots$, from the inner product definition of two functions $f_1,f_2$ defined on the interval $0,T$, that$$\begin{array}{}\text{(92)}& \sin m\omega t,\cos n\omega t=0,\cos m\omega t,\cos n\omega t=\sin m\omega t,\sin n\omega t=\begin{array}{cc}0,& \text{if}\hspace{1em}m\ne n.\\ \frac{\pi }{\omega },& \text{if}\hspace{1em}m=n.\end{array}\end{array}$$

Hence, equation (91) is reduced to$$\begin{array}{c}\text{(93a)}& \dot{a}=\Xi a,b,\dot{b}=\Lambda a,b,\text{(93b)}& \Xi a,b=-\frac{1}{4\omega }{\omega }^2-{\omega }_n^{2}b+\epsilon \lambda \omega a+\frac{1}{4}\epsilon \mu \omega a{a}^2+b^2+\frac{3}{4}\epsilon \nu {\omega }^{3}a{a}^2+b^2-\frac{3}{4}\Omega b{a}^2+b^2+F_2,\\ \text{(93c)}& \Lambda a,b=\frac{1}{4\omega }{\omega }^2-{\omega }_n^{2}a-\epsilon \lambda \omega b-\frac{1}{4}\epsilon \mu \omega b{a}^2+b^2-\frac{3}{4}\epsilon \nu {\omega }^{3}b{a}^2+b^2-\frac{3}{4}\Omega a{a}^2+b^2+F_1.\end{array}$$

Consequently, the initial conditions can be given in terms of those for equation (35) by$$\begin{array}{c}\text{(94)}& a0=x0,b0=\frac{\dot{x}0}{\omega }.\end{array}$$

5.1. Relative Equilibria and Stability

According to equation (93), it can be deduced that all equilibrium solutions are given by $\dot{a}=\dot{b}=0$. In addition to the fact that the other paths corresponding to solutions of equation (35) are nonperiodic in general, let $a_o,b_o$ represent any of the three fixed points, $a_1,b_1,a_2,b_2\text{\hspace{0.17em}and\hspace{0.17em}}{a}_3,b_3$, of equation (93). Hence, the linear analysis in the neighborhood of point $a_o,b_o$ is assumed to be$$\begin{array}{}\text{(95)}& a=a_o+{\eta }_1,b=b_o+{\eta }_2,\end{array}$$

and then the corresponding linearized system is$$\begin{array}{}\text{(96a)}& {\dot{\eta }}_1={\Xi }_1{a}_o,b_o{\eta }_1+{\Xi }_2{a}_o,b_o{\eta }_2,\text{(96b)}& {\dot{\eta }}_2={\Lambda }_1{a}_o,b_o{\eta }_1+{\Lambda }_2{a}_o,b_o{\eta }_2,\end{array}$$where$$\begin{array}{}\text{(97a)}& {\Xi }_1{a}_o,b_o=-\frac{1}{4}\epsilon \lambda +\frac{1}{4}\epsilon \mu 3a_o^2+b_o^2+\frac{3}{4}\epsilon \nu {\omega }^{2}3a_o^2+b_o^2-\frac{3}{2\omega }\Omega a_o{b}_o,\text{(97b)}& {\Xi }_2{a}_o,b_o=-\frac{1}{4\omega }{\omega }^2-{\omega }_n^2+\frac{1}{2}\mathrm{ϵ}\mu \omega a_o{b}_o+\frac{3}{2}\epsilon \nu {\omega }^3{a}_o{b}_o-\frac{3}{4}\Omega a_o^2+3b_o^2,\text{(97c)}& {\Lambda }_1{a}_o,b_o=\frac{1}{4\omega }{\omega }^2-{\omega }_n^2-\frac{1}{2}\epsilon \mu \omega a_o{b}_o-\frac{3}{2}\epsilon \nu {\omega }^3{a}_o{b}_o-\frac{3}{4}\Omega 3a_o^2+b_o^2,\text{(97d)}& {\Lambda }_2{a}_o,b_o=\frac{1}{4\omega }-\epsilon \lambda \omega -\frac{1}{4}\epsilon \mu \omega a_o^2+3b_o^2-\frac{3}{4}\epsilon \nu {\omega }^3{a}_o^2+3b_o^2-\frac{3}{2}\Omega a_o{b}_o.\end{array}$$

The general conditions of stability for every stable fixed point are$$\begin{array}{}\text{(98)}& {\Xi }_1{\Lambda }_2>{\Xi }_2{\Lambda }_1,{\Xi }_1+{\Lambda }_2\le 0.\end{array}$$

Moreover, the boundary between nodes and spirals can be given by the following identity:$$\begin{array}{}\text{(99)}& {{\Xi }_1+{\Lambda }_2}^2=4{\Xi }_1{\Lambda }_2-{\Xi }_2{\Lambda }_1.\end{array}$$

The phase diagram computed for a particular case $\epsilon =0$ is shown in Figure 3 for $F=0$ and $F=0.005$, respectively. In Figure 3(a), point $P_1$ is a center and represents a stable oscillation, but in Figure 3(b) there exist three equilibrium points: $P_1$ and $P_2$ are centers and $P_3$ is a saddle representing an unstable oscillation.

[figures omitted; refer to PDF]

For the particular case of $\epsilon =0.005$ and $\lambda$  =  $\mu$  =  $\nu$  = 1, the phase plane is shown in Figure 4 for $F=0$ and $F=0.005$, respectively. In Figure 4(a), point $P_1$ is a spiral point and represents a stable oscillation, but in Figure 4(b) there exist three equilibrium points: $P_1$ and $P_2$ are stable spiral points and $P_3$ is a saddle representing an unstable oscillation.

[figures omitted; refer to PDF]

The phase diagram computed for a particular case where $\epsilon =0.005$, $\lambda$  = −1, and $\mu$  =  $\nu$  = 1 is shown in Figure 5 for $F=0$ and $F=0.005$, respectively. In Figure 5(a), point $P_1$ is an unstable point as it bifurcates to generate a stable limit cycle, but in Figure 5(b) there exist three equilibrium points: $P_1$ and $P_2$ are centers and $P_3$ is a saddle representing an unstable oscillation.

[figures omitted; refer to PDF]

5.2. The Jump Phenomenon

From equation (93), three are three different amplitudes at the equilibria. Consequently, due to the existence of the unstable node, the amplitude surface $\omega ,F,\rho$ has unstable oscillations. Then, the surface exhibits a fold catastrophe or a cusp which leads to the jump phenomenon, cf. [14]. Illustratively, the cubic amplitude equation in $\mathrm{\varrho }$ representing the surface of equilibrium points $\mathrm{\Xi }a,b=\mathrm{\Lambda }a,b=0$ is given by$$\begin{array}{}\text{(100)}& F^2=\mathrm{\varrho }{{\omega }^2-{\omega }_n^2-\frac{3}{4}\Omega \mathrm{\varrho }}^2+{\epsilon \lambda \omega +\frac{1}{4}\epsilon \omega \mu +3\nu {\omega }^2\mathrm{\varrho }}^2,\end{array}$$where$$\begin{array}{}\text{(101)}& \mathrm{\varrho }={\rho }^2=a^2+b^2.\end{array}$$

Hence, there may be as many as three positive values of $\mathrm{\varrho }$ satisfying the equation of amplitude. This explicates that, for certain ranges of the problem parameters, there might be three distinct periodic solutions of the governing equation. Also, the responses are stable when $dF/d\rho >0$ and unstable when $dF/d\rho <0$.

Consider the following:$$\begin{array}{c}\text{(102)}& G\mathrm{\varrho }=\frac{9}{16}{\Omega }^2+\frac{1}{16}{\epsilon }^2{\omega }^2\mu +3\nu {\omega }^2{\mathrm{\varrho }}^3+\frac{1}{2}{\epsilon }^2\lambda {\omega }^2\mu +3\nu {\omega }^2-\frac{3}{4}\Omega {\omega }^2-{\omega }_n^2{\mathrm{\varrho }}^2+{{\omega }^2-{\omega }_n^2}^2+{\epsilon }^2{\lambda }^2{\omega }^2\mathrm{\varrho }.\end{array}$$

Since $G0=0$ and $G\mathrm{\varrho }⟶\infty$ as $\mathrm{\varrho }⟶\infty$, the amplitude equation must have at least one positive root. For some parameters’ values, the amplitude equation has three different real roots corresponding to its cubic form of $\mathrm{\varrho }$ if$$\begin{array}{c}\text{(103)}& \frac{dG}{d\mathrm{\varrho }}=\frac{3}{16}9{\Omega }^2+{\epsilon }^2{\omega }^2\mu +3\nu {\omega }^2{\mathrm{\varrho }}^2+{\epsilon }^2\lambda {\omega }^2\mu +3\nu {\omega }^2-\frac{3}{2}\Omega {\omega }^2-{\omega }_n^2\mathrm{\varrho }+{{\omega }^2-{\omega }_n^2}^2+{\epsilon }^2{\lambda }^2{\omega }^2=0\end{array}$$

has two distinct solutions; thus,$$\begin{array}{c}\text{(104)}& {\mathrm{\varrho }}_1,{\mathrm{\varrho }}_2=\frac{1}{3/89{\Omega }^2+{\epsilon }^2{\omega }^2\mu +3\nu {\omega }^2}-{\epsilon }^2\lambda {\omega }^2\mu +3\nu {\omega }^2-\frac{3}{2}\Omega {\omega }^2-{\omega }_n^2\pm \sqrt{{{\epsilon }^2\lambda {\omega }^2\mu +3\nu {\omega }^2-\frac{3}{2}\Omega {\omega }^2-{\omega }_n^2}^2-\frac{3}{4}9{\Omega }^2+{\epsilon }^2{\omega }^2\mu +3\nu {\omega }^2{{\omega }^2-{\omega }_n^2}^2+{\epsilon }^2{\lambda }^2{\omega }^2}.\end{array}$$

So the solution of the amplitude equation is interested in ${\mathrm{\varrho }}_1,{\mathrm{\varrho }}_2>0$; then,$$\begin{array}{}\text{(105)}& {\mathrm{\varrho }}_1+{\mathrm{\varrho }}_2>0,{\mathrm{\varrho }}_1{\mathrm{\varrho }}_2>0.\end{array}$$

Thus, the amplitude equation has three distinct real roots if $F$ is controlled by the overlap region.$$\begin{array}{}\text{(106)}& \sqrt{G_2}<F<\sqrt{G_1},G_1=G{\mathrm{\varrho }}_1,G_2=G{\mathrm{\varrho }}_2.\end{array}$$

To explain, we take the case of undamped motion $\epsilon =0$; then, we have$$\begin{array}{}\text{(107)}& {\mathrm{\varrho }}_1=\frac{2}{3}\frac{{\omega }^2-{\omega }_n^2}{\Omega },{\mathrm{\varrho }}_2=\frac{10}{9}\frac{{\omega }^2-{\omega }_n^2}{\Omega },\end{array}$$

and then, as shown in Figure 6, the overlap region is controlled by$$\begin{array}{}\text{(108)}& \sqrt{\frac{155}{162}}\frac{{{\omega }^2-{\omega }_n^2}^{3/2}}{\sqrt{\Omega }}<F<\frac{1}{\sqrt{6}}\frac{{{\omega }^2-{\omega }_n^2}^{3/2}}{\sqrt{\Omega }}.\end{array}$$

[figure omitted; refer to PDF]

Hence, the cusps are exactly determined at$$\begin{array}{}\text{(109)}& \omega ={\omega }_n.\end{array}$$

In Figure 7, it is clear that the surface under the fold gets along with unstable oscillations which will not be obtainable or accessible using equation (100). The surface displays what is known as a fold catastrophe which produces the jump phenomenon controlled by parameter $\epsilon$ at the values of 0, 0.1, and 0.25. It is clear that, with the increase of parameter $\epsilon$ at the constant values of damping parameters $\lambda$, $\mu$, and $\nu$, the jump is going to disappear.

[figures omitted; refer to PDF]

As shown in Figure 7, it can be predicted that, for certain values of $\omega$ and $F$ at the lower values of $\epsilon \approx 0$, there are three responses: the oscillations with greatest and smallest amplitudes are stable but the middle one is unstable. But when $\epsilon$ is increased, then there exists only a single stable response.

The consequences of this fold are shown illustratively using Figure 8 and 9. So, consider that $F$ is held constant at a value which intersects the fold, and the applied frequency $\omega$ is varied in such a way that it crosses the edges of the fold. At the small values of $\epsilon \approx 0$, as $\omega$ increases beyond, after some irregular motion, the amplitude jumps at the critical value; henceforth, it follows the smooth curve. On the way back, the response point moves along the upper curve. Then, here it must go down to the lower curve and then back to the upper one crossing the jumping point of the amplitude.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

6. Homoclinic Bifurcation

The existence of homoclinic bifurcation with the support of the other problem parameters, particularly $\alpha$ and $\beta$, is considered. The homoclinic bifurcation problem in equation (35) can be considered at certain values of the governing problem parameters. In general, Melnikov’s method can be used to predict the critical value to the homoclinic bifurcation, but in most of the cases the perturbation method has a less accuracy when the perturbation parameter increases further, cf. [37, 38]. Some critical values are detected for some special perturbed systems as shown in Table 1 at some specific values of the problem parameters, cf. [2, 39, 40].

Table 1

Homoclinic critical relations for some special type systems.

The system	The homoclinic solution	The critical relation
(i) $\ddot{x}+\epsilon \lambda \dot{x}+\sin x=0$	$x_{0}t=\pm 2{\sin }^{-1}\text{tanh}t-t_0$	$\lambda \le 0$
(ii) $\ddot{x}+\epsilon \lambda \dot{x}-\sin x=0$	$x_{0}t=\pm 2{\tan }^{-1}\text{sinh}t-t_0$	$\lambda \le 0$
(iii) $\ddot{x}+\epsilon \lambda \dot{x}-x+x^3=0$	$x_{0}t=\pm \sqrt{2}\text{sech}t-t_0$	$\lambda \le 0$
(iv) $\ddot{x}+\epsilon \lambda \dot{x}+x-x^3=0$	$x_{0}t=\pm \tanh 1/\sqrt{2}t-t_0$	$\lambda \le 0$
(v) $\ddot{x}+\epsilon \lambda \dot{x}-x+x^5=0$	$x_{0}t=\pm \sqrt[4]{3}\sqrt{\mathrm{sech}2t-t_0}$	$\lambda \le 0$
(vi) $\ddot{x}+\epsilon \lambda \dot{x}-x+3x^5=0$	$x_{0}t=\pm \sqrt{\mathrm{sech}t-t_0}$	$\lambda \le 0$
(vii) $\ddot{x}+\epsilon \mu x^2-\lambda \dot{x}-x+x^3=0$	$x_{0}t=\pm \sqrt{2}\text{sech}t-t_0$	$\lambda \ge 4/5\mu$
(viii) $\ddot{x}+\epsilon \mu x^2-\lambda \dot{x}+x-x^3=0$	$x_{0}t=\pm \tanh 1/\sqrt{2}t-t_0$	$\lambda \ge 1/3\mu$
(ix) $\ddot{x}+\epsilon \lambda -\mu x^2\dot{x}-x+x^5=0$	$x_{0}t=\pm \sqrt[4]{3}\sqrt{\mathrm{sech}2t-t_0}$	$\lambda \le 4\sqrt{3}/3\pi \mu$
(x) $\ddot{x}+\epsilon \mu x^2-\lambda \dot{x}-x+3x^5=0$	$x_{0}t=\pm \sqrt{\text{sech}t-t_0}$	$\lambda \ge 4/3\pi \mu$
(xi) $\ddot{x}+\epsilon \lambda +\mu x^2+\nu {\dot{x}}^2\dot{x}-x+x^3=0$	$x_{0}t=\pm \sqrt{2}\text{sech}t-t_0$	$\lambda \ge 28\hspace{0.17em}\mu +12\nu /35$
(xii) $\ddot{x}+\epsilon \lambda \dot{x}-x+x^3=\epsilon e\cos \omega t$	$x_{0}t=\pm \sqrt{2}\text{sech}t-t_0$	$e\ge 2\sqrt{2}\lambda /3\pi \omega \text{cosh}\pi /2\omega$
(xiii) $\ddot{x}+\epsilon \lambda \dot{x}-x+x^3=\epsilon ex\cos \omega t$	$x_{0}t=\pm \sqrt{2}\text{sech}t-t_0$	$e\ge 4\lambda /3\pi {\omega }^2\text{sinh}\pi /2\omega$
(xiv) $\ddot{x}+\epsilon \lambda \dot{x}-x+x^3=\epsilon e1-x^2\cos \omega t$	$x_{0}t=\pm \sqrt{2}\text{sech}t-t_0$	$e\ge 2\sqrt{2}\lambda /\pi \omega 2-{\omega }^2\text{cosh}\frac{\pi }{2}\omega$
(xv) $\ddot{x}+\epsilon \lambda \dot{x}+\sin x=\epsilon e\sin t\sin x$	$x_{0}t=\pm 2{\sin }^{-1}\text{tanh}t-t_0$	$e\ge 4\lambda \text{sinh}\pi /2$
(xvi) $\ddot{x}+\epsilon \lambda \dot{x}-\sin x=\epsilon e\sin \omega t\sin x$	$x_{0}t=\pm 2{\tan }^{-1}\text{sinh}t-t_0$	$|e|\ge 4\lambda /\pi {\omega }^2\text{sinh}\pi \omega /2$

To verify the theoretical results, system (xi) is considered as depicted in Table 1. System (xi) has the saddle point at $\mathrm{0,0}$ with an associated homoclinic trajectory with solution $x_{0}t$ at $\epsilon =0$. So, using Melnikov’s function, the homoclinic bifurcation occurs approximately where$$\begin{array}{}\text{(110)}& {\lambda }_c=\frac{28\mu +12\nu }{35}.\end{array}$$

From Figure 10, it is shown that there exist three equilibrium points which become stable spirals if $\lambda >{\lambda }_c$ and unstable if $\lambda <{\lambda }_c$. As $\lambda$ touches the critical value ${\lambda }_c=0.81714$ the path changes from heteroclinic spiral saddle connection $\lambda >{\lambda }_c$ to a homoclinic saddle connection at the critical value to heteroclinic saddle connection $\lambda <{\lambda }_c$. This transition is due to the fact that parameter $\lambda$ generates a homoclinic bifurcation at ${\lambda }_c$.

[figures omitted; refer to PDF]

7. Route to Chaos

Considering equation (33), some special cases are looked at to allow the parameter space $\alpha ,\beta ,\kappa ,\gamma$ to vary where the chaotic phenomenon might occur. By using the numerical configuration, the critical value of $\gamma$ at $\lambda =1$, $\mu =1$, $\omega =1$, $\alpha =1$, $\beta =1$, and $\nu$  =  $\kappa$  =  $F$  = 0 might be roughly obtained.

A simple numerical search reveals a stable periodic solution which is shown in Figure 11. Hence, if $\gamma$ is increased further, the solution bifurcates at a specific value of $\gamma$ into a solution which wraps around the origin twice. Further increase of $\gamma$ leads to further bifurcation that takes place giving more periods of solution as vividly shown. The process continues with further period-doubling to have a chaotic attracting set which is known as a strange attractor. At each bifurcation, an unstable periodic orbit will be always resided in, so that the strange attractor must contain a bounded set of unstable limit cycles. For its final stage, the periodic solution becomes closer to the destabilizing effect of the homoclinic paths which appears as $\gamma$ increases.

[figures omitted; refer to PDF]

This analysis is confirmed by computing Poincare’s map and Lyapunov exponents to detect period-doubling and chaos. From Figure 12, Poincare’s map becomes bounded sets of returns without any obvious repetitions implying the existence of strange attractor. In this case, at least one of Lyapunov exponents should be positive and the total sum of Lyapunov exponents should be negative, which are typically satisfied in Figure 12. The solutions’ figures as depicted show wandering solutions of an irregular oscillations type without any uniform pattern. Then, it is inferred that such solutions are said to display a chaotic behavior.

[figures omitted; refer to PDF]

8. Conclusion

This work is concerned with modeling of a motion of an axisymmetric rigid body under an external sinusoidal excitation and a Newtonian attractive field with a general van der Pol damping based on the nutational angle of Euler’s interpretation. Two approaches are established to study the stability analysis: The first presents the linear analysis to obtain the transition curves of the system considering the fractional derivative of the damping term using Liouville’s definition. The effect of the fractional order and the system parameters on the stability of the solution are presented using Ince’s diagram. The second approach discussed the nonlinear analysis of the system. Sufficient conditions that guarantee stable motion of the system are established for both autonomous and nonautonomous cases. Consequently, a sufficient set of conditions that ensure the asymptotic stability of the solution are given. Also, other sets of conditions that ensure the existence of at least one periodic solution and a limit cycle are derived. In addition, averaging analysis to clear the relative equilibria with their stability and the description of the jump phenomenon are shown. The existence of homoclinic bifurcation based on system parameters using Melnikov’s function is illustrated. The transition from homoclinic bifurcation to chaos is identified analytically and numerically using phase trajectories, Poincare’s map, and Lyapunov exponents.

Disclosure

This research was performed as part of the employment of the authors at University of Tanta and University of Kafrelsheikh, Egypt.

Authors’ Contributions

All authors have read and agreed to the published version of the manuscript.

References

[1] F. C. Nelson, "Vibration isolation: a review, I Sinusoidal and random excitations," Shock and Vibration, vol. 1 no. 5, pp. 485-493, DOI: 10.1155/1994/576705, 1994.

[2] M. El-Borhamy, N. Mosalam, "On the existence of periodic solution and the transition to chaos of Rayleigh-Duffing equation with application of gyro dynamic," Appl. Math. Nonl. Sci, vol. 5 no. 1, pp. 43-58, DOI: 10.2478/amns.2020.1.00010, 2020.

[3] J. M. Longuski, "Real solutions for the attitude motion of a self-excited rigid body," Acta Astronautica, vol. 25 no. 3, pp. 131-139, DOI: 10.1016/0094-5765(91)90140-z, 1991.

[4] J. M. Longuski, P. Tsiotras, "Analytical solutions for a spinning rigid body subject to time-varying body-fixed torques, Part I: constant axial torque," Journal of Applied Mechanics, vol. 60 no. 4, pp. 970-975, DOI: 10.1115/1.2901010, 1993.

[5] P. Tsiotras, J. M. Longuski, "Analytic solutions for a spinning rigid body subject to time-varying body-fixed torques, Part II: time-varying axial torque," Journal of Applied Mechanics, vol. 60 no. 4, pp. 976-981, DOI: 10.1115/1.2901011, 1993.

[6] J. A. Vera, "The gyrostat with a fixed point in a Newtonian force field: relative equilibria and stability," Journal of Mathematical Analysis and Applications, vol. 401 no. 2, pp. 836-849, DOI: 10.1016/j.jmaa.2012.11.003, 2013.

[7] L. D. Akulenko, Y. S. Zinkevich, T. A. Kozachenko, D. D. Leshchenko, "The evolution of the motions of a rigid body close to the Lagrange case under the action of an unsteady torque," Journal of Applied Mathematics and Mechanics, vol. 81 no. 2, pp. 79-84, DOI: 10.1016/j.jappmathmech.2017.08.001, 2017.

[8] M. El-Borhamy, Perturbed Rotational Motion of a Rigid Body, 2005. M.Sc. thesis

[9] M. El-Borhamy, "On the existence of new integrable cases for Euler-Poisson equations in Newtonian fields," Alexandria Engineering Journal, vol. 58 no. 2, pp. 733-744, DOI: 10.1016/j.aej.2019.06.004, 2019.

[10] Z. Wei, Y. XiaoLi, "Transverse nonlinear vibrations of a circular spinning disk with a varying rotating speed," Science China, vol. 53 no. 8, pp. 1536-1553, 2010.

[11] M. El-Borhamy, "Chaos transition of the generalized fractional duffing oscillator with a generalized time delayed position feedback," Nonlinear Dynamics, vol. 101 no. 4, pp. 2471-2487, DOI: 10.1007/s11071-020-05840-y, 2020.

[12] M. El-Borhamy, T. Medhat, M. E. Ali, "Chaos prediction in fractional delayed energy-based models of capital accumulation," Complexity, vol. 2021,DOI: 10.1155/2021/8751963, 2021.

[13] Z. Wei, P. Yu, W. Zhang, M. Yao, "Study of hidden attractors, multiple limit cycles from Hopf bifurcation and boundedness of motion in the generalized hyperchaotic Rabinovich system," Nonlinear Dynamics, vol. 82 no. 1-2, pp. 131-141, DOI: 10.1007/s11071-015-2144-8, 2015.

[14] D. W. Jordan, P. Smith, Nonlinear Ordinary Differential Equations: An Introduction for Scientists and Engineers, 2007.

[15] Z. Wei, W. Zhang, "Hidden hyperchaotic attractors in a modified Lorenz-Stenflo system with only one stable equilibrium," International Journal of Bifurcation and Chaos, vol. 24 no. 10,DOI: 10.1142/s0218127414501272, 2014.

[16] M. Yao, W. Zhang, J. W. Zu, "Multi-pulse chaotic dynamics in non-planar motion of parametrically excited viscoelastic moving belt," Journal of Sound and Vibration, vol. 331 no. 11,DOI: 10.1016/j.jsv.2012.01.027, 2012.

[17] J. H. Zhang, W. Zhang, "Multi-pulse chaotic dynamics of non-autonomous nonlinear system for a honeycomb sandwich plate," Acta Mechanica, vol. 223 no. 5, pp. 1047-1066, DOI: 10.1007/s00707-012-0618-z, 2012.

[18] W. Zhang, S. B. Li, "Resonant chaotic motions of a buckled rectangular thin plate with parametrically and externally excitations," Nonlinear Dynamics, vol. 62 no. 3, pp. 673-686, DOI: 10.1007/s11071-010-9753-z, 2010.

[19] A. Wolf, J. B. Swift, H. L. Swinney, J. A. Vastano, "Determining Lyapunov exponents from a time series," Physica D: Nonlinear Phenomena, vol. 16 no. 3, pp. 285-317, DOI: 10.1016/0167-2789(85)90011-9, 1985.

[20] J.-P. Eckmann, D. Ruelle, "Ergodic theory of chaos and strange attractors," Reviews of Modern Physics, vol. 57 no. 3, pp. 617-656, DOI: 10.1103/revmodphys.57.617, 1985.

[21] J.-P. Eckmann, S. O. Kamphorst, D. Ruelle, S. Ciliberto, "Liapunov exponents from time series," Physical Review A, vol. 34 no. 6, pp. 4971-4979, DOI: 10.1103/physreva.34.4971, 1986.

[22] M. Lara, "Complex variables approach to the short-axis-mode rotation of a rigid body," Applied Mathematics and Nonlinear Sciences, vol. 3 no. 2, pp. 537-552, DOI: 10.2478/amns.2018.2.00042, 2018.

[23] F. E. Udwadia, B. Han, "Synchronization of multiple chaotic gyroscopes using the fundamental equation of mechanics," Journal of Applied Mechanics, vol. 75 no. 2,DOI: 10.1115/1.2793132, 2008.

[24] J. G. Alaba, B. S. Ogundare, "On stability and boundedness properties of solutions of certain second order non-autonomous nonlinear ordinary differential equation," Kragujevac Journal of Mathematics, vol. 39 no. 2, pp. 255-266, DOI: 10.5937/kgjmath1502255a, 2015.

[25] E. A. Coddington, N. Levinson, Theory of Ordinary Differential Equations, 1955.

[26] D. R. Merkin, Introduction to the Theory of Stability, 1997.

[27] G. Morosanu, C. Vladimirescu, "Stability for a damped nonlinear oscillator," Nonlinear Analysis, vol. 60 no. 2, pp. 303-310, DOI: 10.1016/s0362-546x(04)00345-1, 2005.

[28] R. Ortega, "Periodic solutions of a Newtonian equation: stability by the third approximation," Journal of Differential Equations, vol. 128 no. 2, pp. 491-518, DOI: 10.1006/jdeq.1996.0103, 1996.

[29] R. Seydel, "New methods for calculating the stability of periodic solutions," Computers & Mathematics with Applications, vol. 14 no. 7, pp. 505-510, DOI: 10.1016/0898-1221(87)90045-9, 1987.

[30] H. K. Wilson, Ordinary differential equations, vol. III, 1970.

[31] I. Podlubny, Fractional Differential Equations, 1999.

[32] M. S. Tavazoei, "A note on fractional-order derivatives of periodic functions," Automatica, vol. 46 no. 5, pp. 945-948, DOI: 10.1016/j.automatica.2010.02.023, 2010.

[33] M. El-Borhamy, E. E. M. Rashad, I. Sobhy, "Floquet analysis of linear dynamic RLC circuits," Open Physics, vol. 18 no. 1, pp. 264-277, DOI: 10.1515/phys-2020-0136, 2020.

[34] M. El-Borhamy, E. M. Rashad, I. Sobhy, M. El-sayed, "Modeling and semi-analytic stability analysis for dynamics of AC machines," Mathematics, vol. 9 no. 1-13, 2021.

[35] R. H. Rand, S. M. Sah, M. K. Suchorsky, "Fractional Mathieu equation," Communications in Nonlinear Science and Numerical Simulation, vol. 15 no. 11, pp. 3254-3262, DOI: 10.1016/j.cnsns.2009.12.009, 2010.

[36] A. Parra-Hinojosa, J. C. Gutiérrez-Vega, "Fractional Ince equation with a Riemann-Liouville fractional derivative," Applied Mathematics and Computation, vol. 219 no. 22, pp. 10695-10705, DOI: 10.1016/j.amc.2013.04.044, 2013.

[37] M. Sun, W. Zhang, J. E. Chen, "Subharmonic Melnikov method of six-dimensional nonlinear systems and application to a laminated composite piezoelectric rectangular plate," Nonlinear Dynamics, vol. 75 no. 1-2, pp. 289-310, DOI: 10.1007/s11071-013-1066-6, 2014.

[38] W. Zhang, M. H. Yao, J. H. Zhang, "Using the extended Melnikov method to study the multi-pulse global bifurcations and chaos of a cantilever beam," Journal of Sound and Vibration, vol. 319 no. 1-2, pp. 541-569, DOI: 10.1016/j.jsv.2008.06.015, 2009.

[39] V. S. Aslanov, "Chaotic behavior of the biharmonic dynamics system," International Journal of Mathematics and Mathematical Sciences, vol. 2009, 2009.

[40] Z.-M. Ge, H.-K. Chen, H.-H. Chen, "The regular and chaotic motions of a symmetric heavy gyroscope with harmonic excitation," Journal of Sound and Vibration, vol. 198 no. 2, pp. 131-147, DOI: 10.1006/jsvi.1996.0561, 1996.