# Existence of Solution and Self-Exciting Attractor in the Fractional-Order Gyrostat Dynamical System

Authors: Muhammad Marwan; Gauhar Ali; Ramla Khan
Journal: Complexity
 (2022)

---

## Abstract

This work identifies the influence of chaos theory on fractional calculus by providing a theorem for the existence and stability of solution in fractional-order gyrostat model with the help of a fixed-point theorem. We modified an integer order gyrostat model consisting of three rotors into fractional order by attaching rotatory fuel-filled tank and provided an iterative scheme for our proposed model as a working rule of obtained analytical results. Moreover, this iterative scheme is injected into algorithms for a system of integer order dynamical systems to observe Lyapunov exponents and a bifurcation diagram for our proposed fractional-order dynamical model. Furthermore, we obtained five equilibrium points, including four unstable spirals and one saddle node, using local dynamical analysis which acted as self-exciting attractors and a separatrix in a global domain.


---

## Body

## 1. Introduction

System of ordinary differential equations [1](1)x˙=fx,β.Is called the dynamical system, and a parameterβ∈ℝ in the velocity vector field is termed as bifurcation parameter if system (1) changes its topological structure with the variation in parameter values, whereas the process of changing in qualitative structures is known as bifurcation. There are several types of bifurcation including saddle node [2], Hopf [3–7], and zero-Hopf [8–11]. The bifurcation diagram [12] for the parameter makes it easy for predicting the type of bifurcation and existence of chaos in system (1). Chaos has a vital role in engineering [13–17], medical [18–20], aeronautics [14, 21] and fluid dynamics [22–24]. Apart from the above cited applications, its great influence can also be found in fractional calculus [25–27] and reference therein. Dynamical systems based on ordinary differential equations with an integer order, γ=1, describe velocity vectors, but for fractional order, γ∈0,1, researchers aim to target velocity vectors and replace it by differential equation with order between 0 and 1. Several discretization techniques such as fractional linear multistep [28], Adam [29], predictor-corrector [30], and Adam–Bashforth/Moulton [31] are used to solve fractional-order dynamical systems since decade, but the most flexible scheme with fast convergence in solving nonlinear problems is the variation iteration method (VIM). This technique was used for integer order dynamical systems, but later on modified for fractional-order systems by introducing Lagrangian multiplier [32] into it. Many researchers have enhanced its importance by using it in several engineering-based complex problems, such as in 2006, the variation iteration scheme was utilized for fractional-order systems by Odibat and Momani [33], whereas new development in the VIM was carried by Wu and Baleanu [32] in 2013 to overcome its limitation. Recently in 2021, Kumar and Gupta [34] worked on application of the VIM in a fuzzy-based system.It has been observed from the above-cited work and our knowledge from the literature that dynamical systems related to spacecrafts or its attached devices such as beam and gyrostat have never been considered for the existence of solution and self-exciting attractors in a fractional-order form. Therefore, we have restructured the gyrostat chaotic system [35] into a fractional order along with the addition of a rotatory liquid-filled tank to discuss its unique solution, bounds, and stability using the fixed point theory. Moreover, for bringing novelty into our work, a variation iteration scheme has been used in our proposed fractional-order system to observe chaos into it. For this purpose, several algorithms such as by Wolf et al. [36] and the bifurcation diagram [37] were modified by injecting the VIM iteration scheme into these algorithms. Finally, analyzing local dynamics of our proposed model, trajectories around five equilibrium points with four unstable spirals and a single saddle node motivated us to search for self-exciting attractors with a separatrix in a global domain.The following pattern can be followed for understanding the rest of the paper. In Section2, the gyrostat chaotic system is remodeled by adding rotatory liquid-filled tank and modified into fractional order. Several theorems have been proved in Section 3 for the existence of solution and stability. An iteration scheme for our proposed model has been introduced in Section 4, while several applications of this scheme related to dynamical analysis are discussed in Section 5. Finally, Section 6 comprises concluding remarks and future target.


## 2. Modeling of Gyro Chaotic System Attached with Fuel-Filled Tank

Gyrostat is a device consisting of rotors, used as an attachment in larger objects for bringing stability in their dynamics with the passage of time. The system of three-dimensional ordinary differential equations for the gyrostat model is designed by Qi et al. [35]:(2)Ixx˙=Iy−Izyz−yhz+zhy−Tx+Lx,Iyy˙=Iz−Ixzx−zhx+xhz−Ty+Ly,Izz˙=Ix−Iyxy−xhy+yhx−Tz+Lz,where X=xyzT is the angular velocity vector, I=Ix,Iy,Iz are the principal moments of inertia of the gyrostat in the body axis frame, H=hx,hy,hz are constants of total angular momentum, whereas L=Lx,Ly,Lz and T=Tx,Ty,Tz are external and disturbed torques applied on the gyrostat, respectively.A tank, rotating about an angleθ at desired point x¯,y¯,z¯ (shown in Figure 1), is attached with an originally disturbed gyrostat system given in (2). It is observed that L=Lx,Ly,Lz is a vector of external forces applied on the gyrostat. Therefore, we have attached a tank containing fuel which exert external forces on the gyrostat due to rotation of the attached tank with respect to z− axis about angle θ at desired point x¯,y¯,z¯. Hence, we replaced L with(3)R=Rzθd¯x,where R=Rx,Ry,Rz, Rzθ is the rotation matrix for z− axis and d¯x is a desired point about which one can rotate the attached tank. Therefore, using vector T in ((2) [38]) and R given in (3) into (2), we obtain the following equation:(4)IxCDγx=Iy−Izyz−yhz+zhy−μxx+x¯cosθ−y¯sinθ,IyCDγy=Iz−Ixzx+xhz+μyy+x¯sinθ+y¯cosθ,IzCDγz=Ix−Iyxy−xhy−μzz+z¯,Figure 1
Spacecraft attached with partially filled fuel tank and three rotors. This figure is reproduced from the work of Sabir et al. [7] with the assumption that the third rotor is rotating about z-axis.System (4) is a fractional-order mathematical representation of the model given in Figure 1 in which γ is the fractional number between 0 and 1 exclusive, u = ux,uy,uz is a damping constant vector, while X, I, and H are defined in equation (2). Moreover, system (4) shows chaotic behavior for X=0.1,0.1,0.1, I=0.85,0.45,0.2, u=6,6.42,5.8, H=0,0.57416,2,38, and R=−20,2,20. The phase portrait of system (4) with given initial and parameter values can be seen in region 10 of Figure 2.Figure 2
Existence of chaos in clockwise and anticlockwise directions.


## 3. Existence and Stability of Solution

In this part of our paper, we determined results based on the existence theory for system (4) using the fixed point theorem with Banach space. Therefore, basic definitions and important lemmas are considered for the understanding of this work.Definition 1 (see [39]).
The integral of fractional orderγ0 for a function ϖ is given by(5)ℐγϖt=1Γγ∫0tt−ηγ−1ϖηdη.Definition 2 (see [39]).
The Caputo fractional derivative of orderγ0 of a continuous function ϖ is given by(6)DCγϖt=1Γn−γ∫0tt−ηn−γ−1ϖnηdη,where n=γ+1.
Following two lemmas have importance in achieving the solutions of the systems consisting of fractional differential equations.Lemma 1 (see [39]).
Assumeϖ∈C0,1, then the solution of fractional differential equation(7)DCγϖt=0,
Of orderγ0 is(8)ϖt=c0+c1t+⋯+cn−1tn−1,ci∈ℛ,i=0,2,…,n−1.Lemma 2 (see [39]).
Let us considerϖ∈C0,1, with a derivative of fractional order γ, then(9)ℐγDCγϖt=ϖt+c0+c1t+⋯+cn−1tn−1,ci∈ℛ,i=0,2,…,n−1.We begin our work by introducingϖ1, ϖ2, and ϖ3 on the right side of equation (4) and for convenience, we use the following notions:(10)χt=xt,yt,zt,χ0=x0,y0,z0,Ψt,χt=ϖ1t,x,y,z,ϖ2t,x,y,z,ϖ3t,x,y,z,and(11)Ψ0=ϖ10,x0,y0,z0,ϖ20,x0,y0,z0,ϖ30,x0,y0,z0,System (4) can be rewritten, using (5) as(12)DCtγχt=Ψt,χt,t∈0,τ,χ0=χ0.According to Lemma2, problem (12) can be converted into an integral equation(13)χt=χ0+1Γγ∫0tt−ηγ−1Ψη,χηdη.Definition 3.
Let us considerX=C0,τ a Banach space under the suitable norm(14)χ=supt∈Jχ:χ∈X,
and the operator is defined as(15)τχt=χ0+1Γγ∫0tt−ηγ−1Ψη,χηdη,where 0≤t≤τ〈∞ and J=0,τ. Then, the following assumptions are true:
A1 There exists a positive constant ðΨ〉0 such that(16)Ψt,χt−Ψt,χ¯t≤ðΨχt−χ¯t.
A2 The following inequality holds for positive constants MΨ,GΨ〉0:(17)Ψt,χt≤MΨχ+GΨ.Theorem 1.
Let us consider thatðΨτγ〈Γγ+1 and assumption A2 is satisfied. Then, there exists a unique solution of system (4) with the contraction of operator τ.Proof.
Letχ,χ¯∈X, then one has(18)τχ−τχ¯=supt∈J|χ0+1Γγ∫0tt−ηγ−1Ψη,χηdη−χ0+1Γγ∫0tt−ηγ−1Ψη,χ¯ηdη|≤1Γγ∫0τt−ηγ−1Ψη,χη−Ψη,χ¯ηdη≤1Γγ∫0τt−ηγ−1dηðΨχη−χ¯η≤ðΨτγΓγ+1χ−χ¯.
This shows thatτ is a contraction. Hence, our desired result is obtained, that system (4) has a unique solution.Theorem 2.
The integral equation (8) has at least one solution if ðΨ〈Gγ under the assumptions of A1 and A2.Proof.
For existence of a solution for operatorτ, it is enough to show that τ is completely continuous, and there exists an element χ∈X such that χ=δτχ for δ∈0,1. Therefore, our proof will pass through three steps for achieving our desired results.Step 1.
Let us consider a sequenceχn⟶χ in X and for each t∈J, we have(19)τχn−τχ=supt∈J|χ0+1Γγ∫0tt−ηγ−1Ψη,χnηdη−χ0+1Γγ∫0tt−ηγ−1Ψη,χηdη|≤1Γγ∫0τt−ηγ−1Ψη,χnη−Ψη,χηdη.
Hence,τχn approaches τχ as time t tends to infinity(20)τχn−τχ⟶0asn⟶∞.
Equation (20) identifies continuity of an operator τ.Step 2.
Let us consider a bounded setBr=χ∈X:χ≤r, where r is a positive real number. Then, for any χ∈Br, we have(21)τχ=supt∈Jχ0+1Γγ∫0tt−ηγ−1Ψη,χηdη≤χ0+1Γγ∫0ττ−ηγ−1Ψη,χηdη≤χ0+MΨχ+GΨΓγ∫0ττ−ηγ−1dη≤χ0+MΨr+GΨτγΓγ+1=l.
Hence,τ maps a bounded set into a bounded set.Step 3.
The image of a bounded set underτ is equicontinuous in X.
Lett1〈t2 in J and χ∈Br, we have(22)τχt2−τχt1=|χ0+1Γγ∫0t2t2−ηγ−1Ψη,χηdη−χ0+1Γγ∫0t1t1−ηγ−1Ψη,χηdη|=|1Γγ∫0t1t2−ηγ−1Ψη,χηdη+∫t1t2t2−ηγ−1Ψη,χηdη−∫0t1t1−ηγ−1Ψη,χηdη|≤1Γγ∫0t1t2−ηγ−1+t1−ηγ−1Ψη,χηdη+∫t1t2t2−ηγ−1Ψη,χηdη≤MΨr+GΨΓγ+1∫0t2t2−ηγ−1dη−∫0t1t1−ηγ−1dη≤MΨr+GΨΓγ+1t2γ−t1γ.
Ast2⟶t1, then τχt2−τχt1⟶0, and thus, τ is continuous and bounded. Hence, τχt2−τχt1⟶0 shows uniform continuity of τ. Therefore, steps 1−3 show that τ is completely continuous.Step 4.
Finally, we have to show thatℬ=χ∈X:χ=δτχ for some δ∈0,1, is bounded. Let χ∈ℬ and for any t, we have(23)χ=supt∈JδTχ=supt∈J|δχ0+δΓγ∫0tt−ηγ−1Ψη,χηdη|≤χ0+1Γγ∫0ττ−ηγ−1Ψη,χηdη,asδ≤1,≤χ0+MΨχ+GΨΓγ∫0ττ−ηγ−1dη,χ≤χ0+MΨχ+GΨτγΓγ+1.
Simplifying inequality (19) yields(24)χ≤χ0Γγ+1+GΨτγΓγ+1−MΨτγ.
This shows that the defined setℬ is bounded. Hence using Schaefer’s theorem [40], system (4) has at least one solution.
For achieving stability, a negligible perturbation parameterθt can be included in D0CFtγχt such that(i)
D0CFtγχt=Ψt,χt+θt(ii)
θt〈ε for ε0Lemma 3.
Solution of the perturbed problem(25)D0CFtγχt=Ψt,χt+θtχ0=χ0,
Satisfies the following relation:(26)χt−χ0+Ψt,χt−Ψ01−γGγ+γGγ∫0tΨη,χηdη≤LΨε.Theorem 3.
Gyrostat system (4) achieves Ullam–Hyers stability if LΨ〈1 and assumption A2, together with Lemma 3 is satisfied.Proof.
Letχ∈X be any solution and χ¯∈X is a unique solution, then(27)χt−Y¯t=χt−χ0+Ψt,χ¯t−Ψ01−γGγ+γGγ∫0tΨη,χ¯ηdη≤χt−χ0+Ψt,χt−Ψ01−γGγ+γGγ∫0tΨη,χηdη+|χ0+Ψt,χt−Ψ01−γGγ+γGγ∫0tΨη,χηdη−χ0+Ψt,χ¯t−Ψ01−γGγ+γGγ∫0tΨη,χ¯ηdη|≤LΨε+LΨχ−χ¯.
This implies that(28)χ−χ¯≤LΨ1−LΨε.
Hence, solution of the proposed system (4) is Ullam–Hyers stable.


## 4. Variational Iterative Scheme for System (4)

An iterative scheme, variational iterative method (VIM), is introduced in this section using successive approximations of the solution for rapid convergence and analytical results discussed in Section2.


### 4.1. Working Rule

To express the VIM, we consider general nonlinear differential equation as(29)ℒχt+Nχt=ht,where ℒ, N, and h are linear, nonlinear, and source functions, while the corrector function for (29) is considered as(30)χn+1t=χnt+∫0tΛηℒχη+N˜χη−hηdη.Λ in (30) is defined as −1nt−ηn−1/Γn, whereas N˜ is used as a restricted value with σN˜=0. Then, the exact solution can be obtained as(31)χt=limn⟶∞χnt.System (4) can be discretized using VIM as(32)xn+1=xn−ℐγIxCDγxn−Iy−Izynzn+ynhz−znhy+μxxn−x¯cosθ+y¯sinθ,yn+1=yn−ℐγIyCDγyn−Iz−Ixznxn−xnhz−μyyn−x¯sinθ−y¯cosθ,zn+1=zn−ℐγIzCDγz−Ix−Iyxnyn+xnhy+μzzn+z¯,where x0=x0,y0=y0, and z0=z0.Forn=0,1(33)x1=x0+xγκx1Γγ+1,y1=y0+yγκy1Γγ+1,z1=z0+zγκz1Γγ+1,x2=x0+κx2Γγ+1xγ−μxκx1Γ2γ+1x2γ,y2=y0+κy2Γγ+1yγ+κy1Γ2γ+1y2γ,z2=z0+κz2Γγ+1zγ−μz−hyκz1Γ2γ+1z2γ,where(34)κx1=−Iy−Izy0z0+y0hz−z0hy+μxx0−x¯cosθ+y¯sinθ,κx2=κx11+Ix−Iy−Izy1z1+y1hz−z1hy+μxx0−x¯cosθ+y¯sinθ,κy1=−Iz−Ixz0x0−x0hz−μyy0−x¯sinθ−y¯cosθ,κy2=Ly1+Iy−Iz−Ixz1x1−x1hz−μyy0−x¯sinθ−y¯cosθ,κz1=−Ix−Iyx0y0+x0hy+μzz0−z¯,κz2=Lz1+Iz−Ix−Iyx1y1+x1hy+μzz1−z¯.Forn=2, we have(35)x3=x0+ρ1/xΓγ+1xγ−ρ2/xΓ2γ+1x2γ+ρ3/xΓ3γ+1x3γ,y3=y0+1/yΓγ+1yγ+2/yΓ2γ+1y2γ+3/yΓ3γ+1y3γ,z3=z0+ς1/zΓγ+1zγ−ς2/zΓ2γ+1z2γ+ς3/zΓ3γ+1z3γ.The values ofρxi, ϱyi and ςzi, i=1,2,3 in (35) are given in Appendix A. In the next section, we have discussed system (4) analytically and qualitatively. For numerical simulations, our designed algorithm is used to plot Lyapunov exponents and bifurcation diagram in integer order as well as fractional-order chaotic systems.


## 4.1. Working Rule

To express the VIM, we consider general nonlinear differential equation as(29)ℒχt+Nχt=ht,where ℒ, N, and h are linear, nonlinear, and source functions, while the corrector function for (29) is considered as(30)χn+1t=χnt+∫0tΛηℒχη+N˜χη−hηdη.Λ in (30) is defined as −1nt−ηn−1/Γn, whereas N˜ is used as a restricted value with σN˜=0. Then, the exact solution can be obtained as(31)χt=limn⟶∞χnt.System (4) can be discretized using VIM as(32)xn+1=xn−ℐγIxCDγxn−Iy−Izynzn+ynhz−znhy+μxxn−x¯cosθ+y¯sinθ,yn+1=yn−ℐγIyCDγyn−Iz−Ixznxn−xnhz−μyyn−x¯sinθ−y¯cosθ,zn+1=zn−ℐγIzCDγz−Ix−Iyxnyn+xnhy+μzzn+z¯,where x0=x0,y0=y0, and z0=z0.Forn=0,1(33)x1=x0+xγκx1Γγ+1,y1=y0+yγκy1Γγ+1,z1=z0+zγκz1Γγ+1,x2=x0+κx2Γγ+1xγ−μxκx1Γ2γ+1x2γ,y2=y0+κy2Γγ+1yγ+κy1Γ2γ+1y2γ,z2=z0+κz2Γγ+1zγ−μz−hyκz1Γ2γ+1z2γ,where(34)κx1=−Iy−Izy0z0+y0hz−z0hy+μxx0−x¯cosθ+y¯sinθ,κx2=κx11+Ix−Iy−Izy1z1+y1hz−z1hy+μxx0−x¯cosθ+y¯sinθ,κy1=−Iz−Ixz0x0−x0hz−μyy0−x¯sinθ−y¯cosθ,κy2=Ly1+Iy−Iz−Ixz1x1−x1hz−μyy0−x¯sinθ−y¯cosθ,κz1=−Ix−Iyx0y0+x0hy+μzz0−z¯,κz2=Lz1+Iz−Ix−Iyx1y1+x1hy+μzz1−z¯.Forn=2, we have(35)x3=x0+ρ1/xΓγ+1xγ−ρ2/xΓ2γ+1x2γ+ρ3/xΓ3γ+1x3γ,y3=y0+1/yΓγ+1yγ+2/yΓ2γ+1y2γ+3/yΓ3γ+1y3γ,z3=z0+ς1/zΓγ+1zγ−ς2/zΓ2γ+1z2γ+ς3/zΓ3γ+1z3γ.The values ofρxi, ϱyi and ςzi, i=1,2,3 in (35) are given in Appendix A. In the next section, we have discussed system (4) analytically and qualitatively. For numerical simulations, our designed algorithm is used to plot Lyapunov exponents and bifurcation diagram in integer order as well as fractional-order chaotic systems.


## 5. Dynamical Analysis

The fractional-order dynamical system exhibits chaos for some values of fractional term,γ, but using a hit and trial method for such purpose is difficult to investigate chaos in dynamical systems. Therefore, we plotted the bifurcation diagram for system (4) with respect to fractional term, γ. A noisy dense area is observed in Figure 3 that illustrates occurrence of chaos in the fractional-order gyrostat system beginning with γ=0.86.Figure 3
Bifurcation diagram of system (4) with respect to γ.In section2, we used the concept of the fixed point theory to obtain at least one solution of system (4). Hence, for fixed points, we consider the following function ϖi, i=1,2,3 equals to zero:(36)ϖ1=Iy−Izyz−yhz+zhy−μxx+x¯cosθ−y¯sinθ=0,ϖ2=Iz−Ixzx+xhz+μyy+x¯sinθ+y¯cosθ=0,ϖ3=Ix−Iyxy−xhy−μzz+z¯=0.After fixing all parameters given in section (1), then we solve equation (36) to get the following five equilibrium points:(37)E1=−0.819,3.141,3.352T,E2=−11.078,14.8464,−6.79773T,E3=12.445,−15.8945,−11.4255T,E4=11.5523,22.4277,20.1731T,E5=−12.4328,−23.1278,24.5096T.In Theorem4 local dynamical analysis of system (4) is used for observing trajectories around equilibrium points (27).Theorem 4.
A gyrostat chaotic system (4) is composed of five equilibrium points, in which E1 is the saddle node and E2,3,4,5 are all unstable saddle spirals. Moreover, these spirals lead to four attractors and one saddle node E1 that act as a separatrix as t extends.Proof.
Five equilibrium points are calculated in equation (27). The Jacobian matrix plays a vital role in the system of differential equations for local dynamical analysis. Therefore, the Jacobian matrix of system (4) is(38)J=−120175z17−1455y17+71771062523845−13z921415−13x92y−717725002x−29.
And the Jacobian matrix for fixed parameter values atE1 is(39)J|E1=−7.0588−1.81421.59930.447314.26671.18423.4111−1.6397−29.0000.
The characteristic equations of the Jacobian matrix (29) is(40)λ13+21.79λ12−312.44λ1−2796.9.
Solution of (40) results into single positive and two negative eigenvalues:(41)λ11=−29.21,λ12=−6.75,λ13=14.17.
Equation (41) illustrates that two states will move away from E1, while a single state will move inward towards equilibria: E1, and such information shows occurrence of the saddle. The Jacobian matrix at E2 is(42)J|E2=−7.0588−4.79935.042115.107814.266716.001626.8220−22.1560−29.0000.
And the corresponding characteristic equation is(43)λ23+21.79λ22−17.93λ2+7361.79.
Solution of (43) gives three eigenvalues with one negative real and two complex numbers with positive real part:(44)λ21=−30.37,λ223=4.2871±14.9684ι.
Equation (44) describes occurrence of the unstable spiral. In a similar fashion, one can achieve(45)λ31=−30.4682,λ323=4.3380±19.4641ι,λ41=−35.8534,λ423=7.0306±17.3188ι,λ51=−35.8109,λ523=7.0094±20.5301ι.
Eigenvalues ofE3, E4, and E5. In view of (45), equilibrium points, E3,4,5 are also unstable spirals.
Analytical results (29–35) are explained in Figure 4, which illustrate trajectories of system (4) around their equilibrium points. Five different colors are used for each equilibrium point, which are also highlighted as a legend in Figure 4. It is observed that E2,3,4,5 are unstable spirals plotted in green, brown, blue, and black colors, while red color shows the saddle node. In detail, we can see that the red trajectory starts from E1 and passes through the regions of E2, E3, E4, and E5 with the passage of time. The trajectory for E2 shows a spiral emerging from its equilibrium point and is moving away from it. After some time it has been observed that the green trajectory is acting as a heteroclinic orbit: E5. The same theory can be observed between E2 and E3, when a brown colored orbit starts with high unstable oscillations and approaches to a region occupied by E2. Apart from these four unstable spirals, one can also locate saddle node equilibria E1 in the red color, in which its trajectory passes through regions acquired by unstable spirals and act as a separatrix between them. For further analysis, we have extended time for observing the trajectories around five equilibrium points in the greater domain. It has been analyzed that four unstable equilibrium points are self-exciting attractors and occupy four basins. Moreover, the combination of all these four regions leads to the concept of a strange attractor in system (4). Studying in more depth, it has been also observed that the saddle node in the global domain is busy in separating regions of self-exciting attractors. For getting more knowledge about chaoticity in the fractional-order gyrostat system (4), some basic results are used for the possibility and detection of chaos.Figure 4
Trajectories of system (4) around equilibrium points.


### 5.1. Lyapunov Exponents

The Lyapunov exponent is one of the fundamental results, which help researchers in pointing out existence of unpredictability in trajectories of their corresponding systems. Moreover, in a three-dimensional autonomous system of ordinary differential equations, there exist three Lyapunov exponentsλi, i=1,2,3. Now, if λ1,λ2,λ3 = +ive,0,−ive, then it shows existence of chaos, whereas −ive,−ive,+ive illustrates the existence of periodic solutions. In Figure 5, three Lyapunov exponents can be observed, emerging from −9,11,−30 and leading to −26.51,0,4.87, which motivated us to work further on it and find out chaotic trajectories in it. For further investigation, we have used the concept of the bifurcation diagram [12].Figure 5
Lyapunov exponents of system (4).


### 5.2. Bifurcation Leading to Chaos

For confirmation of existence of chaos in our proposed model (4), we fixed all other parameter values except for μy. For damping coefficient, μy, it is observed in Figure 6 that there seems no bifurcation in system (4) for μy∈1,3.6. The single bifurcation emerges at μy=3.6 and continues till μy=5, which changes into period doubling bifurcation (PDB) for 5≤μy≤5.2. Trajectories of our proposed system jump into the chaotic region for μy lying in interval 5.2,6.5. One can observe symmetric behavior in sense of bifurcation leading to chaos in Figure 6. If we start from μy=11, two lines can be observed that are converted into period doubling bifurcations, then period 2⟶ period 4⟶8⟶ chaotic region. This concept is also explained with the aid of a series of phase portraits, which confirms chaotic behavior in our proposed system. Therefore, we have divided bifurcation diagram 5 in nine regions and plotted phase portraits to their corresponding values.Figure 6
Bifurcation diagram for damping coefficient,μy in system (4).Figure2 is validation of Figure 6, which explains existence of chaos in detail by moving clockwise or anticlockwise. Therefore, we have an indexed sequence of phase portraits for μy. If we start from region 1, a spiral trajectory can be observed and is expanding in regions 2 and 3. This trajectory is converted into period doubling and period 4 bifurcations in region 4 and 5 for μy=4 to 5.1, respectively. In region 6, chaotic movement of trajectories can be observed, which gradually declines to period 8⟶ period 4⟶ period 2 bifurcations by moving in the anticlockwise direction from region 6 to 9. Similarly, if we begin in the clockwise direction, region 9 to 1, one can see trajectory starts with period doubling bifurcation for μy=9.6 is gradually increasing to period 4⟶ period 8⟶ chaos from region 9 to 6, then decline in a symmetric way is observed from chaos to the period doubling bifurcation till region 4 which finally shrinks into spiral and bifurcation disappearing in region 3 to 1.Figure7 is the series of Lyapunov exponents corresponding to each subregion of the bifurcation diagram (plotted in Figure 6). In Figure 2, the existence of chaos in a symmetrical way is thoroughly discussed, but in Figures 7(a)–7(h) the same concept is explained in more detail where for each value of the damping coefficient μy; there exist different values of Lyapunov exponents. Moreover, it is also observed that the Lyapunov exponent of system (4) tends to −ive,0,+ive as the damping coefficient μy approaches to 6.43.Figure 7
Lyapunov exponents corresponding to each subregion in Figure6.
(a)(b)(c)(d)(e)(f)(g)(h)


## 5.1. Lyapunov Exponents

The Lyapunov exponent is one of the fundamental results, which help researchers in pointing out existence of unpredictability in trajectories of their corresponding systems. Moreover, in a three-dimensional autonomous system of ordinary differential equations, there exist three Lyapunov exponentsλi, i=1,2,3. Now, if λ1,λ2,λ3 = +ive,0,−ive, then it shows existence of chaos, whereas −ive,−ive,+ive illustrates the existence of periodic solutions. In Figure 5, three Lyapunov exponents can be observed, emerging from −9,11,−30 and leading to −26.51,0,4.87, which motivated us to work further on it and find out chaotic trajectories in it. For further investigation, we have used the concept of the bifurcation diagram [12].Figure 5
Lyapunov exponents of system (4).


## 5.2. Bifurcation Leading to Chaos

For confirmation of existence of chaos in our proposed model (4), we fixed all other parameter values except for μy. For damping coefficient, μy, it is observed in Figure 6 that there seems no bifurcation in system (4) for μy∈1,3.6. The single bifurcation emerges at μy=3.6 and continues till μy=5, which changes into period doubling bifurcation (PDB) for 5≤μy≤5.2. Trajectories of our proposed system jump into the chaotic region for μy lying in interval 5.2,6.5. One can observe symmetric behavior in sense of bifurcation leading to chaos in Figure 6. If we start from μy=11, two lines can be observed that are converted into period doubling bifurcations, then period 2⟶ period 4⟶8⟶ chaotic region. This concept is also explained with the aid of a series of phase portraits, which confirms chaotic behavior in our proposed system. Therefore, we have divided bifurcation diagram 5 in nine regions and plotted phase portraits to their corresponding values.Figure 6
Bifurcation diagram for damping coefficient,μy in system (4).Figure2 is validation of Figure 6, which explains existence of chaos in detail by moving clockwise or anticlockwise. Therefore, we have an indexed sequence of phase portraits for μy. If we start from region 1, a spiral trajectory can be observed and is expanding in regions 2 and 3. This trajectory is converted into period doubling and period 4 bifurcations in region 4 and 5 for μy=4 to 5.1, respectively. In region 6, chaotic movement of trajectories can be observed, which gradually declines to period 8⟶ period 4⟶ period 2 bifurcations by moving in the anticlockwise direction from region 6 to 9. Similarly, if we begin in the clockwise direction, region 9 to 1, one can see trajectory starts with period doubling bifurcation for μy=9.6 is gradually increasing to period 4⟶ period 8⟶ chaos from region 9 to 6, then decline in a symmetric way is observed from chaos to the period doubling bifurcation till region 4 which finally shrinks into spiral and bifurcation disappearing in region 3 to 1.Figure7 is the series of Lyapunov exponents corresponding to each subregion of the bifurcation diagram (plotted in Figure 6). In Figure 2, the existence of chaos in a symmetrical way is thoroughly discussed, but in Figures 7(a)–7(h) the same concept is explained in more detail where for each value of the damping coefficient μy; there exist different values of Lyapunov exponents. Moreover, it is also observed that the Lyapunov exponent of system (4) tends to −ive,0,+ive as the damping coefficient μy approaches to 6.43.Figure 7
Lyapunov exponents corresponding to each subregion in Figure6.
(a)(b)(c)(d)(e)(f)(g)(h)


## 6. Conclusion and Future Work

An integer ordered dynamical system of the gyrostat was considered by researchers since decade, and a variety of work related to chaos was achieved with the help of sensitivity in its initial conditions. But we have analyzed the gyrostat model with modification by attaching a rotatory cylinder and conversion into fractional order for the first time. Several theorems were proved in this work for the existence of solution and Ullam–Hyers stability. Moreover, dealing with the fractional-order system does not work onODE45; therefore, an iterative scheme was designed for system (4) to attain chaos in the fractional order. Studying local dynamics of system (4) leads to five solutions with four unstable spirals and one saddle node, but observing trajectories around these equilibrium points in global domain acted as a self-exciting attractor and separatrix. In future, we aim to target fractional-order dynamical systems for codimension 2 bifurcations, which itself is a tedious task due to a large number of involved parameters. Apart from bifurcation, our future aim also involves application of (integer and fractional) ordered chaotic systems in strategy-based mobile gaming.



---

*Source: 3505634-2022-07-18.xml*