Malika Belouerghi
email: belouerghimalika@yahoo.fr

Tidjani Menacer
email: tidjanimenacer@yahoo.fr

René Lozi
email: rene.lozi@univ-cotedazur.fr

Integration Duration-Based Method for Unveiling Hidden Patterns of Even Number of Spirals of Chua Chaotic Attractor

Keywords: hidden bifurcation, Chua's circuit, pattern of even number of spirals

The purpose of this article is to introduce a novel method for unveiling hidden patterns of even number of spirals in the multi-spiral Chua Chaotic attractor. This method is based on the duration of integration of this system and allows displaying every pattern containing from 1 to c+1 spirals of this attractor. After having given the equation of the multispiral Chua Chaotic attractor, the Menacer-Lozi-Chua (MLC) method for uncovering hidden bifurcations is explained again and a numerical example of a route of bifurcation is given. In the chosen example, the attractors found along this hidden bifurcation route display an odd number of spirals. With the novel method, during the integration process, before reaching the asymptotical attractor which possesses an odd the number of spirals, the number of spirals increases one by one until it reaches the maximum number corresponding to the value fixed by ε, unveiling both patterns of even and odd number of spirals.

INTRODUCTION

The name 'bifurcation' was first introduced by Henri Poincaré in 1885 [START_REF] Poincaré | L'Équilibre d'une masse fluide animée d'un mouvement de rotation[END_REF] in the frame of the qualitative theory of dynamical systems. A bifurcation occurs when, for a fixed value of a parameter, the number of solutions of a family of differential equations is increased from one to two, like the bifurcation of a branch of a tree [START_REF] Sattinger | Topics in stability and bifurcation theory[END_REF][START_REF] Lozi | A computing method for bifurcation boughs of nonlinear eigen value problems[END_REF], or when the topological nature of the solution changes from steady state to periodic function (Hopf bifurcation [START_REF] Marsden | The Hopf bifurcation and its applications[END_REF] ).

In 1963, the discovery of the first chaotic strange attractor by Lorenz [START_REF] Lorenz | Deterministic non-periodic flow[END_REF] has broadened the field of the bifurcation theory. The first system of differential equation modelling real electronic device for which the asymptotic attractor is chaotic (hence a strange attractor) was invented by Chua [START_REF] Chua | The genesis of Chua's circuit[END_REF] during a stay in the laboratory of professor Matsumoto, at Waseda University, Tokyo.

Chua attractor is nowadays incredibly used because both of its realisations: electronic circuit or system of differential equations can be combined for multiple purposes [START_REF] Duan | Multiinput and multi-output nonlinear systems interconnected Chua's circuits[END_REF] . Chua circuit is the simplest electronic circuit exhibiting chaos and possesses a very rich dynamical behaviour which was verified from numerous laboratory experiments [START_REF] Zhong | Experimental confirmation of chaos from Chua's circuit[END_REF] , computer simulations [START_REF] Matsumoto | A chaotic attractor from Chua's Circuit[END_REF] and rigorous mathematical analysis [START_REF] Chua | The double scroll family[END_REF][START_REF] Lozi | The theory of confinors in Chua's circuit: accurate analysis of bifurcations and attractors[END_REF][START_REF] Boughaba | Fitting trapping regions for Chua's attractor-a novel method based on isochronic lines[END_REF] .

Since that discovery, the Chua circuit has been generalised to several versions. One generalisation substitutes the continuous piecewise linear function by a smooth function, such as a cubic polynomial [START_REF] Khibnik | On periodic orbits and homoclinic bifurcations in Chua's circuit with a smooth nonlinearity[END_REF][START_REF] Shilnikov | Chua's circuit: rigorous results and future problems[END_REF][START_REF] Huang | Chua's equation with cubic nonlinearity[END_REF][START_REF] Hirsch | Differential equations, dynamical systems and an introduction to chaos[END_REF] , etc. Complex attractors with n-double spirals were reported in Suykens and Vandewalle [START_REF] Suykens | Quasilinear approach to nonlinear systems and the design of ndouble scroll (n = 1[END_REF] by introducing additional break points in the non-linear element in the Chua's circuit or using cellular neural networks with a piecewise linear output function [START_REF] Arena | Generation of n-double scrolls via cellular neural networks[END_REF] .

It is also demonstrated in Suykens et al. and Yalçin et al. [START_REF] Suykens | A family of n-scroll attractors from a generalized Chua's circuit[END_REF][START_REF] Yalçin | Experimental confirmation of 3-and 5-scroll theory attractors from a generalized Chua's circuit[END_REF] ; odd number of scrolls can be observed by a similar modification. In Tang et al. [START_REF] Tang | Generation of n-scroll attractors via sine function[END_REF], it is shown that the n-spirals attractors can be obtained with a simple sine or a cosine function.

In 2011, Leonov et al. introduced a new classification of attractors of dynamical systems with the concept of hidden attractor and self-excited attractor [START_REF] Leonov | Localization of hidden Chua's attractors[END_REF][START_REF] Leonov | Analytical-numerical methods for investigation of hidden oscillations in nonlinear control systems[END_REF] . They showed that the basin of attraction of a self-excited attractor overlaps with the neighbourhood of an equilibrium point; therefore, these attractors are very easy to be found, because it suffices to give the equilibrium point as an initial condition, if this equilibrium point is unstable. On the contrary, a hidden attractor has a basin of attraction that does not intersect with small neighbourhoods of any equilibrium points, thereby making it very difficult to find. The name hidden comes from this property.

Hidden attractors are important in engineering applications because they allow unexpected and potentially disastrous responses to perturbations in a structure like aircrafts control systems (windup and anti-windup) and electrical machines. An effective method for the numerical localisation of hidden attractors in multi-dimensional dynamical systems has been proposed by these Leonov and Leonov et al. [START_REF] Leonov | On harmonic linearization method[END_REF][START_REF] Leonov | Algorithm for localizing Chua attractors based on the harmonic linearization method[END_REF] . This method is based on homotopy and numerical continuation. They construct a sequence of similar systems such that for the first (starting) system the initial condition for numerical computation of oscillating solution (starting oscillation) can be obtained analytically. Then the transformation of this starting solution is tracked numerically in passing from one system to another. The first example of a hidden chaotic strange attractor was found in the Chua attractor [START_REF] Leonov | Hidden attractor in smooth Chua systems[END_REF] .

In their study of the modified multi-spiral Chua system with sine function [START_REF] Tang | Generation of n-scroll attractors via sine function[END_REF] , Menacer et al. [START_REF] Menacer | Hidden bifurcations in the multispiral Chua attractor[END_REF] , using the Kuznetsov and Leonov method in another way, found a sequence of hidden bifurcations.

In this system, the parameter c governing the number of spirals is an integer; hence, it is not possible to vary it continuously and therefore it is not possible to observe bifurcations of attractors from n to n+2 spirals when the parameter c changes. Moreover, it is not possible to use non-integer real values for c. To overcome this obstacle, they introduced a new bifurcation parameter ε provided by the method of researching hidden attractor.

Next, in 'Multiple-spiral attractors in Chua's sine function system' section, the multispiral Chua system with sine function is recalled. In 'The Menacer-Lozi-Chua (MLC) method for method for uncovering hidden bifurcations' section, the Menacer-Lozi-Chua (MLC) method for uncovering hidden bifurcations is explained again. In the chosen example, it shown that the attractors along the hidden bifurcation route display an odd number of spirals. In 'Integration duration method for unveiling hidden patterns of even number of spirals' section, we introduce a novel method for unveiling hidden patterns of even number of spirals. This method based on the duration of integration of the considered system. During the integration process, before reaching the asymptotical attractor which possesses and odd the number of spirals, the number of spirals increases until it reaches the maximum number corresponding to the value fixed by ε. Finally, in 'Conclusion' section, a brief conclusion is drawn.

MULTIPLE-SPIRAL ATTRACTORS IN CHUA'S SINE FUNCTION SYSTEM

Complex chaotic attractors with n-double scrolls using cellular neural networks with a piecewise linear output function were reported in Arena et al. [START_REF] Arena | Generation of n-double scrolls via cellular neural networks[END_REF] . In Tang et al. [START_REF] Tang | Generation of n-scroll attractors via sine function[END_REF] , a new family of continuous functions for generating n-scroll attractors is proposed. It is shown that n-scroll attractors can be obtained with a simple sine or cosine function. Since then in the last decades, multi-scroll chaotic attractors generation has been extensively studied due to their promising applications in various real-word chaos-based technologies including secure and digital communications, random bit generation and so on.

The system of differential equations, describing the behaviour of Chua's circuits, is threedimensional with a piecewise linear non-linearity [START_REF] Chua | The genesis of Chua's circuit[END_REF][START_REF] Madan | Chua's circuit: a paradigm for chaos[END_REF] . In this article, rather than considering the fractal structure of the attractors, we highlight their global geometric feature, that is, we only consider their general shape described in terms of the number of spirals. We fix the real parameter values 𝛼 = 11, 𝛽 = 15, 𝑎 = 2, 𝑏 = 0.2 for which one can find chaotic attractors topologically equivalent to the attractors described in Ref. [21].

{ 𝑥(𝑡) = 𝛼(𝑦(𝑡) -𝑓(𝑥(𝑡))) 𝑦̇(𝑡) = 𝑥(𝑡) -𝑦(𝑡) -𝑧(𝑡) 𝑧(𝑡) = -𝛽𝑦(𝑡) (1)
The parameter c is an integer which governs the number n of spirals. In order to compare the new integration duration based method versus the method previously introduced by Menacer et al. [START_REF] Menacer | Hidden bifurcations in the multispiral Chua attractor[END_REF] we need a benchmark sequence of hidden bifurcations. In this goal, we limit our study to the case where parameter c is set to the value 12 which gives a 13 spirals attractor (Figure 2). This number of scrolls is enough to obtain a significant comparison. For a complete study of this model, with other values of c see Ref. [27].

𝑛 = 𝑐 + 1 (

THE MLC METHOD FOR METHOD FOR UNCOVERING HIDDEN BIFURCATIONS

In the modified Chua circuit with sine function ( 1)-( 2), the parameter c governing the number of spirals is an integer; hence, it is not possible to vary it continuously and therefore it is not possible to observe bifurcations of attractors from n to n + 2 spirals when the parameter c changes. Moreover, it is not possible to use non-integer real values for c. To overcome this obstacle, [27] introduced a new method for uncovering hidden bifurcations based on the core idea of the Leonov and Kuznetsov for searching hidden attractors (i.e. homotopy and numerical continuation, see Appendix). While keeping c constant, a new bifurcation parameter ε is introduced. We recall briefly this method in this section in which the value of parameters is fixed at 𝛼 = 11, 𝛽 = 15, 𝑎 = 2, 𝑏 = 0.2, 𝑐 = 12, 𝑑 = 𝜋.

The Method

For this purpose, we rewrite System (1)-( 2) in the form Lur'e's form [29] 𝑑𝑋 𝑑𝑡 = 𝑀𝑋 + 𝑞𝐻(𝑟 𝑇 𝑋) 𝑋 ∈ 𝑅 3 , (

where 𝑀 = (

0 𝛼 0 1 -1 1 0 -𝛽 0 ) , 𝑋 = ( 𝑥 𝑦 𝑧 ) , 𝑞 = ( -𝛼 0 0
) , 𝑟 = ( 1 0 0

), and 𝐻(𝜎) = 𝑓(𝜎).

We introduce a coefficient  and a small parameter  in order to transform System (5) into the form similar to the system (V) as follows 03 0 ( ),

T dX M X q g r X X R dt  = +  (6)
where

𝑀 0 = 𝑀 + 𝜒𝑞𝑟 𝑇 = ( -𝛼𝜒 𝛼 0 1 -1 1 0 -𝛽 0 ),
and 𝑔(𝜎) = 𝐻(𝜎) -𝜒𝜎 = -𝛼𝑓(𝜎) -𝜒𝜎.

Then considering the transfer function 𝑊(𝑚) (XI) solving the equations 𝐼𝑚 𝑊 (𝑖𝜔 0 ) = 0 and 𝜒 = -𝑅𝑒 𝑊 (𝑖𝜔 0 ) -1 , which obtains 𝜔 0 = 2.1018 and 𝜒 = 0.03796.

By the non-singular linear transformation, 𝑋 = 𝑆𝑍 the System ( 6) is reduced to the form similar to Equation (VI) 𝑑𝑍 𝑑𝑡 = 𝐴𝑍 + 𝐵𝜀𝑔 0 (𝐶 𝑇 𝑍), 𝑍 ∈ 𝑅 3 (7) where

𝐴 = ( 0 -𝜔 0 0 𝜔 0 0 0 0 0 -𝑑 1 ) , 𝑍 = ( 𝑧 1 𝑧 2 𝑧 3 ) , 𝐵 = ( 𝑏 1 𝑏 2 1
) , 𝐶 = (

1 0 -ℎ ) . (8)
Remark: The System (1) being 3-dimensional, 𝑍 3 of Equation (VI) is a scalar denoted here 𝑧 3 This implies that ).

𝐴 = 𝑆 -1 𝑀 0 𝑆, 𝐵 = 𝑆 -1
The transfer function 𝑊 𝐴 (𝑚) of System (7) reads

𝑊 𝐴 (𝑚) = -𝑏 1 𝑚 + 𝑏 2 𝜔 0 𝑚 2 + 𝜔 0 2 + ℎ 𝑚 + 𝑑 1 .
Then, using the equality of transfer of both functions 𝑊 𝐴 (𝑚) and 𝑊 𝑀 0 (𝑚) = 𝑟 𝑇 (𝑀 0 -𝑚𝐼) -1 𝑞 of Systems ( 6) and ( 7), we obtain the following relations

𝜒 = 𝑎+𝜔 0 2 -𝛽 𝛼 , 𝑑 1 = 𝛼 + 𝜔 0 2 -𝛽 + 1, ℎ = 𝛼(𝛽 -𝑑 1 + 𝑑 1 2 ) 𝜔 0 2 + 𝑑 1 2 , 𝑏 1 = 𝛼(𝛽 -𝜔 0 2 -𝑑 1 ) 𝜔 0 2 + 𝑑 1 2 , 𝑏 2 = 𝛼(1 -𝑑 1 )(𝜔 0 -𝛽𝑑 1 ) 𝜔 0 2 (𝜔 0 2 + 𝑑 1 2
) .

Numerical Example of Hidden Bifurcation Route

For the value of parameters fixed in this section, one obtains 𝜒 = 0.03796, 𝑑 ).

Using Theorem 1 in the Appendix, for 𝜀 small enough, one obtains the following initial condition for the first step of the multi-stage localisation procedure:

𝑋 0 (0) = 𝑆𝑍(0) = 𝑆 ( 𝜏 0 0 0 ) = ( 𝜏 0 𝑠 11 𝜏 0 𝑠 21 𝜏 0 𝑠 31 ), (10)
Let us apply the localisation procedure described in the Appendix. First, we compute the starting frequency 𝜔 0 and the coefficient of harmonic linearisation 𝜒 (already computed in 'The method' section). Using Equation (10), we obtain the initial conditions for the first step. Then, we compute the solutions of System ( 6) with the non-linearity 𝜀𝑔(𝑥) = 𝜀(𝐻(𝑥) -𝜒𝑥), by increasing sequentially 𝜀 from the value 𝜀 = 0.1 to 𝜀 = 1 by starting at the beginning with the Step 0.1 and decreasing progressively it, down to 0.001 between 𝜀 = 0.8 and 𝜀 = 1.

If the stable periodic solution 𝑋 1 (𝑡) (corresponding to the small value 𝜀 = 0.1 near the harmonic one is determined, all the points of this stable periodic solution are located in the domain of attraction of the stable periodic solution 𝑋 2 (𝑡) of the system corresponding to 𝜀 = 0.2). The solution 𝑋 2 (𝑡) can be found numerically by searching one trajectory of System (5) with 𝜀 = 0.2 taking as an initial point 𝑋 1 (𝑡 𝑚𝑎𝑥 ) where tmax is the last value of the integration time. For the value of the parameter 𝜀 = 0.2 after a transient regime, the computational procedure reaches the periodic solution 𝑋 2 (𝑡). We continue by increasing the parameter  using the same numerical procedure to obtain 𝑋 3 (𝑡), 𝑋 [START_REF] Marsden | The Hopf bifurcation and its applications[END_REF] (𝑡), ⋯ , 𝑋 𝑗 (𝑡), ⋯ which are solutions of Systems ( 1)-( 2) for peculiar initial conditions. When 𝜀 = 0.8, the first chaotic solution possessing one scroll is found.

We apply the above procedure to obtain the initial conditions for the solutions for increasing values of 𝜀 as shown in Table 1. Using these initial conditions, we get the solutions 𝑋 8 (𝑡) (Figure 3a) with one scroll to 𝑋 [START_REF] Shilnikov | Chua's circuit: rigorous results and future problems[END_REF] (𝑡) (Figure 5g). In each figure, there are a different odd number of spirals in the attractor. The number of spiral increases by 2 at each step as displayed in Table 2 from 1 to 13 spirals. The values of 𝜀 in this table are exactly the values of bifurcation points. In Figure 3, we display the hidden bifurcations of the multi-spiral Chua's attractor. 𝑋 [START_REF] Matsumoto | A chaotic attractor from Chua's Circuit[END_REF] (0) 𝑋 10 (0) 𝑋 [START_REF] Lozi | The theory of confinors in Chua's circuit: accurate analysis of bifurcations and attractors[END_REF] (0) 𝑋 [START_REF] Boughaba | Fitting trapping regions for Chua's attractor-a novel method based on isochronic lines[END_REF] (0) 𝑋 [START_REF] Khibnik | On periodic orbits and homoclinic bifurcations in Chua's circuit with a smooth nonlinearity[END_REF] (0) 𝑋 [START_REF] Shilnikov | Chua's circuit: rigorous results and future problems[END_REF] (0) 1 spiral 3 spirals 5 spirals 7 spirals 9 spirals 11 spirals 13 spirals

INTEGRATION DURATION METHOD FOR UNVEILING HIDDEN PATTERNS OF EVEN NUMBER OF SPIRALS

In the previous section, it is shown that the attractors along the hidden bifurcation route display an odd number of spirals. In this section, we introduce a novel method for unveiling hidden patterns of even number of spirals. This method based on the duration of integration of Systems ( 1)- (2). For this novel method, we fix ε and we increase the duration of integration tmax. During the integration process, before reaching the asymptotical attractor which possesses and odd the number of spirals, the number of spirals increases until it reaches the maximum number corresponding to the value fixed by ε.

As a numerical example, we set the value of ε to 0.9953. For this value, 11 spirals belong to the asymptotic attractor (Figure3f). All the patterns are displayed in Figure 4 and the value of tstepmax is given in Table 3. In each figure, the number of spiral increases by 1 (instead of 2 in the MLC method). The MATLAB's standard solver for ordinary differential equations ode45 is used in order to integrate the Chua's system. This function implements a Runge-Kutta method with a variable time step. Therefore, it is difficult to know tmax. The value tstepmax (shown in Table 3) which is the maximum number of steps used (i.e. corresponding to the duration of integration). In order to highlight the efficiency of this integration duration based method, we repeat the same procedure for all values of ε of Table 2, we are going to observe the change of scroll number and we resume our results in Table 4 and Figure 5.

CONCLUSION

In this article, we have introduced a novel method for unveiling hidden patterns of even number of spirals in the multi-spiral Chua Chaotic attractor. This method is based on the duration of integration of this system and allows displaying every pattern containing from 1 to c + 1 spirals of this attractor. After having given the equation of the multi-spiral Chua Chaotic attractor, the MLC method for uncovering hidden bifurcations is explained again and a numerical example of a route of bifurcation is given. In the chosen example, it is shown that the attractors found along this hidden bifurcation route display an odd number of spirals. With the novel method, during the integration process, before reaching the asymptotical attractor which possesses and odd the number of spirals, the number of spirals increases one by one until it reach the maximum number corresponding to the value fixed by ε, unveiling both patterns of even and odd number of spirals.

APPENDIX ANALYTICAL-NUMERICAL METHOD FOR HIDDEN ATTRACTOR LOCALISATION

Ten years ago, Kuznetsov and Leonov [START_REF] Leonov | Algorithm for localizing Chua attractors based on the harmonic linearization method[END_REF][START_REF] Leonov | Hidden attractor in smooth Chua systems[END_REF] proposed an effective method for the numerical localisation of hidden attractors in multi-dimensional dynamical systems.

Their method is based on homotopy and numerical continuation. They construct a sequence of similar systems such that for the first (starting) system the initial data for numerical computation of an oscillating solution (starting oscillation) can be obtained analytically. Then the transformation of this starting solution is tracked numerically in passing from one system to another. This proposed approach is generalised in the work [23] to the systems of the form

𝑑𝑋 𝑑𝑡 = 𝑀𝑋 + 𝐻̱ (𝑋) (I)
where 𝑀 is a constant 𝑛 × 𝑛 matrix, 𝐻̱ (𝑋) is a continuous vector-function and 𝐻̱ (0) = 0. They consider a system having one scalar non-linearity, which takes the following form:

𝑑𝑋 𝑑𝑡 = 𝑀𝑋 + 𝑞𝐻(𝑟 𝑇 𝑋), 𝑋 ∈ ℝ 𝑛 (II)
Here 𝑞, 𝑟 are constant n-dimensional vectors, T is the transposition operation, 𝐻(𝜎) is a continuous piecewise differentiable scalar function satisfying the condition 𝐻(0) = 0.

They then define a coefficient of harmonic linearisation 𝜒 in such a way that the matrix

𝑀 0 = 𝑀 + 𝜒𝑞𝑟 𝑇 (III)
has a pair of purely imaginary eigenvalues ±𝜔 0 , (𝜔 0 > 0) and the rest of its eigenvalues have negative real parts. They assume that such 𝜒 exists and they can rewrite the System (II) as

𝑑𝑋 𝑑𝑡 = 𝑀 0 𝑋 + 𝑞𝑔(𝑟 𝑇 𝑋) (IV)
where 𝑔(𝜏) = 𝐻(𝜏) -𝜒𝜏.

After that, they introduce a finite sequence of functions 𝑔 0 (𝜏), 𝑔 1 (𝜏), … , 𝑔 𝑚 (𝜏) such that the graph of the functions 𝑔 𝑗 (𝜏) and 𝑔 𝑗+1 (𝜏), (𝑗 = 0, … , 𝑚 -1) differs slightly from one to another, where the function 𝑔 0 (𝜏) is the small and 𝑔 𝑚 (𝜏) = 𝑔(𝜏). Using the smallness of function, they can apply the method of harmonic linearisation (describing function method) for the system

𝑑𝑋 𝑑𝑡 = 𝑀 0 𝑋 + 𝑞𝑔 0 (𝑟 𝑇 𝑋) (V)
in order to determine the non-trivial stable periodic solution 𝑋 0 (𝑡). For the localisation of the attractor of the original System (IV), they follow numerically the transformation of this periodic solution (a starting oscillating attractor, not including equilibria, denoted further by 𝐴 0 with increasing j).

To determine the initial condition 𝑋 0 (0) of the periodic solution, the System (V) can be transformed by a linear non-singular transformation S (X = SZ) to the form { 𝑧̇1(𝑡) = - 𝜔 0 𝑧 2 (𝑡) + 𝑏 1 𝑔 0 (𝑧 1 (𝑡) + 𝑐 3 𝑇 𝑍 3 (𝑡)) 𝑧̇2(𝑡) = 𝜔 0 𝑧 1 (𝑡) + 𝑏 2 𝑔 0 (𝑧 1 (𝑡) + 𝑐 3 𝑇 𝑍 3 (𝑡))

𝑍 ̇3(𝑡) = 𝐴 3 𝑍 3 (𝑡) + 𝐵 3 𝑔 0 (𝑧 1 (𝑡) + 𝑐 3 𝑇 𝑍 3 (𝑡))

(VI)
Here, 𝑧 1 (𝑡), 𝑧 2 (𝑡) are scalar values, 𝑍 3 (𝑡) is a(𝑛 -2) dimensional vector; 𝐵 3 and 𝑐 3 are (𝑛 -2) dimensional vectors, 𝑏 1 and𝑏 2 are real numbers; 𝐴 3 is (𝑛 -2) × (𝑛 -2) matrix, whose all of its eigenvalues have negative real parts. Without loss of generality, it can be assumed that for the matrix 𝐴 3 there exists a positive number 𝑑 1 > 0 such that then for the initial condition of the periodic solution 𝑋 0 (0) = 𝑆(𝑧 1 (0), 𝑧 2 (0), 𝑍 3 (0)) 𝑇 . At the initial step of the algorithm, one has 𝑧 1 (0) = 𝜏 0 + 𝑂(𝜀), 𝑧 2 (0) = 0, 𝑍 3 (0) = 𝑂 𝑛-2 (𝜀), where 𝑂 𝑛-2 (𝜀) is an (𝑛 -2) dimensional vector such that all its components are 𝑂(𝜀).

For the stability of 𝑋 0 (𝑡) (where stability is defined in the sense that for all solutions with the initial data sufficiently close to 𝑋 0 (0) the modulus of their difference with 𝑋 0 (𝑡) is uniformly bounded for all 𝑡 > 0), it is sufficient to require the following condition is true:

𝑏 1 𝑑𝐺(𝜏) 𝑑𝜏 | 𝜏=𝜏 0 < 0 (X)
In practice, to determine 𝜒 and 𝜔 0 , one uses the transfer function 𝑊(𝑚) of the System (IV) 𝑊(𝑚) = 𝑟(𝑀 -𝑚𝐼) -1 𝑞 (XI) where m is a complex variable. The number 𝜔 0 is determined from the equation 𝐼𝑚 𝑊 (𝑖𝜔 0 ) = 0 and 𝜒 is calculated then by the formula 𝜒 = -𝑅𝑒 𝑊 (𝑖𝜔 0 ) -1 .

2 )Figure 1 .

 21 Figure 1. Proposed Sine Function () fxwith Parameters Values 𝑎 = 2, 𝑏 = 0.2, 𝑐 = 12, 𝑑 = 𝜋.

Figure 2 .

 2 Figure 2. The 13-Spiral Attractor Generated by Equations (1) and (2) for 𝑐 = 12. (a) 3-Dimensional Figure, (b) Projection onto the Plane (x, y).

Figure 3 .Figure 3 .

 33 Figure 3. (Continued)

Figure 4 .Figure 4 .

 44 Figure 4. (Continued)

Theorem 1 [ 25 ] 1 𝑑𝐺

 1251 )(𝐴 3 + 𝐴 3 𝑇 )𝐴𝑍 3 (𝑡) ≤ -2𝑑 1 |(𝑍 3 (𝑡)) 2 |, ∀𝑍 3 ∈ 𝑅 𝑛-2 , (VII)In the scalar case, they define the following describing function 𝐺 of a real variable 𝜏 and assume the existence of its derivative: 𝐺(𝜏) = ∫ 𝑔(𝑐𝑜𝑠( 𝜔 0 𝑡)𝜏) 𝑐𝑜𝑠( 𝜔 0 𝑡)If it can found a positive 𝜏 0 such that 𝐺(𝜏 0 ) = 0, 𝑏

  𝑞, 𝐶 𝑇 = 𝑟 𝑇 𝑆, and 𝑆 = ( 𝑠 11 𝑠 12 𝑠 13 𝑠 21 𝑠 22 𝑠 23 𝑠 31 𝑠 32 𝑠 33

  1 = 1.4176, ℎ = 26.686, 𝑏 1 = 15.686, 𝑏 2 = 3.1003,

		1	0	-26.686
	And solving the matrix Equation (9) one gets 𝑆 = (	0.03795 -0.19107 -2.4264
		-1.3636 -0.27084 -25.674

Table 1 .

 1 Initial Condition According to the Values of 𝜀

	𝜺	𝑿 𝒋 (𝟎)	𝒙 𝟎	𝒚 𝟎	𝒛 𝟎
	0.1	𝑋 1 (0) = 𝑋 0 (𝑡 𝑚𝑎𝑥 )	74.72	2.8356	-101.89
	0.2	𝑋 2 (0) = 𝑋 1 (𝑡 𝑚𝑎𝑥 )	-3.5603	0.0978	5.1616
	0.3	𝑋 3 (0) = 𝑋 2 (𝑡 𝑚𝑎𝑥 )	3.7396	0.1166	-5.1258
	0.4	𝑋 4 (0) = 𝑋 3 (𝑡 𝑚𝑎𝑥 )	1.6312	0.7179	-1.3015
	0.5	𝑋 5 (0) = 𝑋 4 (𝑡 𝑚𝑎𝑥 )	2.9797	0.5334	-3.5896
	0.6	𝑋 6 (0) = 𝑋 5 (𝑡 𝑚𝑎𝑥 )	3.1137	-0.2868	-4.7176
	0.7	𝑋 7 (0) = 𝑋 6 (𝑡 𝑚𝑎𝑥 )	-3.2894	-0.5358	3.9445
	0.8	𝑋 8 (0) = 𝑋 7 (𝑡 𝑚𝑎𝑥 )	-3.5390	0.2601	4.9005
	0.86	𝑋 9 (0) = 𝑋 8 (𝑡 𝑚𝑎𝑥 )	-3.5920	0.1733	5.0637
	0.9785	𝑋 10 (0) = 𝑋 9 (𝑡 𝑚𝑎𝑥 )	10.5751	0.1758	-11.8656
	0.989	𝑋 11 (0) = 𝑋 10 (𝑡 𝑚𝑎𝑥 )	-12.1911	0.3828	11.2562
	0.993	𝑋 12 (0) = 𝑋 11 (𝑡 𝑚𝑎𝑥 )	-7.4755	-0.6595	5.6561
	0.9953	𝑋 13 (0) = 𝑋 12 (𝑡 𝑚𝑎𝑥 )	8.9249	0.7895	-7.6282
	0.9994	𝑋 14 (0) = 𝑋 13 (𝑡 𝑚𝑎𝑥 )	-2.8006	-0.2619	3.9268
	1	𝑋 15 (0) = 𝑋 14 (𝑡 𝑚𝑎𝑥 )	-2.6705	0.2845	2.2206

Table 2 .

 2 Values of Parameter 𝜀 at the Bifurcation Points for c = 12 (13 Spirals)

	𝜺

Table 3 .

 3 Values of tstepmax for  = 0.9953 and c= 12

	tstepmax	144	2050	2070	2700	6090	6700
	Number of spirals 1	2	3	4	5	6
	Figure	4a	4b	4c	4d	4e	4f
	tstepmax	9590	14,000	70,000	88,210	200,000
	Number of spirals 7	8	9	10	11
	Figure	4g	4h	4i	4j	4k

Table 4 .

 4 Values of tstepmax for 𝜀 = 0.80 to 𝜀 = 0.9994 and c = 12

	ε = 0.80	tstepmax Nb of spirals	50,000 1
	ε = 0.86	tstepmax Nb of spirals	7200 1	19,000 41,430 2 3
	ε = 0.97885	tstepmax Nb of spirals	195 1	680 2	20,000 90,000 100,000 3 4 5
		tstepmax	95	350	1000	1380	6000
	ε = 0.989	Nb of spirals tstepmax	1 80,000 90,110 2	3	4	5
		Nb of spirals	6	7
		tstepmax	85	2260	2300	2340	9310
	ε = 0.993	Nb of spirals tstepmax	1 9410	2 36,090 36,322 90,000 3 4	5
		Nb of spirals	6	7	8	9
		tstepmax	144	2050	2070	2700	6090
		Nb of spirals	1	2	3	4	5
	ε = 0.9953	tstepmax Nb of spirals	6700 6	9590 7	14,000 60,000 88,210 8 9 10
		tstepmax	200,000
		Nb of spirals	11
		tstepmax	91	910	1435	2250	3590
		Nb of spirals	1	2	3	4	5
	ε = 0.9994	tstepmax Nb of spirals	4450 6	5610 7	8070 8	8773 9	10,125 10
		tstepmax	15,000 20,000 70,000
		Nb of spirals	11	12	13

Figure 3. (Continued)

Figure 4. (Continued)

Values of ε