**1. Introduction**

In the nineteenth century, fractional calculus had its origin in the generalization of integer order differentiation and integration to non-integer order (fractional order) ones [1–5]. The first treatise concerning fractional calculus was conducted by K.B Oldham and J. Spanier (1974) [6]. In the early stages of fractional calculus, the focus was on continuous time fractional dynamical systems [7]. A few decades later, interest shifted to discrete fractional calculus. Diaz and Olser (1974) gave the first definition of a discrete fractional operator. It was introduced by discretization of the continuous time fractional operator, and this led to representations by recurrence relations. The first definition of fractional order operators was pioneered by Miler and Ross [8]. Following these definitions, and in a series of papers, Atici and Eloe [1,9–11] defined and proved many properties of these operators. Balanu et al. provided significant contributions for the subject concerning stability, and chaos of some discrete fractional systems [4,12–14]. Discrete fractional calculus is a sub-discipline of dynamical system that has been of significant importance since its first emergence. It first appeared in relation to real-world phenomena, and in a broad range of disciplines such as biology, economics, physics, engineering [15], nonlinear optics [16], finance [17], demography [18], medicine [19], and so on. It has become an interesting field of research in the last decades.

The concept of chaos in dynamical systems has attracted researchers due to the unpredictability of the dynamic behaviors of dynamical systems. This means that two systems may start very close, but they move very far for apart. In fact, Li and York [20] were the first researchers who gave a mathematical definition of chaos. More studies and investigations for a variety of chaotic maps have been presented in [21]. In addition, many considerable results have been proposed in [14,21–24]. Those maps possess sensitivity to initial conditions and random-like trajectories. One of the characteristics

of chaos is the existence of at least one positive Lyapunov exponent. Another efficient method for detecting chaos is the 0–1 test. That is, if the output is zero, this means that the dynamic is regular; it is chaotic if the test yields one [25]. The fact that research in discrete fractional chaotic systems is still in development [26–34], and very few difference equations have been considered, was the motivation of our work.

The aim of this paper is to study two new chaotic maps with specific types of fixed points by the application of a new test approach (0–1 test) in order to find out whether or not these systems are chaotic. The dynamic behaviors of the considered chaotic systems are investigated numerically using bifurcation diagrams and Lyapunov exponents. These systems possess an interesting property: symmetry. The rest of the paper is organized as follows. In Section 2, we present certain essential definitions and results from difference fractional calculus. Then, chaotic dynamics of the two systems will be discussed in Section 3. Section 4 is about the chaotic analysis based on the 0–1 test approach. Finally, we conclude with a summary.

## **2. Preliminaries**

We give in this section some basic definitions, and one major theorem from discrete fractional calculus. Let *X* : N*a* → R, with N*a* = {*<sup>a</sup>*, *a* + 1, *a* + 2, ... }, the general *nth* order difference can be written as:

$$
\Delta^n x(t) = \Delta^{n-1} X(t+1) - \Delta^{n-1} X(t) = \sum\_{k=0}^n C\_k^n (-1)^k X(t+n-k). \tag{1}
$$

Prolonging this concept for fractional-order difference, the fractional sum of order *υ*, using the fact that *t*(*υ*) is defined by *t*(*υ*) = <sup>Γ</sup>(*t*+<sup>1</sup>) <sup>Γ</sup>(*<sup>t</sup>*+1−*<sup>υ</sup>*), is given in the following definition.

**Definition 1** ([35])**.** *The υth fractional sum of X can be defined as:*

$$
\Delta\_a^{-\upsilon} X(t) = \frac{1}{\Gamma(\upsilon)} \sum\_{s=a}^{t-\upsilon} (t - \sigma'(s))^{(\upsilon - 1)} X(s), \tag{2}
$$

*where υ* ∈/ N*, with υ* > 0 *is a fractional order, t* ∈ N*a*+*n*−*υ, n* = [*υ*] + 1; *σ*(*t*) = *t* + 1.

**Definition 2** ([35])**.** *The υ-Caputo-like difference, denoted by <sup>C</sup>*Δ*<sup>υ</sup>a<sup>X</sup>*(*t*)*, of X*(*t*) *is defined as:*

$$\, \_a^C\Delta\_a^v X(t) = \Delta\_a^{-(n-v)} \Delta^n X(t) = \frac{1}{\Gamma(n-v)} \sum\_{s=a}^{t-(n-v)} (t-\sigma'(s))^{(n-v-1)} \Delta\_a^n X(s). \tag{3}$$

**Theorem 1** ([35])**.** *The corresponding discrete integral equation of the delta fractional difference equation:*

$$\begin{cases} ^C\Delta\_a^v u(t) &= f(t+\nu-1, u(t+\nu-1)) \\ ^C\Delta^k u(a) &= u\_k \quad n = [\nu]+1, \; k = 0, \ldots, n-1, \end{cases} \tag{4}$$

*is given by:*

$$u(t) = u\_0(t) + \frac{1}{\Gamma(\nu)} \sum\_{s=a+n-\nu}^{t-\nu} (t - \sigma(s))^{(\upsilon - 1)} f(s + \upsilon - 1, u(s + \upsilon - 1)), \ t \in \mathbb{N}\_{a+n\nu} \tag{5}$$

*where:*

$$
\Delta u\_0(t) = \sum\_{k=0}^{n-1} \frac{(t-a)^k}{k} \Delta^k u(a). \tag{6}
$$

#### **3. Fractional–Order Maps with Closed Curve Fixed Points**

In this section, we study the chaotic behavior of two new two-dimensional fractional systems with closed curve fixed points. In [36], a new class of two-dimensional chaotic maps with different types of closed curve fixed points are constructed. Examples of these systems are given by:

$$\begin{cases} \mathbf{x}\_{k+1} &= \mathbf{x}\_k + \mathfrak{a}f(\mathbf{x}\_k, \mathbf{y}\_k), \\ \mathbf{y}\_{k+1} &= \mathbf{y}\_k + f(\mathbf{x}\_k, \mathbf{y}\_k)(\beta \mathbf{x}\_k + \gamma y\_k + \delta \mathbf{x}\_k^2 + \varepsilon \mathbf{y}\_k^2 + \varepsilon \mathbf{x}\_k y\_k + \theta), \end{cases} \tag{7}$$

where *xk*, *yk* are system states, {*<sup>α</sup>*, *β*, *γ*, *δ*, ,*ε*, *θ*} are system parameters and *f*(*<sup>x</sup>*, *y*) is a nonlinear function which can be written as

$$f(x,y) = (\frac{x}{m})^p + (\frac{y}{n})^p - r^2,\tag{8}$$

where *m* > 0, *n* > 0, *r* > 0, and *p* is an integer order such that *p* ≤ 2. The authors showed that the fixed points of these systems are nonhyperbolic for different values of *m*, *n*, *p*, and *r*, and they presented the phase-basin portrait of the chaotic attractors. Also, they explored them numerically by Lyapunov exponents and Kaplan–Yorke dimension.

#### *3.1. Fractional–Order Map with Square-Shaped Fixed Points*

Here we consider the system (7) when *m* = 1, *n* = 1, *p* = 12, and *r* = 1. In this case, the map is described by the following equation:

$$\begin{cases} \mathbf{x}\_{k+1} &= \mathbf{x}\_k - \mathfrak{a} (\mathbf{x}\_k^{12} + y\_k^{12} - 1), \\ y\_{k+1} &= y\_k + \beta \mathbf{x}\_k y\_k (\mathbf{x}\_k^{12} + y\_k^{12} - 1), \end{cases} \tag{9}$$

where *α* and *β* are bifurcation parameters. Based on the above equation, we construct a new fractional-order map by introducing the *υ*-Caputo like difference operator as follows:

$$\begin{cases} ^C\Delta\_a^{v\_1}\mathbf{x}(t) &= -\mathfrak{a}(\mathbf{x}^{12}(t-1+\upsilon\_1) + \mathfrak{y}^{12}(t-1+\upsilon\_1)-1), \\ ^C\Delta\_a^{v\_2}\mathbf{y}(t) &= \mathfrak{beta}\mathfrak{x}(t-1+\upsilon\_2)\mathfrak{y}(t-1+\upsilon\_2)(\mathbf{x}^{12}(t-1+\upsilon\_2) + \mathfrak{y}^{12}(t-1+\upsilon\_2)-1), \end{cases} \tag{10}$$

where *υ* = (*<sup>υ</sup>*1, *<sup>υ</sup>*2) is the fractional order. The vector field *<sup>F</sup>*(*<sup>x</sup>*, *y*) associated to the system (10) possesses a symmetry with respect to the *x*-axis. That is, if we let *<sup>G</sup>*(*<sup>x</sup>*, *y*)=(*<sup>x</sup>*, <sup>−</sup>*y*), then *<sup>F</sup>*(*G*(*<sup>x</sup>*, *y*)) = *<sup>G</sup>*(*F*(*<sup>x</sup>*, *y*)). This property will be reflected in the fixed points and attractors of system (10). The fixed points of the new fractional-order system (10) are obtained by taking *<sup>C</sup>*Δ*<sup>υ</sup>*1 *a x* = *<sup>C</sup>*Δ*<sup>υ</sup>*2 *a y* = 0, as follows:

$$\begin{cases} -a(\mathbf{x}^{12} + y^{12} - 1) &= 0, \\ \beta xy(\mathbf{x}^{12} + y^{12} - 1) &= 0. \end{cases} \tag{11}$$

From Equation (11) we must have

$$x^{12} + y^{12} - 1 = 0.\tag{12}$$

Consequently, the system (10) has square-shaped fixed points. By choosing suitable system parameters the fractional-order map (10), we can generate hidden attractors with square-shaped fixed points.

*Symmetry* **2020**, *12*, 756

In order to further investigate the dynamical behavior of the new fractional-order map (10), we need to determine the corresponding numerical formula. Using Theorem 1, the equivalent discrete integral equation of the fractional-order map (10) can be written as:

$$\begin{cases} \mathbf{x}(t) = \mathbf{x}(a) - \frac{\mathbf{a}}{\Gamma(v\_1)} \sum\_{s=a+1-v\_1}^{t-v\_1} (t-s-1)^{(v\_1-1)} (\mathbf{x}^{12}(s-1+v\_1) + \mathbf{y}^{12}(s-1+v\_1) - 1), \\\ y(t) = y(a) + \frac{\theta}{\Gamma(v\_2)} \sum\_{s=a+1-v2}^{t-v\_2} (t-s-1)^{(v\_2-1)} \\\ \qquad\qquad\qquad\qquad\qquad \mathbf{x}(s-1+v\_2)y(s-1+v\_2) (\mathbf{x}^{12}(s-1+v\_2) + \mathbf{y}^{12}(s-1+v\_2) - 1), \end{cases} \tag{13}$$

in which *t* ∈ N*a*+1−*<sup>υ</sup>* and *a* is the starting point. Replacing the discrete kernel function (*<sup>t</sup>*−*<sup>σ</sup>*(*s*))(*<sup>υ</sup>*−<sup>1</sup>) <sup>Γ</sup>(*υ*) by <sup>Γ</sup>(*<sup>t</sup>*−*<sup>s</sup>*) <sup>Γ</sup>(*υ*)Γ(*<sup>t</sup>*−*s*−*υ*+<sup>1</sup>) and assuming that *a* = 0, the above equation is converted to:

$$\begin{cases} \mathbf{x}(n) = \mathbf{x}(0) - \frac{\mathbf{a}}{\Gamma(v\_1)} \sum\_{j=1}^{n} \frac{\Gamma(n-j+v\_1)}{\Gamma(n-j+1)} (\mathbf{x}^{12}(j-1) + \mathbf{y}^{12}(j-1) - 1), \\\mathbf{y}(n) = \mathbf{y}(0) + \frac{\beta}{\Gamma(v\_2)} \sum\_{j=1}^{n} \frac{\Gamma(n-j+v\_2)}{\Gamma(n-j+1)} \left(\mathbf{x}(j-1)\mathbf{y}(j-1) \left(\mathbf{x}^{12}(j-1) + \mathbf{y}^{12}(j-1)\right)\right). \end{cases} \tag{14}$$

where *x*(0) and *y*(0) are the initial condition. This numerical formula will allow us to examine the sensitivity of the fractional-order map throughout the next subsection.

Bifurcation Diagrams and Largest Lyapunov Exponents

In this section, we will study the dynamical behavior of the novel fractional-order map (10) for some specific parameters and initial values. Figure 1 shows the chaotic attractor along with its bifurcation plots with respect to *α* and *β*, and the time evolution of states *x* for *υ* = 1. The chaotic attractor of the fractional-order map (10) under system parameters *α* = 0.2, *β* = 2.4 and initial value (*x*(0), *y*(0)) = (0.26, −0.14) is shown in Figure 1a, where the black points represent the square-shaped fixed points. Figure 1b gives the time evolutions of the variables *x* and *y* which illustrates that the system (10) has two symmetrical states. Figure 1c,d shows the bifurcation plots produced when fixing one of the parameters and varying the second. These bifurcation diagrams are constructed for 50 initial points. We see that the maximum chaotic range is observed for *α* ≈ 0.35 and *β* ≈ 2.455. Also, we notice that the range of *β* for which we obtain a chaotic behavior is very short. Hence, it is more interesting to visualize the effect of *α* on the map's dynamics.

**Figure 1.** Numerical analysis of the fractional-order map (10) for *ν* = 1. (**a**) Square-shaped fixed point of the fractional-order map are illustrated with black dots and the hidden chaotic attractor in red. (**b**) Evolution states of the fractional-order map (10) of *x* and *y*. (**c**) Bifurcation diagram of the fractional-order map (10) versus *α* for *β* = 2.4. (**d**) Bifurcation diagram of the fractional-order map (10) versus *β* for *α* = 0.2.

To investigate the sensitivity of the fractional-order map (10) with respect to the bifurcation parameter *α*, we fix *β* = 2.4 and we vary *α* in the range [0, 0.35]. The bifurcation diagrams and largest Lyapunov exponents under the fractional order *υ*1 = *υ*2 = 0.95 and *υ*1 = *υ*2 = 0.9 are given in Figure 2. By comparing Figure 2a,b with Figure 2c,d, we find that the interval where chaos exists shrinks as *υ* decreases. In particular, when the value of *υ* = (*<sup>υ</sup>*1, *<sup>υ</sup>*2) changes from 0.95 to 0.9, the areas of chaotic motion shrinks from [1.022, 1.116] ∪ [1.164, 1.216] to [0.2405, 0.2429] ∪ [0.3126, 0.3423]. Also, the transient states region increases.

**Figure 2.** (**a**) Bifurcation diagram of the fractional-order map (10) versus *α* for *υ*1 = *υ*2 = 0.95. (**b**) Largest Lyapunov exponents corresponding to (**a**). (**c**) Bifurcation diagram of the fractional-order map (10) versus *α* for *υ*1 = *υ*2 = 0.9. (**d**) Largest Lyapunov exponents corresponding to (**c**).

Now, reset the parameter *β* = 2.2 and keep the parameter *α* as *α* = 0.2, Figure 3 shows the bifurcation diagram of the fractional-order map (10) with the decreasing of *υ*2 where *υ*1 = 1. When *υ*2 decreases from 1, the fractional-order map (10) presents a periodic state until it enters into chaos at 0.9035. It is noticed that the periodic state of the integer order map becomes chaotic as *υ*2 decreases. Another interesting dynamic is observed if we increase the value of *α* to 0.23 and decrease the value of *β* to 0.5. In this case, two chaotic attractors are observed in the fractional-order map with respect to the fractional orders *υ*1 = 1 and *υ*2 = 0.01, as shown in Figure 4, where the red color attractor is yielded from initial value (0.1, 0.3) and blue color attractor is yielded from initial value (0.1, <sup>−</sup>0.3). Therefore, we conclude that the fractional order improves the complexity of the chaotic motion of the original system (7) with square-shaped fixed points.

**Figure 3.** Bifurcation diagram for system (10) in the *x* − *υ*2 plane for *α* = 0.2 and *β* = 2.2.

**Figure 4.** Coexisting chaotic attractor of the fractional-order map (16) for *α* = 0.2 and *β* = 2.2.

#### *3.2. A New Fractional Map with Rectangle-Shaped Fixed Points*

Now, we consider system (7) when *m* = 5, *n* = 1, *p* = 12, *r* = 1. In this case we have the following system:

$$\begin{cases} \mathbf{x}\_{k+1} &= \mathbf{x}\_k + ((\frac{\mathbf{x}\_k}{5})^{12} + y\_k^{12} - 1), \\ y\_{k+1} &= y\_k + a y\_k ((\frac{\mathbf{x}\_k}{5})^{12} + y\_k^{12} - 1), \end{cases} \tag{15}$$

where *x* and *y* are the system variables. Using the same argumen<sup>t</sup> as with the previous system, we obtain the following fractional-order map:

$$\begin{cases} ^C\Delta\_a^{v\_1}x(t) = \left(\frac{x(t-1+v\_1)}{5}\right)^{12} + y^{12}(t-1+v\_1) - 1, \\ ^C\Delta\_a^{v\_2}y(t) = ay(t-1+v\_2)((\frac{x(t-1+v\_2)}{5})^{12} + y^{12}(t-1+v\_2) - 1), \end{cases} \tag{16}$$

where *υ* ∈ (0, 1) is the fractional order. The fixed points of the new fractional-order map (16) are obtained by solving the following system of equation:

$$\begin{cases} \left(\frac{x}{5}\right)^{12} + y^{12} - 1 &= 0, \\ ay((\frac{x}{5})^{12} + y^{12} - 1) &= 0. \end{cases} \tag{17}$$

Thus, we must have

$$(\frac{\chi}{\mathfrak{F}})^{12} + y^{12} - 1 = 0.\tag{18}$$

Therefore, the novel fractional-order map (16) has rectangle-shaped fixed points.

To study the dynamic behavior of the novel fractional-order map (16), we define the numerical formula as:

$$\begin{cases} x(n) = x(0) + \frac{1}{\Gamma(v\_1)} \sum\_{j=1}^{n} \frac{\Gamma(n-j+v\_1)}{\Gamma(n-j+1)} \left( \left( \frac{x(j-1)}{5} \right)^{12} + \left( \frac{y(j-1)}{12} \right)^{12} - 1 \right), \\ y(n) = y(0) + \frac{n}{\Gamma(v\_2)} \sum\_{j=1}^{n} \frac{\Gamma(n-j+v\_2)}{\Gamma(n-j+1)} \left( y(j-1) \left( \left( \frac{x(j-1)}{5} \right)^{12} + y^{12}(j-1) - 1 \right) \right). \end{cases} \tag{19}$$

where *x*(0) and *y*(0) are the initial conditions.

#### Bifurcation and Largest Lyapunov Exponents

In this section we investigate the dynamic behavior of the new fractional-order map (16) using numerical simulations. When *α* = 2.2 the fractional-order map (16) shows the hidden chaotic attractor as depicted in Figure 5. It can be shown that it also has a symmetry with respect to the *x*-axis. Figure 6 shows the bifurcation diagrams and largest Lyapunov exponents of the fractional-order map (16). According to Figure 6, it is very clear that the dynamic behavior of system (16) goes from periodic state to chaotic state with the increase of the control parameter *α*. Obviously, as the fractional order *υ* = (*<sup>υ</sup>*1, *<sup>υ</sup>*2) decreases from 0.998 to 0.995, the bifurcation diagram gradually shrinks to the left.

**Figure 5.** Chaotic attractor of the fractional-order map (16) for *α* = 2.2.

**Figure 6.** (**a**) Bifurcation diagram of the fractional-order map (16) versus *α* for *υ*1 = *υ*2 = 0.988. (**b**) Largest Lyapunov exponents corresponding to (**a**). (**c**) Bifurcation diagram of the fractional-order map (16) versus *α* for *υ*1= *υ*2 = 0.995. (**d**) Largest Lyapunov exponents corresponding to (**c**).

#### *3.3. 0–1 Test*

In order to further confirm the influence of the fractional order *υ* = (*<sup>υ</sup>*1, *<sup>υ</sup>*2) on the properties of the fractional-order maps, we apply a new approach proposed by Gottwald and Melbourne called the 0–1 test method [25]. Unlike other methods, the test does not involve phase space reconstruction, and we can applied it directly to the series data *<sup>x</sup>*(*n*). The output is either 0 or 1 depending on whether the nature of the system is regular or chaotic. Naturally, we do not obtain the values 1 and 0 exactly, so the test remains valid when *K* is close enough to these values.

Here, we applied the 0–1 test method directly to the state *x*(*n*) that was obtained from the numerical formulas (14) and (19), respectively. For an arbitrary constant *c* in (0, *<sup>π</sup>*), we calculate the translation variables

$$p\_c(n) = \sum\_{j=1}^{n} \mathbf{x}(j)\cos(j\mathbf{c}), \quad q\_c(n) = \sum\_{j=1}^{n} \mathbf{x}(j)\sin(j\mathbf{c}), \ n = 1, 2, \ldots, N. \tag{20}$$

*Symmetry* **2020**, *12*, 756

The dynamics of the translation components (*pc*, *qc*) provide a visual test. Basically, if the dynamic is regular then the behavior of trajectories in the (*pc* − *qc*) plane is bounded, whereas if the dynamic is chaotic then the (*pc* − *qc*) trajectories depict Brownian like behavior. In order to examine the diffusive behavior of *pc* and *qc*, we define the mean square displacement as:

$$M\_{\varepsilon}(n) = \lim\_{N \to \infty} \frac{1}{N} \sum\_{j=1}^{N} \left( \left( p\_{\varepsilon}(j+n) - p\_{\varepsilon}(j) \right)^{2} + \left( q\_{\varepsilon}(j+n) - q\_{\varepsilon}(j) \right)^{2} \right), \quad n \le \frac{N}{10} \tag{21}$$

with the asymptotic growth rate of

$$K\_{\mathfrak{c}} = \lim\_{n \to \infty} \frac{\log M\_{\mathfrak{c}}(n)}{\log n}. \tag{22}$$

In practice, the final result *K* can be determined numerically by computing the median of *Kc*. When *K* 1, the behavior is classified as chaotic and when *K* is close to 0 the behavior is regular.

In Figure 7, we present the 0–1 test of the fractional-order map (10), when *υ* = 0.95 and *α* = 0.2. It is clear that the output *K* converges to the values 1, which implies that system (10) is chaotic. Moreover, when 0.98 ≤ *υ* < 0.998 the fractional-order map (10) has a transient state as shown in Figure 8. As it can be seen that the state first converges into a bounded attractor and then it becomes divergent.

Similarly, for the second fractional-order map (16), we use Equation (19). The 0–1 test method is applied for *υ* = 0.998, and *υ* = 0.995, and the resulting diagrams are plotted in Figure 9. It shows that the asymptotic growth rate *K* approaches 1 as *α* increases, which confirms very well the results in Figure 6.

**Figure 7.** The 0–1 test method of the fractional-order map (10).

**Figure 8.** Transient state of the fractional-order map for *υ* = 0.98 (**a**) when *n* = 570, (**b**) when *n* = 578.

**Figure 9.** The 0–1 test of the fractional-order map (16). (**a**) The asymptotic growth rate versus *α* for *υ*1 = *υ*2 = 0.998. (**b**) The asymptotic growth rate versus *α* for *υ*1 = *υ*2 = 0.995.

## **4. Conclusions**

In this article, we investigated the chaotic behavior of new two-dimensional fractional chaotic maps with closed curve fixed points. Numerical methods, including computations of Lyapunov exponents, bifurcation diagrams, and the 0–1 test have been used to illustrate the complex dynamics of these systems. Through this paper, we have shown that chaos exists in these fractional-order maps

and that the type and range of chaotic behavior is dependent on the fractional order. Results shows that the fractional-order maps are more complex then their integer order counterparts. Also, numerical experiments have shown that the system exhibits the property of coexisting attractors. Future studies will explore more dynamic behaviors of those systems, such as, control and synchronization.

**Author Contributions:** Conceptualization, A.O., F.H.; N.S., and A.-A.K.; Data curation, A.-A.K. and F.H.; A.O. Investigation, F.H., A.O., and N.S.; Methodology, A.O., F.H., Supervision, N.S., A.O., and G.G.; Validation, A.O. and N.S. All authors have read and agreed to the published version of the manuscript.

**Acknowledgments:** The author Adel Ouannas was supported by the Directorate General for Scientific Research and Technological Development of Algeria.

**Conflicts of Interest:** The authors declare no conflict of interest.
