Approximating hidden chaotic attractors via parameter switching

In this paper, the problem of approximating hidden chaotic attractors of a general class of nonlinear systems is investigated. The Parameter Switching (PS) algorithm is utilized, which switches the control parameter within a given set of values with the initial value problem numerically solved. The PS-generated generated attractor approximates the attractor obtained by averaging the control parameter with the switched values, which represents the hidden chaotic attractor. The hidden chaotic attractors of a generalized Lorenz system and the Rabinovich-Fabrikant system are simulated for illustration

In [10,11,12] it is proved that the attractors of a chaotic system, considered as the unique numerical approximations of the underlying ω-limit sets (see e.g.[31]), after neglecting sufficiently long transients, can be numerically approximated by switching the control parameter in some deterministic or random manner while the underlying initial value problem (IVP) is numerically integrated with the Parameter Switching (PS) algorithm.The attractors, whose basins of attractions are not connected with equilibria are called hidden attractors, while the attractors for which the trajectories starting from a point in a neighborhood of an unstable equilibrium are attracted by some attractor, are called self-excited attractors [4,5,6].In this paper we prove analytically and verified numerically that the PS algorithm can be used to approximate any desired hidden attractors of a class of general systems which model systems like Lorenz, Chen, Rössler, etc.

Introduction
One main task in the investigation of a dynamical model is to study the limiting behavior of the system states after the transient processes, i.e., the problem of localization and analysis of attractors (limiting sets of system's states).Here, one of the challenging problems is to study models with multistability, whose states can alternate between some mutually exclusive attractors over time [1].In such models, particularly in the case of the existence of attractors with very small basins or unidentified attractors, one can observe sudden switching to unexpected (unpredictable or unknown) attractors, since such systems are sensitive to noise, initial conditions, and system parameters [2,3].While trivial attractors (stable equilibria) can be found either analytically or numerically of any dynamical system, the search for nontrivial attractors could be a very challenging task (e.g. the famous Hilbert 16th problem on periodic attractors in two-dimensional polynomial dynamical systems is still far from being solved).The structures of many classical physical dynamical models guarantee that attractors exist because the trajectories cannot tend to infinity and the oscillations are excited by an unstable equilibrium.Such attractors are called self-excited attractors, which can be easily found by constructing a solution using initial data from a small neighborhood of the equilibrium, observing how it is attracted thereby visualizing the attractor.However, there are attractors of another type, called hidden attractors [4,5,6], whose basins of attractions are not connected with equilibria and, thus, the search and study of such attractors are very challenging [3,7].For example, hidden attractors can be in systems with no equilibria [8] or in a multistable system with only one stable equilibrium [9].
Self-excited attractors can be numerically visualized through a standard computational procedure, in which after the transient process a trajectory starting from a point in a neighborhood of an unstable equilibrium is attracted to an attractor.In contrast, the basin of attraction for a hidden attractor is not connected with any small neighborhood of any equilibrium and, thus for the numerical localization of a hidden attractor it is necessary to develop a special analytical-numerical procedure, in which an initial point is chosen from the basin of attraction.To numerically verify that a chaotic attractor is hidden, one has to check that all trajectories starting in small neighborhoods of unstable equilibria are either attracted by stable attractors or diverging to infinity.
The known autonomous chaotic dynamical systems depending on a single real control parameter p ∈ R, such as the Lorenz system, Rössler system, Chen system, Lotka-Volterra system, Rabinovich-Fabrikant system, Hindmarsh-Rose system, Lü system, classes of minimal networks and many others, are modeled by the following Initial Value Problem (IVP): where t ∈ I = [0, T ], x 0 ∈ R n , p ∈ R the control parameter, A ∈ R n×n a constant matrix, and f : R n → R n a continuous nonlinear function.
For example, for the Lorenz system with n = 3 and the standard parameter values a = 10 and c = 8/3, if one considers p = ρ then system (1) has The PS algorithm approximates numerically any solution of the IVP (1) [10,11,12].If one chooses a finite set of values of the underlying control parameter, P N = {p 1 , p 2 , ..., p N }, N ≥ 2, and then switches p within P N for a relatively short period of time, while the underlying IVP is numerically integrated, then the resultant "switched" numerical solution will converge to the "averaged" solution of the system.Consequently, any attractor of the underlying system, obtained for p being replaced with the average of switched values, can be approximated by the attractor generated from the switching operations.
The PS algorithm is useful e.g. when one intends to obtain an attractor but for some reason the underlying parameter of the attractor cannot be set.Also, the PS algorithm could explain why, in some natural systems, alternations between different dynamics could lead to unexpected behavior.
In this paper, the PS algorithm is used to approximate some hidden chaotic attractors which, as mentioned above, is a challenging task.
The paper is organized as follows: Section 2 presents the PS algorithm, its convergence and its numerical implementation.In Section 3, the PS algorithm is used to approximate hidden attractors in a generalized Lorenz system and the Rabinovich-Fabrikant system, respectively.The short conclusion section ends the investigation.

Description and convergence of the PS algorithm
Let P N = {p 1 , p 2 , ..., p N } ⊂ R, a set of N values of parameter-p, N ≥ 2. Consider the IVP (1), numerically integrated it over I with p switching periodically its values within P N for relatively short periods of time.The PS algorithm is associated with the "switching" equation in the following form: with p h : I → P N being a T p -periodic piece-wise constant function, depending on a small h > 0, which switches periodically its values p h (t) = p h (t + T p ) = p i ∈ P N , i ∈ {1, 2, ..., N }, for t ∈ I i,j , and j = 1, 2, ..., where I i,j are subintervals of the time interval I = j N i=1 I i,j (see the sketch in Fig. 1 for the case of N = 3 and P 3 = {p 1 , p 2 , p 3 }).Denote the average of the switched values by p * , which is a constant having the same value for all t ∈ I = [0, T ]: where p h (•) is one of the parameter values p i , i ∈ {1, 2, ..., N }.Then, the "averaged" equation of ( 1), obtained for p being replaced with p * , reads It can be proved that switching p within P N in (3) while the switching equation ( 3) is integrated, the obtained solution approximates the solution of the averaged equation ( 5).Before proceeding, the following assumption is needed.Assumption H1.In the IVP (1), f satisfies the Lipschitz condition for some L > 0. Denote p h (t) := P (t/h) and let • 0 be the maximum norm on C([0, T ], R n ), i.e., x 0 := max t∈[0,T ] |x(t)|.Then, under Assumption H1, on [0, T ], the following theorem holds1 . where Proof.From ( 3) and ( 5), one has On the other hand, one has Hence, By the Gronwall inequality [28], one obtains (7).
i) The above proof is more general than the proof presented in [12], where the convergence is obtained via the averaging method [29] and the initial conditions of (3) and (5) are equal.In [11], the proof is made numerically on the basis of the global error of Runge-Kutta.In [30], beside the convergence of the PS algorithm, numerical approximation estimation and Lyapunov method are presented, and moreover the PS convergence for any utilized Runge-Kutta method is proved.
ii) The periodicity assumption on p in (3) is too strong.Actually, the convergence proof in [11] and the numerical experiments in [21,13,14] (or experimental applications in [26]) show that the PS algorithm can be implemented in some random way as well.
For instance, once P N is set, the order in which p switches its values, p = p i ∈ P N , can be random.Therefore, one may assume that random or periodic switches of parameters in natural systems have a real meaning, such as in ecological systems or circuitry.Also, random parameter switches in some systems explain why chaotic (hidden) attractors could appear unexpectedly.
iii) If the averaged system has a hyperbolic invariant compact set, then the switching equation ( 3) has also a near hyperbolic invariant compact set.
A global attractor is a compact and invariant set composed of all bounded global trajectories and contains all the dynamics evolving from all possible initial conditions.In other words, it contains all solutions, including stationary solutions, periodic solutions, as well as chaotic attractors, relevant to the asymptotic behaviors of the system.On the contrary, a local attractor is a compact invariant set, which attracts its neighboring trajectories.A global attractor is hence composed of the set of all local attractors, where each local attractor only attracts trajectories from a subset of initial conditions, specified by its basin of attraction.In some cases, a unique local attractor may also be the global one.When a global attractor is composed of several local attractors, the initial conditions are essential for the numerical approximations of these attractors, respectively.Therefore, the following assumption is made.Assumption H2 Suppose that x 0 and x0 belongs to the same attraction basin of solutions to the IVP (3).
The ω-limit set of a trajectory through x ∈ R n is given as ω(x) = s≥0 t≥s Φ(t, x), where Φ(t, x) is the flow of the system.
As common in numerical investigations of nonlinear systems, for every p and x 0 , by an attractor one considers the unique numerical approximation of the underlying ω-limit sets (see e.g.[31]), neglecting sufficiently long transients.
By Theorem 1, which characterizes the PS algorithm, the following result can be obtained.
Corollary 2. Every attractor of the system modeled by the IVP (3) can be numerically approximated using the PS algorithm.
In other words, using the PS algorithm, the attractor A * (switched attractor ) obtained from equation (3) by switching p within P N , will approximate numerically the attractor denoted A p * (averaged attractor ) obtained from (5).

Implementation of the PS algorithm
To implement numerically the PS algorithm, let h be the step-size of the utilized explicit numerical method for integrating the corresponding IVP (such as the standard Runge-Kutta method, used in this paper).
Symbolically, for a given h > 0, the PS algorithm can be denoted as where m i ∈ N * , i = 1, 2, ..., N , are some positive integers, called "weights" of the p values.Then, p * can be expressed as The scheme (8) reads as follows: while the IVP (3) is numerically integrated, for the first m 1 integration steps, p = p 1 ; for the next m 2 steps, p = p 2 ; and so on, till the last m N step, when p = p N .So, the first set of subintervals I i,1 , for i = 1, ..., N , is covered.
Next, the algorithm repeats on the next set of subintervals I i,2 , and so on, until the entire time interval is covered.Time subintervals I i,j have lengths m i h, for i = 1, 2, ..., N and j = 1, 2, ..., and the switching period is Remark 2. The main characteristic of the PS algorithm relies on the linear dependence on p of the right-hand side of system (1) and on the convexity of the relation (9).By denoting α j = m j / N i=1 m i , j = 1, 2, ..., N , relation (9) becomes p * = N i=1 α i p i with N i=1 α i = 1.Thus, for any set P N , N > 1, and any weights m i , i = 1, 2, ..., N , p * is always inside the interval (p min , p max ), with p min ≡ min{P N } and p max ≡ max{P N }.Therefore, to approximate some attractor A p * using the PS algorithm, the set P N has to be chosen such as Consider a dynamical system modeled by the IVP (3), with the set of its attractors A and some set of admissible parameter values P N , N ≥ 2. Based on (10) and on the convexity of the relation (9), beside the fact that every attractor can be approximated with the PS algorithm (Corollary 2), the following important result can be proved.
Corollary 3. Given a set P N , with N ≥ 2 and weights m i , i = 1, 2, • • • , N , the attractor A * obtained with the PS algorithm belongs to A.
To visualize the results, i.e., to underline the match between the averaged attractor, A p * , which is to be approximated and the approximating attractor, A * , a computer-graphic criterion is now introduced.Criterion 4. Two attractors are considered to be almost-identical if a. their geometrical forms in the phase space (almost) coincide; b. the orientation of the motion is preserved.
The above criterion is a suitable modification and adaptation of the known concept of topological equivalence (see e.g.[41]), for practical use rather than for theoretical rigor.
Criterion 4.b is easy to implement computationally (e.g. with Matlab comet3 function) and has been verified for all examples studied later in this paper.
To apply Criterion 4.a (the match between the two attractors), A p * (blue or green plot) and A * (red plot) are overplotted in the phase space and also for their Poincaré sections.Visually, the histograms reveal the match between attractors.
Also, the match between the two attractors can be verified by calculating the Hausdorff distance D H (A * , A p * ) between them.The Hausdorff distance between two sets A and B in the metric space R 3 , D H (A, B), is given by [32, p.
for i, j = 1, 2, ..., M .Therefore, the Hausdorff distance ( 11) can be calculated numerically by A study of numerical limitations of the PS algorithm is presented in [20] and [11].

Hidden chaotic attractors approximated with the PS algorithm
In this section, under Assumptions H1 and H2, hidden chaotic attractors of a generalized Lorenz system and the Rabinovitch-Fabrikant system are computed with the PS algorithm.The approximated hidden chaotic attractor is the averaged attractor A p * .To approximate A p * , one chooses a set of N parameter values for the switching process, P N , such that condition (10) holds and, with weights m i , i = 1, 2, ..., N , relation ( 9) is verified.Then, scheme ( 8) is applied to obtain the switched attractor A * which approximates the desired hidden attractor A p * .
The numerical and simulations results of the PS algorithm were realized with the standard RK algorithm, which allows to implement easily the switches imposed by the algorithm for every m i steps, i = 1, 2, ..., N .The integration step-size was taken as h = 0.0002−0.001, the histograms for the x 1 component use 512 bars, and the integration time interval is For the case of stable cycles, D H (A * , A p * ) is of order 10 −3 − 10 −2 .In the case of chaotic attractors, D H (A * , A p * ) is larger, e.g., for I = [0, 300], D H (A * , A p * ) is of order 10 −1 , while for I = [0, 500], D H (A * , A p * ) diminishes at 10 −2 .
In the case of chaotic attractors, in order to reduce numerical errors, Assumption H2 is strengthened by using identical initial conditions x 0 and x0 .
Since the eigenvalues of X * 0 are (−7.355,−1, 2.855), the equilibrium X * 0 is a saddle.X * 1,2 are stable (focus node) equilibria since their eigenvalues are (−5.497,−0.002 + 3.900i, −0.002 − 3.900i).The zoomed vicinity V X * 0 of the unstable equilibrium X * 0 reveals the fact that all trajectories started from V X * 0 are attracted either by the stable equilibrium X * 1 (red trajectories) or by the stable equilibrium X * 2 (blue trajectories).Therefore, the chaotic attractor is a hidden attractor.
Example 1.A stable cycle corresponding to p * = 25.5 is obtained, situated in a relative large periodic window (Fig. 2).The attractor, considered as the averaged attractor A p * , can be obtained using e.g.scheme (8), with e.g.
Example 2. Other sets of P N with appropriate weights can be used to obtain the same stable cycle, using for example the scheme where m 1 = m 2 = 1, which gives p * = (p 1 + p 2 )/2 = 25.5.In this case, A p * is obtained by alternating every integration step, p, within P 2 = {21, 30}.
Example 3.Not only stable cycles can be approximated by the PS algorithm, but also chaotic (self-excited) attractors, for example the one corresponding to p = 34.2(see Fig. 2), can be approximated by using e.g.
[2p 1 , 3p 2 ], with P 2 = {25.5,40}, which gives p * = 34.2(Figs. 5 (a)-(d); where the Poincaré sections are obtained for x 3 = 38).Because of the infinite integration time needed to generate a chaotic attractor, between the two attractors there are some small differences as can be seen from the histograms (Figs.

(c), (d)).
Example 4. To generate the hidden attractor H 1 , one needs to find a scheme (8) which gives p * = 7.One of the possible choices is The match between the switched attractor A * and the averaged attractor A p * is underlined by the phase overplot of the averaged (hidden) attractor (green plot) and switched attractor (red plot) in Fig. 6
The equilibrium X * 0 is unstable (saddle) since its eigenvalues are (−0.5752,0.1 − i, 0.1 + i).Equilibria X 3,4 are also unstable (saddle) with eigenvalues (0.1869, −0.281+5.397i,−0.281− 5.397i), while equilibria X * 1,2 are stable (focus nodes), with eigenvalues (−0.2561, −0.060 − 1.473i, −0.060 − 1.473i).As can bee seen, trajectories from a small vicinity of the unstable equilibrium X * 0 or X 3,4 are attracted either by infinity or by the stable equilibria X * 1,2 .Compared with the hidden chaotic attractor H 1 of the Lorenz system (12), due to the presence of complex eigenvalues, in the case of the hidden chaotic attractor H 2 here, trajectories starting from X * 0 (grey and black) and also from X * 3,4 (light brown and blue, for the case of X * 4 ) exit the vicinities by spiralling routes (see detail in Fig. 3 and Fig. 7  (b) and (c)).
To easily choose the set P N , the bifurcation diagram for p ∈ [0.24, 0.295] may be utilized (Fig. 8).
Example 5. To generate H 2 , one can use e.g. the scheme [1p 1 , 2p 2 , 2p 3 ], with P 3 = {0.28,0.289, 0.29}, for which, by (9), p * = 0.2876.This generates the hidden attractor H 2 .The match between the obtained switched attractor A * (red plot) and the approximated hidden attractor H 2 (green plot) can be observed from the match in the phase plot (Fig. 9  Example 6.Similarly, the other hidden attractor H 3 (Fig. 8)3 , which corresponds to p = 0.2715, can be be obtained with the PS algorithm, by alternating e.g. the values of the set P 2 = {0.265,0.278}, using the scheme In this case, the Poincaré section is set at x 3 = 0.3.Again, there exists a perfect match between the obtained attractor A * and the hidden attractor H 2 (see the phase overplots, overploted Poincaré sections and histograms in Figs. 10 (a)-(d), respectively).

Conclusion
In this paper, it has been proved and verified numerically that the hidden chaotic attractors of dynamical systems modeled by a general initial value problem can be approximated by switching the control parameter, while the problem is integrated.The approximations is verified with numerical tools by means of phase portraits, histograms, Poincaré sections and Hausdorff distance.In order to facilitate the choice of the switching parameters, the bifurcation diagrams are also utilized.The algorithm has been applied successfully to a generalized Lorenz system and the Rabinovich-Fabrikant system.(12).Trajectories starting from the vicinity V X * 0 of the unstable equilibrium X * 0 are attracted either to the stable equilibrium X * 1 (red plot) or to the stable equilibrium X * 2 (blue plot).
a, b) is the Euclidean distance between two points a = (x 1 , x 2 , x 3 ) and b = (y 1 , y 2 , y 3 ) in R 3 .Since the two numerically generated attractors A * and A p * are curves with the same number of M ordered pairs of coordinates A * = {a 1 , a 2 , ..., a M } and A p * = {b 1 , b 2 , ..., b M }, the distance between a point a i ∈ A * to the set A p * is given by (a), overplots of Poincaré sections with x 3 = 6 (Fig. 6 (b)), and the histograms (Figs. 6 (c), (d)).

Figure 3 :
Figure 3: Hidden chaotic attractor H1 (green) of the generalized Lorenz system(12).Trajectories starting from the vicinity V X * 0 of the unstable equilibrium X * 0 are attracted either to the stable equilibrium X * 1 (red plot) or to the stable equilibrium X * 2 (blue plot).

Figure 7 :
Figure 7: Hidden chaotic attractor H2 (green plot) of the RF system (13) for p * = 0.2876.(a) Phase portrait.(b) Zoomed vicinity V X * 4 of the unstable equilibrium X * 4 .(c) Zoomed vicinity V X * 0 of the unstable equilibrium X * 0 .Trajectories starting from the unstable points X * 0 and X * 3,4 are either attracted to the stable equilibria X1,2 (dotted blue and red respectively) or to infinity (light and dark brown, and grey and black).