Complex dynamics, hidden attractors and continuous approximation of a fractional-order hyperchaotic PWC system

In this paper, a continuous approximation to studying a class of PWC systems of fractional-order is presented. Some known results of set-valued analysis and differential inclusions are utilized. The example of a hyperchaotic PWC system of fractional order is analyzed. It is found that without equilibria, the system has hidden attractors.


Introduction
This paper studies a new class of piecewise continuous (PWC) fractional-order (FO) systems modeled by the following general initial value problem (IVP): where the PWC function f : R n → R n has the form of f (x(t)) = g(x(t)) + A(x(t))s(x(t)), (2) in which q ∈ (0, 1), g : R n → R n a scalar vectorvalued function, at least continuous, with s : R n → R n , s(x) = (s 1 (x 1 ), s 2 (x 2 ), . . . , s n (x n )) T a vector-valued piecewise function, with s i : R → R, i = 1, 2, . . . , n, real piecewise constant functions (e.g., sgn functions), and A n×n a square matrix of real functions. Let M be the discontinuity set. Moreover, let D q * denote the Caputo differential operator of order q with starting point 0 [1]: One of the reasons to use Caputo differential operator is that it has a physically meaningful interpretation for the initial conditions just like in the integerorder problems [a unique condition x(0) in the case of q ∈ (0, 1)].
Remark 1 Fractional-order differential equations (FDEs) do not define dynamical systems in the usual sense: by denoting the solution of (1) as Φ(t, x 0 ), one does not have the flow property Φ s • Φ t = Φ t+s [2]. However, in this paper, by numerical calculation of solutions, the definition of an integer-order dynamical system is adopted, which states that if the underlying IVP admits a solution, the problem defines a dynamical system [3, Definition 2.1.2].
Because the systems modeled by the IVP (1) are autonomous, hereafter, unless otherwise indicated, the time variable will be dropped in writing.
For numerical integration of discontinuous ordinary differential equations (ODEs) of integer order, there exist dedicated numerical methods (see, e.g., the survey [4]) and, there are two main strategies to approach discontinuous systems: one strategy for treating discontinuities is simply to ignore them (time stepping methods) and to rely on a local error estimator such that the error remains acceptably small; the other strategy is to use a scalar event function h : R n → R, which defines the discontinuity Σ = {x ∈ R n |h(x) = 0}, to determine the intersection point as the new starting point for the numerical solution (event-driven methods).
The following aspects of discontinuous dynamical systems should be mentioned: A numerical method for discontinuous systems may become either inaccurate or inefficient, or both, in a region where discontinuities of the solution or its derivatives occur and the local truncation error analysis, which forms the basis of most step-size control techniques, fails if there is not sufficient (local) smoothness.
Several known numerical methods assume that trajectories will cross the discontinuity surface as they reach it, or this will never happen. But there will always be errors in finding discontinuities.
Actually, systems of PWC systems are mostly ideal, since switch-like functions like sgn are used, where the hysteresis or delay of the real switching operation is not considered, or the regularization represents a good approach in these cases.
Although there are numerical methods for FDEs (see, e.g., [5][6][7]) and also for DEs with discontinuous right-hand sides (see, e.g., [8][9][10] or the survey [4]), to the best of our knowledge, there are no numerical methods for FDEs with discontinuous right-hand sides. Consequently, modeling continuously or smoothly the underlying systems is of real important.
To approximate the PWC problem (1)-(2), one has to regularize the right-hand side first using, e.g., the Filippov regularization, transforming the discontinuous problem to a set-valued IVP, i.e., an FO ordinary differential inclusion (DI) of FO. Then, Cellina's Theorem ensures the existence of continuous approximations. 1 In this way, the obtained continuous FO problem can be numerically integrated using, e.g., the multistep predictor-corrector Adams-Bashforth-Moulton (ABM) scheme for FDEs.
Two kinds of continuous approximations are proposed and utilized in this paper: global approximation and local approximation.
As an example of the PWC system modeled by (1)-(2), consider a fractional variant of one PWC system, proposed in [11] for the integer case, 2 as follows: where a, b > 0 are real positive parameters. In this paper, let a = 1, b be the bifurcation parameter and, unless otherwise specified, q = 0.98. Comparing (3) with (1), one has 1 Graphically, by approximation one understands that not all graphic points of the PWC function to be approximated need to be located on the created figure, compared with interpolation where all graphic points of the PWC function have to be located on the created figure. 2 See also [12], where several PWC and non-smooth jerk systems are proposed. Fig. 1a. Thus, R 4 is splitted by the discontinuity surface, x 1 = 0, to two open half spaces, Ω ± = {x ∈ R 4 : ±x 1 > 0}. Besides rich dynamics (similar to those presented in [11] for the integer-order case), it will be shown that the FO setting reveals some new behavior. It will be shown numerically that the local approximation to the PWC system by a continuous one is preferable instead of the classical global sigmoid approximation. Finite-time local Lyapunov exponents are determined, and chaotic as well as hidden hyperchaotic attractors are found. Sliding motion is numerically investigated. Moreover, using existing results on the periodicity of solution of FDEs, it will be shown that such a system cannot admit stable cycles.
The paper is organized as follows. Section 2 presents the continuous approximation of PWC systems (1). In Sect. 3, the PWC system (3) is approximated. In Sect. 4, the dynamics of the obtained continuous system are numerically investigated. In Conclusion section, the obtained results are enumerated, while the three appendices present some mathematical notions and known results utilized by the paper.

Regularization
The PWC IVP (1) will be transformed to the following FO DI: The set-valued function F : R n ⇒ 2 R n can be defined in several ways. A simple (convex) expression of F, obtained by the Filippov regularization, is given by where F(x) is the convex hull of f (x), μ the Lebesgue measure and ε the radius of the ball centered at x. At those points where f is continuous, F(x) consists of one single point, which coincides with the value of f at this point (i.e., F(x) = { f (x)}). At the points belonging to M, F(x) is given by (5) ( [13, p. 85]; see also [14]).
Set-valued function F, given by (5), is upper semicontinuous (see Definition A.1 in Appendix), with compact and convex values (see also Remark A.1 in Appendix). Then, applying [2, Theorem 3.2] (see also [15]), one obtains the following result about the existence of generalized (Filippov) solutions (see Definition A.2 in Appendix) to the DI (4).
Even after regularization, the DI (4) admits generalized solutions, because single-valued IVPs offer numerical opportunities. The interest in this paper is to transform the set-valued problem (4) to a single-valued continuous one.
In order to justify the use of the Filippov regularization to some physical systems, one must choose small values for ε, so that the motion of the physical system is arbitrarily close to a certain solution of the underlying DI (it tends to the solution, as ε → 0). However, extremely small values of ε can lead to large values of derivatives and, consequently, to stiff systems. Therefore, as one can see in the following sections, special attention has to be focussed on ε.
If the piecewise constant functions s i are sgn, their set-valued forms obtained with Filippov regularization and denoted by Sgn : R ⇒ R are defined as follows: Here, the conventional sgn(0) is replaced with the whole interval [− 1, 1] connecting the points − 1 and + 1 (see Fig. 2a, b). By applying the Filippov regularization to the function f defined by (2), one arrives at the following problem: in which After regularization. c After local continuous approximation for a large ε value in the case of sgn function).

Continuous approximation
The continuous approximation of the function f defined by (2) is realized next, as described in [16]. For brevity, only the most important steps are presented. Let C 0 ε (R) be the class of real continuous approximations (see Definition A.3 in Appendix) s : R → R of the set-valued function F, which satisfy Here, B(x, ε) is the disk of radius ε centered at x. Figure  2c presents the case of the set-valued function The approximation of the set-valued function S, with the single-valued functions, will be denoted as The set-valued functions S i , i = 1, 2, . . . , n, can be approximated due to Approximate Theorem, called Cellina's Theorem (Theorem A.1 in Appendix), which states that a set-valued function F, with closed graph and convex values (Remark A.1 in Appendix), admits C 0 ε approximations. By global approximation (GA) of the PWC functions S in (6), denoteds, one understands a function defined on the entire axis R, while by local approximation (LA), denotedS ε , a function defined on some ε-neighborhood [− ε, ε] of discontinuity, with ε being a small positive number [16].
Theorem 2 [16] The PWC system (1)-(2), with g continuous, can be locally or globally continuously approximated with the following problem: wheref is either a local or a global approximation of f . If g is smooth, then the approximation is smooth.

Global approximation
Any single-valued function on R, with the graph located in the ε-neighborhood, can be considered as a global approximation of S by Cellina's Theorem (see the sketch in Fig. 2c, or Weiestrass Theorem A.2). However, some of the best candidates for s are the sigmoid functions, which provide the required flexibility and to which the abruptness of the discontinuity can be easily smoothened. If S(x) = Sgn(x), one of the mostly utilized sigmoid approximations is the following function sgn: 3 3 Sigmoid functions include the ordinary arctangent such as 2 π arctan x ε , the hyperbolic tangent used especially in modeling neural networks (see, e.g., [17]), the logistic function, some algebraic functions like x √ +x 2 , and so on. where δ is a positive parameter which controls the slope of the curve in the neighborhood of the discontinuity x = 0. In Fig. 3a, sgn is plotted for three distinct values of δ.
The smallest ε values, necessarily to embed sgn within an ε-neighborhood of Sgn (as stated by Cellina's Theorem), depend proportionally on δ. Note that for x = 0, sgn is identical to the single-valued branches of Sgn (the horizontal lines ± 1) only asymptotically, as x → ± ∞. For example, for δ = 0.01, at the point x = 0.06, the difference is of order 10 −3 , even the two graphs look apparently identical at the underlying points A or B (Fig. 3b). To reduce the ε value, e.g., to 10 −4 , δ should be of order 10 −5 . On the other hand, as one can see in Sect. 3, lower values of δ does not necessarily imply substantial reduced values of ε.
In order not to significantly affect the physical characteristics of the underlying system, it is desirable to approximate S only on some tight ε-neighborhoods of the discontinuity x = 0, not on the entire real axis, since the difference between S ands persists along the entire real axis R.

Local approximation
A better approximation is LA,s ε : [− ε, ε] → R, wheres ε is some continuous function satisfying the gluing conditions Obviously,s ε ands are both C 0 ε functions. If g is continuous, then for every ε > 0, there exists an LA of f ,f ε : R n → R n , such that [16] (see Fig. 4a) Approximations ε can also be continuously extended on R, yielding a GA,s, as follows: Among the simplest functionss ε which, compared to the exponential function in the GA (7), are the cubic Fig. 3 a Graph of sigmoid function sgn (7) for three values of δ. b Overplots of set-valued function Sgn and its approximate function sgn. Detail reveals the difference between the two curves. ε-neighborhood (see Fig. 2c) is not shown here polynomials (called cubic splines) s ε : R → R defined by Remark 2 Polynomials have the advantage to be directly evaluated by computers with the four arithmetic operations of adding, subtracting, multiplication and division.
By imposing, near the gluing conditions (8), the supplementary differentiability conditions at the boundary of the discontinuity neighborhood 4 for the particular case of the Sgn function, the smooth LA function, denoted by sgn ε , is 4 Since s i are PW constant functions on x = 0, they are differentiable on x = 0.
Using (9), Sgn can be continuously approximated on R by the following piecewise function (see Fig. 4b): Remark 3 Both GA (7) and LA (10) are also smooth approximations.
Concluding, in order to obtain numerical solutions to (4), the simplest way is to replace the discontinuous problem with the continuous (smooth) one, using one of the locally or globally approximations (see the sketch in Fig. 5).

Approximations and numerical integration of the investigated FO PWC system
Following the way the problem was transformed to DI, the PWC problem (3) can be transformed to the following set-valued problem: or Note that only the second equation is a DI. For x 4 = 0, the underlying set-valued function is present in Fig.  1b.
By applying Theorem 2, via relation (7), GA leads to the following problem: while with LA, (11), one has the following problem: Remark 4 Although both GA and LA are smooth, the approximating functionf is only continuous due to the modulus (third equation). However, the existence of solution to (12) is not affected by the non-smoothness (see, e.g., [18] for FO DIs).
The graph of approximating surface of the second component on the right-hand side of the function F, with GA for a large value of δ, δ = 1, is present in Fig.  1c. As can be seen, no significant graphical differences are indicated between both approximations.
The numerical integration of the approximate system is obtained in this paper with the predictorcorrector multi-step ABM method [5], implemented using the MATLAB subroutine FDE12.m [19] with default tolerance 1E −6, for fixed q = 0.98 and, unless otherwise specified, integration step size h = 0.002.
The following numerical analysis is performed for b = 1.25 and q = 0.98 with default double precision (roughly 15-16 decimal digits), sufficient for the current purpose.
To compare the approximation results, Fig. 6 shows the time series from component x 1 for both GA and LA (with δ = 1E − 6 and ε = 1E − 6, respectively), and also the time series obtained without approximation (WA) by using the ABM method on the nonapproximated system, which are overplotted together.
The tests have been repeated with different initial conditions and time-step sizes.
Results are summarized as follows: -All the results (for GA, LA and WA) are accurate for all t ∈ [0, t * 1 ), with t * 1 ≈ 74.95. Thus, for t ∈ [0, t * 1 ), the underlying numerical solutions coincide, because the differences between them are 0; -GA "escapes" early from this coincidence, near t * 1 (Fig. 6b, while LA and WA remain further coinciding, till t * 2 ≈ 75.60; -For decreasing values of δ and ε, 1E − 10 < δ < 1E − 6 and 1E − 8 < ε < 1E − 6, the time superior limit slightly increases for both approximations (t * max 76); -For LA, with step size 0.002 and ε < 1E − 8, the ε-neighborhood (necessary for the LA algorithm) can no longer be identified; -For δ < 1E −10 in GA, numerically sgn is the same as sgn, therefore, there exists no more approximation, for which the code works further as FDE12 applied to the original WA system; -With a smaller step size of FDE12, no significant accuracy can increase. 5 Therefore, up to t ≈ t * 1 , both approximations are time-step independent 5 The predictor-corrector ABM method has an error which is roughly proportional to h 2 . Thus, to obtain an error of, e.g., 1.0e− 6, which is sufficient in applications, a step size close to h = 1.0E − 3 must be considered.
(invariant) in the sense that the error between the two computed solutions, starting from the same initial conditions, remains zero, which are also close to WA; -GA takes longer computational time, despite its simpler form, because it is calculated over the entire axis x 1 , compared to LA which is calculated only on [− ε, ε] (see also Remark 2).
Note that, for a smaller step-size value (e.g., h = 0.0002), in the crossing and sliding surface (i.e., x 1 = 0; see Sect. 4.1), the two approximations are identical and also with the WA trajectory (Fig. 6c).
Summarizing, by comparing LA and GA and also with WA, it can be concluded that LA is better than GA in using with ABM for FDEs.
It is important to note that even FDE12.m can be applied directly to PWC problems, although the routine was not designed for this kind of problems.

Remark 5
The above results are in concordance with those obtained in [20][21][22]. Beyond numerical artifacts that might occur when numerically integrating a system of DEs, notions such as "shadowing time" and "maximally effective computational time" reveal that it is possible to have reliable numerical simulations only for some chaotic systems on a finite-time interval (see, e.g., [20][21][22]). For example, in the classical Lorenz system, obtaining a precise solution for, e.g., t ∈ [0, 100] represents a real challenge. Therefore, the case of FO systems is even more delicate.
On the other hand, larger intervals must be considered, so that possible phenomena like transient behaviors could be studied.
In this paper, after intensive numerical tests it is concluded that the obtained numerical results could be acceptable even on relatively larger intervals, as large as [0, 800].

Dynamics of the investigated FO PWC system
Next, because of the above-mentioned advantages, the LA will be utilized.

Sliding motion
In order to check that there exist crossing, sliding and grazing phenomena in system (3), it is first considered without approximation.
Denote the open half spaces by Ω ± = {x ∈ R 4 : ± x 1 > 0}. For x 0 ∈ Ω ± , the switching time t s is given by the following equation (see "Appendix B"): Since φ ± (t) are analytic functions, they may only have isolated zeros. Consequently, crossing, sliding and grazing at x(t s ) are possible, which are determined by the sign behavior of functions φ + (t) and φ − (t) near t s .
However, one cannot use the local behavior of the function f given by (2). Therefore, one does not know if the solution is really sliding on x 1 = 0, except the obtained numerical results. Moreover, one does not know how the solution is crossing the discontinuous surface. When a possible solution is crossing the discontinuous surface at some points, there is no theoretical proof for the existence of a solution. Therefore, the graphical approach is adopted to study the trajectories of the continuous approximated system, to see what happens near the discontinuity surface x 1 = 0.
As can be observed from Fig. 7a, b, d, where the case of b = 2.2 is considered, the trajectory seems to slide along the plane x 1 = 0. However, the tubular representations shown in Fig. 7 (especially the detail in Fig. 7d) would indicate that the trajectory actually crosses the plane x 1 = 0 for several (but finite number of) times. It is remarked that this phenomenon happens for other values of b too.
Note that, after approximation, the system still remains in class C 0 (see Remark 4). Therefore, beside the "smoothed" corners caused by the approximated discontinuity, which appear along the plane x 1 = 0, 6 the trajectories have some other corners due to the modulus |x 1 | (see the 3D phase projections in Fig. 7a, c, d).
An interesting behavior has been found, especially when the, respectively, initial conditions are situated in a relatively larger distance from the origin of the phase space, where the trajectories scroll toward a  (x 2 , x 3 , x 4 ) and zoomed view. Possible sliding phenomenon can be observed in a, b, d direction parallel with the axis x 3 in the 3D phase projection (x 1 , x 2 , x 3 ), or with some axis parallel to the axis x 1 in the plane x 3 = 0 in the 3D phase projection (x 1 , x 3 , x 4 ) (see Fig. 7b). This also happens, on the plane projection (x 4 , x 3 ), for the hidden hyperchaotic attractor corresponding to b = 0.5 in Fig.  4d.

Periodicity
The bifurcation diagram of the approximated system, with bifurcation parameter b (Fig. 8), shows that there are some ranges for b where the system would have stable periodic cycles (see also Fig. 9 for b = 2). In reality, the following result on the periodicity of solutions of FDEs should be noted.
Theorem 3 [24,25] (see also [26]) The fractionalorder system (1) cannot have any exact non-constant periodic solution. In [27,28], it is proved that a long-time non-constant periodic solution may have a steady-state behavior, but with − ∞ as the lower limit in the Caputo operator. In these cases, the underlying FDEs may have asymptotically T -periodic solutions for which lim t→∞ x(t + T ) − x(t) = 0 (see also [29]), for a certain T > 0.
For example, consider the following simple scalar linear FDE using Caputo derivative with the lower limit at − ∞: where α, β, γ , Ω ∈ R and [2] A periodic solution of (15) is (see the proof in "Appendix C") On the other hand, the scalar linear FDE using Caputo derivative and the lower limit at 0 (D q * ), has not periodic solutions (Theorem 3). However, continuous dynamical systems of FO, modeled by Caputo's derivative D q * , could have non- The discontinuity plane x 1 = 0 reveals the corners, typical to continuous non-smooth systems, and also the sliding phenomena along the plane x 1 constant periodic trajectories, if the system variables are impulsed periodically (impulsive fractional-order systems) [30].
Also, the corresponding torus in the bifurcation diagram for b > 2.1 in Fig. 8b could contain numerically periodic oscillations (see Fig. 10). As remarked in the previous section, tori also present scrolls before they are reached by system trajectories (Fig. 10d).

Hidden Chaotic and hyperchaotic attractors
From a computational point of view, it is natural to suggest the following classification of attractors, based on the complexity in finding basins of attraction in the phase space [32][33][34][35]: an attractor is called a self-excited attractor if its basin of attraction intersects with any open neighborhood of a stationary state (equilibrium); otherwise, it is called a hidden attractor.
For a hidden attractor, chaotic or not, its basin of attraction is not connected with equilibria. Thus, the On the other hand, hidden attractors can be attractors in systems without equilibria (the case of our system (3), see [36]), or in systems with only one stable equilibrium [37]. Hidden hyperchaotic attractors have been reported, e.g., in [38][39][40], and for a FO system, e.g., in [41].
In order to numerically find finite-time local Lyapunov exponents (LEs), one can locally approximate the non-smoothness, caused by the modulus function, by a quadratic polynomial p(x) = 1 2ε x 2 + ε 2 within the neighborhood [− ε, ε].
In general, sustained chaos is numerically indistinguishable from transient chaos, which can persist for a long time (see [17,42]). For example, for q = 0.9725 and b = 1.77, because one of the LEs {0.0058, − 0.0000, − 0.0042, − 0.0447} (measured with a precision of 1E − 5 and from the initial conditions (1, 2, 0.0, 0.1)) is 0, the system evolves first for a relatively long time (t ∈ [0, t * ], with t * ≈ 220, with hidden transient chaos, after which the trajectory is attracted by a hidden torus which, as discussed in Sect. 4.2, is characterized by a numerical periodic trajectory (see phase projections (x 1 , x 2 , x 3 ) and time series in Fig. 11a, b). On the other hand, if one uses q = 0.936 and b = 1.19, the spectrum of the LEs (with the same initial value conditions and step size) is Λ = {0.0034, 0.0022, 0.0000, − 0.0592)}, which means that the existing transient is hidden hyperchaotic, since two LEs are positive (see phase projections (x 1 , x 2 , x 3 ) and the time series in Fig. 11c, d).
Because, as indicated by the bifurcation diagram (Fig. 8), the system presents multistability, so to obtain different hidden hyperchaotic attractors, one needs to choose the initial points in their respective basins of attraction (see, e.g., the case present in Fig.  10c, where the initial conditions are (1, 2, 0, 1) and (11, − 1, 0, .1)).

Conclusion
In this paper, it has been shown numerically that local approximations of the discontinuities in systems (1)-(2) are more useful than global approximations. Therefore, the considered system is locally continuously and smoothly approximated, for which its dynamics are numerically analyzed. Without equilibria, the system admits only hidden attractors. It has been found for what values of the fractional order q and parameter b, the system admits a single positive finite-time Lyapunov exponent, with which the system behaves chaotically. It has also been found when the system admits two positive time-finite local Lyapunov exponents, with which the system behaves hyperchaotically. Finally, it has been shown that the system cannot have exact periodic oscillations; therefore, for correctness the seeming oscillations are referred to as numerically periodic oscillations. An important and yet challenging future research topic would be to implement the methodology and algorithms by physical means such as electronic circuits and mechanical structures.

A Basic notions and results
Because the set-valued property of F in (6) is generated by S i , which are real functions, the notions and results presented here are considered in R, for the case of n = 1, but they are also valid in the general cases of n > 1.
The graph of a set-valued function F is defined as follows: Remark A.1 Due to the symmetric interpretation of a set-valued function as a graph (see, e.g., [43]), a setvalued function satisfies a property if and only if its graph satisfies it. For instance, a set-valued function is closed or convex if and only if its graph is closed or convex. F is u.s.c. if it is so at every x 0 ∈ R, which means that the graph of F is closed. Generally, a set-valued function admits (infinitely) many approximations.  The last formula of (B.3) gives explicit solutions of (12) on each Ω ± , respectively. where α, β, γ , Ω ∈ R. Note that [2] To find a solution of (C.1) in the form of one can use formulas [28, (24), (26)] to derive