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

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.

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 Section 2, 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 number 𝑞, the Riemann-Liouville and Caputo fractional derivatives are defined, respectively, as _𝑎𝐷^𝑞_𝑡𝑓(𝑡)=1Γ𝑑𝑚−𝑞𝑑𝑡^𝑚𝑚𝑡_𝑎(𝑡−𝜏)^{𝑚−𝑞−1}𝑓(𝜏)𝑞,𝑑𝜏𝑚=+1_𝐶𝑎𝐷^𝑞_𝑡𝑓(𝑡)=1Γ𝑚−𝑞^𝑡_𝑎𝑓^(𝑚)(𝜏)(𝑡−𝜏)^{𝑞+1−𝑚}𝑞,𝑑𝜏𝑚=+1 (1)(2) where[𝑞]stands for the integer part of𝑞; that is,𝑚is the first integer which is not less than𝑞.

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 𝑑^𝑞𝑥𝑑𝑡^𝑞()𝑥=𝑓𝑡,𝑥,0≤𝑡≤𝑇,_(𝑘)(0)=𝑥₀^(𝑘),𝑘=0,1,2,…,𝑚−1, (3) is equivalent to the Volterra integral equation (see [18, Lemma  2.1]) 𝑥=^{⌈𝑞⌉−1}_𝑘=0𝑥₀^(𝑘)𝑡^𝑘+1𝑘!Γ𝑞^𝑡₀𝑓(𝜏,𝑥)(𝑡−𝜏)^1−𝑞𝑑𝜏. (4) Setℎ=𝑇/𝑁, 𝑡_𝑛=𝑛ℎ, 𝑛=0,1,2,…,𝑁, and system (4) can be discretized as 𝑥_ℎ𝑡_𝑛+1=^{⌈𝑞⌉−1}_𝑘=0𝑥₀^(𝑘)𝑡^𝑘_𝑛+1+ℎ𝑘!Γ^𝑞𝑓𝑡𝑞+2^𝑛+1,𝑥^𝑝_ℎ𝑡_𝑛+1+ℎ^𝑞Γ𝑞+2^𝑛_𝑗=0𝑎_𝑗,𝑛+1𝑓𝑡_𝑗,𝑥_ℎ𝑡_𝑗, (5) where 𝑎_𝑗,𝑛+1=𝑛^𝑞+1−𝑛−𝑞(𝑛+1)^𝑞,𝑗=0,𝑛−𝑗+2^𝑞+1+𝑛−𝑗^𝑞+1−2𝑛−𝑗+1^𝑞+1,𝑥1≤𝑗≤𝑛,1,𝑗=𝑛+1,_𝑝ℎ𝑡_𝑛+1=^{⌈𝑞⌉−1}_𝑘=0𝑥₀^(𝑘)𝑡^𝑘_𝑛+1+1𝑘!Γ𝑞^𝑛_𝑗=0𝑏_𝑗,𝑛+1𝑓𝑡_𝑗,𝑥_ℎ𝑡_𝑗,𝑏_𝑗,𝑛+1=ℎ^𝑞𝑞𝑛+1−𝑗^𝑞−𝑛−𝑗^𝑞. (6) The error estimate is𝑒=max|𝑥(𝑡_𝑗)−𝑥_ℎ(𝑡_𝑗)|=𝑂(ℎ^𝑝) (𝑗=0,1,2,…,𝑁), in which𝑝=min(2,1+𝑞).

3. Fractional Rössler Toroidal System

The Rössler toroidal system, which was firstly introduced in [20], is defined by 𝑑𝑥𝑑𝑡=−𝑦−𝑧,𝑑𝑦𝑑𝑡=𝑥,𝑑𝑧𝑑𝑡=𝑎𝑦−𝑎𝑦²−𝑏𝑧, (7) where 𝑎>0 and 𝑏>0 are real parameters.

For what follows, it is useful to transform this system in the equivalent third-order differential equation: 𝑑³𝑦𝑑𝑡³𝑑+𝑏²𝑦𝑑𝑡²+𝑑𝑦+()𝑑t𝑎+𝑏𝑦−𝑎𝑦²=0. (8) Now we introduce fractional order into the system (8) and obtain the corresponding fractional-order Rössler toroidal model: 𝑑^2+𝑞𝑦𝑑𝑡^2+𝑞𝑑+𝑏²𝑦𝑑𝑡²+𝑑𝑦+()𝑑𝑡𝑎+𝑏𝑦−𝑎𝑦²=0, (9) where 𝑎>0, 𝑏>0, 0<𝑞<1, and 𝑑^𝑞/𝑑𝑡^𝑞 is in sense of the Caputo fractional derivative defined in (2) with the initial time 𝑡=0.

According to the Poincare-Bendixson theorem [21], if𝑞is very close to zero, then system (9) behaves like a two-dimensional system and cannot be chaotic. However, if𝑞tends 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) when𝑞∈(0,1).

First, system (9) can be represented as a nonlinear feedback structure, made of a linear dynamic part𝐺and a nonlinear dynamic part𝑁(shown in Figure 1).

Figure 1 

Nonlinear feedback structure of the system (9).

The linear subsystem𝐺is modeled by the following transfer function: 𝐺(𝑠)=1𝑠^2+𝑞+𝑏𝑠²,+𝑠+𝑎+𝑏 (10) and the nonlinear subsystem𝑁is described by a real function 𝑓𝑦=−𝑎𝑦². (11) Suppose that the hypothetical solution of system (9) is represented by 𝑦₀(𝑡)=𝐴+𝐵sin𝜔𝑡, (12) where𝐵>0and𝜔>0. Also, the nonlinear subsystem𝑓(𝑦₀)can be approximated by 𝑓𝑦₀≈𝑁₀(𝐴,𝐵)𝐴+𝑁₁(𝐴,𝐵)𝐵sin𝜔𝑡, (13) where𝑁₀and𝑁₁are represented by 𝑁₀()=1𝐴,𝐵2𝜋𝐴^𝜋_−𝜋𝑓(𝐴+𝐵sin𝜔𝑡)𝑑𝜔𝑡=−𝑎𝐴−𝑎𝐵²,𝑁2𝐴₁()=1𝐴,𝐵𝜋𝐵^𝜋_−𝜋𝑓(𝐴+𝐵sin𝜔𝑡)sin𝜔𝑡𝑑𝜔𝑡=−2𝑎𝐴. (14)(15) Note that this approximation can be satisfied, if the linear part𝐺(𝑗𝜔)has the property of a low pass filter. From Figure 1, it follows −𝐺(0)𝑁₀(𝐴,𝐵)𝐴+𝐺𝑗𝜔𝑁₁(𝐴,𝐵)𝐵sin𝜔𝑡=𝐴+𝐵sin𝜔𝑡. (16) Thus, the following equations are obtained: 𝐴1+𝐺(0)𝑁₀()𝑁𝐴,𝐵=0,1+𝐺𝑗𝜔₁(𝐴,𝐵)=0. (17)(18)

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 cycle𝑦₀(𝑡)=𝐴+𝐵sin𝜔𝑡, (ii) existence of an unstable equilibrium point𝐸which is different from the fixed point that originates the predicted limit cycle in (i), (iii) an interaction between the observed limit cycle𝑦₀(𝑡)and equilibrium point𝐸; that is,𝑦₀(𝑡)=𝐸for some𝑡, (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 by𝑦₀(𝑡)=𝐴₀+𝐵₀sin𝜔₀𝑡. Substituting (15) into (18) and letting the real and imaginary parts be equal to zero, respectively, their solutions can be expressed as follows: 𝜔₀=1sin𝑞𝜋/2^1/(1+𝑞),𝐴₀=1−12𝑎sin𝑞𝜋/2^{(2+𝑞)/(1+𝑞)}cos𝑞𝜋2−𝑏sin𝑞𝜋2^2/(1+𝑞)+𝑎+𝑏. (19)(20) From (14) and (17), it follows 𝐵₀=2(𝑎+𝑏)𝐴₀−2𝑎𝐴²₀𝑎, (21) which implies that 𝐴₀𝑎+𝑏−𝑎𝐴₀>0. (22) Since𝑎>0and𝑏>0, the necessary condition [17] for existence of a limit cycle is −1sin𝑞𝜋/2^{(2+𝑞)/(1+𝑞)}cos𝑞𝜋21−𝑏sin𝑞𝜋/2^2/(1+𝑞)+𝑎+𝑏>0. (23) Based on the approximate Loeb criterion [23], the stability of the predicted limit cycle depends on the relative location of𝐺(𝑗𝜔)with respect to𝜔and the locus of1/𝑁₁(𝐴(𝐵),𝐵)with respect to𝐵, where𝐴is considered as a function of𝐵, that is,𝐴(𝐵). The predicted limit cycle is stable if the locus of1/𝑁₁(𝐴(𝐵),𝐵)enters in a stable interval of the real axis by increasing𝐵from𝐵₀ and, vice versa, by decreasing𝐵(see Figure6in [24]). From [17], it is sufficient to check𝑑𝑁₁(𝐴(𝐵),𝐵)/𝑑𝐵∣_𝐵=𝐵0<0.

According to (14)–(18), it follows 𝑁₁(𝐴(𝐵),𝐵)=−2𝑎²𝐵²(2/𝐺(0))−1/𝐺𝑗𝜔=−2𝑎²𝐵²×2(𝑎+𝑏)−−𝜔^2+𝑞cos𝑞𝜋2+𝑗sin𝑞𝜋2−𝑏𝜔²+𝑗𝜔+𝑎+𝑏⁻¹. (24) Thus, from inequalities (22) and (23), it can be calculated that 𝑑𝑁₁(𝐴(𝐵),𝐵)𝑑𝐵_𝐵=𝐵0=−4𝑎²𝐵₀×2()−−1𝑎+𝑏sin𝑞𝜋/2^{(2+𝑞)/(1+𝑞)}cos𝑞𝜋2−𝑏sin𝑞𝜋2^2/(1+𝑞)+𝑎+𝑏⁻¹=−2𝑎²𝐵₀𝑎+𝑏−𝑎𝐴₀<0. (25) Therefore, the predicted limit cycle is stable according to the Loeb criterion [23].

4.2. An Unstable Equilibrium Point

By solving the equation(𝑎+𝑏)𝑦−𝑎𝑦²=0, the equilibrium points of system (9) can be derived as follows: 𝑦₁=0,𝑦₂=(𝑎+𝑏)𝑎. (26) The equilibrium point𝑦₂=(𝑎+𝑏)/𝑎is the equilibrium point of interest in this step of the described function method because the equilibrium point𝑦₁=0originates the limit cycle described by (19)–(21). Moreover, the characteristic equation of linearized system of the system (9) around the equilibrium point𝑦₂=(𝑎+𝑏)/𝑎is given by 𝑃(𝜆)=𝜆^2+𝑞+𝑏𝜆²+𝜆−𝑎−𝑏. (27) Since𝑎>0and𝑏>0, it is easy to check that𝑃(0)<0and𝑃(𝑎+𝑏)>0. This means that there at least exists a positive real number𝜆^∗between0and𝑎+𝑏, which satisfies𝑃(𝜆^∗)=0. Hence, the equilibrium point𝑦₂is 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 cycle𝑦₀(𝑡)and equilibrium point𝑦₂should be verified. This condition can be restated by the following inequality: 𝐵₀≥𝑦₂−𝐴₀. (28) Because of heuristic nature of the method, condition (28) should be considered in the more indefinite form such as 𝐵₀≥𝜃𝑦₂−𝐴₀, (29) where𝜃is a number close to1[22]. According to (20) and (21), condition (29) can be verified if 2+𝜃²2−𝜃²−1sin𝑞𝜋/2^{(2+𝑞)/(1+𝑞)}cos𝑞𝜋21−𝑏sin𝑞𝜋/2^2/(1+𝑞)+𝑎+𝑏>0, (30) which is stronger than the necessary condition (23). Hence, an interaction between the observed limit cycle𝑦₀(𝑡)and equilibrium point𝑦₂is certain.

4.4. Filter Effect

Since the reliability of the stable limit cycle is based on the filtering hypothesis, the higher harmonic terms of𝑓(𝐴+𝐵sin𝜔𝑡)are discarded. The amount of the neglected higher harmonics concerning the predicted limit cycle with𝑦₀(𝑡)of frequency𝜔can be expressed by the so-called distortion index 𝑦Δ=⁰(𝑡)−𝑦₀(𝑡)₂𝑦₀(𝑡)₂, (31) where the symbol‖⋅‖₂stands for the𝐿₂norm with the period2𝜋/𝜔 and𝑦₀(𝑡)denotes the steady-state periodic output of the system [26]. But, it is too difficult to calculate the undistorted signal𝑦₀(𝑡)in the fractional-order system. As suggested in [16], the limitslim_𝑡→∞𝑦₀(𝑡)andlim_𝑡→∞𝑦₀(𝑡)are used to replace𝑦₀(𝑡)and𝑦₀(𝑡). Based on the final value theorem of Laplace transform, it follows lim_𝑡→∞𝑦₀(𝑡)=lim_𝑠→0𝑠ℒ𝑦₀(𝑡)=lim_𝑠→0𝑠𝐴₀𝑠+𝐵₀𝜔𝑠²+𝜔²=𝐴₀,lim_𝑡→∞𝑦₀(𝑡)=lim_𝑠→0𝑠ℒ𝑦₀(𝑡)²𝐺(𝑠)=lim_𝑠→0𝑠𝐴²₀𝑠+2𝐴₀𝐵₀𝜔𝑠²+𝜔²+2𝐵²₀𝜔²𝑠³+4𝑠𝜔²𝐺(𝑠)=2𝐴²₀+𝐵²₀.2𝑎+2𝑏 (32)(33) Substituting relations (32) and (33) into (31), the distortion indexΔcan be obtained for different values𝑎,𝑏, and𝑞. Indeed, small value ofΔimplies that the predicted limit cycle is reliable.

5. Numerical Simulations

By fixing𝑏=0.2, the possible chaotic parameter region for fractional-order Rössler toroidal system based on the Genesio-Tesi conjecture in the𝑞−𝑎plane 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 with𝑏=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 values𝑏and𝑞cannot affect the value of distortion indexΔ. Moreover, the values of distortion indexΔfor different values𝑎are shown in Table 1.

Table 1 

Distortion indices Δ for different values of the parameter 𝑎.

According to Table 1, the filtering effect of the system decreases by increasing𝑎. It means that the possibility of existence of chaos increases by choosing larger 𝑎 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(𝑞,𝑎)=(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 (𝑎,𝑏,𝑞)=(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.

Acknowledgment

This work was supported by the Science and Technology Innovation Projects of Education Bureau of Guangdong Province (no. 2012KJCX0073).