On the stability domain of a class of linear systems of fractional order

In this paper, the shape of the stability domain S^q for a class of difference systems defined by the Caputo forward difference operator D^q of order q\in (0, 1) is numerically analyzed. It is shown numerically that due to of power of the negative base in the expression of the stability domain, in addition to the known cardioid-like shapes, S^q could present supplementary regions where the stability is not verified. The Mandelbrot map of fractional order is considered as an illustrative example. In addition, it is conjectured that for $q<0.5$, the shape of S^q does not cover the main body of the underlying Mandelbrot set of fractional order as in the case of integer order.


Introduction
In the last decade, fractional calculus has been considered as an important tool in many scientific and engineering fields. The basic theory of fractional calculus modeling and control systems can be found in, e.g., [1][2][3][4][5][6][7]. Studies of their applications are presented in [8][9][10]. Systems with fractional variable orders are analyzed in [11,12]. In terms of the stability analysis of fractional differential equations, one of the important properties that is analysed in order to study the behavior of the considered systems is presented in, e.g., [13][14][15][16], while applications and stability studies of discrete fractional difference equations can be found in [17][18][19].
For commensurate fractional-order systems, several powerful criteria are established. The most well-known Matignon's stability theorem [20] (see also [21]) determines the system stability by searching for the location of the eigenvalues in the complex plane, which represents a starting point for most research into the stability of fractional-order systems.
In [22], a stability criterion for a fractional difference linear system is presented: ∆ q y(n + 1 − q) = Ay(n), n = 0, 1, ..., where A ∈ R d×d and ∆ q is the Caputo forward difference operator (see, e.g., [23,24]). As shown in [22], the stability criterion for (1) is proved to be fully explicit, also involving the decay rate of the solutions.
The asymptotic stability property of (1) for both scalar and vector cases is stated by the following result: Theorem 1 ( [25]). Let q ∈ (0, 1) and A ∈ R d×d . Then, (1) is asymptotically stable if and only if the isolated zeros, off the non-negative real axis, of det( In [22], the aim is to formulate an alternative stability criterion for the integer difference system: ∆y(n) = Ay(n), n = 0, 1, ..., where ∆y(n) = y(n + 1) − y(n) is the standard operator here for q = 1. As is well known, the system (2) is asymptotically stable if and only if all the eigenvalues of I + A, where I is the identity matrix, are located inside the unit circle. If one considers the stability set for Equation (2) in the polar form S = z ∈ C : |z|< − cos(arg z) and |arg z|> π 2 , the stability result can be reformulated in the following form: 22]). The linear difference system (2) is asymptotically stable if and only if λ ∈ S for all the eigenvalues λ of A. In this case, the solutions of (2) decay towards zero exponentially as n → ∞.
In order to introduce an alternative stability criterion for (1), as a direct extension of Theorem (1), consider the following set: Note that the second inequality in (3) represents the Matignon criterion.
The main result in [22] is: 22]). Let q ∈ (0, 1) and A ∈ R d×d . If λ ∈ S q for all the eigenvalues λ of A, then system (1) is asymptotically stable. In this case, the solutions of (1) decay towards zero in such a way that y(n) = O(n −q ) as n → ∞, for any solution y of (1). Furthermore, if λ ∈ C \ cl(S q ) for an eigenvalue λ of A, then (1) is not stable.
Because of the explicit and convenient form of S q , drawn in the complex plane for eigenvalues, Theorem 3 became widely used in practical applications.
In this paper, it is shown that for some empirically found values of q ∈ (0, 1), the form of S q might present other unexpected shapes. In addition, a conjecture related to the shape of S q in the case of the Mandelbrot set of fractional order is introduced.

On the Shape of the Stability Area S q
The parametric equations of the frontier of the stability region S q for a fixed point, z = x+ıy, x, y ∈ R, can be easily drawn in the following form [22]: On the other hand, to analyze the shape of S q in more detail, consider the form (3). As is known, the argument of a complex number z = x + iy, arg z, used in (3), is the angle between the positive real axis and the line joining the origin and the image in the complex plane of z. Usually, the value of the argument is numerically computed as argz = atan2(y, x), where the two-argument function atan2 is an available implemented function in the math libraries of many programming languages.
In contrast to the inverse tangent function, tan −1 (arctan, or atan), the function atan2 computes the principal value of the arctangent of y/x, determining the quadrant of the returned value by using the signs of both arguments. Note that a domain error may occur if both arguments are zero. (Regarding the 1999 ISO C, 1978 ANSI Fortran, or 1982 ISO Pascal Standard standard descriptions and other issues related to atan2 for the case of x = 0 and y = 0, they are referred to on p. 70 [26]). However, as will be shown below, the problem of arg in (3) can be quite complicated, even if both arguments are not zero, as shown next with simple reasoning.
In order to determine in which quadrant the values of E a (q) are situated, let us analyze and sketch the graph of the restriction of the function E a to (0, 1). The derivative of E a is E a (q) = (a − π)/(2 − q) 2 < 0. Therefore, from the monotonicity, one deduces that E a decreases from (a − π)/2, for q ↓ 0, but to a − π, for q ↑ 1. For example, for the limit case of a = 0, E 0 (q) takes values within the range (−π, −π/2), i.e., E 0 (q) are situated in quadrant III (Figure 1a). For the case of a = 0.7π/2, the values of E 0.7 π 2 (q) can be situated in both quadrants III and IV (Figure 1b). For a = π/2, E π 2 moves within in the quadrant IV ( Figure  1c). The case of a = 0.8π is presented in Figure 1d. Therefore, for q ∈ (0, 1) and a ∈ [0, π), the argument E a can take values situated in quadrant III and/or IV, where cos(E) could be either negative or positive.
The relation between q and a on the boundary between the two quadrants III and IV is To obtain the shape of S q via relation (3), consider a variable fixed point with the underlying variable eigenvalue, z = x + iy, to be replaced in (3). Denote hereafter by Γ the curve surrounding S q . The two relations in (3) which define S q represent the implicit inequalities of two variables x and y that can be drawn with implicit functions plotting commands available in software such as Matlab (e.g., fimplicit), Mathematica (e.g., ContourPlot), or by using the free software Desmos [27]. In Figure 2a, the case of S q for q = 0.5 is presented.

Remark 4.
As is known, if in the expression x y , x is negative and y is not an integer, the mathematical situation is somewhat ambiguous. With infinite numeric precision, the correct result of x y is mathematically well-defined without ambiguity. Certain values of y yield an imaginary number as a result, while other values of y result in a real-valued result. Specifically, when x < 0, the result of x y is real-valued exactly when y can be written as a fraction, m/n, where m is an integer and n is an odd integer. Furthermore, the result is positive when m is even but negative when m is odd. When y cannot be written as such, the result would be an imaginary number (see, e.g., [28]). For example, in Matlab, for negative base x and non-integer y, the power function returns complex results. A solution could be to use the function nthroot. However, note that in IEEE floating-point computations x y = exp(y log(x)).
Moreover, the log function domain includes negative and complex numbers, which can lead to unexpected results if used unintentionally. Therefore, x y may be a complex number, if x < 0 and y is noninteger. In these cases, in some software, such as Matlab, imaginary parts of complex arguments are ignored to the detriment of the shape of S q .
To summarize, for those values x and y for which E a < −π/2, they are situated in quadrant III, cos(E) < 0 and, therefore, cos q (E) is not well-defined . The consequence is Table 1: Value of A q at points X i , i = 1, ..., 4, for q = 0.5 (see also Figure 2a). Table 2: Value of A q at points Y i , i = 1, 2, ..., 7, for q = 0.8 (see also Figure 2b).
that, in addition to the expected area, for some values of q, S q could present unexpected additional parts. Thus, for q = 0.8, using Desmos, the domain S 0.8 is presented in Figure  2b.
In these cases, the supplementary domain can be determined via the condition |arg z|> qπ/2 (Matignon condition), for q ∈ (0, 1), and the lines |arg z|= qπ/2 being tangent to S q . However, this is not possible in some more complicated cases (Figure 3d, Section 3).
To verify if S q in the considered cases of q = 0.5 and q = 0.8 contains the stability points, consider the points X i , i = 1, ..., 4 (images of the underlying complex numbers) for the case of q = 0.5 and, Y i , i = 1, 2, ..., 7 for the case of q = 0.8, respectively (see Figure 2), and denote The graph of A q (z) is the boundary curve Γ and S q = z ∈ C : A q (z) < 0 and |arg z|> qπ 2 .
Therefore, the images of z for which A q (z) > 0, or A q (z) are complex values are not S q points. The values of A q at the points X i , and Y i are presented in Table 1 and Table 2, respectively.
Regarding the numerical implementation, if in the atan2 function, x and y have different signs, then atan2 and arctan (atan) could have different values. Thus, if in some software the atan2 function is not implemented, it can be emulated with atan2(y, x) = arctan y x + π 2 sign(y) (1 − sign(x)).
An animation which shows the variation in the stability domain for q ∈ (0, 1), where several such cases appear, is presented as a supplementary video.

Stability of the Mandelbrot Map of Fractional Order
Next, consider the fractional discretization of the Mandelbrot map ∆ q z(t) = f c (z(t + q − 1), z(0) = 0, q ∈ (0, 1), t = 0, 1, 2, ..., where f c (z(t + q − 1) = z 2 (t + q − 1) + c, and z, c ∈ C, with c being the complex parameter. For the initial value problems for real fractional discrete systems, see, e.g., [29,30] (in [31], several properties of the complex Mandelbrot map of fractional order and the matlab code to draw the fractal are presented). The numerical integral of the initial value problem (6) is (see [30] for the solution of real fractional-order systems) For q = 0.85, one obtains the fractional-order Mandelbrot (FOM) set, presented in Figure  3b. The fractal is obtained with the FO mandelbrot.m code [32].
To obtain the stability domain of the stable fixed points (period 1) of the FOM map, consider next, for computational reasons, the problem in the Cartesian plane, with c = c x + ic y , c x , c y ∈ R. Contrary to the integer-order case, where the fixed points are found from the equation f c (z) = z, here, the fixed points are obtained by solving the equation f c (z) = 0, i.e., z 2 + c = 0, with the solutions: z * 1,2 = ±i √ c. For the fixed point z * 1 = i √ c, after some calculations, the eigenvalues are obtained as while, for z * 2 = −i √ c, Consider the eigenvalue λ 1 , the reasoning for λ 2 being similar. Then, the useful parametric representation of Γ in coordinates c x and c y , defining Γ(c), are determined from the following equations (see (4)): x + c 2 y + 2c x = −2 q cos q θ sin((2 − q)θ), |θ|≤ π/2. The solutions, which represent the parametric equations of the curve Γ, are c x = −2 2q−2 cos 2q θ cos(2θ(q − 2)), c y = 2 2q−1 sin(θ(2 − q)) cos 2q θ cos(θ(q − 2)), |θ|≤ π/2.
In Figure 3a, the integer-order Mandelbrot (IOM) set is drawn together with the main cardioid Γ which has the equation (see, e.g., [33][34][35]) or, in parametric form, 4c x = 2 cos θ − cos(2θ), 4c y = 2 sin θ − sin(2θ), |θ|≤ π. 6 The interior of the curve Γ corresponds to the stable fixed points, obtained from the equation f c (z) = z, i.e. z * 1,2 = (1 ± √ 1 − 4c)/2. The curve Γ for the FOM set for q = 0.85, which surrounds the stability domain of the fixed point z * 1 , S 0.85 , is presented in Figure 3c. Note that in both the IO case ( Figure 3a) and the FO case (Figure 3b), the Mandelbrot sets (defined on the parametric plane of c in coordinates (c x , c y )) and the stability domains (defined on the eigenvalues λ in coordinates (λ x , λ y )) are overplotted.
As detailed in Section 2, to a fixed point z * , the eigenvalues λ given by (8) or (9) correspond. For a point λ ∈ Γ, z * loses its stability, while if λ / ∈ S 0.85 , z * is unstable. If z * is asymptotically stable, λ belongs to S 0.85 , and reversely, a point λ within S 0.85 corresponds to a point c for which z * is asymptotically stable.
Let us numerically verify this property. Consider λ = −0.5701+i0.3019 ∈ S 0.85 (magenta point 1 of coordinates (−0.5701, 0.3019) in Figure 3b). Because λ < 0 and λ > 0, one chooses (9) with right-hand side having signs − and +, respectively. To find the corresponding point c of this value of λ, one has to solve the problem (8)  Among several differences and resemblances between the FOM set and the IOM set [36], one can see that the origin of the so-called "Elephant Valley" (see, e.g., [33]) in the case of the IOM set is located at (1/4, 0), while for the FOM set, it is located at (0, 0) (see dotted line in Figure 3).
While for q ≥ 0.5, S q fills the main body of the FOM set well enough, for q < 0.5, the shape of S q does not cover the main body. In Figure 4a-c, for q = 0.5, q = 0.1 and q = 1 × 10 −15 , the numerically near 0 value, q = 1 × 10 −15 , is chosen instead of Γ(0), which is not defined. In this case, in order to obtain more fractal details (which look similar to the IOM set), only 30 iterations were used, compared to 1000 iterations for the other cases). One can see that q = 0.5 is the ultimate value for q. Thus, for q values below this limit, S q is shrinking with respect to the main body of the FOM set. Another characteristic is that the shape of Γ obtained with Desmos ( Figure 4d) presents an additional part (see also Figure 1b). This image was obtained using both representations: the representation (3) and the parametric form (10).
Therefore, the following conjecture is numerically sustained (see also [37]): Conjecture 5. For q < 0.5, the surface of the stability domain S q of the fixed points of the FOM shrinks compared to the main body of the fractal set.
Remark 6. After the first result on the non-periodicity of non-constant solutions of FO, continuous-time systems appeared [38]; this intriguing result has been extended to FO discrete systems too (see, e.g., [39]). Thus, contrary to the IOM set, the head contains points generating two-period stable cycles, and where the bulbs contain the points, generating multiple stable cycles of period 3.4, and so on (Figure 3a). In the case of the FOM set, these apparent periodic orbits do not exist! For example, points in the head (yellow fill in Figure  3b) do not generate periodic orbits. The same situations happen to all bulbs.