Coupled discrete fractional-order logistic maps

: This paper studies a system of coupled discrete fractional-order logistic maps, modeled by Caputo’s delta fractional difference, regarding its numerical integration and chaotic dynamics. Some interesting new dynamical properties and unusual phenomena from this coupled chaotic-map system are revealed. Moreover, the coexistence of attractors, a necessary ingredient of the existence of hidden attractors, is proved and analyzed.


Introduction
Nonlinear phenomena are difficult to describe by models analysis based only on smoothness; thereby, fractional calculus has been used to model many such processes for which the standard integer-order derivatives cannot be applied adequately.The generalization of the concept of derivatives of non-integer values dates back to the beginning of the theory of differential calculus, while the rapid development of the theory of fractional calculus started from the work of Euler, Liouville, Riemann, Letnikov, and so on [1,2].In the past, new results of fractional modeling and applications were reported every year.The fractional derivatives and integrals are useful in engineering and mathematics, being helpful for scientists and researchers working with real-life applications (see, e.g., [3,4]).
It is well-known that the classical derivative of a continuous-time periodic function is a periodic function with the same period.However, with respect to derivatives of fractional order, this is not necessarily the case [5][6][7][8][9][10][11][12].The non-periodicity of solutions in fractional-order (FO) systems was first discovered by engineers (see, e.g., [7]), and then proved by mathematicians (see, e.g., [5,12]).Generally, FO systems have no non-constant periodic solutions by their nature (verified, e.g., using Laplace or Z transformations).Nevertheless, both continuous and discrete FO systems may have asymptotic periodic solutions.However, just like continuous FO systems, the periodicity of solutions in discrete FO systems is a delicate issue [10, [13][14][15][16][17].As a consequence, all reported results based on the "periodicity" of continuous-or discrete-time autonomous FO systems became questionable.
In this paper, orbits apparently indicating some regular behavior will be called "periodic-like" orbits.Recall that for some values of the bifurcation parameter of the integral-order (IO) logistic map, unstable periodic orbits (UPOs) will emerge, leading to chaos.For FO logistic maps, this will be referred to as "chaotic-like" behavior.In [18], these kinds of periodic-like orbits are called "numerically periodic orbits".It is also wellknown that in the theory of dynamical systems, every emerging abrupt period-doubling phenomenon is considered as bifurcation.Therefore, in this paper, the term bifurcation or bifurcation diagram is understood in the above sense of a periodic-like phenomenon.
On the other hand, from a computational point of view, based on the complexity or simplicity in finding a basin of attraction in the phase space, it is natural to consider the following classification of attractors: self-excited attractors, which can be revealed numerically by integrating the systems with initial conditions within small neighborhoods of unstable equilibria, and hidden attractors, which have the basins of attraction not connected with any equilibria [19][20][21][22][23]. Examples of hidden attractors in continuous-time FO systems exist in some classical systems, such as the Rabinovich-Fabrikant system [24][25][26], Hopfield neuronal system [27], economic system [28], hyperchaotic discontinuous system [29], and so on [30].
In a numerical approach, the need for previous history in numerical integration requires a trade-off between the calculation time and the approximation precision.This is a basic principle for some classical numerical methods for fractional differential equations (FDEs), such as the Adams-Bashforth-Moulton method [31], for which a tutorial can be found in [32].Finally, it should also be noted that the mathematical theory of FDEs is still quite limited today, although the subject has been studied since 1956 [33].
With all the above motivations and background, the present paper is devoted to studying a system of discrete coupled FO logistic maps, with respect to its numerical bifurcation analysis and hidden attractor search, which reveals some very interesting new dynamical properties and unusual phenomena.This paper will also discuss the stability of discrete FO equations, where [34][35][36] can be referred to for more details.

A Model of Coupled FO Logistic Maps
In [37], the bistability of some "aggregates" of logistic maps with excitation-type coupling is studied.One such model is composed by two functional units, a neuron or a group of neurons (voxels), as a discrete nonlinear oscillator with two possible states: active (meaning one type of activity) or not (meaning other type of activity).A reasonable modality to take the most elemental local nonlinearity is, for instance, the logistic evolution: The first equation with the coupled functional units (1) refers to excitation coupling, while the second equation, to inhibition coupling.These kinds of models are proposed in [37] to mimic the waking-sleeping bistability and even multistability found in brain systems (see [37] for details).
In this paper, the Fractional Order (FO) variant of the system (1) is considered.
To obtain the FO form, let N a = {a, a + 1, a + 2, . ..}.Then, for q > 0 and q ∈ N, the q-th Caputo-like discrete fractional difference of a function u : for t ∈ N a+n−q and n = [q] + 1.In (4), ∆ n is the n-th order forward difference operator, and ∆ −q a represents the fractional sum of order q of u, namely, with the falling factorial t (q) in the following form: .
Note that the fractional operator ∆ −q a maps functions defined on N a to functions on N a+q .
For the case considered in this paper, q ∈ (0, 1), when n = 1, ∆u(s) = u(s + 1) − u(s), and Caputo's fractional difference, denoted hereafter ∆ q * , becomes Assuming the starting point of the fractional sum (3), a = 0, one can consider the following discrete autonomous Initial Value Problem (IVP) of FO (the non-autonomous case can be found in [40,41]): for q ∈ (0, 1) and f is a continuous map.
The solution of the IVP ( 5) is given by [40,41] A convenient iterative form of the integral ( 6) is [35] Now, one can introduce the FO variant of the system (1) modeled by the Caputo delta fractional difference equation: where p ∈ R is a parameter.In this case, and the integral (7) has the following form, which will be used to integrate and simulate the system dynamics later: (10)

Bounds and Global Dynamics
The theory of orbits for discrete FO systems is technically rather sophisticated in general and is still under development.Additionally, finding either bounds of solutions or qualitative properties of global dynamics are difficult tasks.
From (10), the degree k of polynomials (11) is The formulas a nj (x(0), y(0)) and b nj (x(0), y(0)) are difficult to find explicitly for general n, but nevertheless, one can find that Using [42] (Proposition 1), one obtains the following result.
Concerning global dynamics, the following result holds (for the dynamical behavior of discrete-time linear FO systems, see [43]).
for some n 0 ∈ N, where S is one of the following subsets of R 2 : x with δ > 0. Proof.
The whole proof is thus completed.
Remark 1. Theorem 2 asserts the non-existence of an attractor of (8) within any subsets shown in (13).

Non-Existence of Hidden Attractors
To numerically find hidden attractors, it is necessary to determine the stability of the system equilibria.
The Jacobian of the function ( 9) evaluated at equilibria E i , i = 1, 2, . . ., 5, gives the spectrum of eigenvalues σ as shown in Table 1.

E σ(J)
The stability of the linearized FO system implies the stability of the nonlinear FO system (8), which conforms to [44] (Theorem 1.4).The system (8) is asymptotically stable if all the eigenvalues belong to the following set S q : where λ denotes the argument of the eigenvalue.Obviously, if there exist eigenvalues not belonging to S q , then the system is unstable.
Proof.Consider E 1 , for which the arguments of eigenvalues {e 1 , e 2 } are λ 1,2 = ± π 2 and |z| = p 4  3 .For both eigenvalues {e 1 , e 2 }, the first inequality in S q becomes (see the first column in Table 2) The second inequality in S q is verified for both arguments λ 1,2 .
For E 5 , λ 1,2 = −π, and the first inequality reads 4p < 2 cos 0 2 − q q = 2 q , which is the result in the last column of Table 2.The second inequality is verified for all q ∈ (0, 1).
Because for equilibria E 2,3,4 , one argument is zero, the second inequality is not verified, so the eigenvalues do not belong to S q , e / ∈ S q .Therefore, E 2,3,4 are unstable and the system is unstable for all q ∈ (0, 1) and p (the second column in Table 2).
Unstable for all p and q ∈ (0, 1) stable for p < 2 q−2 In the following simulations, unless otherwise mentioned, 3500 iterations were performed for all examples, and the last 600 points are plotted.
The stability regions in the plane (q, p) for equilibria E 1,5 are plotted in Figure 1a,b.The region where all equilibria are unstable is shown in Figure 1c (light brawn).Figure 1d,e presents orbits starting from initial conditions near equilibria E 1 and E 5 .The spiral orbit of E 1 corresponds to complex eigenvalues of the system at this equilibrium, while the orbit of E 5 corresponds to real eigenvalues., with p = 0.5 and q = 0.5 (Figure 1a) from initial condition close to E 1 ; (e) phase plot of a representative orbit from a point within the stability region of E 1 , with p = 0.1 and q = 0.5 (Figure 1b) from initial condition close to E 5 .
It is observed that, sometimes, possible system dynamics are richer than what can be revealed through examining bifurcation diagrams (BDs).As experienced, the BDs are suggestive and they show what parameter values the system can take on, and therefore helps to identify potential hidden attractors coexisting with other attractors.Consider, for example, the BDs of the component x vs. parameter q for fixed p = 0.55 with initial condition (x 0 , y 0 ) = (0.1, 0.3) (Figure 2a), and for q = 0.4 vs. parameter p ∈ (0.3, 0.6) (the largest range of p for q = 0.4) (Figure 2b).As expected, Figure 2 reveals the influence of the IO logistic map.However, there are some significant differences, which will be illustrated and discussed next.The first interesting phenomenon observed is related to the apparent "bifurcation points" (see the zoomed in rectangular region in Figure 2a).This "explosion" takes place for a relatively large range of q values, which was also found from some other discrete FO systems (e.g., [45,46]).
Note that this system cannot be numerically characterized in some regions of the (p, q)-space.There exist some values of p and q for which the orbits are unbounded for whatever initial conditions (x 0 , y 0 ), as shown by the BD in Figure 2b, where for q = 0.4, the upper bound of the range of p values could only be chosen to be p = 0.6 (see also Theorem 2).Now, consider the BDs with respect to q ∈ [0.1, 0.5] and p = 0.4, for two different initial conditions (0.1, 0.3) (blue) and (0.5, 0.5) (red) (Figure 3a), and with respect to p ∈ [0.3, 0.5] and q = 0.295, for initial conditions (0.1, 0.3) (blue) and (0.9, 0.6) (red) (Figure 3b).
Within the BD with respect to q, we denote the bifurcation sets of points corresponding to a single initial condition (Poincaré vertical slices through BD) as bifurcative sets (BSs).For better visualization, within the BD in Figure 3a, the BS corresponding to (0.1, 0.3) is colored blue, while the BD corresponding to (0.5, 0.5) is colored red.Similarly, the BD with respect to p, shown in Figure 3b, is composed of two BSs.To each point on the por q-axis, there is a vertical line of points, colored red or blue on the corresponding BS.
The dotted line in Figure 3a indicates a chosen representative case, with q = 0.28 for p = 0.55.Two different attractors can be seen in the two red and blue BSs: a four-periodlike attractor (red bullets numbered 1, 2, 3, and 4; see also Figure 4a, where the light red lines indicate the way in which the points 1, 2, 3, 4 are visited) and a two-band, chaotic-like attractor (dark blue tick lines; see also Figure 4b).The zoomed region around the point 1 (Figure 4c) shows the last 600 points, which reveals the slow convergence of this orbit towards a regular-like state.The schematic arrows marked on the time series in Figure 4d indicate the order in which the points 1, 2, 3, and 4 are visited by the orbit.The BD for p = 0.55 vs. q with two BSs: one obtained for (x 0 , y 0 ) = (0.1, 0.3), in Figure 2a (red plot), and one for a second initial condition (x 0 , y 0 ) = (0.5, 0.5) (blue plot); (b) the BD for q = 0.4 vs. p ∈ (0.3, 0.6) with two BSs: one for (x 0 , y 0 ) = (0.1, 0.3), in Figure 2b (red plot), and a new one for (x 0 , y 0 ) = (0.9, 0.6) (blue plot).Dotted line in the BD in Figure 2b indicates the existence of two different attractors: a four-periodic-like orbit (red bullets) and a two-period chaotic band orbit (dark blue segments).From experience, this "coexistence" of attractors or "multistability" suggests the possible existence of hidden attractor(s).However, instead of the two initial conditions (0.1, 0.3) (red BS) and (0.5, 0.5) (blue BS) (Figure 5a), if one considers three initial conditions (0.1, 0.3) (red BS), (0.5, 0.5) (blue BS) and (0.01, 0.7) (green plot) (Figure 5b), then one can see that there exist three BSs, which allow to find three possibly different attractors (see dotted lines I, II, and III within the three attractors) and this phenomenon seems to continue indefinitely.To verify the influence of the maximum number of iterations on this phenomenon, consider two BDs generated with the above same three initial conditions, but with different numbers of iterations, 3500 and 5000 (Figure 6a,b), respectively.The iteration number affects the shape of BSs only slightly (compare the green two-period chaotic bands generated after 3500 iterations in Figure 6a and the one-period chaotic band in Figure 6b generated after 5000 iterations).However, the existence of different BSs is not affected by the maximum number of iterations.
To this end, one can conclude that the BSs are non-invariant with respect to the initial conditions and, in fact, their positions change significantly with initial conditions in the fractional order space.Similarly, this happens also in the parameter space.While in the parameter space, the non-invariance is evident for all parameter values p, in the fractional order space, this property is conceived only for small values of q (once q grows over q ≈ 0.75, the phenomenon vanishes, see the right vertical dotted line in Figure 3a).The above-described numerical simulations on the complex dynamics of system (8) are summarized as follows.
Proposition 2. The BDs of the FO system (8), obtained with the discretization (7), presents non-invariance with respect to the initial conditions in both the fractional order space and the parameter space.Remark 2. Some interesting observations are worth highlighting.
(i) The shapes of BSs approximately preserve the shapes for different initial conditions, but move along the p-axis.
(ii) This delay-like phenomenon with respect to the initial conditions (the BSs seem to move "forward" or "backward" on the BDs (see the dotted lines I, I I, and I I I in Figure 5, as the initial conditions are changing) was already found in a continuous-time FO system [47], where the "delay" was observed with respect to the integration step-size of the numerical method used.(iii) It is interesting to compare the above results with the case of the Feigenbaum attractor of the IO logistic map x(n + 1) = px(n)(1 − x(n)), for the limiting value p ∞ = 3.569946 . . .[48], which, however, is not an attracting set and for which there is no sensitive dependence on initial conditions.
By summarizing the investigation in this work, it is concluded that, because of the mentioned dependence on initial conditions, it is impossible to find hidden attractors in the FO system (8) by numerically searching the paths of system orbits by testing initial conditions within neighborhoods of equilibria.

Discussion
In this paper, the system of coupled logistic maps modeled by Caputo's delta fractional difference was studied, both analytically and numerically.The system boundedness and global dynamics were analyzed in detail.Extensive numerical simulations were performed on the system dynamics, revealing the impossibility of finding hidden attractors by numerically testing the orbits starting from initial conditions within neighborhoods of equilibria.The main reason appears to be that, at least for the considered system, the BSs forming the BDs for small values of q (about q < 0.75) are different for different initial conditions; thereby, the existence of hidden attractors cannot be realized numerically, as typical for continuous integer-order systems.It seems that the coexistence of attractors as necessary for the existence of hidden attractors cannot be well-defined for such a discrete fractional-order system, perhaps also for other discrete fractional-order systems.In general, the phenomenon of "coexistence of attractors" for discrete fractional-order systems needs further investigation.
A possible future research direction is to consider k-periodic problems given by the condition u(k) = u(0) for k ≥ 1, regarding the existence, uniqueness, and bifurcation of solutions, similarly to the periodic boundary value problem x(0) = x(T) for functional differential equations [49][50][51].

Figure 1 .
Figure 1.Stability regions of equilibria.(a) Stability region of E 1 (grey plot); (b) Stability region of E 5 (grey plot); (c) Instability region of all equilibria (light red plot); (d) Phase plot of a representative orbit from the point within the stability region of E 1, with p = 0.5 and q = 0.5 (Figure1a) from initial condition close to E 1 ; (e) phase plot of a representative orbit from a point within the stability region of E 1 , with p = 0.1 and q = 0.5 (Figure1b) from initial condition close to E 5 .

Figure 6 .
Figure 6.The BDs of system (8) for p = 0.55, obtained with the same three different initial conditions but with different maximum iteration numbers: (a) 3500 iterations; (b) 5000 iterations.