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

Authors: Mohamed El-Borhamy; Zeinab El-Sheikh; Manal E. Ali
Journal: Mathematical Problems in Engineering
 (2022)

---

## 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.


---

## Body

## 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.(1)hx,y=ελ+μx2+νy2y,y=x˙.where the parameters λ, μ, and ν are real constants; 0≤ε≪1.Using the dimensionless small parameterε, 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: Section2 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 massm is mounted on a vibrating base of a sinusoidal excitation frequency ω 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 I^, J^, and 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 i^, j^, and k^, respectively.Figure 1
A rigid body mounted on a vibrating base in a Newtonian field of attraction.Hence, the translational displacement of the center of mass (C.G.) with respect to the fixed axesX,Y,Z reads(2)rc=xcI^+ycJ^+zcK^+dt,dt=d1tI^+d2tJ^+d3tK^.here, dt is the excitation displacement vector. Then,(3)rc=d1t+rcsinϕsinθI^+d2t−rcsinϕcosθJ^+d3t+rccosθK^.Then, we obtain the velocity of the center of mass:(4)vc=drcdt.We assume that the principal axes rotate in space with the same angular velocity of the bodyϖ defined by (in terms of Euler’s angles)(5)ϖt=ψ˙sinϕsinθ+θ.cosϕi^+ψ˙cosϕsinθ−θ.sinϕj^+ϕ˙+ψ˙cosθk^,where θt, ϕt, and ψt are the nutation, spin, and precession angles, respectively, and ⋅ refers to the derivative with respect to time.LetK^t be the Poisson vector with the following direction cosines with respect to the moving axes:(6)K^t=γ1ti^+γ2tj^+γ3tk^.The relations of the direction cosines with Euler’s angles read(7)γ1=sinθsinϕ,γ2=sinθcosϕ,γ3=cosθ.The potential energy of a Newtonian force field is considered by the Newtonian attraction of a pointP, of mass M, on the body mass m. If P is located on the negative part of the axis OZ, then we have(8)R^=d3t−RK^,do≪Rg=GMR2.here, 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(9)UN=U0+U1+U2+U3+…where(10)U0=mgd3t−R,(11)U1=mgrc⋅K^,(12)U2=−3g2RAγ12+Bγ22+Cγ32,(13)U3=−3g2R2Jxxx+Jxyy+Jxzzγ1+Jxxy+Jyyy+Jyzzγ2+Jxxz+Jyyz+Jzzzγ3+5g2R2Jxxxγ13+Jyyyγ23+Jzzzγ33+3Jxxyγ12γ2+3Jxxzγ12γ3+3Jyyxγ22γ1+3Jyyzγ22γ3+3Jxzzγ1γ32+3Jyzzγ2γ32+6Jxyzγ1γ2γ3.here, A, B, and C are the principal moments of inertia along the moving axes ox, oy, and oz, respectively, and J… 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(14)U=UN.Consider the following approximation:(15)U≈U0+U1+U2,and then(16)U=−mgR+mgxcsinθsinϕ+ycsinθcosϕ+zccosθ−3g2RAsin2θsin2ϕ+Bsin2θcos2ϕ+Ccos2θ+mgd3t.The kinetic energy reads(17)K.E.=12mvc2+12ϖIcϖ.Then, the system admits the following LagrangianL on the Lie group SO (3):(18)L=12mvc2+12ϖIcϖ−U,where Ic=I1,I2,I3 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:(19)I1=I2,xc=yc=0,d1t=d2t=0,d3t=dt=dosinωt,(20)A=B=I,I3=C,I=I1+mzc2.Then,(21)U=−mgR+mgzccosθ−3g2R2Isin2θ+I3cos2θ+Mgdt,where dt=d0sinωt is the time-varying amplitude of the vertical support motion along the vertical fixed axis OZ with constant amplitude do and forced frequency ω.Using Routhian’s functionℜ,(22)ℜ=L−pψψ˙pψ,pϕ,θ−pϕϕ˙pψ,pϕ,θ,where pϕ and pψ are the corresponding generalized momenta to the cyclic coordinates ϕ and ψ, respectively, cf. [8, 9, 22, 23]:(23)pϕ=I3ϕ˙+ψ˙cosθ=const.,pψ=Iψ˙sin2θ+pϕcosθ=const.,and then the equation of motion reads(24)ddt∂ℜ∂θ.−∂ℜ∂θ=−Dθ,θ.+F˜1cosωt+F˜2sinωt.The right terms in equation (24) represent the damping term D and the harmonic forcing terms of amplitudes F˜1 and F˜2>0 with circular frequency (ω) applied to the body on such nonlinear spring.The damping termD is assumed to have the following generalized van der Pol’s damping form:(25)Dθ,θ.=ελ˜+μ˜θ2+ν˜θ.2θ.,λ˜,μ˜,ν˜∈ℝ,0≤ε<1.Using equation (24) and if pψ=pϕ=po, the governing equation becomes(26)θ¨+1I+po21−cosθ2I2sin3θ−mgzcIsinθ−3g2RI4I−2I3sinθcosθ−mzcIsinθd¨t=F˜1Icosωt+F˜2Isinωt.Equation (26) can be generalized as follows:(27)θ¨+hθ,θ.+fθ,t=et,where(28)fθ,t=po2I21−cosθ2sin3θ−mgzcIsinθ−3g2RI2I−I3sin2θ+1Iω2mzcdosinωtsinθ,(29)hθ,θ.=ελ+μθ2+νθ.2θ.,et=F1cosωt+F2sinωt.Consider the following notations:(30)λ=λ˜I,μ=μ˜I,ν=ν˜I,F1=F˜1I,F2=F˜2I,F=F12+F22.For simplicity, let us assume that(31)α=po2I2,β=mgzcI,γ=1Iω2mzcdo,κ=3g2RI2I−I3.Then, equation (27) reads(32)θ¨+ελ+μθ2+νθ.2θ.+α1−cosθ2sin3θ−βsinθ−κsin2θ+γsinωtsinθ=et.Ifθ=x and θ.=y, then equation (32) is converted to the following system:(33a)x˙=y,(33b)y˙=−ελ+μx2+νy2y−α1−cosx2sin3x+βsinx+κsin2x−γsinωtsinx+et.Use the following approximations:(34)sinθ=θ−θ36,1−cosθ2sin3θ=θ4+θ312,ωn2=α4−β−2κ,Ω=α12+β6+4κ3,and then equation (32) reads(35)θ¨+ελ+μθ2+νθ.2θ.+ωn2+γsinωtθ+Ω−γ6sinωtθ3=et,and the approximate system of equation (35) reads(36a)x˙=y,(36b)y˙=−ελ+μx2+νy2y−ωn2+γsinωtx−Ω−γ6sinωtx3+et.Ifθ≪1, then equation (35) has the following linear form under a proposed fractional damping:(37)θ¨+ελLDqθ+ωn2+γsinωtθ=et,0<q<2,where DLqθ is Liouville’s fractional derivative defined as(38)DqLθt=1Γn−q∫totθnξt−ξ−n+q+1dξ,to⟶−∞,where Γ. is Euler’s Gamma function and n−1<q<n, n=1,2,3,…. Liouville’s fractional derivative possesses the following fractional derivative properties:(39a)DLqcosωt=ωqcosωt+πq2,(39b)DLqsinωt=ωqsinωt+πq2,(39c)DLqexpωt=ωqexpωt.


## 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∗,y∗=0,0, which describes a motion in which one of the body principal axes coincides with the local vertical axis(ii)
x∗,y∗=±π,0 which describes an inverted motion of the rigid bodyThus, the linearized homogenous form of equation (37) at a general equilibrium point θ∗  = 0 π or -π for the variable ξt=θt−θ∗ with the time scale 2τ=ωt; then it reads(40)ξ″+2q−2ελωq−2DLqξ+4ω2ωn2+γsin2τξ=0.Then, consider that(41)ω=ιωn,ωq−2=ιq−2ωnq−2.Hence,(42)ξ″+2q−2ελιq−2ωnq−2DLqξ+4ω2ωn2+γsin2τξ=0,(43)ξ″+2q−2ελιq−2ωnq−2DLqξ+4ι2+4ωn2γι2sin2τξ=0.Then,(44)ξ″+ζDLqξ+δ+εsin2τξ=0.whereζ=2q−2ελιq−2ωnq−2,δ=4ι2,ε=4ωn2γι2,′=ddτ.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 to and the point of interest (t), in this case, we are interested in looking for periodic solutions of equation (44) in δ,ϵ 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 π or 2π, cf. [33, 34].Typically, for given values of the parametersδ and ϵ, all solutions of this equation should be either stable, unstable, or periodic. The δ,ϵ 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:(45)ξτ≈Ao+∑ℓ=1∞Aℓcosℓτ+Bℓsinℓτ.Then, inserting equation (45) into equation (44) and collecting terms, we obtain the following recurrence relations:(46)ℓ=0,−ε2B2=δAo,(47)ℓ=1,1−ζcosπ2qA1−ε2+ζsinπ2qB1−ε2B3=δA1,(48)−ϵ2+ζsinπ2qA1+1−ζcosπ2qB1+ϵ2A3=δB1,(49)ℓ=2,22−2qζcosπq2A2−2qζsinπq2B2−ε2B4=δA2,(50)−εAo+2qζsinπq2A2+22−2qζcosπq2B2+ε2A4=δB2,(51)ℓ≥3,ε2Bℓ−2+ℓ2−ζℓqcosπ2qAℓ−ζℓqsinπ2qBℓ−ε2Bℓ+2=δAℓ,(52)−ε2Aℓ−2+ζℓqsinπ2qAℓ+ℓ2−ζℓqcosπ2qBℓ+ε2Aℓ+2=δBℓ,As shown in Figure2, 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 ζ 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].Figure 2
Transition curves at different values ofq and ζ using equation (44). (a) q = 1. (b) ζ = 0.5.
(a)(b)


### 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).(53)θ¨+hθ,θ˙+fθ,0=0,,θ∈0,π.
Then, the sufficient condition for a stable motion is(54)α4>β+2κ.Proof.
The autonomous form of equation (32), when γ=0, is equivalent to(55)x˙=y,y˙=−ελ+μx2+νy2y−α1−cosx2sin3x+βsinx+κsin2x.
We assume that(56)Hx,y=−ελ+μx2+νy2y−α1−cosx2sin3x+βsinx+κsin2x.
By taking Lyapunov function in the form(57)Vx,y=12y2−∫0xHx,0dx,
we obtain(58)V.=yHx,y−Hx,0.
Hence, we have the following conditions:(59a)i.Hx,0x<0,forx≠0,(59b)ii.Hx,y−Hx,0y<0fory≠0,(59c)iii.−∫0xHu,0du⟶∞asx⟶∞.
Apply the conditions onHx,y; then, the conclusion holds.Theorem 2.
Consider the following form of equation (32):(60)θ¨+h˜θ,θ˙θ˙+fθ,0=et,,θ∈0,π.
Then, the sufficient conditions for a stable motion are(61)λ,μ,ν>0,α4>β+2κ.Proof.
The form of equation (60) is equivalent to the system(62)x˙=y,y˙=−h˜x,yy−fx,0+et,where,(63)h˜x,y=λ+μx2+νy2,fx,0=α1−cosx2sin3x−βsinx−κsin2x.
By taking Lyapunov function in the form(64)Vt,x,y=e−2Et12y2+∫0xfu,0du+1,
we obtain(65)V.<−e−2E∞h˜y2,E∞=∫0∞eudu.
So, every solutionxt exists in the future and xt⟶0 and yt⟶0 as t⟶∞ if the following conditions are satisfied:(66a)i.h˜x,yis continuous and nonnegative on,(66b)ℝ×ℝand is bounded whenx2+y2is bounded,ii.fx,0is continuous inℝ,xfx,0>0,forx≠0,(66c)and∫0xfu,0du⟶∞asx⟶0,iii.etis continuous on0,∞,andEt=∫0teudu<∞,
Apply the conditions; then, the conclusion holds.Theorem 3.
The governing equation, equation (27),(67)θ¨+hθ,θ˙+fθ,t=et,,θ∈0,π,has an asymptotically stable solution if(68)λ,μ,ν>0,α4>β+2κ.Proof.
By considering that(69a)hx,y=ελ+μx2+νy2y,(69b)fx,t=12αtanx2sec2x2−βsinx−κsin2x+γsinωtsinx.
Take that Lyapunov functionV is(70)Vx,y,t=expy2+f∗x,t+k−∫0t|eu|du−kt,in the domain k∈ℝ+, 0≤t≤T, |x|<∞, |y|<∞, and x2+y2≥k2. Hence, we have(71)V.≤0.
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π/ω(iii)
f∗x,t=∫0xfu,tdu>−k for all t and x and ft∗/f∗x,t+k is bounded
Apply the conditions; then, the conclusion holds.


## 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∗,y∗=0,0, which describes a motion in which one of the body principal axes coincides with the local vertical axis(ii)
x∗,y∗=±π,0 which describes an inverted motion of the rigid bodyThus, the linearized homogenous form of equation (37) at a general equilibrium point θ∗  = 0 π or -π for the variable ξt=θt−θ∗ with the time scale 2τ=ωt; then it reads(40)ξ″+2q−2ελωq−2DLqξ+4ω2ωn2+γsin2τξ=0.Then, consider that(41)ω=ιωn,ωq−2=ιq−2ωnq−2.Hence,(42)ξ″+2q−2ελιq−2ωnq−2DLqξ+4ω2ωn2+γsin2τξ=0,(43)ξ″+2q−2ελιq−2ωnq−2DLqξ+4ι2+4ωn2γι2sin2τξ=0.Then,(44)ξ″+ζDLqξ+δ+εsin2τξ=0.whereζ=2q−2ελιq−2ωnq−2,δ=4ι2,ε=4ωn2γι2,′=ddτ.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 to and the point of interest (t), in this case, we are interested in looking for periodic solutions of equation (44) in δ,ϵ 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 π or 2π, cf. [33, 34].Typically, for given values of the parametersδ and ϵ, all solutions of this equation should be either stable, unstable, or periodic. The δ,ϵ 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:(45)ξτ≈Ao+∑ℓ=1∞Aℓcosℓτ+Bℓsinℓτ.Then, inserting equation (45) into equation (44) and collecting terms, we obtain the following recurrence relations:(46)ℓ=0,−ε2B2=δAo,(47)ℓ=1,1−ζcosπ2qA1−ε2+ζsinπ2qB1−ε2B3=δA1,(48)−ϵ2+ζsinπ2qA1+1−ζcosπ2qB1+ϵ2A3=δB1,(49)ℓ=2,22−2qζcosπq2A2−2qζsinπq2B2−ε2B4=δA2,(50)−εAo+2qζsinπq2A2+22−2qζcosπq2B2+ε2A4=δB2,(51)ℓ≥3,ε2Bℓ−2+ℓ2−ζℓqcosπ2qAℓ−ζℓqsinπ2qBℓ−ε2Bℓ+2=δAℓ,(52)−ε2Aℓ−2+ζℓqsinπ2qAℓ+ℓ2−ζℓqcosπ2qBℓ+ε2Aℓ+2=δBℓ,As shown in Figure2, 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 ζ 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].Figure 2
Transition curves at different values ofq and ζ using equation (44). (a) q = 1. (b) ζ = 0.5.
(a)(b)


## 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).(53)θ¨+hθ,θ˙+fθ,0=0,,θ∈0,π.
Then, the sufficient condition for a stable motion is(54)α4>β+2κ.Proof.
The autonomous form of equation (32), when γ=0, is equivalent to(55)x˙=y,y˙=−ελ+μx2+νy2y−α1−cosx2sin3x+βsinx+κsin2x.
We assume that(56)Hx,y=−ελ+μx2+νy2y−α1−cosx2sin3x+βsinx+κsin2x.
By taking Lyapunov function in the form(57)Vx,y=12y2−∫0xHx,0dx,
we obtain(58)V.=yHx,y−Hx,0.
Hence, we have the following conditions:(59a)i.Hx,0x<0,forx≠0,(59b)ii.Hx,y−Hx,0y<0fory≠0,(59c)iii.−∫0xHu,0du⟶∞asx⟶∞.
Apply the conditions onHx,y; then, the conclusion holds.Theorem 2.
Consider the following form of equation (32):(60)θ¨+h˜θ,θ˙θ˙+fθ,0=et,,θ∈0,π.
Then, the sufficient conditions for a stable motion are(61)λ,μ,ν>0,α4>β+2κ.Proof.
The form of equation (60) is equivalent to the system(62)x˙=y,y˙=−h˜x,yy−fx,0+et,where,(63)h˜x,y=λ+μx2+νy2,fx,0=α1−cosx2sin3x−βsinx−κsin2x.
By taking Lyapunov function in the form(64)Vt,x,y=e−2Et12y2+∫0xfu,0du+1,
we obtain(65)V.<−e−2E∞h˜y2,E∞=∫0∞eudu.
So, every solutionxt exists in the future and xt⟶0 and yt⟶0 as t⟶∞ if the following conditions are satisfied:(66a)i.h˜x,yis continuous and nonnegative on,(66b)ℝ×ℝand is bounded whenx2+y2is bounded,ii.fx,0is continuous inℝ,xfx,0>0,forx≠0,(66c)and∫0xfu,0du⟶∞asx⟶0,iii.etis continuous on0,∞,andEt=∫0teudu<∞,
Apply the conditions; then, the conclusion holds.Theorem 3.
The governing equation, equation (27),(67)θ¨+hθ,θ˙+fθ,t=et,,θ∈0,π,has an asymptotically stable solution if(68)λ,μ,ν>0,α4>β+2κ.Proof.
By considering that(69a)hx,y=ελ+μx2+νy2y,(69b)fx,t=12αtanx2sec2x2−βsinx−κsin2x+γsinωtsinx.
Take that Lyapunov functionV is(70)Vx,y,t=expy2+f∗x,t+k−∫0t|eu|du−kt,in the domain k∈ℝ+, 0≤t≤T, |x|<∞, |y|<∞, and x2+y2≥k2. Hence, we have(71)V.≤0.
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π/ω(iii)
f∗x,t=∫0xfu,tdu>−k for all t and x and ft∗/f∗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π/ω. Also, it is considered that, at the equilibrium point θ∗=0, the following assumption holds:(72)h0,0+f0,t=0.Hence, the only solution of equation (72) is θ=θ∗. Equation (32) is rewritten as follows (θ=x, θ.=y):(73)x¨+fx,t=et−hx,y,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π/ω 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π/ω(iii)
f∗x,t=∫0xfu,tdu>−k for all t and ft∗/f∗x,t+k is bounded(iv)
hx,yy≥0(v)
There exists1,s2, s1<s2, hs1,0+fs1,t−et≥0, and hs2,0+fs2,t−et≤0(vi)
fx,t+et≤koProof.
For the domaint∈0,2π/ω, let |x|<∞, |y|<∞, and x2+y2≥k2 and k∈ℝ+. Consider that there exists a continuous function such that(74)|hx,y|≤hoy.
Assume that the following identities hold:(75)∫−∞uhou+kodu⟶∞,∫∞uhou+kodu⟶∞,
and then let the two functionsV1 and V2 be as follows:(76)V1=ex+∫kyuhou+kodu,V2=ex+∫−kyuhou+kodu.
Hence, it can be concluded that there exist in the interior of two domains(77)V.1≥0,V.2≤0.
Then, the solutionsxt and yt exist and are bounded for all t≥0. Then, the conclusion holds.Theorem 5.
In the autonomous case of the approximate governing equationγ=0, there exists a limit cycle under the following conditions:(78)ωn=1,Ω=0,λ<0.Proof.
By using the polar coordinates(79)ρo2=x2+y2,tanσ=yx,
one gets the following:(80)ρ˙o=xx˙+yy˙ρo,σ˙=xy˙−x˙yρo2.
After substitutions and using the conditionsωn=1 and Ω=0, we obtain(81a)ρ˙o=−ερosin2σλ+ρo2μcos2σ+νsin2σ,(81b)σ˙=−1+ε2sin2σλ+r2μcos2σ+νsin2σ,
and one particular solution is(82)ρo2=−λμcos2σ+νsin2σ,σ˙=−1.
Corresponding to the limit cycle,(83)x=−λμcos2σ+νsin2σcosσ,y=−λμcos2σ+νsin2σsinσ.
It is noticed that(84)ρ˙o>0if0<ρo<λμ2+ν2,andρ˙o<0ifρo>λμ2+ν2.
This means that all paths approach the shown radius of the limit cycleρo especially when μ=ν. 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(85)xt=atcosωt+btsinωt,With time-varying coefficientsat 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 x˙=yt have the same form as the linear case,(86)yt=ωbcosωt−aωsinωt,and the following acceleration term is obtained:(87)y˙t=−ω2a+2ωb˙cosωt−bω2+2ωa˙sinωt.Consider the following condition:(88)a˙cosωt+b˙sinωt=0.Then, we have to transform the system equation (85) to equation (88) with equation (35) into the two unknowns at and bt as follows:(89)a˙=−12ωfosinωt,b˙=12ωfocosωt,where(90)foa,b,t=348γba2+b2−12γb+cosωtω2−ωn2a−ελωb−12εμωba2+b2−34ενbω3a2+b2−34aΩa2+b2+F1+sinωtω2−ωn2b+ελωa+14εμωaa2+b2+34ενω3aa2+b2−34Ωba2+b2+F2+cos2ωt12Γb−348Γba2+b2+148Γb3a2−b2+sin2ωt−12Γb+348Γaa2+b2−148Γaa2−3b2+higher harmonics.Consequently, by considering that the major parts ofa and b are slowly varying functions of time, they change very little during the time period T=2π/ω 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(91)a˙=−14π∫0Tfosinωtdt,b˙=14π∫0Tfocosωtdt.Sincea and b can be considered constant in the interval 0,T, a and b (and so a˙ and b˙) can be taken outside the integral signs in equation (91). Typically, the set of functions 1,cosnωt,sinnωt forms an orthogonal set over the interval 0,T. Then, we have, for all m,n=0,1,2,…, from the inner product definition of two functions f1,f2 defined on the interval 0,T, that(92)sinmωt,cosnωt=0,cosmωt,cosnωt=sinmωt,sinnωt=0,ifm≠n.πω,ifm=n.Hence, equation (91) is reduced to(93a)a˙=Ξa,b,b˙=Λa,b,(93b)Ξa,b=−14ωω2−ωn2b+ελωa+14εμωaa2+b2+34ενω3aa2+b2−34Ωba2+b2+F2,(93c)Λa,b=14ωω2−ωn2a−ελωb−14εμωba2+b2−34ενω3ba2+b2−34Ωaa2+b2+F1.Consequently, the initial conditions can be given in terms of those for equation (35) by(94)a0=x0,b0=x˙0ω.


### 5.1. Relative Equilibria and Stability

According to equation (93), it can be deduced that all equilibrium solutions are given by a˙=b˙=0. In addition to the fact that the other paths corresponding to solutions of equation (35) are nonperiodic in general, let ao,bo represent any of the three fixed points, a1,b1,a2,b2anda3,b3, of equation (93). Hence, the linear analysis in the neighborhood of point ao,bo is assumed to be(95)a=ao+η1,b=bo+η2,and then the corresponding linearized system is(96a)η˙1=Ξ1ao,boη1+Ξ2ao,boη2,(96b)η˙2=Λ1ao,boη1+Λ2ao,boη2,where(97a)Ξ1ao,bo=−14ελ+14εμ3ao2+bo2+34ενω23ao2+bo2−32ωΩaobo,(97b)Ξ2ao,bo=−14ωω2−ωn2+12ϵμωaobo+32ενω3aobo−34Ωao2+3bo2,(97c)Λ1ao,bo=14ωω2−ωn2−12εμωaobo−32ενω3aobo−34Ω3ao2+bo2,(97d)Λ2ao,bo=14ω−ελω−14εμωao2+3bo2−34ενω3ao2+3bo2−32Ωaobo.The general conditions of stability for every stable fixed point are(98)Ξ1Λ2>Ξ2Λ1,Ξ1+Λ2≤0.Moreover, the boundary between nodes and spirals can be given by the following identity:(99)Ξ1+Λ22=4Ξ1Λ2−Ξ2Λ1.The phase diagram computed for a particular caseε=0 is shown in Figure 3 for F=0 and F=0.005, respectively. In Figure 3(a), point P1 is a center and represents a stable oscillation, but in Figure 3(b) there exist three equilibrium points: P1 and P2 are centers and P3 is a saddle representing an unstable oscillation.Figure 3
Phase plane trajectories atε  = 0, ω  = 0.975, ωn  = 1, and Ω  = −0.167. (a) F = 1, (b) F = 0.005.
(a)(b)For the particular case ofε=0.005 and λ  =  μ  =  ν  = 1, the phase plane is shown in Figure 4 for F=0 and F=0.005, respectively. In Figure 4(a), point P1 is a spiral point and represents a stable oscillation, but in Figure 4(b) there exist three equilibrium points: P1 and P2 are stable spiral points and P3 is a saddle representing an unstable oscillation.Figure 4
Phase plane trajectories atε  = 0.005, λ  =  μ  =  ν  = 1, ω  = 0.975, ωn  = 1, and Ω  = −0.167. (a) F = 1, (b) F = 0.005.
(a)(b)The phase diagram computed for a particular case whereε=0.005, λ  = −1, and μ  =  ν  = 1 is shown in Figure 5 for F=0 and F=0.005, respectively. In Figure 5(a), point P1 is an unstable point as it bifurcates to generate a stable limit cycle, but in Figure 5(b) there exist three equilibrium points: P1 and P2 are centers and P3 is a saddle representing an unstable oscillation.Figure 5
Phase plane trajectories atε  = 0.005, λ  = −1, μ  =  ν  = 1, ω  = 0.975, ωn  = 1, and Ω  = −0.167. (a) F = 1, (b) F = 0.005.
(a)(b)


### 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 ω,F,ρ 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 ϱ representing the surface of equilibrium points Ξa,b=Λa,b=0 is given by(100)F2=ϱω2−ωn2−34Ωϱ2+ελω+14εωμ+3νω2ϱ2,where(101)ϱ=ρ2=a2+b2.Hence, there may be as many as three positive values ofϱ 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ρ>0 and unstable when dF/dρ<0.Consider the following:(102)Gϱ=916Ω2+116ε2ω2μ+3νω2ϱ3+12ε2λω2μ+3νω2−34Ωω2−ωn2ϱ2+ω2−ωn22+ε2λ2ω2ϱ.SinceG0=0 and Gϱ⟶∞ as ϱ⟶∞, 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 ϱ if(103)dGdϱ=3169Ω2+ε2ω2μ+3νω2ϱ2+ε2λω2μ+3νω2−32Ωω2−ωn2ϱ+ω2−ωn22+ε2λ2ω2=0has two distinct solutions; thus,(104)ϱ1,ϱ2=13/89Ω2+ε2ω2μ+3νω2−ε2λω2μ+3νω2−32Ωω2−ωn2±ε2λω2μ+3νω2−32Ωω2−ωn22−349Ω2+ε2ω2μ+3νω2ω2−ωn22+ε2λ2ω2.So the solution of the amplitude equation is interested inϱ1,ϱ2>0; then,(105)ϱ1+ϱ2>0,ϱ1ϱ2>0.Thus, the amplitude equation has three distinct real roots ifF is controlled by the overlap region.(106)G2<F<G1,G1=Gϱ1,G2=Gϱ2.To explain, we take the case of undamped motionε=0; then, we have(107)ϱ1=23ω2−ωn2Ω,ϱ2=109ω2−ωn2Ω,and then, as shown in Figure6, the overlap region is controlled by(108)155162ω2−ωn23/2Ω<F<16ω2−ωn23/2Ω.Figure 6
Overlap regions for the undamped motion atΩ=1.Hence, the cusps are exactly determined at(109)ω=ωn.In Figure7, 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 ε at the values of 0, 0.1, and 0.25. It is clear that, with the increase of parameter ε at the constant values of damping parameters λ, μ, and ν, the jump is going to disappear.Figure 7
The amplitude surface atωn=1, Ω=1, and λ=μ=ν=1. (a) ε = 0. (b) ε = 0.1. (c) ε = 0.25.
(a)(b)(c)As shown in Figure7, it can be predicted that, for certain values of ω and F at the lower values of ε≈0, there are three responses: the oscillations with greatest and smallest amplitudes are stable but the middle one is unstable. But when ε is increased, then there exists only a single stable response.The consequences of this fold are shown illustratively using Figure8 and 9. So, consider that F is held constant at a value which intersects the fold, and the applied frequency ω is varied in such a way that it crosses the edges of the fold. At the small values of ε≈0, as ω 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.Figure 8
The amplitude versus the frequency atωn=1, Ω=1, and λ=μ=ν=1. (a) ε = 0. (b) ε = 0.1. (c) ε = 0.25.
(a)(b)(c)Figure 9
The amplitude versusF at ωn=1, Ω=1, and λ=μ=ν=1. (a) ε = 0. (b) ε = 0.1. (c) ε = 0.25.
(a)(b)(c)


## 5.1. Relative Equilibria and Stability

According to equation (93), it can be deduced that all equilibrium solutions are given by a˙=b˙=0. In addition to the fact that the other paths corresponding to solutions of equation (35) are nonperiodic in general, let ao,bo represent any of the three fixed points, a1,b1,a2,b2anda3,b3, of equation (93). Hence, the linear analysis in the neighborhood of point ao,bo is assumed to be(95)a=ao+η1,b=bo+η2,and then the corresponding linearized system is(96a)η˙1=Ξ1ao,boη1+Ξ2ao,boη2,(96b)η˙2=Λ1ao,boη1+Λ2ao,boη2,where(97a)Ξ1ao,bo=−14ελ+14εμ3ao2+bo2+34ενω23ao2+bo2−32ωΩaobo,(97b)Ξ2ao,bo=−14ωω2−ωn2+12ϵμωaobo+32ενω3aobo−34Ωao2+3bo2,(97c)Λ1ao,bo=14ωω2−ωn2−12εμωaobo−32ενω3aobo−34Ω3ao2+bo2,(97d)Λ2ao,bo=14ω−ελω−14εμωao2+3bo2−34ενω3ao2+3bo2−32Ωaobo.The general conditions of stability for every stable fixed point are(98)Ξ1Λ2>Ξ2Λ1,Ξ1+Λ2≤0.Moreover, the boundary between nodes and spirals can be given by the following identity:(99)Ξ1+Λ22=4Ξ1Λ2−Ξ2Λ1.The phase diagram computed for a particular caseε=0 is shown in Figure 3 for F=0 and F=0.005, respectively. In Figure 3(a), point P1 is a center and represents a stable oscillation, but in Figure 3(b) there exist three equilibrium points: P1 and P2 are centers and P3 is a saddle representing an unstable oscillation.Figure 3
Phase plane trajectories atε  = 0, ω  = 0.975, ωn  = 1, and Ω  = −0.167. (a) F = 1, (b) F = 0.005.
(a)(b)For the particular case ofε=0.005 and λ  =  μ  =  ν  = 1, the phase plane is shown in Figure 4 for F=0 and F=0.005, respectively. In Figure 4(a), point P1 is a spiral point and represents a stable oscillation, but in Figure 4(b) there exist three equilibrium points: P1 and P2 are stable spiral points and P3 is a saddle representing an unstable oscillation.Figure 4
Phase plane trajectories atε  = 0.005, λ  =  μ  =  ν  = 1, ω  = 0.975, ωn  = 1, and Ω  = −0.167. (a) F = 1, (b) F = 0.005.
(a)(b)The phase diagram computed for a particular case whereε=0.005, λ  = −1, and μ  =  ν  = 1 is shown in Figure 5 for F=0 and F=0.005, respectively. In Figure 5(a), point P1 is an unstable point as it bifurcates to generate a stable limit cycle, but in Figure 5(b) there exist three equilibrium points: P1 and P2 are centers and P3 is a saddle representing an unstable oscillation.Figure 5
Phase plane trajectories atε  = 0.005, λ  = −1, μ  =  ν  = 1, ω  = 0.975, ωn  = 1, and Ω  = −0.167. (a) F = 1, (b) F = 0.005.
(a)(b)


## 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 ω,F,ρ 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 ϱ representing the surface of equilibrium points Ξa,b=Λa,b=0 is given by(100)F2=ϱω2−ωn2−34Ωϱ2+ελω+14εωμ+3νω2ϱ2,where(101)ϱ=ρ2=a2+b2.Hence, there may be as many as three positive values ofϱ 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ρ>0 and unstable when dF/dρ<0.Consider the following:(102)Gϱ=916Ω2+116ε2ω2μ+3νω2ϱ3+12ε2λω2μ+3νω2−34Ωω2−ωn2ϱ2+ω2−ωn22+ε2λ2ω2ϱ.SinceG0=0 and Gϱ⟶∞ as ϱ⟶∞, 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 ϱ if(103)dGdϱ=3169Ω2+ε2ω2μ+3νω2ϱ2+ε2λω2μ+3νω2−32Ωω2−ωn2ϱ+ω2−ωn22+ε2λ2ω2=0has two distinct solutions; thus,(104)ϱ1,ϱ2=13/89Ω2+ε2ω2μ+3νω2−ε2λω2μ+3νω2−32Ωω2−ωn2±ε2λω2μ+3νω2−32Ωω2−ωn22−349Ω2+ε2ω2μ+3νω2ω2−ωn22+ε2λ2ω2.So the solution of the amplitude equation is interested inϱ1,ϱ2>0; then,(105)ϱ1+ϱ2>0,ϱ1ϱ2>0.Thus, the amplitude equation has three distinct real roots ifF is controlled by the overlap region.(106)G2<F<G1,G1=Gϱ1,G2=Gϱ2.To explain, we take the case of undamped motionε=0; then, we have(107)ϱ1=23ω2−ωn2Ω,ϱ2=109ω2−ωn2Ω,and then, as shown in Figure6, the overlap region is controlled by(108)155162ω2−ωn23/2Ω<F<16ω2−ωn23/2Ω.Figure 6
Overlap regions for the undamped motion atΩ=1.Hence, the cusps are exactly determined at(109)ω=ωn.In Figure7, 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 ε at the values of 0, 0.1, and 0.25. It is clear that, with the increase of parameter ε at the constant values of damping parameters λ, μ, and ν, the jump is going to disappear.Figure 7
The amplitude surface atωn=1, Ω=1, and λ=μ=ν=1. (a) ε = 0. (b) ε = 0.1. (c) ε = 0.25.
(a)(b)(c)As shown in Figure7, it can be predicted that, for certain values of ω and F at the lower values of ε≈0, there are three responses: the oscillations with greatest and smallest amplitudes are stable but the middle one is unstable. But when ε is increased, then there exists only a single stable response.The consequences of this fold are shown illustratively using Figure8 and 9. So, consider that F is held constant at a value which intersects the fold, and the applied frequency ω is varied in such a way that it crosses the edges of the fold. At the small values of ε≈0, as ω 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.Figure 8
The amplitude versus the frequency atωn=1, Ω=1, and λ=μ=ν=1. (a) ε = 0. (b) ε = 0.1. (c) ε = 0.25.
(a)(b)(c)Figure 9
The amplitude versusF at ωn=1, Ω=1, and λ=μ=ν=1. (a) ε = 0. (b) ε = 0.1. (c) ε = 0.25.
(a)(b)(c)


## 6. Homoclinic Bifurcation

The existence of homoclinic bifurcation with the support of the other problem parameters, particularlyα and β, 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 systemThe homoclinic solutionThe critical relation(i)x¨+ελx˙+sinx=0x0t=±2sin−1tanht−t0λ≤0(ii)x¨+ελx˙−sinx=0x0t=±2tan−1sinht−t0λ≤0(iii)x¨+ελx˙−x+x3=0x0t=±2secht−t0λ≤0(iv)x¨+ελx˙+x−x3=0x0t=±tanh1/2t−t0λ≤0(v)x¨+ελx˙−x+x5=0x0t=±34sech2t−t0λ≤0(vi)x¨+ελx˙−x+3x5=0x0t=±secht−t0λ≤0(vii)x¨+εμx2−λx˙−x+x3=0x0t=±2secht−t0λ≥4/5μ(viii)x¨+εμx2−λx˙+x−x3=0x0t=±tanh1/2t−t0λ≥1/3μ(ix)x¨+ελ−μx2x˙−x+x5=0x0t=±34sech2t−t0λ≤43/3πμ(x)x¨+εμx2−λx˙−x+3x5=0x0t=±secht−t0λ≥4/3πμ(xi)x¨+ελ+μx2+νx˙2x˙−x+x3=0x0t=±2secht−t0λ≥28μ+12ν/35(xii)x¨+ελx˙−x+x3=εecosωtx0t=±2secht−t0e≥22λ/3πωcoshπ/2ω(xiii)x¨+ελx˙−x+x3=εexcosωtx0t=±2secht−t0e≥4λ/3πω2sinhπ/2ω(xiv)x¨+ελx˙−x+x3=εe1−x2cosωtx0t=±2secht−t0e≥22λ/πω2−ω2coshπ2ω(xv)x¨+ελx˙+sinx=εesintsinxx0t=±2sin−1tanht−t0e≥4λsinhπ/2(xvi)x¨+ελx˙−sinx=εesinωtsinxx0t=±2tan−1sinht−t0|e|≥4λ/πω2sinhπω/2To verify the theoretical results, system (xi) is considered as depicted in Table1. System (xi) has the saddle point at 0,0 with an associated homoclinic trajectory with solution x0t at ε=0. So, using Melnikov’s function, the homoclinic bifurcation occurs approximately where(110)λc=28μ+12ν35.From Figure10, it is shown that there exist three equilibrium points which become stable spirals if λ>λc and unstable if λ<λc. As λ touches the critical value λc=0.81714 the path changes from heteroclinic spiral saddle connection λ>λc to a homoclinic saddle connection at the critical value to heteroclinic saddle connection λ<λc. This transition is due to the fact that parameter λ generates a homoclinic bifurcation at λc.Figure 10
Phase plane for system (xi) at different values ofλ showing the occurrence of homoclinic bifurcation due to parameter λ at λc  = 0.81714. (a) ε = 0.25, λ = 2, μ = −1, v = −0.05. (b) ε = 0.25, λ = 1, μ = 1, v = 0.05. (c) ε = 0.25, λ = 0.817, μ = 1, v = 0.05. (d) ε = 0.25, λ = 0.801, μ = 1, v = 0.05.
(a)(b)(c)(d)


## 7. Route to Chaos

Considering equation (33), some special cases are looked at to allow the parameter space α,β,κ,γ to vary where the chaotic phenomenon might occur. By using the numerical configuration, the critical value of γ at λ=1, μ=1, ω=1, α=1, β=1, and ν  =  κ  =  F  = 0 might be roughly obtained.A simple numerical search reveals a stable periodic solution which is shown in Figure11. Hence, if γ is increased further, the solution bifurcates at a specific value of γ into a solution which wraps around the origin twice. Further increase of γ 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 γ increases.Figure 11
Homoclinic bifurcation and transition to chaos atα  = 1, β  = 1, and κ  =  F  = 0. (a) ε = 0.25, λ = −1, μ = 1, v = 0 and γ = 0.25. (b) ε = 0.25, λ = −1, μ = 1, v = 0 and γ = 0.35. (c) ε = 0.25, λ = −1, μ = 1, v = −0 and γ = 0.43. (d) ε = 0.25, λ = −1, μ = 1, v = −0 and γ = 0.5.
(a)(b)(c)(d)This analysis is confirmed by computing Poincare’s map and Lyapunov exponents to detect period-doubling and chaos. From Figure12, 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.Figure 12
Poincare’s map and Lyapunov exponents atα  = 1, β  = 1, and κ  =  F  = 0. (a) ε = 0.25, λ = −1, μ = 1, v = 0, γ = 0.45. (b) ε = 0.25, λ = −1, μ = 1, v = 0, γ = 0.45.
(a)(b)


## 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.



---

*Source: 4329906-2022-03-20.xml*