The Routh Hurwitz Criteria for Studying The Stability and Bifurcation in Multispiral Chua Chaotic Attractor

MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX ISSN: 1803-3814 (Printed), 2571-3701 (Online) https://doi.org/10.13164/mendel.2023.1.071 The Routh Hurwitz Criteria for Studying The Stability and Bifurcation in Multispiral Chua Chaotic Attractor Malika Belouerghi1, Tidjani Menacer1,  , Ren´e Lozi2 1Laboratory of Applied Mathematics, University Mohamed Khider of Biskra, P.O. Box 145, Biskra 07000, Algeria 2University Cˆote d’Azur, CNRS, LJAD, France, Parc Valrose 06108, Nice C´edex 02 m.b[email protected], [email protected]  , [email protected] Abstract This article discusses the multispiral Chua Chaotic attractor’s hidden bifurcations that are generated by the sine function. The number of spirals (also known as a multiscroll attractor) that are controlled by the integer parameter ccan be used to describe the basic shape of chaotic attractors. Since this parameter is an integer, increasing it by one does not allow the observation of bifurcations from nto n+ 2 spirals. The method of hidden bifurcations, however, enables the observation of such bifurcations by adding a real parameter ε. Chaotic attractors with either an even or an odd number of spirals are visible along the marked paths of bifurcation. Moreover, this additional hidden parameter allows finding the bifurcation of the multispiral Chua attractor from a stable state to a chaotic state. Furthermore, the Routh-Hurwitz criteria are used to study the stability of the original equilibrium point of the Chua attractor. Keywords: Bifurcation, Chua’s Circuit, Hidden Attractor, Hidden Bifurcation, Stability. Received: 17 April 2023 Accepted: 25 June 2023 Online: 30 June 2023 Published: 30 June 2023 1 Introduction In the framework of the qualitative theory of dynamic systems, the name ”bifurcation” was coined by Henri Poincar´e in 1885 [24]. For a given value of a parameter, a bifurcation is seen when the number of solutions to a family of differential equations increases from one to two, like the pitchfork of a branch of a tree [25], [18], or when the topological structure of the solution is changed from steady state to periodic function (Hopf bifurcation [21]), or from periodic to quasi-periodic function (secondary Hopf bifurcation). Lorenz [17], extended the scope of the bifurcation theory by identifying the first chaotic strange attractor in 1963. Chua invented the first differential equation modeling a real electronic device system with a chaotic asymptotic attractor (hence the strange attractor) while he was invited Professor at Matsumoto’s laboratory, Waseda University, Tokyo [4]. Nowadays, Chua’s attractor is widely used, due to both its realizations: electronic circuit or its mathematical model. The electronic circuit and the system of differential equations may be combined to reach multiple objectives by Duan et al. [7]. One can consider Chua’s circuit as the simplest electronic circuit presenting chaos and possessing a highly interesting dynamical behavior. This was checked in many laboratory experiments for Zhong and Ayrom [31], computer simulations for Matsumoto [22] and rigorously done mathematics [5], [19], [3]. Chua’s circuit is, surely the most extensively studied chaotic electrical system. It has the extraordinary feature to be able to generate a large variety of dynamical behaviours with just a few self-electronic components. In particular, it consists in two linear capacitors (a resistor and an inductor) and one non-linear resistor known as Chua’s diode. By appropriately choosing these components, the circuit becomes chaotic, and the trajectories tend to a limit set called strange attractor. The best-known attractor generated by Chua’s circuit is called double scroll [28], [2]. It had been generalized in several ways, replacing continuous piecewise-linear functions with some smooth functions, such as polynomials [10], [26], [9], and [8], etc. The modified family defined in [27], [30], generates attractors with even or odd number of scrolls. The n-double scrolls correspond to the even case (2n-scroll attractors in the new family) with Chua’s double scroll being a 2-scroll attractor. In the work of Tang et al. [29], it is demonstrated that the n-scroll attractors can be generated using a simple sine or a cosine function. In 2011, Leonov et al. [12], [11] introduced a new classification of attractors of non-linear dynamical systems; in which attractors are dispatched in self-excited or hidden attractors. The first ones can be localized numerically via a standard computational procedure. 71 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX After a transient process, a trajectory, starts from a point of unstable manifold in a neighborhood of equilibrium to reach a state of oscillation. Hence, one can easily identify it. In contrast, for hidden attractors, a basin attraction does not intersect with small neighborhoods of equilibria. To localize them, it is necessary to develop special procedures. Therefore, there are no transient processes leading to such attractors. The name hidden comes from this property. Hidden attractors play a key role in engineering applications (such as air crafts control systems, electrical machines, etc.) as they allow unexpected and potentially catastrophic responses to perturbations in a structure. An effective method for the numerical localization of hidden attractors in multidimensional dynamical systems has been proposed by the previously cited authors [15], [14]. Their method is based on numerical continuation: a sequence of similar systems is constructed such that for the first (starting) system the initial data for numerical computation of possible oscillating solution (starting oscillation) can be obtained analytically and its transformation from one system to another is followed numerically. One can find the first example of a hidden chaotic strange attractor in Chua’s attractor [13]. The modified multispiral Chua system with a sine function is presented in [29]. In [23], Menacer et al. using the Kuznetsov and Leonov method in another way, found a sequence of hidden bifurcations. In this system, the parameter cgoverning the number of spirals is an integer; hence it is not possible to vary it continuously. Therefore, it is not possible to observe classical bifurcations of attractors from nto n+ 2 spirals when the parameter cchanges, because the theory of bifurcation concerns the variation of solutions with respect to a continuous parameter. Also, it is not possible to use non-integer real values for c. To overcome this obstacle, a new bifurcation parameter ε is introduced by the hidden attractor search method. In this method, the parameter εcontrols the birth of an attractor and the increase of the number of spirals until the final number corresponding to the value of the integer c(in this article cis set to 11 and 12) is reached. In this paper, a second role has been assigned to the parameter ε: the control of the nature of solutions (transition from stability to chaos). The Routh-Hurwitz criterion given in [1] is used for the study of the stability of the origin equilibrium point with respect to the parameter εfor the multispiral chaotic Chua system considered in this paper. The organization of the paper is as follows. Attractors with multiSpiral in Chua’s sine function system are recalled in Section 2. In Section 3, the hidden bifurcations uncovering method is explained. In the chosen example every attractor located along the hidden bifurcation path have an odd number of spirals. In Section 4, one studies of the stability of the origin equilibrium point E0with respect to εusing the Routh-Hurwitz Criteria. Some numerical results are presented in Section 5. Finally, in Section 6, a brief conclusion is drawn. 2 Attractors With Multi-Spiral in Chua’s Sine Function System In [8], the use of cellular neural networks with a piecewise-linear output function to create Complex chaotic attractors with n-double scrolls was highlighted. Another simpler mechanism for generating n-scroll attractors was introduced by Tang in [29], using sine or cosine functions. Since then, in the last decades, multi-scroll chaotic attractors generation has been extensively studied due to their promising applications in various real-world chaos-based technologies including secure and digital communications, random bit generation, etc. The system of differential equations, describing the behavior of Chua’s circuits, that we consider in this article is three-dimensional with a combination of piecewise-linear and sinusoidal nonlinearity [4], [29], [6](see Fig. 1)    ˙x(t) = α(y(t)−f(x(t))), ˙y(t) = x(t)−y(t) + z(t), ˙z(t) = −βy(t), (1) where ˙x(t) = dx(t) dt ,˙y(t) = dy(t) dt ,˙z(t) = dz(t) dt . f(x(t)) =    bπ 2a(x(t)−2ac) if x(t)≥2ac, −bsin(πx(t) 2a+d) if −2ac < x(t)<2ac, bπ 2a(x(t)+2ac) if x(t)≤ −2ac, (2) Figure 1: Proposed sine function f(x) with parameters values a= 2, b = 0.2, c = 12, d =π. The goal of this article is to analyze the general shape of the attractors and their global geometric features, 72 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX Belouerghi-g2igHX,g The Routh Hurwitz Criteria for Studying The Stability and Bifurcation in ... which can be described in terms of the number of spirals. In this article the real parameter values have been fixed to α= 11, β = 15, a = 2, b = 0.2 for which topologically chaotic attractors have been found. They are equivalent to the attractors found in electronic device in [29]. The parameter c is an integer that governs the number nof spirals according to the following formula n=c+ 1 (3) and dis chosen such that d=πif cis even. 0 if cis odd. (4) A straightforward computation gives the equilibrium points of (1) which are (−xeq,0, xeq) with xeq = 2ak, k= 0,±1,±2, ..., ±c[29]. In the case c= 12, one obtains 13 spiral attractors (see Fig. 2). For a complete study of this model with other values of csee [23]. Figure 2: The 13-spiral attractor generated by Equation (1) and (2) for c= 12.(a) 3−dimensional figure, (b) Projection into the plane (x, y). 3 Hidden Bifurcations Uncovering Method Menacer et al. [23] technique for finding hidden bifurcations is based on Leonov and Kuznetsov’s fundamental concept for searching hidden attractors (i.e. homotopy and numerical continuation, see Appendix 6). While keeping cconstant, a new bifurcation parameter εis introduced. This method is briefly recalled and applied to (1)-(2) in this section. The value of parameters are fixed to α= 11, β = 15, a = 2, b = 0.2, c = 10, d=π. 3.1 Hidden Bifurcations The system(1)-(2) is rewritten in the Lure’s form [20] ˙ U=MU +qH(rTU), U = (x, y, z)∈R3.(5) where M=  0α0 1−1 1 0−β0  , q =  −α 0 0  , r =  1 0 0   and H(σ)=f(σ). In order to transform system (5) into the form similar to the system (22, see AppendixA), a new coefficient χ, and a small parameter ε, are introduced as follows ˙ U=M0U+qεg0(rTU), U = (x, y, z)∈R3.(6) where M0=M+χqrT=  −αχ α 0 1−1 1 0−β0  , and g(σ) = H(σ)−χσ =−α(f(σ)−χσ). Therefore, (6) is written as    ˙x(t) = −α(χx(t)−y(t)) + εg(x(t)) ˙y(t) = x(t)−y(t) + z(t) ˙z(t) = −βy(t) (7) Then considering the transfer function W(m) (21, see AppendixA) and solving the equations Im W(iω0) = 0 and χ=−Re W(iω0)−1, one obtains ω0= 2.1018 and χ= 0.03796. Using the nonsingular linear transformation U=SZ defined in the AppendixA, the system (7) is reduced to the form similar to equation (23) ˙ Z=AZ +Bεg0(CTZ), Z = (z1, z2, z3)∈R3(8) where A=  0−ω00 −ω00 0 0 0 −d1  , B =  b1 b2 1  , C =  1 0 −h   The transfer function WA(m) of system (8) reads WA(m) = −b1m+b2ω0 m2+ω2 0 +h m+d1 . Then, using the equality of transfer of both functions WA(m) and WM0(m) = rT(M0−mI)−1qto systems (7) and (8) one obtains the following relations χ=a+ω2 0−β α; d1=α+ω2 0−β+ 1; h=α(β−d1+d2 1) ω2 0+d2 1 ; b1=α(β−ω2 0−d1) ω2 0+d2 1 ; b2=α(1 −d1)(ω0−βd1) ω2 0(ω2 0+d2 1). Remark 1. The matrix S is defined by the following relations easily computed A=S−1M0S, B =S−1q, CT=rTS, S=   1 0 −h χ−ω0 −α hd−hαχ α −β α χβ ω0 βh(d−αχ) dα   .(9) 73 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX 3.2 Hidden Bifurcations Route In this section, we present a numerical example of a hidden bifurcations route. For the values of parameters fixed at the beginning of this section, we obtain χ= 0.03796, d1= 1.4176, h = 26.686, b1= 15.686, b2= 3.1003, therefore, the matrix (9), is S=  1 0 −26.686 0.03795 −0.19107 −2.4264 −1.3636 −0.27084 −25.674   Using theorem 1 of the AppendixA, for εsmall enough, one obtains the initial condition U0(0) = SZ(0) = S  τ0 0 0  =  τ0s11 τ0s21 τ0s31  ,(10) Using the notation of Section 2, one obtains for the determination of the initial condition of starting solution for the multistage localization procedure, Chua system x(0) = τ0, y(0) = τ0χ, z(0) = −τ0 β α,(11) Then, the localization procedure described in Appendix A is applied to the system (5)-(7). The starting frequency ω0and the coefficient of harmonic linearization χhave been already computed in Sec. 3.1. Equation(11), allows to obtain the initial conditions for the first step. Then, the solutions of system (7) with the nonlinearity εg(x) = ε(H(x)−χx) are computed, by increasing sequentially εfrom the value ε= 0.1 to ε= 1, with step size 0.1, except between ε= 0.9 and ε= 1 , where one uses 0.001 as increasing step. . All the points of the stable periodic solution U1(t) corresponding-to ε= 0.1 belong to the domain of attraction of the stable periodic solution U2(t) corresponding to ε= 0.2. This stable periodic solution U2(t) is then simply obtained numerically by solving system (5) with ε= 0.2 and U1(tmax) as initial point, where tmax represents the last value of time of integration after discarding the transitory regime. The same numerical procedure is reitered by increasing the ε-value, to obtain the next periodic solutions U3(t), U4(t), ..., Uj(t), ... corresponding to the next values of ε. When ε= 0.8 the first chaotic solution one scroll is found. The initial conditions for recovering the solutions for increasing values of εas shown in the Table 1. Table 1: Initials conditions according to the values of ε. ε Xj(0) x0y0z0 0.1 U1(0) = U0(tmax) 38.3269 -0.0257 -37.9339 0.2 U2(0) = U1(tmax) -3.7487 0.1619 -5.0818 0.3 U3(0) = U2(tmax) 3.3406 -0.2062 -4.9772 0.4 U4(0) = U3(tmax) -3.6009 -0.2981 4.7285 0.5 U5(0) = U4(tmax) 3.4797 -0.1484 -5.0577 0.6 U6(0) = U5(tmax) 3.6919 0.0202 -5.1617 0.7 U7(0) = U6(tmax) 4.1512 -0.0232 -5.5316 0.8 U8(0) = U7(tmax) 1.3915 -0.5758 -2.6930 0.86 U9(0) = U8(tmax) -0.4213 -0.3265 0.2085 0.978 U10(0) = U9(tmax) 1.8233 0.2554 -2.6074 0.989 U11(0) = U10(tmax) 1.4667 0.7801 -0.6740 0.993 U12(0) = U11(tmax) 3.2919 0.5697 -3.9526 0.995 U13(0) = U12(tmax) 15.5612 -0.3315 -16.0715 1U14(0) = U13(tmax) -33.7244 -0.0092 35.1636 Using these initial conditions, we get the solutions U8(t) (Fig. 3(a)) with one spiral to U13(t) (Fig. 4(c)) with 11 spirals. In each figure, there are a different odd number of spirals in the attractor. The number of spirals increases by 2 at each step as displayed in Table 2from 1 to 11 spirals. The values of εin this table are exactly the values of bifurcation points with respect to ε. Therefore Fig. 3and Fig. 4display the hidden bifurcations of the multi-spiral Chua’s attractor. Table 2: Values of the parameter εat the bifurcation points for c= 10 (11 spirals). ε 0.80 0.86 0.978 0.989 0.993 0.995 U8(0) U9(0) U10(0) U11(0) U12(0) U13(0) 1 spiral 3 spirals 5 spirals 7 spirals 9 spirals 11 spirals 74 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX Belouerghi-g2igHX,g The Routh Hurwitz Criteria for Studying The Stability and Bifurcation in ... Figure 3: The increasing number of spirals of system (6) with respect to increasing εvalues. (a) 1-spiral for ε= 0.80, (b) 3-spirals for ε= 0.86, (c) 5-spirals for ε= 0.978. Figure 4: The increasing number of spirals of system (6) with respect to increasing εvalues. (a) 7-spirals for ε= 0.989, (b) 9-spirals for ε= 0.993, (c) 11-spirals for ε= 0.995. 4 Stability of the Origin Equilibrium Point E0With Respect to ε In this section, we study the stability of the equilibrium point E0with respect to epsilon of the system (7) using Routh-Hurwitz’s conditions [1]. In [23] Menacer et al. introduce the concept of hidden bifurcation in the system of Chua adding a new parameter epsilon, which controls the number of spirals. When the value of εincreases from 0 to 1 the number of scrolls decreases. Let E(xe, ye, ze) be an equilibrium solution of the general three-dimensional system:    ˙x(t) = Q(x(t), y(t), z(t)) ˙y(t) = R(x(t), y(t), z(t)) ˙z(t) = V(x(t), y(t), z(t)) (12) The eigenvalues equation corresponding to this equilibrium point is given by the following polynomial: P(λ) = λ3+a1λ2+a2λ+a3.(13) Using the result of the Routh-Hurwitz conditions, where the necessary and sufficient condition for the equilibrium point Eto be locally asymptotically stable is a1>0, a3>0 and a1×a2−a3>0. 75 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX In this section, the parameter cis constant, and the bifurcation is studied with respect to parameter ε, and the values of parameters are α= 11, β = 15, a = 2, b= 0.2. 4.1 Stability of the Origin Equilibrium Point E0 The origin E0(0,0,0) is an equilibrium point of System (7) independently to epsilon. In the rest of this article, we consider both cases c= 11 with d= 0, and c= 12 with d=π. For c= 11 and d= 0, the Jacobian matrix for evaluated at the equilibrium point E0(0,0,0) is: JE0=  −αχ +αε(χ+πb 2a)α0 1−1 1 0−β0   =  −0.41756 + ε(0.41756 + 6.908 4) 11 0 1−1 1 0−15 0   Its characteristic polynomial is: P(λ) = λ3+ (1.4176 −2.1446ε)λ2 +(4.4176 −2.1446ε)λ+ (6.2634 −32.168ε). According to Routh-Hurwitz conditions, the necessary and sufficient condition, for the equilibrium point E0to be stable is 0.00005139 <ε<0.1947. Proof. a1= 1.4176 −2.1446ε > 0 =⇒ε < 0.66101; a3= 6.2634 −32.168ε > 0 =⇒ε < 0.19472; a1×a2−a3= 4.5993ε2+ 19.655ε−1.6102 ×10−3>0 =⇒ε < −4.2733 or ε > 5.1399 ×10−5. For c= 12 and d=π, the Jacobian matrix evaluated at the equilibrium point E0is: JE0=  −αχ +αε(χ−πb 2a)α0 1−1 1 0−β0   =  −0.41756 + ε(0.41756 −6.908 4) 11 0 1−1 1 0−15 0   Its characteristic polynomial is: P(λ) = λ3+ (1.4176 + 1.3094ε)λ2 +(1.3094ε+ 4.4176)λ+ (19.642ε+ 6.2634). According to Routh-Hurwitz conditions, the equilibrium point E0is unstable. Proof. a1= 1.4176 + 1.3094ε > 0 =⇒ε > −1.0826; a3= 19.642ε+ 6.2634 >0 =⇒ε > −0.3188; a1×a2−a3= 1.7145ε2−12.001ε−1.012 ×10−3>0 =⇒ε < −8.4175 ×10−5or ε > 6.998. because there is a contradiction with 0 <ε<1.Therefore, the limit cycle is unstable. Special case ε= 0 : the system (7) becomes linear:    ˙x(t) = −αχx(t) + αy(t) ˙y(t) = x(t)−y(t) + z(t) ˙z(t) = −βy(t) (14) The Jacobian matrix is: JE0=  −αχ +αε(χ+πb 2a)α0 1−1 1 0−β0   =  −0.41756 + 11ε(0.03796 + 0.628 4) 11 0 1−1 1 0−15 0  . The characteristic polynomial is given by: P(λ) = λ3+ (1.4176 −2.1446ε)λ2 +(4.4176 −2.1446ε)λ+ (6.2634 −32.168ε). One has a1= 1.4176 −2.1446ε > 0 =⇒ε < 0.66101; a3= 6.2634 −32.168ε > 0 =⇒ε < 0.19472; a1×a2−a3= 4.5993ε2+ 19.655ε−1.6102 ×10−3>0 =⇒ε < −4.2733 or ε > 5.1399 ×10−5. Therefore, the Routh-Hurwitz conditions are not verified, because ε= 0, and then the system (7) is unstable. Fig. 5displays the corresponding unstable limit cycle . Figure 5: The attractor of the system (7) where ε= 0: limit cycle unstable. 76 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX Belouerghi-g2igHX,g The Routh Hurwitz Criteria for Studying The Stability and Bifurcation in ... 4.2 Other Equilibrium Points The equilibrium point Eof system (7) is obtained solving,    ˙x(t)=0 ˙y(t) = 0 ˙z(t) = 0 ⇔   −α(χx(t)−y(t)) + εg(x(t)) = 0; x(t)−y(t) + z(t) = 0; −βy(t) = 0; (15) In addition to the origin equilibrium point E0(0,0,0), there are several other equilibrium points: Ek+(xeq,0,−xeq) and Ek−(−xeq,0, xeq). The solution of (15) are : •If x≥2ac or x≤ −2ac xeq =2εbπac 2a(εχ −χ)−εbπ . For the values of parameters above, with c= 11; one finds: xeq =±27.6460ε 0.1518(ε−1) −0.6283ε. •If −2ac < x < 2ac; one obtains: −α(χx(t)−y(t)) + εg(x(t)) = 0; (16) that is −α(χx(t)−εb sin(πx(t) 2a)−εχx(t)) = 0. Case ε= 1 : the system (7) is the original system (1). In this case, in addition to the origin E0(0,0,0) the other equilibrium points of this system are (xeq,0,−xeq), with xeq = 2ak and k=±1,±2, ..., ±c [29]. Case 0 < ε < 1 : in addition to the origin E0(0,0,0) the other equilibrium points cannot be obtained using closed formula. It is possible to compute them numerically solving ˙x(t) = 0 in equation (16). The number of equilibrium points set are shown in Table. 3 and Fig. 6, Fig. 7replacing the values of parameters above and c= 11 in the interval [−44,44] Table 3: Number of equilibrium points in the equation (16) for different values of εand c= 11. Values of ε0.30 0.50 0.70 0.80 0.90 0.95 1 Number of equilibrium points 2 2 6 10 22 22 22 77 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX Figure 6: Number of equilibrium points in the equation (16) for different values of εand c= 11. (a) for ε= 0.30, (b) for ε= 0.50, (c) for ε= 0.70, (d) for ε= 0.80. Figure 7: Number of equilibrium points in the equation (16) for different values of εand c= 11. (a) for ε= 0.90, (b) for ε= 0.95, (c) for ε= 1. 5 Numerical Analysis of Bifurcations 5.1 Case c= 11 In this section, a numerical analysis of the bifurcations of this system is done for c = 11. In this case one can observe only an even number of scrolls. The value of εis increased from 0 to 1 . For each value of εthe same initial conditions [0.001,0,0.001] are chosen. The bifurcation diagram with respect to εis given in Fig. 8, Fig. 9and Fig. 10.Figure 8: Bifurcation diagram with respect to εof the component x, with c=11. 78 MENDEL — Soft Computing Journal, Volume 29, No.g1, June 2023, Brno, Czech RepublicX Belouerghi,g2igHX,g The Routh Hurwitz Criteria for Studying The Stability and Bifurcation in ... Figure 9: Bifurcation diagram with respect to εof the component z, with c=11. Figure 10: Bifurcation diagram with respect to εof the component y, with c=11. In order to highlight the symmetry of the bifurcation diagrams versus the components xand z, Fig. 11 displays the superimposition of both Fig. 8and Fig. 9. Figure 11: Bifurcation diagram with respect to εof the superimposed components xand zwith c=11. In Figs. 12 and 13 one displays some corresponding attractors. •For ε < 0.195, the equilibrium point E0is a locally asymptotically stable focus. •For 0.195 < ε ≤0.43 the fixed point E0becomes unstable, a period-one limit cycle appears, as shown in Fig. 12a. •When ε≈0.46, a new bifurcation occurs, the period-one limit cycle becomes unstable and a period-two limit cycle appears (Fig. 12b). •For ε≈0.495, a period-4 limit cycle appears through a new bifurcation, as shown in Fig. 12c, followed by a bifurcation to a period-8 limit cycles at ε≈0.501 Fig. 12d). This doubling period bifurcation process continues up to the critical value ε= 0.56, where one chaotic attractor appears (see Fig. 12e). For ε= 0.57 a two-scrolls chaotic attractor appears (see Fig. 12f). Figure 12: Phase portrait of Chua system for different values εand c= 11(attractors). (a) a period-one limit cycle for ε= 0.43, (b) a period-two limit cycle for ε= 0.46, (c) a period-4 limit cycle for ε= 0.495, (d) a period-8 limit cycle for ε= 0.501, (e) a chaotic attractor chaotic attractor for ε= 0.56, (f) a 2-spirals for ε= 0.57 . 79