Asymptotic behaviour of solutions to non-commensurate fractional-order planar systems

This paper is devoted to studying non-commensurate fractional order planar systems. Our contributions are to derive sufficient conditions for the global attractivity of non-trivial solutions to fractional-order inhomogeneous linear planar systems and for the Mittag-Leffler stability of an equilibrium point to fractional order nonlinear planar systems. To achieve these goals, our approach is as follows. Firstly, based on Cauchy's argument principle in complex analysis, we obtain various explicit sufficient conditions for the asymptotic stability of linear systems whose coefficient matrices are constant. Secondly, by using Hankel type contours, we derive some important estimates of special functions arising from a variation of constants formula of solutions to inhomogeneous linear systems. Then, by proposing new weighted norms combined with the Banach fixed point theorem for appropriate Banach spaces, we get the desired conclusions. Finally, numerical examples are provided to illustrate the effect of the main theoretical results.


Introduction
Fractional calculus and fractional order differential equations are research topics that have generated a great amount of interest in recent years. For details on their various applications in in science and engineering, we refer the interested reader to the collections [2,3,14,18,19] and the references therein.
To our knowledge, the first contribution in the qualitative study of the fractional order autonomous linear systems was published by Matignon [13]. In that paper, using Laplace transform and the final value theorem, the author has obtained an algebraic criterion to ensure the attractiveness of solutions. The BIBO (bounded input, bounded output) stability for non-commensurate fractional order systems, i.e. for systems whose differential equations are not all of the same order, was investigated by Bonnet and Partington [4], and their result shows that the systems are stable if and only if their transfer function has no pole in the closed right hand side of the complex plane.
Starting from [4], a new difficult task appears: finding the conditions to ensure that the poles of the characteristic polynomial of the system lie on the open left side of the complex plane. Trigeassou et al. [20] have proposed a method based on Nyquist's theorem. In particular, they have derived Routh-like stability conditions for fractional order systems involving at most two fractional derivations. Unfortunately, for higher numbers of differential operators, this approach seems to be unsuitable by its numerical implementation. After that, Sabatier et al. [16] have introduced another realization of the fractional system. This realization is recursively defined and involves nested closed-loops. Based on this realization, they have obtained a recursive algorithm that involves, at each step, Cauchy's argument principle on a frequency range and removes the numerical limitation mentioned in [20] above.
In addition to the algorithmic approach as in [16], a number of analytic approaches have been used to investigate the zeros of characteristic polynomials of systems of fractional order systems. In [12], the stability and resonance conditions are established for fractional systems of second order in terms of a pseudo-damping factor and a fractional differentiation order. The method in [12] has been successfully extended in [26] for a wide class of second kind non-commensurate elementary systems. By the substitution method, a variation of constants formula and the properties of the Mittag-Leffler function in the stable domain, in [10], the authors have shown the asymptotic stability for fractional order systems with (block) triangular coefficient matrices. By combining a variation of constants formula, properties of Mittag-Leffler functions, a special weighted norm type and Banach's fixed-point theorem, Tuan and Trinh [23] have proved the global attractivity and asymptotic stability for a class of mixed-order linear fractional systems when the coefficient matrices are strictly diagonally dominant and the elements on the main diagonal of these matrices are negative. Using the positivity of the system and developing a novel comparison principle, Shen and Lam [17] have considered the stability and performance analysis of positive mixed fractional order linear systems with bounded delays. Tuan et al. [24] have established a necessary and sufficient condition for the asymptotic stability of positive mixed fractional-order linear systems with bounded or unbounded time-varying delays.
Although there have been some articles on mixed fractional order systems as listed above, in our view, the qualitative theory of non-commensurate fractional order systems is still a challenging topic whose development is in its infancy. Even in the simplest case when the coefficient matrix is constant, the current results seem to be far away from a complete characterization of the stability of these systems. In particular, the entire theory for non commensurate systems is far less well developed than the corresponding theory for commensurate systems (i.e. systems all of whose associated differential equations are of the same order) that have been extensively discussed, e.g., in the papers mentioned above or in [7,8] and the references cited therein.
For these reasons, we study in this paper the fractional-order planar system with Caputo fractional derivatives C D α 0 + x(t) = Ax(t) + f (t, x(t)), t > 0, (1) where α = (α 1 , α 2 ) ∈ (0, 1] 2 is a multi-index, A ∈ R 2×2 is a square matrix and f : [0, ∞) × R 2 → R 2 is vector valued continuous function. It is worth noting that for the case f = 0, in [6], by constructing a smooth parameter curve and using Rouché's theorem, Brandibur and Kaslik have provided criteria for the asymptotic stability and for the instability of solutions, respectively. However, these conditions are not explicit and are quite difficult to verify. Motivated by [6], our aim is as follows. First, we want to give sufficient simple and clear conditions that can guarantee the Mittag-Leffler stability of the system (1) in the homogeneous case. Then, by establishing a variation of constants formula, estimates for general Mittag-Leffler type functions, and proposing new weighted norms, we show the asymptotic behavior of the system when the vector field f is inhomogeneous or represents small nonlinear noise around its equilibrium point.
The paper is organized as follows. Section 2 contains a brief summary of existence and uniqueness results for solutions to multi-order fractional differential systems and a variation of constants formula for solutions to fractional order inhomogeneous linear planar systems. Section 3 deals with some properties of the characteristic function to a general fractional order homogeneous linear planar system whose coefficient matrix is constant. Section 4 is devoted to studying important estimates for special functions arising from the variation of constants formula for the solutions. Our main contributions are presented in Section 5 where we show the asymptotic behaviour of solutions to fractional-order linear planar systems and the Mittag-Leffler stability of an equilibrium point to fractional nonlinear planar systems. Numerical examples are provided in Section 6 to illustrate the main theoretical results.
To conclude the introduction, we present some notations that will be used throughout the rest of the paper. In R 2 , we define the norm · by x := max{|x 1 |, |x 2 |} for every x ∈ R 2 . For any r > 0, the closed ball of radius r centered at the origin 0 in R 2 is given by B(0, r) := {x ∈ R 2 : x ≤ r}. The space of all continuous functions ξ : For α ∈ (0, 1) and J = [0, T ] or J = [0, ∞), we define the Riemann-Liouville fractional integral of a function f : J → R as and the Caputo fractional derivative of the order α ∈ (0, 1) of a function f : J → R as where Γ(·) is the Gamma function and d dt is the usual derivative. Letting α = (α 1 , α 2 ) ∈ (0, 1] × (0, 1] be a multi-index and f = (f 1 , f 2 ) with f i : J → R, i = 1, 2, be a vector valued function, we write . See, e.g., [9, Chapter III] and [25] for more details on the Caputo fractional derivative.

Existence and uniqueness of global solutions and exponential boundedness of solutions
Consider the two-component incommensurate fractional-order initial value problem with Caputo fractional derivatives where α = (α 1 , α 2 ) ∈ (0, 1] 2 is a multi-index and f : [0, ∞) × R 2 → R 2 is a continuous function.
with respect to its second variable. Then, for any initial value x 0 ∈ R 2 , the two-component incommensurate fractional-order system (3) has a unique global solution ϕ(·, x 0 ) on the interval [0, ∞).
Then, for any initial value x 0 ∈ R 2 , the two-component incommensurate fractional-order system (3) has a unique global solution ϕ(·, x 0 ) ∈ C ([0, ∞), R 2 ) and where M is some positive constant which depends on x 0 .

The variation of constants formula for the solutions
Consider the non-homogeneous two-component incommensurate fractional-order linear system with initial condition where α = (α 1 , α 2 ) ∈ (0, 1] 2 , A = (a ij ) ∈ R 2×2 is a square real matrix and f = ( for some M > 0 and some γ > 0. Then, we have Due to Theorems 2.1 and 2.2, for any initial condition x 0 ∈ R 2 , the system (4) has a unique exponentially bounded solution in C ([0, ∞), R 2 ). Taking Laplace transform on both sides of the system (4), we obtain the algebraic system where X i (s) and F i (s), i = 1, 2, are the Laplace transforms of x i (t) and f i (t), respectively. By Cramer's rule, we see that and where Q(s) := s α 1 +α 2 − a 11 s α 2 − a 22 s α 1 + det A. Put with l(α) := α 1 + α 2 . Then, with each i ∈ {1, 2}, we obtain where " * " is the Laplace convolution operator.
From the arguments above, the unique solution to the initial value problem (4) has the following form.
3 Some properties of the characteristic function In this paper, we only focus on incommensurate systems, i.e. on systems of the form (1) with α 1 = α 2 , because the case α 1 = α 2 has already been discussed in detail elsewhere [13]. Thus, without loss of generality, we assume 0 < α 1 < α 2 ≤ 1. Our first auxiliary statement in this context deals with functions of the form the characteristic functions of the problems under consideration will be of precisely this structure.
(i) If c < 0 then Q has at least one positive real zero.
(ii) If s ∈ C is a zero of Q then its complex conjugate is also a zero of Q.
(iv) If c > 0, then s = iω with ω > 0 is a zero of Q if and only if where Proof. (i) and (ii) are obvious.
(iii) First, we assume that c = 0. Then Q(0) = 0. Due to the continuity of Q at 0, we can find ε which is small enough such that Q has no zero in {z ∈ C : |z| < ε} . Moreover, because |Q(s)| ≥ |s| α 1 +α 2 − |a| · |s| α 2 − |b| · |s| α 1 − |c|, we have that lim |s|→∞ |Q(s)| = ∞ uniformly for all arg s. This implies that there is a positive real number R such that Q has no zero in the domain {z ∈ C : |z| > R} . Hence, all zeros of Q in {z ∈ C : | arg (z)| ≤ ω} (if they exist) belong to the set Ω := {z ∈ C : ε ≤ |z| ≤ R, | arg (z)| ≤ ω} . Notice that Ω is a compact set and Q is analytic on this domain. If now Q has infinitely many zeros in Ω then, because of the compactness of Ω, the set of zeros has a cluster point. This implies, in view of the analyticity of Q, that Q(s) = 0 for all s which contradicts the definition of Q. Hence, Q has only a finite number of zeros in Ω. This shows that Q has only a finite number of zeros in the domain C if c = 0.
To deal with the case c = 0, we write where P (s) = s α 2 − as α 2 −α 1 − b. By repeating the above arguments for P , the proof is complete.
Corollary 3.2. Assume that a, b, c > 0 and that one of conditions is satisfied where ρ 1 , ρ 2 are defined in (14). Then, the function Q defined in (12) has no purely imaginary zero.
Proof. Consider the system (13). Due to the fact that ρ 2 = 0, this system is equivalent to Thus, we obtain Setting X = ω α 2 , equation (16) takes the form The discriminant of the quadratic equation (17) is Hence, the quadratic equation (17) has no real roots. This implies that the system (13) has no root ω > 0. This together with Lemma 3.1(ii) and Lemma 3.1(iv) shows that Q has no purely imaginary zero.
Moreover, a, b, c, ρ 1 > 0, thus the two roots of the quadratic equation (17) are negative. Hence, in view of the relation X = ω α 2 with 0 < α 2 ≤ 1 between the solution X of (17) and the solution ω of (13), the system (13) has no root ω > 0. Using Lemma 3.1(ii) and Lemma 3.1(iv), we see that Q has no purely imaginary zero.
Recall that if c ≤ 0, then Q has at least one non-negative real zero, which precludes any kind of stability. Thus, in this section, we only consider the case c > 0. As shown above, because Q has only a finite number of zeros in the domain C, there exists a constant R > 0 which is large enough such that Q has no zero in {z ∈ C : |z| ≥ R} . On the other hand, Q is continuous at 0 with Q(0) > 0, so we can find a small constant ε > 0 such that Q(z) = 0 in {z ∈ C : |z| ≤ ε}. We define an oriented contour γ formed by four segments: Clearly, if Q has no purely imaginary zero, then all zeros in the closed right hand side of the complex plane {s = r(cos φ + i sin φ) ∈ C : r ≥ 0, φ ∈ (−π, π]} of Q (if they exist) must lie inside the contour γ. Based on Cauchy's argument principle in complex analysis, we see that n(Q(C), 0) = Z − P, where n(Q(C), 0) is the number of encirclements in the positive direction (counter-clockwise) around the origin of the the Nyquist plot Q(γ), Z and P are the number of zeros and number of poles of Q inside the contour γ in the s-plane, respectively. Due to the fact that Q is analytic inside γ, we have P = 0 and thus n(Q(C), 0) = Z. This implies that if Q has no purely imaginary zero, then all roots of the equation Q(s) = 0 lie in the open left-half complex plane if and only if . It is easy to see that n(Q(C), 0) = 0 if Q(iω) > 0 for any ω > 0 that satisfies Q(iω) = 0. Consider ω > 0 and put If there exists some ω > 0 such that h 2 (ω) = 0, then Thus, the variable ω > 0 satisfies the system with It is then clear from our assumptions on α 1 and α 2 that q 1 , q 2 > 0. Based on the analysis above, we obtain some sufficient conditions that ensure that the function Q(·) has no zero lying in the closed right half of the complex plane.
Then, all zeros of Q lie in the open left-half of the complex plane.
Then, all zeros of Q lie in the open left-half complex plane.
4 Estimates for the functions R λ and S β This section is devoted to the study of some important estimates of the functions R λ and S β on (0, ∞). We recall their definitions from (9), viz.
Lemma 4.1. Let α 1 , α 2 ∈ (0, 1] and denote ν = min{α 1 , α 2 }. Assume there are no zeros of the characteristic function Q in the closed right-half complex plane. Then, the following estimates hold for λ ∈ {0, α 1 , α 2 } and β ∈ {α 1 , α 2 , l(α)}: Moreover, The proof of the lemma is quite lengthy and technical. Therefore, in order not to distract the reader and to make it easier to focus on the main results, we provide the proof in the Appendix at the end of the paper.
Our first application of Lemma 4.1 deals with estimates for the convolution of S β and a continuous function.
(i) If g is bounded thenF β g and F β g are also bounded.
Proof. We denote ν = min{α 1 , α 2 }. Since, by definition, |F β g (t)| ≤F β g (t), the claims for F β g immediately follow from those forF β g , and therefore it suffices to explicitly prove the latter.
Statement (i) is merely the special case of η = 0 of part (iii).
To prove (ii), we note thatF β g (t) ≥ 0 by definition. Therefore, it is sufficient to show that for every ε > 0 there exists a constant T = T (ε) such that Since this is trivially fulfilled if g(t) = 0 for all t, we from now on assume that g(t) = 0 for some t, and hence g ∞ > 0.
Our first observation is then that, from (28) and (29), we know that there exists some constant C > 0 such that Given an arbitrary ε > 0, due to our assumption on g we may then find someT > 0 such that |g(t)| < νε/(3C) for all t >T . Using these valuesT and C, we then define For t > T ≥T + 1, we can then writē Our goal now is to show that, under these assumptions, F j (t) ≤ ε/3 for j = 1, 2, 3, which implies (33) and thus suffices to prove part (ii) of the Theorem. In this context, we see that because here t − s ≥ t −T > T −T ≥ 1, so that we may use the first of the bounds given in (34). In the penultimate step, we have bounded the integral by the product of the length of the integration interval and the maximum of the integrand, and in the last step, we have used the fact that t > T and the definition of T . Furthermore, because here s ≥T , so that |g(s)| ≤ νε/(3C), and t − s ≥ 1, so we may once again use the first bound of (34). Finally, where now t and s are such that we may invoke the second bound of (34) but, as in the previous step, s ≥T , so that once again |g(s)| ≤ νε/(3C). This completes the proof of part (ii) of the Theorem.
For the proof of (iii), we note that (34) is valid in this case too. Moreover, since we are interested in the asymptotic behaviour ofF β g (t) for large t, we may assume without loss of generality that t ≥ 2. Then we writē and we need to show that F j (t) = O(t −µ ) for j = 4, 5, 6, 7.
In this connection, we first note that, by assumption, with some C > 0, so that the upper branch of (34) implies as well as For the remaining part we need to invoke the second branch of (34) in combination with (37) to derive thus completing the proof.
As an immediate application of Theorem 4.2(iii), we can conclude that Moreover, assuming ν < 1 and setting we can obtain (using Lemma 4.1 and the classical relation between the incomplete Beta function and the hypergeometric Function 2 F 1 , cf. [1, eq. (6.6.8)]) the following bounds: • If t ≥ 2 then we havē with some C > 0.
with some C > 0.

Asymptotic behaviour of solutions to non-commensurate fractional planar systems
In this section we will study the asymptotic behaviour of solutions to fractional-order linear planar systems and the Mittag-Leffler stability of an equilibrium point to fractional nonlinear planar systems.

Asymptotic behaviour of solutions to fractional linear planar systems
Consider the non-homogeneous linear two-component incommensurate fractional-order system where α = (α 1 , α 2 ) ∈ (0, 1] × (0, 1] is a multi index, A = (a ij ) ∈ R 2×2 is a square real matrix and f = (f 1 , f 2 ) is a continuous vector valued function which is exponentially bounded on [0, ∞). (i) If f is bounded, then for any x 0 ∈ R 2 the solution to (42) is also bounded.
(iii) If f (t) = O(t −η ) as t → ∞ with some η > 0 then every solution x of (42a) behaves as Proof. The proof is straightforward by combining Lemma 2.3, Lemma 4.1 and Theorem 4.2.

Mittag-Leffler stability of fractional nonlinear planar systems
We now look at a different class of systems. Specifically, we now allow the differential equations to contain nonlinearities, but we do require them to have the structure of an autonomous system, i.e., we consider a fractional nonlinear planar system of the form Our aim is to prove the following theorem.
As shown above, we see that Lemmas 3.3, 3.4, 3.5, 3.6 and 3.7 give sufficient conditions which ensure that the characteristic function Q has no zero in the closed right hand side of the complex plane. Thus, by combining these lemmas and Theorem 5.4, we obtain the result below.
Corollary 5.5. Let The statement of Theorem 5.4 is true if one of the following conditions is satisfied.

Numerical examples
To complete this paper, we now give some numerical examples to illustrate the main theoretical results. In all the examples below, we use the functions f 1 and f 2 with For all cases, we have calculated numerical solutions to verify the theoretical findings. These solutions have been computed with Garrappa's MATLAB implementation of the implicit trapezoidal method described in detail in [11]. This algorithm is known to have very favourable stability properties which makes it highly suitable for handling equations like ours over large intervals (which is required in this case to demonstrate the asymptotic behaviour). The step size has always been chosen as h = 1/200.
In this example, the characteristic function is Q(s) = s 5/6 + s 1/2 + 0.5. According to Lemma 3.4, all zeros of Q lie in the open left-half of the complex plane. Furthermore, the function f satisfies the assumption stated in Theorem 5.1. Hence, every solution to (54) tends to the origin as t → ∞ with the rate O(t −1/3 ). This property is illustrated in Figure 1. The left graph shows that the components x 1 (t) and x 2 (t) decay to zero; the right graph visualizes the fact that t 1/3 x j (t) tends to a nonzero constant for t → ∞ and j = 1, 2, thus demonstrating that the decay behaviour of x j (t) is indeed O(t −1/3 ).
Example 6.2. Consider the two-component incommensurate fractional-order nonlinear system Defintion 5.3 states that the boundedness of the solutions cannot be expected for all choices of the initial value any more (as had been the case in Example 6.1) but only for initial values sufficiently close to (0, 0). Indeed we can see this behaviour in Figure 2 for the initial value (0.1, −0.2), whereas Figure 3 shows that this behaviour is not present for initial values farther away from (0, 0) such as, e.g., the initial value (1, −1). In the latter case, the solutions still seem to be bounded, but the decay behaviour appears to be absent. If one moves the initial values even farther away from the equilibrium point, then one cannot even expect this boundedness any more.  Example 6.3. Consider the fractional linear system The characteristic function of the system is Q(s) = s 1.4 + s 0.6 + 2. By Lemma 3.5, all zeros of Q lie in the open left-half of the complex plane and the assumptions of Theorem 5.1 are satisfied. Hence, every solution to this system converges to the origin as t → ∞ with an O(t −0.6 ) convergence rate. As in Example 6.1, we can also reproduce this behaviour numerically. The corresponding graphs are plotted in Figure 4. Example 6.4. Consider the system Based on Lemma 3.5 and Theorem 5.4, we see that the trivial solution of (57) is Mittag-Leffler stable. As in Example 6.2, this is exhibited-together with the decay behaviour predicted by Theorem 5.4-in Figure 5. Example 6.5. Consider the inhomogeneous two-component incommensurate fractionalorder linear system The system (58) has the characteristic function Q(s) = s 0.7 + s 0.4 + s 0.3 + 3. From Lemma 3.6 (i) and Theorem 5.1, it follows that every solution of this system tends to the origin as t → ∞ as O(t −0.3 ). Once again, our numerical results, shown in Figure 6, support this statement. Figure 6: Solution to the differential equation (58) from Example 6.5 with initial conditions x 1 (0) = 1 and x 2 (0) = 2. The left graph shows the components x j (t) of the solutions themselves, the right graph shows the functions t 0.3 x j (t).
Example 6.6. Consider the two-component incommensurate fractional-order nonlinear system Its characteristic function is Q(s) = s 0.7 + 0.2s 0.4 + 0.1s 0.3 + 0.3. It follows from Lemma 3.6(ii) and Theorem 5.4 that the trivial solution is Mittag-Leffler stable. Once again, we can visualize this observarion on the basis of numerical results, cf. Figure 7. Example 6.7. Consider the two-component incommensurate fractional-order linear system The system (60) has the characteristic function Q(s) = s 0.9 + 4s 0.5 − s 0.3 + 6. According to Lemma 3.7(i) and Theorem 5.1, its solution converges to the origin with a rate O(t −0.4 ). As above, the numerical data shown in Figure 8 confirms this theoretical observation.
(i) Because all zeros of Q (if they exist) lie on the left of the contour γ(R, π 2 + δ), using the same argument as in [21,Lemma 4.1], we obtain the representation Choose ε > 0 such that Q has no zero in the ball {s ∈ C : |s| ≤ ε} . From (63), we have where Λ t is the clockwise oriented contour bounding the domain Ω t := s ∈ C : ε t < |s| < R, | arg s| < π 2 + δ , see Figure 10. Notice that (s l(α)−λ−1 e st )/Q(s) is analytic on Ω t ∪ Λ t for all t ≥ 1. Thus, by applying Cauchy's theorem, we obtain I 1 (t) = 0 for all t ≥ 1.
Therefore, for each t ≥ 1, we see that R λ (t) = 1 2πi γ( ε Figure 10: The contours and sets used in the proof of Lemma 4.1: The radii of the green and magenta circular arcs are ε/t and R, respectively. The contour γ(R, π/2+δ) comprises the upper blue ray, the magenta circular arc, and the lower blue ray and is traversed from top to bottom; Ω t is the open set bounded by the magenta and green boundary lines (so these magenta and green lines together form the contour Λ t ). Λ 1 comprises the upper blue ray and the upper green line; Λ 2 denotes the union of the lower blue ray and the lower green line, and Λ 3 is the green circular arc.