1 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 science and engineering, we refer the interested reader to the collections [2, 3, 16, 20, 21] and the references therein.

To our knowledge, the first contribution in the qualitative study of fractional order autonomous linear systems was published by Matignon [15]. In that paper, using Laplace transforms 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. [22] 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. [18] 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 in [22] mentioned above.

In addition to the algorithmic approach as in [18], a number of analytic approaches have been used to investigate the zeros of characteristic polynomials of systems of fractional order systems. In [13], 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 [13] has been successfully extended in [28] 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 [25] 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 [19] have considered the stability and performance analysis of positive mixed fractional order linear systems with bounded delays. Tuan et al. [26] 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

$$\begin{aligned} ^CD^{\alpha }_{0^+}x(t)&=Ax(t)+f(t,x(t)), \quad t>0, \end{aligned}$$
(1.1)
$$\begin{aligned} x(0)&= x^0 \in {\mathbb {R}}^2, \end{aligned}$$
(1.2)

where \(\alpha = (\alpha _1,\alpha _2) \in (0,1]^2\) is a multi-index, \(A\in {\mathbb {R}}^{2\times 2}\) is a square matrix and \(f :[0, \infty ) \times {\mathbb {R}}^2 \rightarrow {\mathbb {R}}^2\) is a 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.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 behaviour 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 \({\mathbb {R}}^2\), we define the norm \(\Vert \cdot \Vert \) by \(\Vert x\Vert := \max \{|x_1|,|x_2|\}\) for every \(x\in {\mathbb {R}}^2\). For any \(r>0\), the closed ball of radius r centered at the origin 0 in \({\mathbb {R}}^2\) is given by \(B(0,r):=\{x\in {\mathbb {R}}^2:\Vert x\Vert \le r\}\). The space of all continuous functions \(\xi :[0,\infty )\rightarrow {\mathbb {R}}^2\) is denoted by \(C([0,\infty );{\mathbb {R}}^2)\). For any \(\xi \in C([0,\infty );{\mathbb {R}}^2)\), let \(\Vert \xi \Vert _\infty :=\sup _{t\ge 0}\Vert \xi (t)\Vert \). Then, we use the notation \( C_\infty ([0,\infty );{\mathbb {R}}^2) := \{\xi \in C([0,\infty );{\mathbb {R}}^2): \Vert \xi \Vert _\infty <\infty \}\) to designate the subspace of \(C([0,\infty );{\mathbb {R}}^2)\) that comprises the bounded continuous functions on \([0, \infty )\).

For \(\alpha \in (0,1]\) and \(J = [0, T]\) or \(J = [0, \infty )\), we define the Riemann-Liouville fractional integral of a function \(f :J \rightarrow {\mathbb {R}}\) as

$$\begin{aligned} I^\alpha _{0^+}f(t) := \frac{1}{\varGamma (\alpha )}\int _{0}^{t}(t-s)^{\alpha -1}f(s)ds,\;t\in J, \end{aligned}$$

and the Caputo fractional derivative of the order \(\alpha \in (0,1]\) of a function \(f : J \rightarrow {\mathbb {R}}\) as

$$\begin{aligned} ^C D^\alpha _{0^+}f(t) := \frac{d}{dt}I^{1-\alpha }_{0^+}(f(t)-f(0)), \;t\in J \setminus \{ 0 \}, \end{aligned}$$

where \(\varGamma (\cdot )\) is the Gamma function and \(\frac{d}{dt}\) is the usual derivative. Letting \(\alpha = (\alpha _1,\alpha _2) \in (0,1]\times (0,1]\) be a multi-index and \(f = (f_1, f_2)\) with \(f_i : J \rightarrow {\mathbb {R}}\), \(i=1,2\), be a vector valued function, we write

$$\begin{aligned} ^C D^\alpha _{0^+} f(t) := \left( ^C D^{\alpha _1}_{0^+}f_1(t),^C D^{\alpha _2}_{0^+}f_2(t) \right) . \end{aligned}$$

See, e.g., [9,  Chapter 3] and [27] for more details on the Caputo fractional derivative.

2 Preliminaries

2.1 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

$$\begin{aligned} ^C D^{\alpha }_{0^+}x(t)&= f(t,x(t)), \quad t>0, \end{aligned}$$
(2.1a)
$$\begin{aligned} x(0)&= x^0 \in {\mathbb {R}}^2, \end{aligned}$$
(2.1b)

where \(\alpha = (\alpha _1,\alpha _2) \in (0,1]^2\) is a multi-index and \(f :[0, \infty ) \times {\mathbb {R}}^2 \rightarrow {\mathbb {R}}^2\) is a continuous function. To make the exposition of this paper self-contained, we recall two important results about the solutions to this two-dimensional initial value problem.

Theorem 1

(Existence and uniqueness of global solutions) Suppose that the function \(f:[0, \infty ) \times {\mathbb {R}}^2 \rightarrow {\mathbb {R}}^2\) is continuous and that, for some constant \(L > 0\), it satisfies the Lipschitz condition

$$\begin{aligned} \Vert f(t,x) - f(t,{\hat{x}})\Vert \le L \Vert x - {\hat{x}}\Vert ,\; \forall t \in [0.\infty ),\; x ,{\hat{x}} \in {\mathbb {R}}^2 \end{aligned}$$

with respect to its second variable. Then, for any initial value \(x^0 \in {\mathbb {R}}^2\), the two-component incommensurate fractional-order system (2.1) has a unique global solution \(\varphi (\cdot ,x^0)\) on the interval \([0,\infty )\).

Proof

See [24, Theorem 2.2 and Remark 2.3]. \(\square \)

Theorem 2

(Exponential boundedness of global solutions) Suppose that the function f satisfies the assumptions of Theorem 1. Moreover, let there exist a constant \(\gamma > 0\) such that

$$\begin{aligned} \sup _{t\ge 0} e^{-\gamma t} \int _0^t(t-s)^{\alpha _i -1}\Vert f(s,0)\Vert ds < \infty . \end{aligned}$$

Then, for any initial value \(x^0 \in {\mathbb {R}}^2,\) the two-component incommensurate fractional-order system (2.1) has a unique global solution \(\varphi (\cdot ,x^0) \in C\left( [0,\infty ),{\mathbb {R}}^2 \right) \) and

$$\begin{aligned} \Vert \varphi (t,x^0)\Vert \le Me^{\gamma t},\; \forall t \ge 0, \end{aligned}$$

where M is some positive constant which depends on \(x^0.\)

Proof

See [24, Theorem 2.4]. \(\square \)

2.2 The variation of constants formula for the solutions

Consider the non-homogeneous two-component incommensurate fractional-order linear system

$$\begin{aligned} ^CD^{\alpha }_{0^+}x(t)=Ax(t) + f(t) , \quad t>0 \end{aligned}$$
(2.2a)

with initial condition

$$\begin{aligned} x(0) = x^0 \in {\mathbb {R}}^2, \end{aligned}$$
(2.2b)

where \(\alpha = (\alpha _1, \alpha _2) \in (0,1]^2\), \(A = \left( a_{ij} \right) \in {\mathbb {R}}^{2 \times 2}\) is a square real matrix and \(f=(f_1,f_2)^{\mathrm{T}} :[0,\infty ) \rightarrow {\mathbb {R}}^2\) is a continuous function such that

$$\begin{aligned} \Vert f(t)\Vert \le Me^{\gamma t},\;\forall t \ge 0 \end{aligned}$$
(2.3)

for some \(M > 0\) and some \(\gamma > 0\). Then, we have

$$\begin{aligned} \int _0^t(t-s)^{\alpha _i -1}|f_i(s)|ds&\le M \int _0^t(t-s)^{\alpha _i -1}e^{\gamma s}ds\\&= -\frac{M e^{\gamma t}}{\gamma ^{\alpha _i}} \int _0^t\left( \gamma (t-s) \right) ^{\alpha _i -1}e^{-\gamma (t-s)}d(\gamma (t-s))\\&= \frac{M e^{\gamma t}}{\gamma ^{\alpha _i}}\int _0^{\gamma t} \tau ^{\alpha _i -1}e^{-\tau }d\tau \\&\le \frac{M \varGamma (\alpha _i)}{\gamma ^{\alpha _i}} e^{\gamma t}. \end{aligned}$$

Due to Theorems 1 and 2, for any initial condition \(x^0 \in {\mathbb {R}}^2\), the system (2.2) has a unique exponentially bounded solution in \(C\left( [0, \infty ), {\mathbb {R}}^2 \right) \). Taking the Laplace transform on both sides of the system (2.2), we obtain the algebraic system

$$\begin{aligned} {\left\{ \begin{array}{ll} (s^{\alpha _1} -a_{11})X_1(s) - a_{12}X_2(s) &{}= s^{\alpha _1 -1 }x_1^0 +F_1(s) \\ -a_{21}X_1(s) +(s^{\alpha _2} -a_{22})X_2(s)&{} = s^{\alpha _2 -1 }x_2^0 +F_2(s) \end{array}\right. } , \end{aligned}$$
(2.4)

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

$$\begin{aligned} X_1(s)&= \frac{ x_1^0(s^{\alpha _1 + \alpha _2-1} \!-\! a_{22}s^{\alpha _1 -1})+ x^0_2a_{12}s^{\alpha _2-1} + F_1(s)(s^{\alpha _2} \!- \!a_{22}) + a_{12}F_2(s) }{Q(s)}\nonumber \\&= \frac{s^{\alpha _1 + \alpha _2} - a_{22}s^{\alpha _1}}{sQ(s)}x^0_1 + \frac{s^{\alpha _2}}{sQ(s)}x^0_2 + \frac{s^{\alpha _2}-a_{22}}{Q(s)}F_1(s ) + \frac{a_{12}F_2(s)}{Q(s)}, \end{aligned}$$
(2.5)

and

$$\begin{aligned} X_2(s)&= \frac{x_2^0(s^{\alpha _1 + \alpha _2-1} \!-\! a_{11}s^{\alpha _2 -1})+ x^0_1a_{21}s^{\alpha _1-1} + a_{12}F_1(s) + F_2(s)(s^{\alpha _1}\!-\!a_{11})}{Q(s)}\nonumber \\&= \frac{s^{\alpha _1 + \alpha _2} - a_{11}s^{\alpha _2 }}{sQ(s)}x_2^0 + \frac{a_{21}s^{\alpha _1}}{sQ(s)}x_1^0 + \frac{a_{21}F_1(s)}{Q(s)} + \frac{s^{\alpha _1}-a_{11}}{Q(s)}F_2(s), \end{aligned}$$
(2.6)

where \(Q(s) := s^{\alpha _1 + \alpha _2} - a_{11}s^{\alpha _2}-a_{22}s^{\alpha _1} +\det A. \) Put

$$\begin{aligned} {\mathcal {R}}^\lambda (t)&= {\mathcal {L}}^{-1}\left\{ \frac{s^{l(\alpha )-\lambda }}{sQ(s)} \right\} (t) , \quad \lambda \in \left\{ 0,\alpha _1, \alpha _2 \right\} , \end{aligned}$$
(2.7a)
$$\begin{aligned} {\mathcal {S}}^\beta (t)&= {\mathcal {L}}^{-1}\left\{ \frac{s^{l(\alpha )-\beta }}{Q(s)} \right\} (t) , \quad \beta \in \left\{ \alpha _1, \alpha _2 ,l(\alpha )\right\} \end{aligned}$$
(2.7b)

with \(l(\alpha ): = \alpha _1 + \alpha _2.\) Then, with each \(i \in \left\{ 1,2 \right\} \), we obtain

$$\begin{aligned} {\mathcal {L}}^{-1}\left\{ \frac{s^{l(\alpha )-\beta }}{Q(s)}F_i(s) \right\} (t)&= {\mathcal {L}}^{-1}\left\{ {\mathcal {L}}\left\{ {\mathcal {S}}^\beta \right\} (s) {\mathcal {L}}\left\{ f_i \right\} (s)\right\} (t)\\&= {\mathcal {L}}^{-1} \left\{ {\mathcal {L}}\left\{ {\mathcal {S}}^\beta *f_i \right\} (s) \right\} (t) \\&={\mathcal {S}}^\beta *f_i(t), \beta \in \left\{ \alpha _1, \alpha _2, l(\alpha ) \right\} , \end{aligned}$$

where “\(*\)” is the Laplace convolution operator.

From the arguments above, we can derive our first new result that shows a precise analytic representation of the unique solution to the initial value problem (2.2).

Lemma 1

On the interval \([0, \infty )\), the non-homogeneous linear two-component incommensurate fractional-order system (2.2) has the unique solution

$$\begin{aligned} \varphi (\cdot ,x^0) = \begin{pmatrix} \varphi _1(\cdot ,x^0)\\ \varphi _2(\cdot ,x^0) \end{pmatrix} \end{aligned}$$

with

$$\begin{aligned} \varphi _1(t,x^0)&= \left( {\mathcal {R}}^0(t)-a_{22}{\mathcal {R}}^{\alpha _2}(t) \right) x_1^0 + a_{12}{\mathcal {R}}^{\alpha _1}(t)x^0_2 \nonumber \\&+ \left( ( {\mathcal {S}}^{\alpha _1} - a_{22}{\mathcal {S}}^{l(\alpha )} )*f_1 \right) (t) + a_{12}\left( {\mathcal {S}}^{l(\alpha )}*f_2 \right) (t), \end{aligned}$$
(2.8)
$$\begin{aligned} \varphi _2(t,x^0)&= a_{21}{\mathcal {R}}^{\alpha _2}(t)x^0_1 + \left( {\mathcal {R}}^0(t) - a_{11}{\mathcal {R}}^{\alpha _1}(t) \right) x_2^0 \nonumber \\&+a_{21}\left( {\mathcal {S}}^{l(\alpha )}*f_1 \right) (t) + \left( ( {\mathcal {S}}^{\alpha _2}- a_{11}{\mathcal {S}}^{l(\alpha )})*f_2 \right) (t). \end{aligned}$$
(2.9)

3 Some properties of the characteristic function

As mentioned in the introduction, in this paper we only focus on incommensurate systems, i.e. on systems of the form (1.1) with \(\alpha _1 \ne \alpha _2\), because the case \(\alpha _1 = \alpha _2\) has already been discussed in detail elsewhere. Thus, without loss of generality, we assume \(0<\alpha _1 < \alpha _2\le 1\). It will turn out, see Section 5, that such systems possess characteristic functions of the form

$$\begin{aligned} Q(s) = s^{\alpha _1 + \alpha _2} - as^{\alpha _2}- bs^{\alpha _1} + c \end{aligned}$$
(3.1)

with certain \(a, b, c \in {\mathbb {R}}\) whose properties play an essential role in the analysis of the asymptotic behaviour of the solutions. In fact, we shall show in Section 5 that the crucial point is the location of the zeros of functions of this type: One needs to know whether or not all the function’s zeros are in the open left half of the complex plane. In this context, one result has already been shown by Brandibur and Kaslik [5]:

Lemma 2

Assume that \(a, b \le 0,\) and \(c > 0\). Then, all zeros of Q are in the open left-half complex plane regardless of \(\alpha _1\) and \(\alpha _2\).

Proof

See [5,  Proposition 1(b)]. \(\square \)

The conditions of Lemma 2 are quite restrictive though, and therefore this statement can only be applied to a very limited class of systems. In the remainder of this section we shall therefore now develop some new alternative sufficient criteria for asserting a “good” distribution of the zeros of the characteristic function that can cover many further cases. Our first auxiliary statement in this context reads as follows.

Lemma 3

Let \(0<\alpha _1 < \alpha _2\le 1\) and \(a,b,c \in {\mathbb {R}}\). Then, the following statements hold for the function Q defined in (3.1).

  1. (i)

    If \( c < 0\) then Q has at least one positive real zero.

  2. (ii)

    If \(s \in {\mathbb {C}}\) is a zero of Q then its complex conjugate is also a zero of Q.

  3. (iii)

    Let \(0< \omega < \pi \). Then, Q has only a finite number of zeros in the set \({\mathcal {C}} = \{z \in {\mathbb {C}} : | \arg {(z)}| \le \omega \}\).

  4. (iv)

    If \(c > 0\), then \(s = i\omega \) with \(\omega > 0\) is a zero of Q if and only if

    $$\begin{aligned}&{\left\{ \begin{array}{ll} &{} a= \rho _2 \omega ^{\alpha _1} -c\rho _1 \omega ^{-\alpha _2}, \\ &{} b = c\rho _2 \omega ^{-\alpha _1} -\rho _1 \omega ^{\alpha _2}, \end{array}\right. } \end{aligned}$$
    (3.2)

    where

    $$\begin{aligned} \rho _1 = \frac{\sin \frac{\alpha _1 \pi }{2}}{\sin \frac{(\alpha _2 - \alpha _1)\pi }{2}} , \rho _2 = \frac{\sin \frac{\alpha _2 \pi }{2}}{\sin \frac{(\alpha _2 - \alpha _1)\pi }{2}}. \end{aligned}$$
    (3.3)

Proof

(i) and (ii) are obvious.

(iii) First, we assume that \(c\ne 0\). Then \(Q(0) \ne 0\). Due to the continuity of Q at 0, we can find \(\varepsilon \) which is small enough such that Q has no zero in \(\left\{ z\in {\mathbb {C}}: |z| < \varepsilon \right\} .\) Moreover, because \(|Q(s)| \ge |s|^{\alpha _1 + \alpha _2} - |a| \cdot |s|^{\alpha _2}-|b| \cdot |s|^{\alpha _1} -|c|,\) we have that \(\lim _{|s| \rightarrow \infty } |Q(s)| = \infty \) uniformly for all \(\arg s\). This implies that there is a positive real number R such that Q has no zero in the domain \(\left\{ z\in {\mathbb {C}}: |z| > R \right\} .\) Hence, all zeros of Q in \(\{z\in {\mathbb {C}}: | \arg {(z)}| \le \omega \}\) (if they exist) belong to the set \(\varOmega := \left\{ z\in {\mathbb {C}}: \varepsilon \le |z| \le R, |\arg {(z)}|\le \omega \right\} .\) Notice that \(\varOmega \) is a compact set and Q is analytic on this domain. If now Q has infinitely many zeros in \(\varOmega \) then, because of the compactness of \(\varOmega \), 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 \(\varOmega \). This shows that Q has only a finite number of zeros in the domain \({\mathcal {C}}\) if \(c \ne 0\).

To deal with the case \(c=0\), we write

$$\begin{aligned} Q(s)&= s^{\alpha _1}\left( s^{\alpha _2} - as^{\alpha _2-\alpha _1}-b \right) =s^{\alpha _1}P(s) \end{aligned}$$

where \(P(s) = s^{\alpha _2} - as^{\alpha _2-\alpha _1}-b.\) By repeating the above arguments for P, the proof is complete.

(iv) See [5,  Proposition 1, Part 3b]. \(\square \)

Corollary 1

Assume that \(a,b,c > 0\) and that one of conditions

  1. (i)

    \(c(\rho ^2_2 - \rho _1^2)< ab < c(\rho _2^2 + \rho _1^2)\),

  2. (ii)

    \(ab \le c\left( \rho _2 -\rho _1 \right) ^2\),

is satisfied where \(\rho _1,\rho _2\) are defined in (3.3). Then, the function Q defined in (3.1) has no purely imaginary zero.

Proof

Consider the system (3.2). Since \(\rho _2 \ne 0\), this system is equivalent to

$$\begin{aligned} \left\{ \begin{array}{rl} \omega ^{\alpha _1} &{} = \displaystyle \frac{a}{\rho _2} + c\frac{\rho _1}{\rho _2}\omega ^{-\alpha _2}, \\ \rho _1 a \omega ^{\alpha _2} + \rho _2 b \omega ^{\alpha _1} &{}= c(\rho _2^2 - \rho _1^2). \end{array} \right. \end{aligned}$$
(3.4)

Thus, we obtain

$$\begin{aligned} a\rho _1\omega ^{\alpha _2} + bc\rho _1\omega ^{-\alpha _2}+ab - c(\rho _2^2-\rho _1^2) = 0. \end{aligned}$$
(3.5)

Setting \(X=\omega ^{\alpha _2}\), equation (3.5) takes the form

$$\begin{aligned} a \rho _1 X^2 + \left[ ab - c (\rho _2^2-\rho _1^2) \right] X + b c \rho _1 = 0. \end{aligned}$$
(3.6)

The discriminant of the quadratic equation (3.6) is

$$\begin{aligned} \varDelta&= \left( ab - c(\rho _2^2-\rho _1^2) \right) ^2 - 4abc\rho _1^2 \\&=a^2b^2 + c^2(\rho _2^2-\rho _1^2)^2-2abc(\rho _1^2 + \rho _2^2) \\&= \left( ab -c(\rho _1^2 + \rho _2^2)\right) ^2 - 4c^2\rho _1^2\rho _2^2 \\&=\left( ab - c(\rho _1+\rho _2)^2 \right) \left( ab - c(\rho _2-\rho _1)^2 \right) . \end{aligned}$$

(i) Clearly, if \(c(\rho ^2_2 - \rho _1^2)< ab < c(\rho _2^2 + \rho _1^2),\) then \(\varDelta < 0.\) Hence, the quadratic equation (3.6) has no real roots. This implies that the system (3.2) has no root \(\omega >0\). This together with Lemma 3(ii) and 3(iv) shows that Q has no purely imaginary zero.

(ii) If \(ab \le c\left( \rho _2 -\rho _1 \right) ^2\), the quadratic equation (3.6) has two (not necessarily distinct) real roots. Because \( 0< \rho _1 < \rho _2,\) we have \(( \rho _2 -\rho _1 )^2 < \rho _2^2-\rho _1^2.\) This implies that \(ab - c(\rho _2^2-\rho _1^2) < 0\). Moreover, \(a,b,c, \rho _1 > 0,\) thus the two roots of the quadratic equation (3.6) are negative. Hence, in view of the relation \(X = \omega ^{\alpha _2}\) with \(0 < \alpha _2 \le 1\) between the solution X of (3.6) and the solution \(\omega \) of (3.2), the system (3.2) has no root \(\omega >0\). Using Lemmas 3(ii) and 3(iv), we see that Q has no purely imaginary zero. \(\square \)

Recall that if \(c \le 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 \({\mathbb {C}}\), there exists a constant \(R > 0\) which is large enough such that Q has no zero in \(\left\{ z\in {\mathbb {C}}: |z|\ge R\right\} .\) On the other hand, Q is continuous at 0 with \(Q(0) > 0\), so we can find a small constant \(\varepsilon > 0\) such that \(Q(z)\ne 0\) in \(\left\{ z\in {\mathbb {C}}: |z| \le \varepsilon \right\} \). We define an oriented contour \(\gamma \) formed by four segments:

$$\begin{aligned} \gamma _1&:=\left\{ s= i\omega : \varepsilon \le \omega \le R \right\} ;&\quad \gamma _2&:=\left\{ s=R e^{i\varphi }: -\frac{\pi }{2}\le \varphi \le \frac{\pi }{2} \right\} ; \\ \gamma _3&:=\left\{ s=\varepsilon e^{i\varphi }: -\frac{\pi }{2}\le \varphi \le \frac{\pi }{2} \right\} ;&\quad \gamma _4&:=\left\{ s= i\omega : -R \le \omega \le -\varepsilon \right\} . \end{aligned}$$

Clearly, if Q has no purely imaginary zero, then all zeros in the closed right hand side of the complex plane \(\{s=r(\cos \phi +i\sin \phi )\in {\mathbb {C}}: r\ge 0, \phi \in (-\pi ,\pi ]\}\) of Q (if they exist) must lie inside the contour \(\gamma \). 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(\gamma )\), Z and P are the number of zeros and number of poles of Q inside the contour \(\gamma \) in the s-plane, respectively. Due to the fact that Q is analytic inside \(\gamma \), 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 \(n(Q(C), 0) = 0.\) Notice that \(Q(0)>0\), \(\lim _{|s|\rightarrow \infty }|Q(s)|=\infty \), \(\Re Q(i\omega ) = \Re Q(-i\omega )\) and \(\Im Q(i\omega ) = -\Im Q(-i\omega )\). It is easy to see that \(n(Q(C), 0) = 0\) if \(\Re Q(i\omega ) > 0\) for any \(\omega >0\) that satisfies \(\Im Q(i\omega ) = 0\). Consider \(\omega >0\) and put

$$\begin{aligned} h_1(\omega )&:= \Re (Q(i\omega )) \nonumber \\&= \omega ^{\alpha _1 + \alpha _2}\cos \frac{(\alpha _1 + \alpha _2)\pi }{2} - a\omega ^{ \alpha _2}\cos \frac{ \alpha _2\pi }{2} - b\omega ^{\alpha _1}\cos \frac{\alpha _1 \pi }{2}+c, \end{aligned}$$
(3.7)
$$\begin{aligned} h_2(\omega )&:= \Im (Q(i\omega )) \nonumber \\&= \omega ^{\alpha _1 + \alpha _2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} - a\omega ^{ \alpha _2}\sin \frac{\alpha _2\pi }{2} - b\omega ^{\alpha _1} \sin \frac{\alpha _1 \pi }{2}. \end{aligned}$$
(3.8)

If there exists some \(\omega > 0\) such that \(h_2(\omega ) = 0,\) then

$$\begin{aligned}&\sin \frac{(\alpha _1 + \alpha _2)\pi }{2}h_1(\omega ) \\&\quad = \sin \frac{(\alpha _1 + \alpha _2)\pi }{2}h_1(\omega ) -\cos \frac{(\alpha _1 + \alpha _2)\pi }{2}h_2(\omega ) \\&\quad = a\omega ^{ \alpha _2}\left( \sin \frac{ \alpha _2\pi }{2}\cos \frac{(\alpha _1 + \alpha _2)\pi }{2} - \cos \frac{ \alpha _2\pi }{2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} \right) \\&\qquad + b\omega ^{\alpha _1}\left( \sin \frac{ \alpha _1\pi }{2}\cos \frac{(\alpha _1 + \alpha _2)\pi }{2} - \cos \frac{ \alpha _1\pi }{2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} \right) \\&\qquad + c\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} \\&\quad = c\sin \frac{(\alpha _1 + \alpha _2)\pi }{2}-a\omega ^{ \alpha _2}\sin \frac{ \alpha _1\pi }{2} - b\omega ^{\alpha _1}\sin \frac{ \alpha _2\pi }{2}. \end{aligned}$$

Thus, the variable \(\omega >0\) satisfies the system

$$\begin{aligned} \left\{ \begin{array}{rl} h_2(\omega ) &{}= 0 \\ h_1(\omega ) &{}= c - a\omega ^{\alpha _2}q_1 - b\omega ^{\alpha _1}q_2 \end{array} \right. \end{aligned}$$
(3.9)

with

$$\begin{aligned} q_1 = \frac{\sin \frac{\alpha _1 \pi }{2}}{\sin \frac{(\alpha _2 + \alpha _1)\pi }{2}}, q_2 = \frac{\sin \frac{\alpha _2 \pi }{2}}{\sin \frac{(\alpha _2 + \alpha _1)\pi }{2}}. \end{aligned}$$
(3.10)

It is then clear from our assumptions on \(\alpha _1\) and \(\alpha _2\) that \(q_1, q_2 > 0\). Based on the analysis above, we obtain some sufficient conditions that ensure that the function \(Q(\cdot )\) has no zero lying in the closed right half of the complex plane.

Lemma 4

Let \(0< \alpha _1 < \alpha _2 \le 1\). Assume that \(a = 0\), \(b > 0\) and

$$\begin{aligned} c > \left( b q_1 \right) ^{\alpha _1/\alpha _2} b q_2. \end{aligned}$$
(3.11)

Then, all zeros of Q lie in the open left-half of the complex plane.

Proof

Because \(a = 0\) and \(b > 0\), we see that \(h_2(\omega _0) = 0\) if and only if \(\omega _0 = \left( bq_1 \right) ^{1/\alpha _2}\). From (3.9), we have \(h_1(\omega _0) = c - \left( bq_1 \right) ^{\alpha _1/\alpha _2} b q_2\). By the assumption (3.11), we obtain \(h_1(\omega _0) > 0\). This implies that all zeros of Q lie in the open left-half of the complex plane. \(\square \)

Lemma 5

Let \(0< \alpha _1 < \alpha _2 \le 1.\) Assume that \(b = 0\), \(a > 0\) and

$$\begin{aligned} c > \left( a q_2 \right) ^{\alpha _2/\alpha _1} a q_1. \end{aligned}$$
(3.12)

Then, all zeros of Q lie in the open left-half complex plane.

Proof

Since \(b = 0\) and \(a > 0\), it is easy to show that \(h_2(\omega _0) = 0\) if and only if \(\omega _0 = \left( aq_2 \right) ^{1/\alpha _1}\). From (3.9), it follows that \(h_1(\omega _0) = c - \left( aq_2 \right) ^{\alpha _2/\alpha _1}aq_1\). By the assumption (3.12), we see that \(h_1(\omega _0) > 0\) which implies that all zeros of Q lie in the open left-half of the complex plane. \(\square \)

Lemma 6

Let \(0< \alpha _1 < \alpha _2 \le 1.\) Assume that \(a, b, c > 0 \). Then, all zeros of Q are in the open left-half of the complex plane if one of the following conditions holds:

  1. (i)

    \( aq_2+bq_1 > 1\) and \(aq_2\left( (a+b)q_2 \right) ^{\alpha _2/\alpha _1} + b(a+b)q_2^2 \le c.\)

  2. (ii)

    \( aq_2+bq_1 \le 1\) and \( aq_1 + bq_2 < c.\)

Proof

We have

$$\begin{aligned} h_2'(\omega )&= (\alpha _1 + \alpha _2)\omega ^{\alpha _1 + \alpha _2-1}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} \\&\qquad \qquad - a\alpha _2\omega ^{ \alpha _2-1}\sin \frac{ \alpha _2\pi }{2} - b \alpha _1\omega ^{\alpha _1-1 }\sin \frac{\alpha _1 \pi }{2} \nonumber \\&= \omega ^{\alpha _1-1} \Bigg ( (\alpha _1 + \alpha _2)\omega ^{ \alpha _2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} \nonumber \\&- a\alpha _2\omega ^{ \alpha _2-\alpha _1}\sin \frac{ \alpha _2\pi }{2} - b\alpha _1 \sin \frac{\alpha _1 \pi }{2} \Bigg )\nonumber \\&= \omega ^{\alpha _1-1}g_2(\omega ) \nonumber \end{aligned}$$
(3.13)

where

$$\begin{aligned} g_2(\omega ) : = (\alpha _1 + \alpha _2)\omega ^{ \alpha _2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} - a \alpha _2\omega ^{ \alpha _2-\alpha _1}\sin \frac{ \alpha _2\pi }{2} - b \alpha _1 \sin \frac{\alpha _1 \pi }{2}. \nonumber \\ \end{aligned}$$
(3.14)

Notice that

$$\begin{aligned} g_2'(\omega )&= \alpha _2(\alpha _1 + \alpha _2) \omega ^{\alpha _2 - 1}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} - (\alpha _2 - \alpha _1)\alpha _2 a\omega ^{\alpha _2 - \alpha _1 -1}\sin \frac{ \alpha _2\pi }{2}. \end{aligned}$$

It is not difficult to check that \(g_2'(\omega ) < 0\) in \((0,\omega _1)\) and \(g_2'(\omega ) > 0\) in \((\omega _1,\infty ),\) where \(\omega _1 = \left( \frac{\alpha _2-\alpha _1}{\alpha _1 + \alpha _2}aq_2 \right) ^{1/\alpha _1}.\) Due to the fact that \(g_2(0) = -\alpha _1 b \sin \frac{\alpha _1\pi }{2} < 0,\) and \(\lim _{\omega \rightarrow +\infty }g_2(\omega ) = +\infty ,\) the equation \(g_2(\omega ) = 0\) has a unique root \(\omega _2 \in (0, \infty ).\) Moreover \(g_2(\omega ) < 0\) in \((0,\omega _2)\) and \(g_2(\omega ) > 0\) in \((\omega _2, \infty ).\) Hence, \(h_2\) is decreasing in \((0,\omega _2)\) and increasing in \((\omega _2, \infty ).\) On the other hand, \(h_2(0) = 0\) and \(\lim _{\omega \rightarrow +\infty }h_2(\omega ) = +\infty \). This shows that the equation \(h_2(\omega ) = 0\) has a unique root \(\omega _3 \in (0,\infty )\) and then \(h_2(\omega ) < 0\) for all \(\omega \in (0,\omega _3)\) and \(h_2(\omega ) > 0\) for all \(\omega \in (\omega _3,\infty ).\)

(i) If \( aq_2+bq_1 > 1\), then \(h_2(1) < 0.\) This implies that \(\omega _3 > 1.\) Moreover, due to \(\alpha _1 < \alpha _2 \le 1\), we have

$$\begin{aligned} h_2(\omega ) > \omega ^{\alpha _1 + \alpha _2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} - (a+b)\omega ^{\alpha _2}\sin \frac{\alpha _2\pi }{2} \end{aligned}$$
(3.15)

for every \(\omega > 1\), and thus \(h_2\left( \left( (a+b)q_2 \right) ^{1/\alpha _1} \right) > 0.\) This implies that \(1< \omega _3 < \left( (a+b)q_2 \right) ^{1/\alpha _1}.\) Hence, if

$$\begin{aligned} c \ge aq_2\left( (a+b)q_2 \right) ^{\alpha _2/\alpha _1} + b(a+b)q_2^2, \end{aligned}$$

due to \(q_2> q_1 > 0,\) we obtain \(c > a\omega _3^{\alpha _2}q_1 + b\omega _3^{\alpha _1}q_2,\) which together with (3.9) leads to \(h_1(\omega _3) > 0.\) The proof of this part is complete.

(ii) If \( aq_2+bq_1 \le 1\), then \(h_2(1) \ge 0.\) This implies that \(0 < \omega _3 \le 1.\) Due to \( c > aq_1 + bq_2,\) we have \(c > a\omega _3^{\alpha _2}q_1 + b\omega _3^{\alpha _1}q_2\). This together with (3.9) shows that \(h_1(\omega _3) > 0.\) The proof is finished. \(\square \)

Lemma 7

Assume that \(a < 0\) and \(b, c > 0\). Then, all zeros of Q are in the open left-half complex plane if one of the following conditions holds:

  1. (i)

    \(a q_2+b q_1 > 1\) and \(\left( b q_1 \right) ^{\alpha _1/\alpha _2} b q_2 \le c.\)

  2. (ii)

    \(a q_2+b q_1 \le 1\) and \(b q_2 \le c\).

Proof

As shown in the proof of Lemma 6, we have \(h_2'(\omega ) = \omega ^{\alpha _1-1} g_2(\omega )\) where \(g_2\) is as in (3.14). Notice that

$$\begin{aligned} g_2'(\omega ) = \alpha _2(\alpha _1 + \alpha _2) \omega ^{\alpha _2 - 1}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} - (\alpha _2 - \alpha _1)\alpha _2 a\omega ^{\alpha _2 - \alpha _1 -1}\sin \frac{ \alpha _2\pi }{2} > 0 \end{aligned}$$

for \(\omega \in (0,\infty )\). Due to the facts that \(g_2(0) = -\alpha _1 b \sin \frac{\alpha _1\pi }{2} < 0\) and \(\lim _{\omega \rightarrow +\infty }g_2(\omega ) = +\infty ,\) the equation \(g_2(\omega ) = 0\) has a unique root \(\omega _1 \in (0, \infty )\). Moreover \(g_2(\omega ) < 0\) in \((0,\omega _1)\) and \(g_2(\omega ) > 0\) in \((\omega _1, \infty ).\) This shows that \(h_2\) is decreasing on \((0,\omega _1)\) and increasing on \((\omega _1, \infty )\). On the other hand, since \(h_2(0) = 0\) and \(\lim _{\omega \rightarrow +\infty }h_2(\omega ) = +\infty \), the equation \(h_2(\omega ) = 0\) has a unique root \(\omega _2 \in (0,\infty )\) and \(h_2(\omega ) < 0\) for all \(\omega \in (0,\omega _2)\) and \(h_2(\omega ) > 0\) for all \(\omega \in (\omega _2,\infty ).\)

(i) If \(aq_2+bq_1 > 1\) then \(h_2(1) < 0.\) Thus \(\omega _2 > 1\). Moreover, since \(a < 0,\) we have

$$\begin{aligned} h_2(\omega ) > \omega ^{\alpha _1 + \alpha _2}\sin \frac{(\alpha _1 + \alpha _2)\pi }{2} - b\omega ^{\alpha _1}\sin \frac{\alpha _1\pi }{2}, \end{aligned}$$

and thus \(h_2\left( (bq_1)^{1/{\alpha _2}} \right) > 0\). This implies that \(1< \omega _2 < (bq_1)^{1/{\alpha _2}}.\) From that if

$$\begin{aligned} \left( bq_1 \right) ^{\alpha _1/\alpha _2}bq_2 \le c, \end{aligned}$$

we obtain \(c > a\omega _2^{\alpha _2}q_1 + b\omega _2^{\alpha _1}q_2,\) which together with (3.9) leads to \(h_1(\omega _2) > 0.\)

(ii) If \( aq_2+bq_1 \le 1\), then \(h_2(1) \ge 0.\) Thus \(0 < \omega _2 \le 1.\) Due to \( c \ge bq_2\), we see that \(c > a \omega _2^{\alpha _2}q_1 + b\omega _2^{\alpha _1}q_2\). The proof is completed. \(\square \)

4 Estimates for the functions \({\mathcal {R}}^\lambda \) and \({\mathcal {S}}^\beta \)

This section is devoted to the derivation of some important and new estimates of the functions \({\mathcal {R}}^\lambda \) and \({\mathcal {S}}^\beta \) on \((0,\infty )\) and to some first applications of these estimates. We recall the definitions of these functions from (2.7), viz.

$$\begin{aligned} {\mathcal {R}}^\lambda (t)&= {\mathcal {L}}^{-1}\left\{ \frac{s^{l(\alpha )-\lambda }}{sQ(s)} \right\} (t) , \quad \lambda \in \left\{ 0,\alpha _1, \alpha _2 \right\} , \\ {\mathcal {S}}^\beta (t)&= {\mathcal {L}}^{-1}\left\{ \frac{s^{l(\alpha )-\beta }}{Q(s)} \right\} (t) , \quad \beta \in \left\{ \alpha _1, \alpha _2 ,l(\alpha )\right\} , \end{aligned}$$

where \(l(\alpha ) = \alpha _1 + \alpha _2\) and \(Q(s) = s^{\alpha _1 + \alpha _2} -a_{11}s^{\alpha _2} - a_{22}s^{\alpha _1} + \det A\).

Lemma 8

Let \(\alpha _1, \alpha _2 \in (0, 1]\) and denote \( \nu = \min \{\alpha _1, \alpha _2 \}\). Assume there are no zeros of the characteristic function Q in the closed right-half complex plane. Then, the following estimates hold for \(\lambda \in \left\{ 0,\alpha _1, \alpha _2 \right\} \) and \(\beta \in \left\{ \alpha _1, \alpha _2 ,l(\alpha )\right\} \):

$$\begin{aligned} {\mathcal {R}}^\lambda (t)&= O(t^{-\nu }) \quad \text{ as } t \rightarrow \infty , \end{aligned}$$
(4.1)
$$\begin{aligned} {\mathcal {S}}^\beta (t)&= O(t^{-\nu -1}) \quad \text{ as } t \rightarrow \infty , \end{aligned}$$
(4.2)
$$\begin{aligned} {\mathcal {S}}^\beta (t)&= O(t^{\nu -1}) \quad \text{ as } t \rightarrow 0 . \end{aligned}$$
(4.3)

Moreover,

$$\begin{aligned} \int _{0}^\infty |{\mathcal {S}}^\beta (t)| dt < \infty . \end{aligned}$$
(4.4)

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 Appendix A at the end of the paper.

Our first application of Lemma 8 deals with estimates for the convolution of \({\mathcal {S}}^\beta \) and a continuous function.

Theorem 3

Let \(\alpha _1, \alpha _2 \in (0,1]\) and \(\beta \in \left\{ \alpha _1,\alpha _2,l(\alpha ) \right\} \). For each continuous function \(g : [0, \infty ) \rightarrow {\mathbb {R}}\), we put

$$\begin{aligned} F_g^\beta (t) := {\mathcal {S}}^\beta *g(t) = \int _0^t{\mathcal {S}}^\beta (t-s)g(s)ds \end{aligned}$$
(4.5)

and

$$\begin{aligned} {\bar{F}}_g^\beta (t) := |{\mathcal {S}}^\beta | *|g| (t) = \int _0^t\mathcal |S^\beta (t-s)| \cdot | g(s) | ds. \end{aligned}$$
(4.6)

Assume that all zeros of the characteristic function Q are in the open left-half complex plane. Then, the following statements hold.

  1. (i)

    If g is bounded then \({\bar{F}}_g^\beta \) and \(F_g^\beta \) are also bounded.

  2. (ii)

    If \(\lim _{t \rightarrow \infty } g(t) = 0\) then \(\lim _{t \rightarrow \infty } {\bar{F}}_g^\beta (t) = \lim _{t \rightarrow \infty } F_g^\beta (t) = 0\).

  3. (iii)

    If there exists some \(\eta \ge 0\) such that \(g(t) = O(t^{-\eta })\) for \(t \rightarrow \infty \) then

    $$\begin{aligned} {\bar{F}}_g^\beta (t) = O(t^{-\mu }) \text{ and } F_g^\beta (t) = O(t^{-\mu }) \text{ for } t \rightarrow \infty \end{aligned}$$

    where \(\mu = \min \left\{ \alpha _1, \alpha _2, \eta \right\} .\)

Proof

We denote \(\nu = \min \{\alpha _1, \alpha _2 \}\). Since, by definition, \(|F_g^\beta (t)| \le {\bar{F}}_g^\beta (t)\), the claims for \(F_g^\beta \) immediately follow from those for \({\bar{F}}_g^\beta \), and therefore it suffices to explicitly prove the latter.

Statement (i) is merely the special case of \(\eta = 0\) of part (iii).

To prove (ii), we note that \({\bar{F}}_g^\beta (t) \ge 0\) by definition. Therefore, it is sufficient to show that for every \(\varepsilon > 0\) there exists a constant \({\widetilde{T}} = {\widetilde{T}}(\varepsilon )\) such that

$$\begin{aligned} {\bar{F}}_g^\beta (t) \le \varepsilon \quad \text{ for } \text{ all } t > {\widetilde{T}}. \end{aligned}$$
(4.7)

Since this is trivially fulfilled if \(g(t) = 0\) for all t, we from now on assume that \(g(t) \ne 0\) for some t, and hence \(\Vert g \Vert _\infty > 0\).

Our first observation is then that, from (4.2) and (4.3), we know that there exists some constant \(C > 0\) such that

$$\begin{aligned} |S^\beta (t)| \le {\left\{ \begin{array}{ll} C t^{-\nu -1} &{} \text{ for } t \ge 1, \\ C t^{\nu -1} &{} \text{ for } t \le 1. \end{array}\right. } \end{aligned}$$
(4.8)

Given an arbitrary \(\varepsilon > 0\), due to our assumption on g we may then find some \({\hat{T}} > 0\) such that \(|g(t)| < \nu \varepsilon / (3 C)\) for all \(t > {\hat{T}}\). Using these values \({\hat{T}}\) and C, we then define

$$\begin{aligned} {\widetilde{T}} = {\hat{T}} + \max \left\{ 1, \left( \frac{3 C \Vert g \Vert _\infty {\hat{T}}}{\varepsilon }\right) ^{1/(\nu +1)} \right\} . \end{aligned}$$

For \(t > {\widetilde{T}} \ge {\hat{T}} + 1\), we can then write

$$\begin{aligned} {\bar{F}}_g^\beta (t)&= \int _0^{{\hat{T}}} |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds + \int _{{\hat{T}}}^{t-1} \! |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds \nonumber \\&\qquad \qquad + \int _{t-1}^{t} \! |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds \nonumber \\&= F_1(t) + F_2(t) +F_3(t). \end{aligned}$$
(4.9)

Our goal now is to show that, under these assumptions, \(F_j(t) \le \varepsilon /3\) for \(j = 1, 2, 3\), which implies (4.7) and thus suffices to prove part (ii) of the theorem. In this context, we see that

$$\begin{aligned} F_1(t)&= \int _0^{{\hat{T}}} |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds \le \Vert g \Vert _\infty C \int _0^{{\hat{T}}} (t - s)^{-\nu -1} d s \\&\le \Vert g \Vert _\infty C {\hat{T}} (t - {\hat{T}})^{-\nu -1} < \frac{\varepsilon }{3}. \end{aligned}$$

because here \(t - s \ge t - {\hat{T}} > {\widetilde{T}} - {\hat{T}} \ge 1\), so that we may use the first of the bounds given in (4.8). 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 > {\widetilde{T}}\) and the definition of \({\widetilde{T}}\).

Furthermore,

$$\begin{aligned} F_2(t)&= \int _{{\hat{T}}}^{t-1} |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds \le \frac{\nu \varepsilon }{3 C} C \int _{{\hat{T}}}^{t-1} (t - s)^{-\nu -1} d s \\&= \frac{\nu \varepsilon }{3} \frac{1}{\nu }\left( 1 - (t - {\hat{T}})^{-\nu } \right) < \frac{\varepsilon }{3} \end{aligned}$$

because here \(s \ge {\hat{T}}\), so that \(|g(s)| \le \nu \varepsilon / (3 C)\), and \(t-s \ge 1\), so we may once again use the first bound of (4.8).

Finally,

$$\begin{aligned} F_3(t)&= \int _{t-1}^t |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds \le \frac{\nu \varepsilon }{3 C} C \int _{t-1}^t (t - s)^{\nu -1} d s = \frac{\nu \varepsilon }{3} \frac{1}{\nu }= \frac{\varepsilon }{3} \end{aligned}$$

where now t and s are such that we may invoke the second bound of (4.8) but, as in the previous step, \(s \ge {\hat{T}}\), so that once again \(|g(s)| \le \nu \varepsilon / (3 C)\). This completes the proof of part (ii) of the theorem.

For the proof of (iii), we note that (4.8) is valid in this case too. Moreover, since we are interested in the asymptotic behaviour of \({\bar{F}}_g^\beta (t)\) for large t, we may assume without loss of generality that \(t \ge 2\). Then we write

$$\begin{aligned} {\bar{F}}_g^\beta (t)&= \int _0^{1}|{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds +\int _1^{t/2} |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)| ds \nonumber \\&{} \quad + \int _{t/2}^{t-1} |{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds + \int _{t-1}^{t}|{\mathcal {S}}^\beta (t-s)| \cdot |g(s)|ds \nonumber \\&= F_4(t) + F_5(t) + F_6(t) + F_7(t), \end{aligned}$$
(4.10)

and we need to show that \(F_j(t) = O(t^{-\mu })\) for \(j = 4, 5, 6, 7\).

In this connection, we first note that, by assumption,

$$\begin{aligned} |g(t)| \le C' t^{-\eta } \quad \forall t \ge 1 \end{aligned}$$
(4.11)

with some \(C' > 0\), so that the upper branch of (4.8) implies

$$\begin{aligned} 0 \le F_4(t)&\le C \Vert g \Vert _\infty \int _0^1 (t-s)^{-\nu -1} ds = \frac{C \Vert g \Vert _\infty }{\nu } \left( (t-1)^{-\nu } - t^{-\nu } \right) \\&< \frac{C \Vert g \Vert _\infty }{\nu } (t-1)^{-\nu } = O(t^{-\nu }) = O(t^{-\mu }) \end{aligned}$$

and

$$\begin{aligned} 0 \le F_5(t)&\le C C' \int _1^{t/2} (t-s)^{-\nu -1} s^{-\eta } ds \le C C' \int _1^{t/2} (t-s)^{-\nu -1} ds \\&= \frac{C C'}{\nu } \left( (t-1)^{-\nu } - (t/2)^{-\nu } \right) < \frac{C C'}{\nu } (t-1)^{-\nu } = O(t^{-\nu }) = O(t^{-\mu }) \end{aligned}$$

as well as

$$\begin{aligned} 0 \le F_6(t)&\le C C' \int _{t/2}^{t-1} (t-s)^{-\nu -1} s^{-\eta } ds \le C C' \left( \frac{t}{2} \right) ^{-\eta } \int _{t/2}^{t-1} (t-s)^{-\nu -1} ds \\&= 2^\eta \frac{C C'}{\nu } t^{-\eta } \left( 1 - (t/2)^{-\nu } \right) < 2^\eta \frac{C C'}{\nu } t^{-\eta } = O(t^{-\eta }) = O(t^{-\mu }). \end{aligned}$$

For the remaining part we need to invoke the second branch of (4.8) in combination with (4.11) to derive

$$\begin{aligned} 0 \le F_7(t)&\le C C' \int _{t-1}^t (t-s)^{\nu -1} s^{-\eta } ds \le C C' (t-1)^{-\eta } \int _{t-1}^t (t-s)^{\nu -1} ds \\&= \frac{C C'}{\nu } (t-1)^{-\eta } = O(t^{-\eta }) = O(t^{-\mu }), \end{aligned}$$

thus completing the proof. \(\square \)

As an immediate application of Theorem 3(iii), we can conclude that

$$\begin{aligned} {\bar{F}}^\beta _{g_1} (t) = O(t^{-\nu }) \quad \text{ as } t \rightarrow \infty \qquad \text{ for } g_1(t) = \min \{ 1, t^{-\nu } \}. \end{aligned}$$
(4.12)

Moreover, assuming \(\nu < 1\) and setting

$$\begin{aligned} g_2(t) = t^{-\nu } - g_1(t) = {\left\{ \begin{array}{ll} t^{-\nu } - 1 &{} \text{ for } t \in [0,1], \\ 0 &{} \text{ for } t > 1, \end{array}\right. } \end{aligned}$$

we can obtain (using Lemma 8 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 \ge 2\) then we have

    $$\begin{aligned} \nonumber {\bar{F}}^\beta _{g_2}(t)&= \int _0^1 | {\mathcal {S}}^\beta (t-s) | (s^{-\nu }-1) ds \le C \int _0^1 (t-s)^{-\nu -1} s^{-\nu } ds \\ \nonumber&= C t^{-2\nu } B_{1/t} (1-\nu , -\nu ) \\&= \frac{C}{1-\nu } t^{-\nu -1} {} {}_2 F_1(1-\nu , 1+\nu ; 2-\nu ; t^{-1}) \le C' t^{-\nu -1} \end{aligned}$$
    (4.13)

    with some \(C' > 0\).

  • If \(t \in [1,2]\) then

    $$\begin{aligned} \nonumber {\bar{F}}^\beta _{g_2}(t)&= \int _0^1 | {\mathcal {S}}^\beta (t-s) | (s^{-\nu }-1) ds \le C \int _0^1 (t-s)^{\nu -1} s^{-\nu } ds \\ \nonumber&= C B_{1/t} (1-\nu , \nu ) \\&= \frac{C}{1-\nu } t^{\nu -1} {} {}_2 F_1(1-\nu , 1-\nu ; 2-\nu ; t^{-1}) \le C'' t^{-\nu } \end{aligned}$$
    (4.14)

    with some \(C'' > 0\).

Since \({\bar{F}}^\beta _{g_1}(t) + {\bar{F}}^\beta _{g_2}(t) = {\bar{F}}^\beta _{g_1 + g_2}(t)\), we can summarize the observations of eqs. (4.12), (4.13) and (4.14) in the following way:

Remark 1

Assuming that \(\nu = \min \{ \alpha _1, \alpha _2\} < 1\), there exists a constant C such that, for all \(t \ge 1\) and \(\beta \in \left\{ \alpha _1, \alpha _2, l(\alpha ) \right\} \),

$$\begin{aligned} t^\nu \int _0^t |{\mathcal {S}}^\beta (t-s)| \frac{1}{s^\nu }ds \le C. \end{aligned}$$
(4.15)

5 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, thus presenting the main new results of our work.

5.1 Asymptotic behaviour of solutions to fractional linear planar systems

Consider the non-homogeneous linear two-component incommensurate fractional-order system

$$\begin{aligned} ^C D^\alpha _{0^+} x(t)&= Ax(t) + f(t),\; t>0, \end{aligned}$$
(5.1a)
$$\begin{aligned} x(0)&= x^0 \in {\mathbb {R}}^2, \end{aligned}$$
(5.1b)

where \(\alpha =(\alpha _1,\alpha _2)\in (0,1]\times (0,1]\) is a multi index, \(A = \left( a_{ij} \right) \in {\mathbb {R}}^{2 \times 2}\) is a square real matrix and \(f = (f_1, f_2)\) is a continuous vector valued function which is exponentially bounded on \([0,\infty )\).

Theorem 4

Suppose that all zeros of the characteristic function

$$\begin{aligned} Q(s) = s^{\alpha _1+\alpha _2}-a_{11}s^{\alpha _2}-a_{22}s^{\alpha _1}+ \det A \end{aligned}$$

of the problem (5.1a) lie in the open left-half of the complex plane. Then, the following statements hold.

  1. (i)

    If f is bounded, then for any \(x^0\in {\mathbb {R}}^2\) the solution to (5.1) is also bounded.

  2. (ii)

    If \(\lim _{t \rightarrow \infty } f(t) = 0\) then the solution to (5.1) tends to 0 when \(t \rightarrow \infty \) for any \(x^0\in {\mathbb {R}}^2\).

  3. (iii)

    If \(\Vert f(t) \Vert = O(t^{-\eta })\) as \(t \rightarrow \infty \) with some \(\eta > 0\) then every solution x of (5.1a) behaves as \(\Vert x(t) \Vert = O(t^{-\mu })\) for \(t \rightarrow \infty \) where \(\mu = \min \left\{ \alpha _1, \alpha _2, \eta \right\} .\)

Proof

The proof is straightforward by combining Lemmas 1 and  8 and Theorem 3. \(\square \)

Based on Theorem 4 and Lemmas 2, 4, 5, 6 and 7, we obtain the following corollary.

Corollary 2

Let

$$\begin{aligned} q_1= \frac{\sin \frac{\alpha _1\pi }{2}}{\sin \frac{(\alpha _1+\alpha _2)\pi }{2}} \quad \text{ and } \quad q_2= \frac{\sin \frac{\alpha _2\pi }{2}}{\sin \frac{(\alpha _1+\alpha _2)\pi }{2}}. \end{aligned}$$

The statements of Theorem 4 (i), (ii) and (iii) are true if one of the following conditions is satisfied.

  1. (i)

    \(a_{11},a_{22}\le 0\) and \(\det A>0\).

  2. (ii)

    \(a_{11}=0\), \(a_{22}>0\), \(\det A>0\) and

    $$\begin{aligned} (a_{22}q_1)^{\alpha _1/\alpha _2}a_{22}q_2< \det A. \end{aligned}$$
  3. (iii)

    \(a_{22}=0\), \(a_{11}>0\), \(\det A>0\) and

    $$\begin{aligned} (a_{11}q_2)^{\alpha _2/\alpha _1}a_{11}q_1< \det A. \end{aligned}$$
  4. (iv)

    \(a_{11}, a_{22}, \det A > 0\) and one of the following conditions holds:

    (iv)\(_1\):

    \(a_{11}q_2+a_{22}q_1 > 1\) and \(a_{11}q_2\left( (a_{11}+a_{22})q_2 \right) ^{\alpha _2/\alpha _1} + a_{22}(a_{11}+a_{22})q_2^2 \le \det A\);

    (iv)\(_2\):

    \( a_{11}q_2+a_{22}q_1 \le 1\) and \( a_{11}q_1 + a_{22}q_2 < \det A.\)

  5. (v)

    \(a_{11} < 0\), \(a_{22}, \det A > 0\) and one of the following conditions holds:

    (v)\(_1\):

    \( a_{11}q_2+a_{22}q_1 > 1\) and \(\left( a_{22}q_1 \right) ^{\alpha _1/\alpha _2}a_{22}q_2 \le \det A\);

    (v)\(_2\):

    \( a_{11}q_2+a_{22}q_1 \le 1\) and \( a_{22}q_2 \le \det A\).

5.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

$$\begin{aligned} ^C D^\alpha _{0^+}x(t)&= Ax(t) +f(x(t)),\quad t > 0, \end{aligned}$$
(5.2a)
$$\begin{aligned} x(0)&= x^0 \in \varOmega \subset {\mathbb {R}}^2, \end{aligned}$$
(5.2b)

where \(\alpha =(\alpha _1,\alpha _2)\in (0,1]\times (0,1]\) is a multi-index, \(A = \left( a_{ij} \right) \in {\mathbb {R}}^{2 \times 2}\) is a square real matrix, \(\varOmega \) is an open subset of \({\mathbb {R}}^2\) containing the origin and \(f:\varOmega \rightarrow {\mathbb {R}}^2\) is locally Lipschitz continuous at the origin such that \(f(0)=0\) and \(\lim _{r\rightarrow 0} l_f (r)=0\) with

$$\begin{aligned} l_f(r) := \sup _{x,y\in B(0,r),\;x\ne y}\frac{\Vert f(x)-f(y)\Vert }{\Vert x-y\Vert }. \end{aligned}$$

In the context of first order differential equations, it is natural to talk about exponential stability; an appropriate generalization of this classical notion to the fractional order setting is the concept of Mittag-Leffler stability. In the case of a commensurate system, this can be defined as follows (cf., e.g., [14,  Definition 4.1]):

Definition 1

The trivial solution of the system (5.2a) with \(\alpha _1 = \alpha _2 \in (0,1)\) is said to be Mittag–Leffler stable if there exist some \(b, \lambda >0\), some \(B \subset {\mathbb {R}}^2\) and some locally Lipschitz function \(m : B \rightarrow [0, \infty )\) with the property \(m(0) = 0\) such that, for all \(x^0 \in B\),

$$\begin{aligned} \Vert \varphi (t, x^0) \Vert \le \left[ m(x^0) E_{\alpha _1} (-\lambda t^{\alpha _1}) \right] ^b \end{aligned}$$
(5.3)

for all \(t \ge 0\), where \(\varphi (\cdot , x^0)\) denotes the solution of the initial value problem (5.2).

In view of the well known asymptotic properties of the Mittag-Leffler function for \(t \rightarrow 0\) and for \(t \rightarrow \infty \), [12,  Section 3.4], the property (5.3) can essentially be reformulated in the form of the requirement that \(\Vert \varphi (t, x^0) \Vert \) remains bounded for \(t \rightarrow 0\) and exhibits an (at least) algebraic decay to zero as \(t \rightarrow \infty \) whenever \(x^0\) is sufficiently close to 0. Therefore, it seems natural to extend this definition to the case of non-commensurate fractional differential equation systems that is relevant in our context in the following way:

Definition 2

The trivial solution of (5.2a) is Mittag-Leffler stable if there exist positive constants \(\gamma ,m\) and \(\delta \) such that for any initial condition \(x^0\in B(0,\delta )\), the solution \(\varphi (\cdot , x^0)\) of the initial value problem (5.2) exists globally on the interval \([0,\infty )\) and

$$\begin{aligned} \displaystyle \max \{\sup _{t\in [0,1]}\Vert \varphi (t, x^0)\Vert ,\sup _{t \ge 1} t^\gamma \Vert \varphi (t, x^0)\Vert \}\le m. \end{aligned}$$

Our aim is to prove the following theorem.

Theorem 5

Suppose that all zeros of the characteristic function \(Q(s) =\) \(s^{\alpha _1+\alpha _2}-a_{11}s^{\alpha _2}-a_{22}s^{\alpha _1}+ \det A\) lie in the open left-half of the complex plane. Then, the trivial solution of differential equation (5.2a) is Mittag-Leffler stable. More precisely, there exist constants \(\delta , \varepsilon > 0\) such that for any \(\Vert x^0\Vert < \delta \), the unique solution \(\varphi (\cdot ,x^0)\) of the initial value problem (5.2) exists globally on \([0,\infty )\) and \(\sup _{t \ge 1} t^\nu \Vert \varphi (t,x^0)\Vert \le \varepsilon \) with \(\nu = \min \{\alpha _1, \alpha _2\}.\)

As shown above, we see that Lemmas 2, 4, 5, 6 and 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, we obtain the result below.

Corollary 3

Let

$$\begin{aligned} q_1= \frac{\sin \frac{\alpha _1\pi }{2}}{\sin \frac{(\alpha _1+\alpha _2)\pi }{2}} \quad \text{ and } \quad q_2= \frac{\sin \frac{\alpha _2\pi }{2}}{\sin \frac{(\alpha _1+\alpha _2)\pi }{2}}. \end{aligned}$$

The statement of Theorem 5 is true if one of the following conditions is satisfied:

  1. (i)

    \(a_{11},a_{22}\le 0\) and \(\det A>0\).

  2. (ii)

    \(a_{11}=0\), \(a_{22}>0\), \(\det A>0\) and

    $$\begin{aligned} (a_{22}q_1)^{\alpha _1/\alpha _2}a_{22}q_2< \det A. \end{aligned}$$
  3. (iii)

    \(a_{22}=0\), \(a_{11}>0\), \(\det A>0\) and

    $$\begin{aligned} (a_{11}q_2)^{\alpha _2/\alpha _1}a_{11}q_1< \det A. \end{aligned}$$
  4. (iv)

    \(a_{11}, a_{22}, det A > 0\) and one of the following conditions holds:

    (iv)\(_1\):

    \(a_{11}q_2+a_{22}q_1 > 1\) and \(a_{11}q_2\left( (a_{11}+a_{22})q_2 \right) ^{\alpha _2/\alpha _1} + a_{22}(a_{11}+a_{22})q_2^2 \le \det A\);

    (iv)\(_2\):

    \( a_{11}q_2+a_{22}q_1 \le 1\) and \( a_{11}q_1 + a_{22}q_2 < \det A.\)

  5. (v)

    \(a_{11} < 0\), \(a_{22}, det A > 0\) and one of the following conditions holds:

    (v)\(_1\):

    \( a_{11}q_2+a_{22}q_1 > 1\) and \(\left( a_{22}q_1 \right) ^{\alpha _1/\alpha _2}a_{22}q_2 \le \det A\);

    (v)\(_2\):

    \( a_{11}q_2+a_{22}q_1 \le 1\) and \( a_{22}q_2 \le \det A\).

Proof of Theorem 5

From the assumption of the theorem that f is locally Lipschitz continuous at the origin, we can find a constant \(\varepsilon _0 > 0\) such that the function f is Lipschitz continuous on \(B(0,\varepsilon _0)\). Denote by \(\hat{f}\) a Lipschitz extension of f to \({\mathbb {R}}^2\). This means that \(\hat{f}\) is globally Lipschitz continuous and \(\hat{f}(x)=f(x)\) on \(B(0,\varepsilon _0)\). We now focus on the system

$$\begin{aligned} ^C D^\alpha _{0^+}x(t)&= Ax(t) +\hat{f}(x(t)), t > 0, \end{aligned}$$
(5.4a)
$$\begin{aligned} x(0)&= x^0. \end{aligned}$$
(5.4b)

Then, for any \(x^0\in B(0,\varepsilon _0)\), its unique solution \({\hat{\varphi }}(\cdot ,x^0)=({\hat{\varphi }}_1(\cdot ,x^0),{\hat{\varphi }}_2(\cdot ,x^0))^{\!\mathrm T}\) on \([0,\infty )\) satisfies the relationships

$$\begin{aligned} {\hat{\varphi }}_1(\cdot ,x^0)&= \left( {\mathcal {R}}^0(t)-a_{22}{\mathcal {R}}^{\alpha _2}(t) \right) x_1^0 + a_{12}{\mathcal {R}}^{\alpha _1}(t)x^0_2 \nonumber \\&\quad + \left( ( {\mathcal {S}}^{\alpha _1} - a_{22}{\mathcal {S}}^{l(\alpha )} ) *\hat{f}_1({\hat{\varphi }}(\cdot ,x^0)) \right) (t) \nonumber \\&\quad + a_{12}\left( {\mathcal {S}}^{l(\alpha )} *\hat{f}_2 ({\hat{\varphi }}(\cdot ,x^0)) \right) (t), \end{aligned}$$
(5.5a)
$$\begin{aligned} {\hat{\varphi }}_2(t,x^0)&= a_{21}{\mathcal {R}}^{\alpha _2}(t)x^0_1 + \left( {\mathcal {R}}^0(t)-a_{11}{\mathcal {R}}^{\alpha _1}(t) \right) x_2^0 \nonumber \\&\quad +a_{21}\left( {\mathcal {S}}^{l(\alpha )} *\hat{f}_1({\hat{\varphi }}(\cdot ,x^0)) \right) (t) \nonumber \\&\quad + \left( ( {\mathcal {S}}^{\alpha _2}- a_{11}{\mathcal {S}}^{l(\alpha )}) *\hat{f}_2({\hat{\varphi }}(\cdot ,x^0)) \right) (t). \end{aligned}$$
(5.5b)

To show the Mittag-Leffler stability of the trivial solution to the original system, we will prove that for any small initial value vector, the unique solution of the system (5.4) is contained in the space \(C_{\infty }([0,\infty );{\mathbb {R}}^2)\) which is equipped with the norm

$$\begin{aligned} \Vert \xi \Vert _w:=\max \{ \sup _{t\in [0,1]}\Vert \xi (t)\Vert , \sup _{t\ge 1}t^\nu \Vert \xi (t)\Vert \}. \end{aligned}$$

It is easy to see that \(C_w([0,\infty );{\mathbb {R}}^2):=\{\xi \in C_{\infty }([0,\infty );{\mathbb {R}}^2):\Vert \xi \Vert _w<\infty \}\) is a Banach space with the norm \(\Vert \cdot \Vert _w\). For \(\varepsilon > 0\), let \(B_{C_w}(0,\varepsilon ) := \{\xi \in C_{\infty }([0,\infty );{\mathbb {R}}^2):\Vert \xi \Vert _w\le \varepsilon \}\).

Based on the representation in eq. (5.5), we establish a Lyapunov–Perron type operator \({\mathcal {T}}_{x^0}\) on the space \(C_w([0,\infty );{\mathbb {R}}^2)\) in the following way: For any \(\xi \in C_w([0,\infty );{\mathbb {R}}^2)\), let

$$\begin{aligned} ({\mathcal {T}}_{x_0}\xi )_1(t)&:= \left( {\mathcal {R}}^0(t)-a_{22}{\mathcal {R}}^{\alpha _2}(t) \right) x_1^0 + a_{12}{\mathcal {R}}^{\alpha _1}(t)x^0_2 \nonumber \\&\quad + \left( ( {\mathcal {S}}^{\alpha _1} - a_{22}{\mathcal {S}}^{l(\alpha )} )*\hat{f}_1(\xi (\cdot )) \right) (t) + a_{12}\left( {\mathcal {S}}^{l(\alpha )}*\hat{f}_2 (\xi (\cdot )) \right) (t), \\ ({\mathcal {T}}_{x_0}\xi )_2(t)&:= a_{21}{\mathcal {R}}^{\alpha _2}(t)x^0_1 + \left( {\mathcal {R}}^0(t)-a_{11}{\mathcal {R}}^{\alpha _1}(t) \right) x_2^0 \nonumber \\&\quad + a_{21}\left( {\mathcal {S}}^{l(\alpha )}*\hat{f}_1(\xi (\cdot )) \right) (t) +\left( ( {\mathcal {S}}^{\alpha _2}- a_{11}{\mathcal {S}}^{l(\alpha )})*\hat{f}_2(\xi (\cdot ))\right) (t). \end{aligned}$$

On the interval [0, 1], we have

$$\begin{aligned} |({\mathcal {T}}_{x^0}\xi )_1(t)|&\le (|{\mathcal {R}}^0(t)|+|a_{22}| \cdot |{\mathcal {R}}^{\alpha _2}(t)| )|x_1^0| + |a_{12}| \cdot |{\mathcal {R}}^{\alpha _1}(t)| \cdot |x^0_2|\nonumber \\&\quad + l_{\hat{f}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \int _0^t \left( |{\mathcal {S}}^{\alpha _1}(s)| + (|a_{22}|+|a_{12}|)|{\mathcal {S}}^{l(\alpha )}(s)|\right) ds\nonumber \\&\le \sup _{t\in [0,1]} \left( \left( |{\mathcal {R}}^0(t)|+|a_{22}| \cdot |{\mathcal {R}}^{\alpha _2}(t)| \right) |x_1^0| + |a_{12}| \cdot |{\mathcal {R}}^{\alpha _1}(t)| \cdot |x^0_2| \right) \nonumber \\&\quad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w C (1+|a_{22}|+|a_{12}|)\int _0^1 \frac{1}{s^{1-\nu }}ds\nonumber \\&\le \sup _{t\in [0,1]}\left( \left( |{\mathcal {R}}^0(t)|+|a_{22}| \cdot |{\mathcal {R}}^{\alpha _2}(t)| \right) |x_1^0| + |a_{12}| \cdot |{\mathcal {R}}^{\alpha _1}(t)| \cdot |x^0_2|\right) \nonumber \\&\quad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \frac{C(1+|a_{22}|+|a_{12}|)}{\nu } \end{aligned}$$
(5.6)

and

$$\begin{aligned} |({\mathcal {T}}_{x^0}\xi )_2(t)|&\le |a_{21}| \cdot |{\mathcal {R}}^{\alpha _2}(t)| \cdot |x^0_1| +(|{\mathcal {R}}^0(t)|+|a_{11}| \cdot |{\mathcal {R}}^{\alpha _1}(t)|) |x_2^0|\nonumber \\&\quad + l_{\hat{f}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \int _0^t \left( |{\mathcal {S}}^{\alpha _2}(s)| + (|a_{12}|+|a_{11}|)|{\mathcal {S}}^{l(\alpha )}(s)|\right) ds\nonumber \\&\le \sup _{t\in [0,1]} \left( |a_{21}| \cdot |{\mathcal {R}}^{\alpha _2}(t)| \cdot |x^0_1|+(|{\mathcal {R}}^0(t)| +|a_{11}| \cdot |{\mathcal {R}}^{\alpha _1}(t)|)|x_2^0| \right) \nonumber \\&\quad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty ) \Vert \xi \Vert _wC(1+|a_{21}|+|a_{11}|)\int _0^1 \frac{1}{s^{1-\nu }}ds\nonumber \\&\le \sup _{t\in [0,1]} \left( |a_{21}| \cdot |{\mathcal {R}}^{\alpha _2}(t)| \cdot |x^0_1| +(|{\mathcal {R}}^0(t)|+|a_{11}| \cdot |{\mathcal {R}}^{\alpha _1}(t)|)|x_2^0| \right) \nonumber \\&\quad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w\frac{C(1+|a_{21}|+|a_{11}|)}{\nu }. \end{aligned}$$
(5.7)

For \(t\in [1,\infty )\), we have

$$\begin{aligned} \nonumber&t^\nu |( {\mathcal {T}}_{x_0}\xi )_1(t)| \\&\quad \le \sup _{t\ge 1}\Big (\left( t^\nu |{\mathcal {R}}^0(t)|+|a_{22}|t^\nu |{\mathcal {R}}^{\alpha _2}(t)| \right) |x_1^0| + |a_{12}|t^\nu |{\mathcal {R}}^{\alpha _1}(t)| \cdot |x^0_2|\Big )\nonumber \\ \nonumber&\qquad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty ) t^\nu \\ \nonumber&\qquad \quad \times \int _0^t \left( |{\mathcal {S}}^{\alpha _1}(t-s)| + (|a_{22}|+|a_{12}|)|{\mathcal {S}}^{l(\alpha )}(t-s)|\right) s^{-\nu }s^{\nu }|\xi (s)|ds\nonumber \\ \nonumber&\quad \le \sup _{t\ge 1}\Big (\left( t^\nu |{\mathcal {R}}^0(t)|+|a_{22}|t^\nu |{\mathcal {R}}^{\alpha _2}(t)| \right) |x_1^0| + |a_{12}|t^\nu |{\mathcal {R}}^{\alpha _1}(t)| \cdot |x^0_2|\Big )\nonumber \\&\qquad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \\&\qquad \quad \times \sup _{t\ge 1} t^\nu \int _0^t \left( |{\mathcal {S}}^{\alpha _1}(t-s)|+ (|a_{22}|+|a_{12}|)|{\mathcal {S}}^{l(\alpha )}(t-s)|\right) s^{-\nu }ds \nonumber \end{aligned}$$
(5.8)

and

$$\begin{aligned} \nonumber&t^\nu |( {\mathcal {T}}_{x_0}\xi )_2(t)|\\&\quad \le \sup _{t\ge 1}\Big (|a_{21}|t^\nu |{\mathcal {R}}^{\alpha _2}(t)| \cdot |x^0_1|+(t^\nu |{\mathcal {R}}^0(t)|+|a_{11}|t^\nu |{\mathcal {R}}^{\alpha _1}(t)|)|x_2^0|\Big )\nonumber \\ \nonumber&\qquad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty )t^\nu \\ \nonumber&\qquad \quad \times \int _0^t \left( |{\mathcal {S}}^{\alpha _2}(t-s)|+ (|a_{21}|+|a_{11}|)|{\mathcal {S}}^{l(\alpha )}(t-s)|\right) s^{-\nu }s^\nu \Vert \xi (s)\Vert ds \nonumber \\ \nonumber&\quad \le \sup _{t\ge 1}\Big (|a_{21}|t^\nu |{\mathcal {R}}^{\alpha _2}(t)| \cdot |x^0_1|+(t^\nu |{\mathcal {R}}^0(t)|+|a_{11}|t^\nu |{\mathcal {R}}^{\alpha _1}(t)|)|x_2^0|\Big )\nonumber \\&\qquad +l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \\&\qquad \quad \times \sup _{t\ge 1}t^\nu \int _0^t \left( |{\mathcal {S}}^{\alpha _2}(t-s)|+ (|a_{21}|+|a_{11}|)|{\mathcal {S}}^{l(\alpha )}(t-s)|\right) s^{-\nu } ds. \nonumber \end{aligned}$$
(5.9)

From (5.8) and (5.9), we obtain the estimates

$$\begin{aligned}&\Vert ({\mathcal {T}}_{x^0}\xi )_1\Vert _{w,1} \\ \nonumber&\quad \le (\Vert {\mathcal {R}}^0\Vert _{w,1}+|a_{22}| \cdot \Vert {\mathcal {R}}^{\alpha _2}\Vert _{w,1} +|a_{12}| \cdot \Vert {\mathcal {R}}^{\alpha _1}\Vert _{w,1})\Vert x^0\Vert \\&\qquad +l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \left( \frac{C}{\nu } (1\!+\!|a_{12}|\!+\!|a_{22}|)+M_{\alpha _1} +(|a_{22}|+|a_{12}|)M_{l(\alpha )}\right) \nonumber \end{aligned}$$
(5.10)

and

$$\begin{aligned}&\Vert ({\mathcal {T}}_{x^0}\xi )_2\Vert _{w,1}\\ \nonumber&\quad \le (\Vert {\mathcal {R}}^0\Vert _{w,1}+|a_{21}| \cdot \Vert {\mathcal {R}}^{\alpha _2}\Vert _{w,1} +|a_{11}| \cdot \Vert {\mathcal {R}}^{\alpha _1}\Vert _{w,1})\Vert x^0\Vert \\&\qquad +l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \left( \frac{C}{\nu } (1\!+\!|a_{21}|\!+\!|a_{11}|) +M_{\alpha _2}+(|a_{21}|+|a_{11}|)M_{l(\alpha )}\right) \nonumber \end{aligned}$$
(5.11)

where \(\Vert \xi \Vert _{w,1}\!:=\!\max \{\sup _{t\in [0,1]}|\xi (t)|,\sup _{t\ge 1}t^\nu |\xi (t)|\}\) for any \(\xi \!\in \! C_{\infty }([0,\infty );{\mathbb {R}})\) and \(M_{\beta }=\sup _{t\ge 1} t^\nu \int _0^t |{\mathcal {S}}^{\beta }(t-s)| s^{-\nu }ds\) for \(\beta \in \{ \alpha _1, \alpha _2, l(\alpha ) \}\). By (5.10) and (5.11), we see that

$$\begin{aligned}&\Vert ({\mathcal {T}}_{x^0}\xi )\Vert _{w} \\&\quad \le \left( 2\Vert {\mathcal {R}}^0\Vert _{w,1}+(|a_{11}|+|a_{12}|)\Vert {\mathcal {R}}^{\alpha _1}\Vert _{w,1} +(|a_{21}|+|a_{22}|) \Vert {\mathcal {R}}^{\alpha _2}\Vert _{w,1} \right) \Vert x^0\Vert \\&\qquad + l_{{\hat{f}}}(\Vert \xi \Vert _\infty )\Vert \xi \Vert _w \\&\qquad \qquad \quad \times \left( \frac{C}{\nu }\left( 2+\sum _{i,j=1}^2|a_{ij}|\right) +M_{\alpha _1} +M_{\alpha _2}+\sum _{i,j=1}^2|a_{ij}|M_{l(\alpha )}\right) . \end{aligned}$$

On the other hand, by virtue of the assumption that \(\lim _{r\rightarrow 0}l(r)=0\), we can choose \(\varepsilon \in (0,\varepsilon _0)\) so that

$$\begin{aligned} r_0 := \left[ \frac{2C}{\nu } + M_{\alpha _1}+M_{\alpha _2} +\sum _{i,j=1}^2 |a_{ij}| \left( M_{l(\alpha )}+\frac{C}{\nu }\right) \right] l_{{\hat{f}}}(\varepsilon )<1. \end{aligned}$$

Take

$$\begin{aligned} \delta = \frac{\varepsilon (1-r)}{2\Vert {\mathcal {R}}^0\Vert _{w,1} +(|a_{11}|+|a_{12}|)\Vert {\mathcal {R}}^{\alpha _1}\Vert _{w,1} +(|a_{21}|+|a_{22}|)\Vert {\mathcal {R}}^{\alpha _2}\Vert _{w,1} }, \end{aligned}$$

then for any initial condition \(x^0\in B(0,\delta )\), we have

$$\begin{aligned} \Vert {\mathcal {T}}_{x^0}\xi \Vert _w\le \varepsilon ,\;\forall \xi \in B_{C_w}(0,\varepsilon ), \end{aligned}$$

that is, \({\mathcal {T}}_{x^0}(B_{C_w}(0,\varepsilon ))\subset B_{C_w}(0,\varepsilon )\). Moreover, for every \(\xi ,{\hat{\xi }}\in B_{C_w}(0,\varepsilon )\),

$$\begin{aligned} \Vert {\mathcal {T}}_{x^0}\xi -{\mathcal {T}}_{x^0}{\hat{\xi }}\Vert _w&\le \left[ \frac{2C}{\nu } \! + \! M_{\alpha _1} \! + \! M_{\alpha _2} \! + \! \sum _{i,j=1}^2 |a_{ij}| \left( M_{l(\alpha )} \! + \! \frac{C}{\nu } \right) \right] l_{{\hat{f}}}(\varepsilon )\Vert \xi -{\hat{\xi }}\Vert _w\\&=r_0 \Vert \xi -{\hat{\xi }}\Vert _w. \end{aligned}$$

Thus, the operator \({\mathcal {T}}_{x^0}\) is contractive on \(B_{C_w}(0,\varepsilon )\), and by Banach’s fixed point theorem, \({\mathcal {T}}_{x^0}\) has a unique fixed point \(\xi ^*\) in this set. Furthermore, this function is the unique solution to the system (5.4) in \(B_{C_w}(0,\varepsilon )\). Notice that if \(\xi ^*\in B_{C_w}(0,\varepsilon )\) then \(f(\xi ^*(t))=\hat{f}(\xi ^*(t))\) for every \(t\in [0,\infty )\), and thus \(\xi ^*\) is also a solution to the system (5.2). This completes the proof. \(\square \)

6 Numerical examples

To complete this paper, we now give some numerical examples to illustrate the main theoretical results. Specifically, Lemmas 2, 4, 5, 6(i), 6(ii), 7(i) and 7(ii) provide different sufficient criteria for the zeros of the characteristic function Q to be located in the desired part of the complex plane. To demonstrate the applicability of these criteria, we have included an example for each of them except for Lemma 2 because this result has already been known for some time now [5].

In all the examples below, we use the functions \(f_1\) and \(f_2\) with

$$\begin{aligned} f_i(t) = {\left\{ \begin{array}{ll} 1&{} \text { if } 0 \le t < 1, \\ t^{-2i}&{} \text { if } t \ge 1 \end{array}\right. } \quad (i = 1, 2). \end{aligned}$$

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\). For each of the examples, we have provided two plots. The plots on the left always show the two components of the respective solution on their own. For these plots, the values of t range between 0 (the initial point) and 100. This interval is sufficiently large to provide a rough impression of the asymptotic behaviour of the functions and yet sufficiently small to still allow a reasonable view of their behaviour in the initial phase. To complement this, we have added (on the right of each figure) a plot of \(t^\beta x_j(t)\), \(j = 1, 2\), where the value of \(\beta \) is chosen such that, according to the theoretical considerations of Theorems 4 or 5 (whichever is applicable in the example under consideration), \(t^\beta x_j(t)\) is bounded as \(t \rightarrow \infty \) for both values of j. The range of t was chosen larger in this case (specifically, the displayed interval is [0, 300]) to clearly exhibit the boundedness (or lack thereof in the case of Fig. 3 where the conditions of Theorem 5 are not fulfilled).

Example 1

Consider the inhomogeneous two-component incommensurate fractional-order linear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{1/3}_{0^+}x_1(t) = 0.25x_2(t) + f_1(t), \\ ^C D^{1/2}_{0^+}x_2(t) = -2x_1(t) + x_2(t)+ f_2(t), \end{array}\right. } \quad t>0. \end{aligned}$$
(6.1)

In this example, the characteristic function is \(Q(s) = s^{5/6} - s^{1/3} + 0.5.\) According to Lemma 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 4. Hence, every solution to (6.1) tends to the origin as \(t\rightarrow \infty \) with the rate \(O(t^{-1/3})\). This property is illustrated in Fig. 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 \rightarrow \infty \) and \(j = 1, 2\), thus demonstrating that the decay behaviour of \(x_j(t)\) is indeed \(O(t^{-1/3})\).

Fig. 1
figure 1

Solution to the differential equation (6.1) from Example 1 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^{1/3} x_j(t)\)

Full size image

Example 2

Consider the two-component incommensurate fractional-order nonlinear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{1/3}_{0^+}x_1(t) = 0.25x_2(t) + x_1^2(t)x_2^2(t), \\ ^C D^{1/2}_{0^+}x_2(t) = -2x_1(t) + x_2(t)+ x_1^2(t) + x_2^2(t), \end{array}\right. } \quad t>0. \end{aligned}$$
(6.2)

It is not difficult to check that all conditions of Lemma 4 and Theorem 5 are verified. Thus, the trivial solution to (6.2) is Mittag-Leffler stable; more precisely, by Theorem 5, we have to expect an \(O(t^{-1/3})\) decay behaviour for nontrivial solutions with initial values sufficiently close to those of the trivial solution.

Definition 2 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 1) but only for initial values sufficiently close to (0, 0). Indeed we can see this behaviour in Fig. 2 for the initial value \((0.1, -0.2)\), whereas Fig. 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.

Fig. 2
figure 2

Solution to the differential equation (6.2) from Example 2 with initial conditions \(x_1(0) = 0.1\) and \(x_2(0) = -0.2\). The left graph shows the components \(x_j(t)\) of the solutions themselves, the right graph shows the functions \(t^{1/3} x_j(t)\)

Full size image
Fig. 3
figure 3

Solution to the differential equation (6.2) from Example 2 with initial conditions \(x_1(0) = 1\) and \(x_2(0) = -1\). The left graph shows the components \(x_j(t)\) of the solutions themselves, the right graph shows the functions \(t^{1/3} x_j(t)\)

Full size image

Example 3

Consider the fractional linear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{0.6}_{0^+}x_1(t) =x_1(t) + 2x_2(t) + f_1(t), \\ ^C D^{0.8}_{0^+}x_2(t) = -x_1(t) + f_2(t), \end{array}\right. } \quad t>0. \end{aligned}$$
(6.3)

The characteristic function of the system is \(Q(s) = s^{1.4} - s^{0.8} + 2.\) By Lemma 5, all zeros of Q lie in the open left-half of the complex plane and the assumptions of Theorem 4 are satisfied. Hence, every solution to this system converges to the origin as \(t\rightarrow \infty \) with an \(O(t^{-0.6})\) convergence rate. As in Example 1, we can also reproduce this behaviour numerically. The corresponding graphs are plotted in Fig. 4.

Fig. 4
figure 4

Solution to the differential equation (6.3) from Example 3 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.6} x_j(t)\)

Full size image

Example 4

Consider the system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{0.6}_{0^+}x_1(t) =x_1(t) + 2x_2(t) + x_1^2(t)x_2^2(t),\\ ^C D^{0.8}_{0^+}x_2(t) = -x_1(t) + x_1^2(t) + x_2^2(t), \end{array}\right. } \quad t > 0. \end{aligned}$$
(6.4)

Based on Lemma 5 and Theorem 5, we see that the trivial solution of (6.4) is Mittag-Leffler stable. As in Example 2, this is exhibited—together with the decay behaviour predicted by Theorem 5—in Fig. 5.

Fig. 5
figure 5

Solution to the differential equation (6.4) from Example 4 with initial conditions \(x_1(0) = 0.1\) and \(x_2(0) = -0.2\). The left graph shows the components \(x_j(t)\) of the solutions themselves, the right graph shows the functions \(t^{0.6} x_j(t)\)

Full size image

Example 5

Consider the inhomogeneous two-component incommensurate fractional-order linear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{0.3}_{0^+}x_1(t) = x_1(t) - x_2(t) + f_1(t) , \\ ^C D^{0.4}_{0^+}x_2(t) = 2x_1(t) +x_2(t) + f_2(t), \end{array}\right. } \quad t > 0. \end{aligned}$$
(6.5)

The system (6.5) has the characteristic function \(Q(s) = s^{0.7} - s^{0.4} - s^{0.3} + 3\). From Lemma 6 (i) and Theorem 4, it follows that every solution of this system tends to the origin as \(t\rightarrow \infty \) as \(O(t^{-0.3})\). Once again, our numerical results, shown in Fig. 6, support this statement.

Fig. 6
figure 6

Solution to the differential equation (6.5) from Example 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)\)

Full size image

Example 6

Consider the two-component incommensurate fractional-order nonlinear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{0.3}_{0^+}x_1(t) =0.1x_1(t) - 0.4x_2(t) + x_1^2(t)x_2^2(t),\\ ^C D^{0.4}_{0^+}x_2(t) = 0.7x_1(t) +0.2x_2(t) + x_1^2(t) + x_2^2(t), \end{array}\right. } \quad t > 0. \end{aligned}$$
(6.6)

Its characteristic function is \(Q(s) = s^{0.7} -0.1 s^{0.4} - 0.2 s^{0.3} + 0.3.\) It follows from Lemma 6(ii) and Theorem 5 that the trivial solution is Mittag-Leffler stable. Once again, we can visualize this observation on the basis of numerical results, cf. Fig. 7.

Fig. 7
figure 7

Solution to the differential equation (6.6) from Example 6 with initial conditions \(x_1(0) = 0.1\) and \(x_2(0) = -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)\)

Full size image

Example 7

Consider the two-component incommensurate fractional-order linear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{0.4}_{0^+}x_1(t) =-x_1(t) + 2x_2(t) + f_1(t),\\ ^C D^{0.5}_{0^+}x_2(t) = -5x_1(t) +4x_2(t)+ f_2(t), \end{array}\right. } \quad t>0. \end{aligned}$$
(6.7)

The system (6.7) has the characteristic function \(Q(s) = s^{0.9} + s^{0.5} - 4 s^{0.4} + 6.\) According to Lemma 7(i) and Theorem 4, its solution converges to the origin with a rate \(O(t^{-0.4})\). As above, the numerical data shown in Fig. 8 confirms this theoretical observation.

Fig. 8
figure 8

Solution to the differential equation (6.7) from Example 7 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.4} x_j(t)\)

Full size image

Example 8

Consider the two-component incommensurate fractional-order nonlinear system

$$\begin{aligned} {\left\{ \begin{array}{ll} ^C D^{0.4}_{0^+}x_1(t) =-x_1(t) - 2x_2(t) + x_1^2(t)x_2^2(t), \\ ^C D^{0.5}_{0^+}x_2(t) = 2x_1(t) +2x_2(t) + x_1^2(t) + x_2^2(t), \end{array}\right. } \quad t > 0. \end{aligned}$$
(6.8)

Its characteristic function \(Q(s) = s^{0.9} + s^{0.5} - 2s^{0.4} + 2\). According to Lemma 7(ii) and Theorem 5, the trivial solution of (6.8) is Mittag-Leffler stable as illustrated graphically in Fig. 9.

Fig. 9
figure 9

Solution to the differential equation (6.8) from Example 8 with initial conditions \(x_1(0) = 0.1\) and \(x_2(0) = -0.2\). The left graph shows the components \(x_j(t)\) of the solutions themselves, the right graph shows the functions \(t^{0.4} x_j(t)\)

Full size image

7 Conclusion

In this paper, we have provided new insight into the behaviour of non-commensurate fractional order planar systems. The main new contributions are Theorem 4, Corollary 2, Theorem 5 and Corollary 3. In Theorem 4 and Corollary 2, we show sufficient conditions for the global attractivity of non-trivial solutions to fractional-order inhomogeneous linear planar systems. In Theorem 5 and Corollary 3, we obtain 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 carefully chosen weighted norms combined with the Banach fixed point theorem for appropriate Banach spaces, we get the desired conclusions. We also provide numerical examples to illustrate the effect of the main theoretical results.