# Harmonic Balance Method for Chaotic Dynamics in Fractional-Order Rössler Toroidal System

Authors: Huijian Zhu
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

This paper deals with the problem of determining the conditions under which fractional order Rössler toroidal system can give rise to chaotic behavior. Based on the harmonic balance method, four detailed steps are presented for predicting the existence and the location of chaotic motions. Numerical simulations are performed to verify the theoretical analysis by straightforward computations.


---

## Body

## 1. Introduction

The concept of fractional calculus has been known since the development of the regular calculus, with the first reference probably being associated with Leibniz and L’Hôspital in 1695. In the past three decades or so, fractional calculus gained considerable popularity and importance, due mainly to its demonstrated applications in numerous seemingly diverse and widespread fields of science and engineering [1–3]. In particular, by utilizing fractional calculus technique, many investigations were devoted to the chaotic behaviors and chaotic control of dynamical systems involving the fractional derivative, called fractional-order chaotic system [4–8]. For example, it has been shown that Chua circuit of the order as low as 2.7 can behave in chaotic manner [4]. In [5], the nonautonomous Duffing systems of the order less than 2 can still produce a chaotic attractor. In [6], chaotic behavior of the fractional-order Lorenz system was further studied. In [7], chaos and hyperchaos in the fractional-order Rössler equations were also studied, in which the authors showed that chaos can exist in the fractional-order Rössler equation with the order as low as 2.4, and hyperchaos can exist in the fractional-order Rössler hyperchaos equation with the order as low as 3.8. Later on, the chaotic behavior and its control in the fractional-order Chen system were investigated in [8]. And recently, more dynamic behaviors of fractional order chaotic systems were analyzed by using different approaches; we refer the readers to [9–12].However, to our knowledge, the conditions for chaos existence in dynamical systems (including integer-order system and fractional-order system) are still incomplete. For a given dynamical system, can we decide (without invoking numerical simulations) whether and in what parameter ranges, chaotic behavior might exist? It is still an open problem [13]. Even though some theorems such as Melnikov’s criteria [14] and Shil’nikov’s theorem [15] may be helpful in some special cases, it seems that a powerful tool is not generally available to determine the accurate parameter ranges for chaos existence in a given dynamical system. Moreover, the fact that fractional-order system possess long memory characteristics makes the system behave in more complicated manner. This means that chaos analysis in fractional-order systems is more complicated than that in integer-order systems. Recently, based on harmonic balance method (or describing function method), predicting chaotic dynamics in fractional-order Genesio-Tesi system was investigated [16]. Later on, two describing function-based techniques were proposed to predict chaos in a class of fractional-order differential equations [17]. Therefore, harmonic balance method can give another way to detect the conditions for chaos existence in fractional-order systems. Following them, the subject of this paper is predicting conditions for chaos existence in fractional-order Rössler toroidal system.The paper is organized as follows. In Section2, some basic concepts of fractional calculus are briefly introduced. Sections 3 and 4 are about predicting chaos in fractional-order Rössler toroidal system by utilizing harmonic balance method. In Section 5, numerical simulations are further presented to confirm the obtained results of Section 4. Finally, the brief conclusions are drawn in Section 6 and close the paper.


## 2. Preliminaries

### 2.1. Basic Definitions

In what follows, the definitions of fractional derivative are introduced in brief. For an arbitrary real numberq, the Riemann-Liouville and Caputo fractional derivatives are defined, respectively, as
(1)aDtqf(t)=1Γ(m-q)dmdtm∫at‍(t-τ)m-q-1f(τ)dτxxxxxxxxxxxxxxxxxxxxxxx(m=[q]+1),(2)aCDtqf(t)=1Γ(m-q)∫at‍f(m)(τ)(t-τ)q+1-mdτ00000000000000000000(m=[q]+1),
where[q]stands for the integer part ofq; that is,mis the first integer which is not less thanq.Note that the main advantage of Caputo approach is that the initial conditions for fractional differential equations with Caputo derivatives take on the same form as for integer-order differential equations. Some other important properties of the fractional derivatives and integrals can be found in [1–3].


### 2.2. Numerical Method

According to [18, 19], the predictor-corrector scheme for fractional-order system is derived. This scheme is the generalization of Adams-Bashforth-Moulton scheme. First, the following fractional-order differential equation involving Caputo derivative
(3)dqxdtq=f(t,x),0≤t≤T,x(k)(0)=x0(k),k=0,1,2,…,m-1,
is equivalent to the Volterra integral equation (see [18, Lemma  2.1])
(4)x=∑k=0⌈q⌉-1‍x0(k)tkk!+1Γ(q)∫0t‍f(τ,x)(t-τ)1-qdτ.
Seth=T/N, tn=nh, n=0,1,2,…,N, and system (4) can be discretized as
(5)xh(tn+1)=∑k=0⌈q⌉-1‍x0(k)tn+1kk!+hqΓ(q+2)f(tn+1,xhp(tn+1))+hqΓ(q+2)∑j=0n‍aj,n+1f(tj,xh(tj)),
where
(6)aj,n+1={nq+1-(n-q)(n+1)q,j=0,(n-j+2)q+1+(n-j)q+1-2(n-j+1)q+1,1≤j≤n,1,j=n+1,xhp(tn+1)=∑k=0⌈q⌉-1‍x0(k)tn+1kk!+1Γ(q)∑j=0n‍bj,n+1f(tj,xh(tj)),bj,n+1=hqq((n+1-j)q-(n-j)q).
The error estimate ise=max|x(tj)-xh(tj)|=O(hp) (j=0,1,2,…,N), in whichp=min(2,1+q).


## 2.1. Basic Definitions

In what follows, the definitions of fractional derivative are introduced in brief. For an arbitrary real numberq, the Riemann-Liouville and Caputo fractional derivatives are defined, respectively, as
(1)aDtqf(t)=1Γ(m-q)dmdtm∫at‍(t-τ)m-q-1f(τ)dτxxxxxxxxxxxxxxxxxxxxxxx(m=[q]+1),(2)aCDtqf(t)=1Γ(m-q)∫at‍f(m)(τ)(t-τ)q+1-mdτ00000000000000000000(m=[q]+1),
where[q]stands for the integer part ofq; that is,mis the first integer which is not less thanq.Note that the main advantage of Caputo approach is that the initial conditions for fractional differential equations with Caputo derivatives take on the same form as for integer-order differential equations. Some other important properties of the fractional derivatives and integrals can be found in [1–3].


## 2.2. Numerical Method

According to [18, 19], the predictor-corrector scheme for fractional-order system is derived. This scheme is the generalization of Adams-Bashforth-Moulton scheme. First, the following fractional-order differential equation involving Caputo derivative
(3)dqxdtq=f(t,x),0≤t≤T,x(k)(0)=x0(k),k=0,1,2,…,m-1,
is equivalent to the Volterra integral equation (see [18, Lemma  2.1])
(4)x=∑k=0⌈q⌉-1‍x0(k)tkk!+1Γ(q)∫0t‍f(τ,x)(t-τ)1-qdτ.
Seth=T/N, tn=nh, n=0,1,2,…,N, and system (4) can be discretized as
(5)xh(tn+1)=∑k=0⌈q⌉-1‍x0(k)tn+1kk!+hqΓ(q+2)f(tn+1,xhp(tn+1))+hqΓ(q+2)∑j=0n‍aj,n+1f(tj,xh(tj)),
where
(6)aj,n+1={nq+1-(n-q)(n+1)q,j=0,(n-j+2)q+1+(n-j)q+1-2(n-j+1)q+1,1≤j≤n,1,j=n+1,xhp(tn+1)=∑k=0⌈q⌉-1‍x0(k)tn+1kk!+1Γ(q)∑j=0n‍bj,n+1f(tj,xh(tj)),bj,n+1=hqq((n+1-j)q-(n-j)q).
The error estimate ise=max|x(tj)-xh(tj)|=O(hp) (j=0,1,2,…,N), in whichp=min(2,1+q).


## 3. Fractional Rössler Toroidal System

The Rössler toroidal system, which was firstly introduced in [20], is defined by
(7)dxdt=-y-z,dydt=x,dzdt=ay-ay2-bz,
where a>0 and b>0 are real parameters.For what follows, it is useful to transform this system in the equivalent third-order differential equation:(8)d3ydt3+bd2ydt2+dydt+(a+b)y-ay2=0.
Now we introduce fractional order into the system (8) and obtain the corresponding fractional-order Rössler toroidal model:
(9)d2+qydt2+q+bd2ydt2+dydt+(a+b)y-ay2=0,
where a>0, b>0, 0<q<1, and dq/dtq is in sense of the Caputo fractional derivative defined in (2) with the initial time t=0.According to the Poincare-Bendixson theorem [21], ifqis very close to zero, then system (9) behaves like a two-dimensional system and cannot be chaotic. However, ifqtends to 1, system (9) dynamically approaches a third-order ordinary differential equation which can be chaotic. So, it is interesting to analyze the conditions under which chaos can occur in system (9) whenq∈(0,1).First, system (9) can be represented as a nonlinear feedback structure, made of a linear dynamic partGand a nonlinear dynamic partN(shown in Figure 1).Figure 1
Nonlinear feedback structure of the system (9).The linear subsystemGis modeled by the following transfer function:
(10)G(s)=1s2+q+bs2+s+a+b,
and the nonlinear subsystemNis described by a real function
(11)f(y)=-ay2.
Suppose that the hypothetical solution of system (9) is represented by
(12)y0(t)=A+Bsinωt,
whereB>0andω>0. Also, the nonlinear subsystemf(y0)can be approximated by
(13)f(y0)≈N0(A,B)A+N1(A,B)Bsinωt,
whereN0andN1are represented by
(14)N0(A,B)=12πA∫-ππ‍f(A+Bsinωt)dωt=-aA-aB22A,(15)N1(A,B)=1πB∫-ππ‍f(A+Bsinωt)sinωtdωt=-2aA.
Note that this approximation can be satisfied, if the linear partG(jω)has the property of a low pass filter. From Figure 1, it follows
(16)-(G(0)N0(A,B)A+G(jω)N1(A,B)Bsinωt)=A+Bsinωt.
Thus, the following equations are obtained:
(17)A(1+G(0)N0(A,B))=0,(18)1+G(jω)N1(A,B)=0.


## 4. The Conditions for Chaos to Occur

Based on the Genesio-Tesi conjecture [22], the predicted existence of chaos can be synthesized, with the following four conditions satisfied: (i) existence of a stable limit cycley0(t)=A+Bsinωt, (ii) existence of an unstable equilibrium pointEwhich is different from the fixed point that originates the predicted limit cycle in (i), (iii) an interaction between the observed limit cycley0(t)and equilibrium pointE; that is,y0(t)=Efor somet, (iv) suitable filtering effect in the Genesio-Tesi conjecture. So, the following four steps must be performed.


### 4.1. A Stable Limit Cycle

Suppose that the predicted limit cycle is denoted byy0(t)=A0+B0sinω0t. Substituting (15) into (18) and letting the real and imaginary parts be equal to zero, respectively, their solutions can be expressed as follows:
(19)ω0=(1sin(qπ/2))1/(1+q),(20)A0=12a[-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)-bsin(qπ2)2/(1+q)+a+b-(1sin(qπ/2))(2+q)/(1+q)].
From (14) and (17), it follows
(21)B0=2(a+b)A0-2aA02a,
which implies that
(22)A0(a+b-aA0)>0.
Sincea>0andb>0, the necessary condition [17] for existence of a limit cycle is
(23)-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)-b(1sin(qπ/2))2/(1+q)+a+b>0.
Based on the approximate Loeb criterion [23], the stability of the predicted limit cycle depends on the relative location ofG(jω)with respect toωand the locus of1/N1(A(B),B)with respect toB, whereAis considered as a function ofB, that is,A(B). The predicted limit cycle is stable if the locus of1/N1(A(B),B)enters in a stable interval of the real axis by increasingBfromB0 and, vice versa, by decreasingB(see Figure6in [24]). From [17], it is sufficient to checkdN1(A(B),B)/dB∣B=B0<0.According to (14)–(18), it follows
(24)N1(A(B),B)=-2a2B2(2/G(0))-(1/G(jω))=-(2a2B2)×(2(a+b)-[-ω2+q(cos(qπ2)+jsin(qπ2))-bω2+jω+a+b(cos(qπ2)+jsin(qπ2))])-1.
Thus, from inequalities (22) and (23), it can be calculated that
(25)dN1(A(B),B)dB|B=B0=-(4a2B0)×(2(a+b)-[-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)00-bsin(qπ2)2/(1+q)+a+b(1sin(qπ/2))(2+q)/(1+q)(1/sin(qπ2))(2+q)/(1+q)(1/sin(qπ/2))(2+q)/(1+q)])-1=-2a2B0a+b-aA0<0.
Therefore, the predicted limit cycle is stable according to the Loeb criterion [23].


### 4.2. An Unstable Equilibrium Point

By solving the equation(a+b)y-ay2=0, the equilibrium points of system (9) can be derived as follows:
(26)y1=0,y2=(a+b)a.
The equilibrium pointy2=(a+b)/ais the equilibrium point of interest in this step of the described function method because the equilibrium pointy1=0originates the limit cycle described by (19)–(21). Moreover, the characteristic equation of linearized system of the system (9) around the equilibrium pointy2=(a+b)/ais given by
(27)P(λ)=λ2+q+bλ2+λ-a-b.
Sincea>0andb>0, it is easy to check thatP(0)<0andP(a+b)>0. This means that there at least exists a positive real numberλ*between0anda+b, which satisfiesP(λ*)=0. Hence, the equilibrium pointy2is unstable according to the stability theorem proved in [25].


### 4.3. Interaction between the Limit Cycle and Equilibria

In this section, an interaction between the observed limit cycley0(t)and equilibrium pointy2should be verified. This condition can be restated by the following inequality:
(28)B0≥|y2-A0|.
Because of heuristic nature of the method, condition (28) should be considered in the more indefinite form such as
(29)B0≥θ|y2-A0|,
whereθis a number close to1[22]. According to (20) and (21), condition (29) can be verified if
(30)(2+θ22-θ2)[-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)-b(1sin(qπ/2))2/(1+q)]+a+b>0,
which is stronger than the necessary condition (23). Hence, an interaction between the observed limit cycley0(t)and equilibrium pointy2is certain.


### 4.4. Filter Effect

Since the reliability of the stable limit cycle is based on the filtering hypothesis, the higher harmonic terms off(A+Bsinωt)are discarded. The amount of the neglected higher harmonics concerning the predicted limit cycle withy0(t)of frequencyωcan be expressed by the so-called distortion index
(31)Δ=∥y~0(t)-y0(t)∥2∥y0(t)∥2,
where the symbol∥·∥2stands for theL2norm with the period2π/ω andy~0(t)denotes the steady-state periodic output of the system [26]. But, it is too difficult to calculate the undistorted signaly~0(t)in the fractional-order system. As suggested in [16], the limitslimt→∞y0(t)andlimt→∞y~0(t)are used to replacey0(t)andy~0(t). Based on the final value theorem of Laplace transform, it follows
(32)limt→∞y0(t)=lims→0sℒ(y0(t))=lims→0s(A0s+B0ωs2+ω2)=A0,(33)limt→∞y~0(t)=lims→0sℒ(y0(t)2)G(s)=lims→0s(A02s+2A0B0ωs2+ω2+2B02ω2s3+4sω2)G(s)=2A02+B022a+2b.
Substituting relations (32) and (33) into (31), the distortion indexΔcan be obtained for different valuesa,b, andq. Indeed, small value ofΔimplies that the predicted limit cycle is reliable.


## 4.1. A Stable Limit Cycle

Suppose that the predicted limit cycle is denoted byy0(t)=A0+B0sinω0t. Substituting (15) into (18) and letting the real and imaginary parts be equal to zero, respectively, their solutions can be expressed as follows:
(19)ω0=(1sin(qπ/2))1/(1+q),(20)A0=12a[-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)-bsin(qπ2)2/(1+q)+a+b-(1sin(qπ/2))(2+q)/(1+q)].
From (14) and (17), it follows
(21)B0=2(a+b)A0-2aA02a,
which implies that
(22)A0(a+b-aA0)>0.
Sincea>0andb>0, the necessary condition [17] for existence of a limit cycle is
(23)-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)-b(1sin(qπ/2))2/(1+q)+a+b>0.
Based on the approximate Loeb criterion [23], the stability of the predicted limit cycle depends on the relative location ofG(jω)with respect toωand the locus of1/N1(A(B),B)with respect toB, whereAis considered as a function ofB, that is,A(B). The predicted limit cycle is stable if the locus of1/N1(A(B),B)enters in a stable interval of the real axis by increasingBfromB0 and, vice versa, by decreasingB(see Figure6in [24]). From [17], it is sufficient to checkdN1(A(B),B)/dB∣B=B0<0.According to (14)–(18), it follows
(24)N1(A(B),B)=-2a2B2(2/G(0))-(1/G(jω))=-(2a2B2)×(2(a+b)-[-ω2+q(cos(qπ2)+jsin(qπ2))-bω2+jω+a+b(cos(qπ2)+jsin(qπ2))])-1.
Thus, from inequalities (22) and (23), it can be calculated that
(25)dN1(A(B),B)dB|B=B0=-(4a2B0)×(2(a+b)-[-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)00-bsin(qπ2)2/(1+q)+a+b(1sin(qπ/2))(2+q)/(1+q)(1/sin(qπ2))(2+q)/(1+q)(1/sin(qπ/2))(2+q)/(1+q)])-1=-2a2B0a+b-aA0<0.
Therefore, the predicted limit cycle is stable according to the Loeb criterion [23].


## 4.2. An Unstable Equilibrium Point

By solving the equation(a+b)y-ay2=0, the equilibrium points of system (9) can be derived as follows:
(26)y1=0,y2=(a+b)a.
The equilibrium pointy2=(a+b)/ais the equilibrium point of interest in this step of the described function method because the equilibrium pointy1=0originates the limit cycle described by (19)–(21). Moreover, the characteristic equation of linearized system of the system (9) around the equilibrium pointy2=(a+b)/ais given by
(27)P(λ)=λ2+q+bλ2+λ-a-b.
Sincea>0andb>0, it is easy to check thatP(0)<0andP(a+b)>0. This means that there at least exists a positive real numberλ*between0anda+b, which satisfiesP(λ*)=0. Hence, the equilibrium pointy2is unstable according to the stability theorem proved in [25].


## 4.3. Interaction between the Limit Cycle and Equilibria

In this section, an interaction between the observed limit cycley0(t)and equilibrium pointy2should be verified. This condition can be restated by the following inequality:
(28)B0≥|y2-A0|.
Because of heuristic nature of the method, condition (28) should be considered in the more indefinite form such as
(29)B0≥θ|y2-A0|,
whereθis a number close to1[22]. According to (20) and (21), condition (29) can be verified if
(30)(2+θ22-θ2)[-(1sin(qπ/2))(2+q)/(1+q)cos(qπ2)-b(1sin(qπ/2))2/(1+q)]+a+b>0,
which is stronger than the necessary condition (23). Hence, an interaction between the observed limit cycley0(t)and equilibrium pointy2is certain.


## 4.4. Filter Effect

Since the reliability of the stable limit cycle is based on the filtering hypothesis, the higher harmonic terms off(A+Bsinωt)are discarded. The amount of the neglected higher harmonics concerning the predicted limit cycle withy0(t)of frequencyωcan be expressed by the so-called distortion index
(31)Δ=∥y~0(t)-y0(t)∥2∥y0(t)∥2,
where the symbol∥·∥2stands for theL2norm with the period2π/ω andy~0(t)denotes the steady-state periodic output of the system [26]. But, it is too difficult to calculate the undistorted signaly~0(t)in the fractional-order system. As suggested in [16], the limitslimt→∞y0(t)andlimt→∞y~0(t)are used to replacey0(t)andy~0(t). Based on the final value theorem of Laplace transform, it follows
(32)limt→∞y0(t)=lims→0sℒ(y0(t))=lims→0s(A0s+B0ωs2+ω2)=A0,(33)limt→∞y~0(t)=lims→0sℒ(y0(t)2)G(s)=lims→0s(A02s+2A0B0ωs2+ω2+2B02ω2s3+4sω2)G(s)=2A02+B022a+2b.
Substituting relations (32) and (33) into (31), the distortion indexΔcan be obtained for different valuesa,b, andq. Indeed, small value ofΔimplies that the predicted limit cycle is reliable.


## 5. Numerical Simulations

By fixingb=0.2, the possible chaotic parameter region for fractional-order Rössler toroidal system based on the Genesio-Tesi conjecture in theq-aplane is shown in Figure 2. In the shaded region of Figure 2, it is possible for chaos to occur.Figure 2
Possible chaotic parameter region for fractional-order Rössler toroidal system withb=0.2.Next consider the filter effect for different parameter values. Based on the above expressions (31)–(33), an interesting thing has been found that different valuesbandqcannot affect the value of distortion indexΔ. Moreover, the values of distortion indexΔfor different valuesaare shown in Table 1.Table 1
Distortion indicesΔ for different values of the parameter a.
a
0.1
0.2
0.3
0.4
0.5
Δ
9
4
2.3333
1.5
1According to Table1, the filtering effect of the system decreases by increasinga. It means that the possibility of existence of chaos increases by choosing larger a in the possible chaotic region (shown in Figure 2). To confirm the existence of chaotic behavior in fractional-order Rössler toroidal system (9), numerical simulations are carried out by applying the numerical method in Section 2.2. For instance, the pair(q,a)=(0.95,0.386)is selected as a paradigm to confirm the existence of chaos (shown in Figure 3).Figure 3
Chaotic attractor of fractional-order Rössler toroidal system with(a,b,q)=(0.386,0.2,0.95).


## 6. Conclusions

In this paper, harmonic balance based method is applied to detect the conditions indicating a chaotic behavior of fractional-order Rössler toroidal system. Based on the four steps of Genesio-Tesi conjecture, the possible chaotic parameter region is obtained. Also, numerical simulations are carried out to verify the accuracy and effectiveness of the obtained results. Finally, the finding in present paper indicates that harmonic balance method can give a useful approach to find rather simple conditions indicating a chaotic behavior of nonlinear fractional-order dynamical systems.



---

*Source: 593856-2013-12-10.xml*