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

Kai Diethelm Faculty of Applied Natural Sciences and Humanities (FANG),
University of Applied Sciences Würzburg-Schweinfurt, Ignaz-Schön-Str. 11,
97421 Schweinfurt, Germany<EMAIL_ADDRESS>Ha Duc Thai111Institute of
Mathematics, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet,
10307 Ha Noi, Vietnam<EMAIL_ADDRESS>Hoang The Tuan222Institute of
Mathematics, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet,
10307 Ha Noi, Vietnam<EMAIL_ADDRESS>
###### Abstract

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.

2020 Mathematics Subject Classification: 34A08, 35A01, 35B20, 35B40, 60H15,
35R60

Key words:  non-commensurate fractional order planar systems; asymptotic
behavior of solutions; global attractivity; Mittag-Leffler stability

##  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 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

$\displaystyle^{C}D^{\alpha}_{0^{+}}x(t)$ $\displaystyle=Ax(t)+f(t,x(t)),\quad
t>0,$ (1) $\displaystyle x(0)$ $\displaystyle=x^{0}\in\mathbb{R}^{2},$ (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 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 $\mathbb{R}^{2}$, we define the norm
$\|\cdot\|$ by $\|x\|:=\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}:\|x\|\leq 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 $\|\xi\|_{\infty}:=\sup_{t\geq 0}\|\xi(t)\|$. Then, we use the notation
$C_{\infty}([0,\infty);\mathbb{R}^{2}):=\\{\xi\in
C([0,\infty);\mathbb{R}^{2}):\|\xi\|_{\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

$I^{\alpha}_{0^{+}}f(t):=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}f(s)ds,\;t\in
J,$

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

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

where $\Gamma(\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

${}^{C}D^{\alpha}_{0^{+}}f(t):=\left({}^{C}D^{\alpha_{1}}_{0^{+}}f_{1}(t),^{C}D^{\alpha_{2}}_{0^{+}}f_{2}(t)\right).$

See, e.g., [9, Chapter III] and [25] 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

${}^{C}D^{\alpha}_{0^{+}}x(t)$ $\displaystyle=f(t,x(t)),\quad t>0,$ (3a)
$\displaystyle x(0)$ $\displaystyle=x^{0}\in\mathbb{R}^{2},$ (3b)

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.

######  Theorem 2.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

$\|f(t,x)-f(t,\hat{x})\|\leq L||x-\hat{x}||,\;\forall
t\in[0.\infty),\;x,\hat{x}\in\mathbb{R}^{2}$

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

###### Proof.

See [22, Theorem 2.2 and Remark 2.3]. ∎

######  Theorem 2.2 (Exponential boundedness of global solutions).

Suppose that the function $f$ satisfies the assumptions of Theorem 2.1.
Moreover, let there exist a constant $\gamma>0$ such that

$\sup_{t\geq 0}e^{-\gamma
t}\int_{0}^{t}(t-s)^{\alpha_{i}-1}\|f(s,0)\|ds<\infty.$

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

$\|\varphi(t,x^{0})\|\leq Me^{\gamma t},\;\forall t\geq 0,$

where $M$ is some positive constant which depends on $x^{0}.$

###### Proof.

See [22, Theorem 2.4]. ∎

###  2.2 The variation of constants formula for the solutions

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

$\displaystyle^{C}D^{\alpha}_{0^{+}}x(t)=Ax(t)+f(t),t>0$ (4a) with initial
condition $\displaystyle x(0)=x^{0}\in\mathbb{R}^{2},$ (4b)

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})^{\rm T}:[0,\infty)\rightarrow\mathbb{R}^{2}$ is a continuous
function such that

$\|f(t)\|\leq Me^{\gamma t},\;\forall t\geq 0$ (5)

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

$\displaystyle\int_{0}^{t}(t-s)^{\alpha_{i}-1}|f_{i}(s)|ds$ $\displaystyle\leq
M\int_{0}^{t}(t-s)^{\alpha_{i}-1}e^{\gamma s}ds$
$\displaystyle=-\frac{Me^{\gamma
t}}{\gamma^{\alpha_{i}}}\int_{0}^{t}\left(\gamma(t-s)\right)^{\alpha_{i}-1}e^{-\gamma(t-s)}d(\gamma(t-s))$
$\displaystyle=\frac{Me^{\gamma t}}{\gamma^{\alpha_{i}}}\int_{0}^{\gamma
t}\tau^{\alpha_{i}-1}e^{-\tau}d\tau$
$\displaystyle\leq\frac{M\Gamma(\alpha_{i})}{\gamma^{\alpha_{i}}}e^{\gamma
t}.$

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

$\displaystyle\begin{cases}(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{cases},$
(6)

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

$\displaystyle X_{1}(s)$
$\displaystyle=\frac{x_{1}^{0}(s^{\alpha_{1}+\alpha_{2}-1}-a_{22}s^{\alpha_{1}-1})+x^{0}_{2}a_{12}s^{\alpha_{2}-1}+F_{1}(s)(s^{\alpha_{2}}-a_{22})+a_{12}F_{2}(s)}{Q(s)}$
$\displaystyle=\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)},$
(7)

and

$\displaystyle X_{2}(s)$
$\displaystyle=\frac{x_{2}^{0}(s^{\alpha_{1}+\alpha_{2}-1}-a_{11}s^{\alpha_{2}-1})+x^{0}_{1}a_{21}s^{\alpha_{1}-1}+a_{12}F_{1}(s)+F_{2}(s)(s^{\alpha_{1}}-a_{11})}{Q(s)}$
$\displaystyle=\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),$
(8)

where
$Q(s):=s^{\alpha_{1}+\alpha_{2}}-a_{11}s^{\alpha_{2}}-a_{22}s^{\alpha_{1}}+\det
A.$ Put

$\displaystyle\mathcal{R}^{\lambda}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{s^{l(\alpha)-\lambda}}{sQ(s)}\right\\}(t),\quad\lambda\in\left\\{0,\alpha_{1},\alpha_{2}\right\\},$
(9a) $\displaystyle\mathcal{S}^{\beta}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{s^{l(\alpha)-\beta}}{Q(s)}\right\\}(t),\quad\beta\in\left\\{\alpha_{1},\alpha_{2},l(\alpha)\right\\}$
(9b)

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

$\displaystyle\mathcal{L}^{-1}\left\\{\frac{s^{l(\alpha)-\beta}}{Q(s)}F_{i}(s)\right\\}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\mathcal{L}\left\\{\mathcal{S}^{\beta}\right\\}(s)\mathcal{L}\left\\{f_{i}\right\\}(s)\right\\}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\mathcal{L}\left\\{\mathcal{S}^{\beta}\ast
f_{i}\right\\}(s)\right\\}(t)$ $\displaystyle=\mathcal{S}^{\beta}\ast
f_{i}(t),\beta\in\left\\{\alpha_{1},\alpha_{2},l(\alpha)\right\\},$

where ”$\ast$” is the Laplace convolution operator.

From the arguments above, the unique solution to the initial value problem (4)
has the following form.

######  Lemma 2.3.

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

$\varphi(\cdot,x^{0})=\begin{pmatrix}\varphi_{1}(\cdot,x^{0})\\\
\varphi_{2}(\cdot,x^{0})\end{pmatrix}$

with

$\displaystyle\varphi_{1}(t,x^{0})$
$\displaystyle=\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}$
$\displaystyle\hskip
28.45274pt+\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),$
(10) $\displaystyle\varphi_{2}(t,x^{0})$
$\displaystyle=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}$
$\displaystyle\hskip
28.45274pt+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).$
(11)

##  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 $\alpha_{1}\neq\alpha_{2}$, because the case
$\alpha_{1}=\alpha_{2}$ has already been discussed in detail elsewhere [13].
Thus, without loss of generality, we assume $0<\alpha_{1}<\alpha_{2}\leq 1$.
Our first auxiliary statement in this context deals with functions of the form

$Q(s)=s^{\alpha_{1}+\alpha_{2}}-as^{\alpha_{2}}-bs^{\alpha_{1}}+c;$ (12)

the characteristic functions of the problems under consideration will be of
precisely this structure.

######  Lemma 3.1.

Let $0<\alpha_{1}<\alpha_{2}\leq 1$ and $a,b,c\in\mathbb{R}$. Then, the
following statements hold for the function $Q$ defined in (12).

  * (i)

If $c<0$ then $Q$ has at least one positive real zero.

  * (ii)

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

  * (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)}|\leq\omega\\}$.

  * (iv)

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

$\displaystyle\begin{cases}&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{cases}$ (13)

where

$\rho_{1}=\frac{\sin\frac{\alpha_{1}\pi}{2}}{\sin\frac{(\alpha_{2}-\alpha_{1})\pi}{2}},\hskip
56.9055pt\rho_{2}=\frac{\sin\frac{\alpha_{2}\pi}{2}}{\sin\frac{(\alpha_{2}-\alpha_{1})\pi}{2}}.$
(14)

###### Proof.

(i) and (ii) are obvious.

(iii) First, we assume that $c\neq 0$. Then $Q(0)\neq 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)|\geq|s|^{\alpha_{1}+\alpha_{2}}-|a|\cdot|s|^{\alpha_{2}}-|b|\cdot|s|^{\alpha_{1}}-|c|,$
we have that $\lim_{|s|\to\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)}|\leq\omega\\}$ (if they exist) belong
to the set $\Omega:=\left\\{z\in\mathbb{C}:\varepsilon\leq|z|\leq
R,|\arg{(z)}|\leq\omega\right\\}.$ Notice that $\Omega$ is a compact set and
$Q$ is analytic on this domain. If now $Q$ has infinitely many zeros in
$\Omega$ then, because of the compactness of $\Omega$, 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 $\Omega$. This shows that $Q$ has only a finite
number of zeros in the domain $\mathcal{C}$ if $c\neq 0$.
To deal with the case $c=0$, we write

$\displaystyle Q(s)$
$\displaystyle=s^{\alpha_{1}}\left(s^{\alpha_{2}}-as^{\alpha_{2}-\alpha_{1}}-b\right)$
$\displaystyle=s^{\alpha_{1}}P(s),$

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]. ∎

######  Corollary 3.2.

Assume that $a,b,c>0$ and that one of conditions

  * (i)

$c(\rho^{2}_{2}-\rho_{1}^{2})<ab<c(\rho_{2}^{2}+\rho_{1}^{2})$,

  * (ii)

$ab\leq c\left(\rho_{2}-\rho_{1}\right)^{2}$,

is satisfied where $\rho_{1},\rho_{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 $\rho_{2}\neq 0$, this system
is equivalent to

$\displaystyle\left\\{\begin{array}[]{rl}\omega^{\alpha_{1}}&=\displaystyle\frac{a}{\rho_{2}}+c\frac{\rho_{1}}{\rho_{2}}\omega^{-\alpha_{2}},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\rho_{1}a\omega^{\alpha_{2}}+\rho_{2}b\omega^{\alpha_{1}}&=c(\rho_{2}^{2}-\rho_{1}^{2}).\end{array}\right.$
(17)

Thus, we obtain

$a\rho_{1}\omega^{\alpha_{2}}+bc\rho_{1}\omega^{-\alpha_{2}}+ab-c(\rho_{2}^{2}-\rho_{1}^{2})=0.$
(18)

Setting $X=\omega^{\alpha_{2}}$, equation (18) takes the form

$\displaystyle
a\rho_{1}X^{2}+\left[ab-c(\rho_{2}^{2}-\rho_{1}^{2})\right]X+bc\rho_{1}=0.$
(19)

The discriminant of the quadratic equation (19) is

$\displaystyle\Delta$
$\displaystyle=\left(ab-c(\rho_{2}^{2}-\rho_{1}^{2})\right)^{2}-4abc\rho_{1}^{2}$
$\displaystyle=a^{2}b^{2}+c^{2}(\rho_{2}^{2}-\rho_{1}^{2})^{2}-2abc(\rho_{1}^{2}+\rho_{2}^{2})$
$\displaystyle=\left(ab-c(\rho_{1}^{2}+\rho_{2}^{2})\right)^{2}-4c^{2}\rho_{1}^{2}\rho_{2}^{2}$
$\displaystyle=\left(ab-c(\rho_{1}+\rho_{2})^{2}\right)\left(ab-c(\rho_{2}-\rho_{1})^{2}\right).$

(i) Clearly, if
$c(\rho^{2}_{2}-\rho_{1}^{2})<ab<c(\rho_{2}^{2}+\rho_{1}^{2}),$ then
$\Delta<0.$ Hence, the quadratic equation (19) has no real roots. This implies
that the system (13) has no root $\omega>0$. This together with Lemma 3.1(ii)
and Lemma 3.1(iv) shows that $Q$ has no purely imaginary zero.

(ii) If $ab\leq c\left(\rho_{2}-\rho_{1}\right)^{2}$, the quadratic equation
(19) 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 (19) are negative. Hence, in view of the
relation $X=\omega^{\alpha_{2}}$ with $0<\alpha_{2}\leq 1$ between the
solution $X$ of (19) and the solution $\omega$ of (13), the system (13) has no
root $\omega>0$. Using Lemma 3.1(ii) and Lemma 3.1(iv), we see that $Q$ has no
purely imaginary zero. ∎

Recall that if $c\leq 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|\geq
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)\neq 0$ in
$\left\\{z\in\mathbb{C}:|z|\leq\varepsilon\right\\}$. We define an oriented
contour $\gamma$ formed by four segments:

$\displaystyle\gamma_{1}$
$\displaystyle:=\left\\{s=i\omega:\varepsilon\leq\omega\leq R\right\\};$
$\displaystyle\gamma_{2}$
$\displaystyle:=\left\\{s=Re^{i\varphi}:-\frac{\pi}{2}\leq\varphi\leq\frac{\pi}{2}\right\\};$
$\displaystyle\gamma_{3}$ $\displaystyle:=\left\\{s=\varepsilon
e^{i\varphi}:-\frac{\pi}{2}\leq\varphi\leq\frac{\pi}{2}\right\\};$
$\displaystyle\gamma_{4}$
$\displaystyle:=\left\\{s=i\omega:-R\leq\omega\leq-\varepsilon\right\\}.$

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\geq 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|\to\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

$\displaystyle h_{1}(\omega)$
$\displaystyle:=\Re(Q(i\omega))=\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,$
(20) $\displaystyle h_{2}(\omega)$
$\displaystyle:=\Im(Q(i\omega))=\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}.$
(21)

If there exists some $\omega>0$ such that $h_{2}(\omega)=0,$ then

$\displaystyle\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}h_{1}(\omega)$
$\displaystyle=\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}h_{1}(\omega)-\cos\frac{(\alpha_{1}+\alpha_{2})\pi}{2}h_{2}(\omega)$
$\displaystyle=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)$
$\displaystyle{}\quad+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)$
$\displaystyle{}\quad+c\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}$
$\displaystyle=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}.$

Thus, the variable $\omega>0$ satisfies the system

$\displaystyle\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.$
(24)

with

$\displaystyle
q_{1}=\frac{\sin\frac{\alpha_{1}\pi}{2}}{\sin\frac{(\alpha_{2}+\alpha_{1})\pi}{2}},\hskip
56.9055ptq_{2}=\frac{\sin\frac{\alpha_{2}\pi}{2}}{\sin\frac{(\alpha_{2}+\alpha_{1})\pi}{2}}.$
(25)

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

Assume that $a,b\leq 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)]. ∎

######  Lemma 3.4.

Let $0<\alpha_{1}<\alpha_{2}\leq 1$. Assume that $a=0$, $b>0$ and

$c>\left(bq_{1}\right)^{\alpha_{1}/\alpha_{2}}bq_{2}.$ (26)

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 (24), we have
$h_{1}(\omega_{0})=c-\left(bq_{1}\right)^{\alpha_{1}/\alpha_{2}}bq_{2}$. By
the assumption (26), we obtain $h_{1}(\omega_{0})>0$. This implies that all
zeros of $Q$ lie in the open left-half of complex plane. ∎

######  Lemma 3.5.

Let $0<\alpha_{1}<\alpha_{2}\leq 1.$ Assume that $b=0$, $a>0$ and

$c>\left(aq_{2}\right)^{\alpha_{2}/\alpha_{1}}aq_{1}.$ (27)

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 (24), it follows
that $h_{1}(\omega_{0})=c-\left(aq_{2}\right)^{\alpha_{2}/\alpha_{1}}aq_{1}$.
By the assumption (27), 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. ∎

######  Lemma 3.6.

Let $0<\alpha_{1}<\alpha_{2}\leq 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:

  * (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}\leq c.$

  * (ii)

$aq_{2}+bq_{1}\leq 1$ and $aq_{1}+bq_{2}<c.$

###### Proof.

We have

$\displaystyle h_{2}^{\prime}(\omega)$
$\displaystyle=(\alpha_{1}+\alpha_{2})\omega^{\alpha_{1}+\alpha_{2}-1}\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}-a\alpha_{2}\omega^{\alpha_{2}-1}\sin\frac{\alpha_{2}\pi}{2}-\alpha_{1}\omega^{\alpha_{1}-1}\sin\frac{\alpha_{1}\pi}{2}$
$\displaystyle=\omega^{\alpha_{1}-1}\left((\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}\right)$
$\displaystyle=\omega^{\alpha_{1}-1}g_{2}(\omega)$ (28)

where

$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}.$
(29)

Notice that

$\displaystyle g_{2}^{\prime}(\omega)$
$\displaystyle=\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}.$

It is not difficult to check that $g_{2}^{\prime}(\omega)<0$ in
$(0,\omega_{1})$ and $g_{2}^{\prime}(\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\to+\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\to+\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}\leq 1$, we have

$\displaystyle
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}$
(30)

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

$\displaystyle c\geq
aq_{2}\left((a+b)q_{2}\right)^{\alpha_{2}/\alpha_{1}}+b(a+b)q_{2}^{2},$

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 (24) leads to $h_{1}(\omega_{3})>0.$ The proof of this part is
complete.

(ii) If $aq_{2}+bq_{1}\leq 1$, then $h_{2}(1)\geq 0.$ This implies that
$0<\omega_{3}\leq 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 (24) shows that $h_{1}(\omega_{3})>0.$ The proof is finished. ∎

######  Lemma 3.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:

  * (i)

$aq_{2}+bq_{1}>1$ and $\left(bq_{1}\right)^{\alpha_{1}/\alpha_{2}}bq_{2}\leq
c.$

  * (ii)

$aq_{2}+bq_{1}\leq 1$ and $bq_{2}\leq c$.

###### Proof.

As shown in the proof of Lemma 3.6, we have
$h_{2}^{\prime}(\omega)=\omega^{\alpha_{1}-1}g_{2}(\omega)$ where $g_{2}$ is
as in (29). Notice that

$g_{2}^{\prime}(\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$

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\to+\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\to+\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

$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},$

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

$\displaystyle\left(bq_{1}\right)^{\alpha_{1}/\alpha_{2}}bq_{2}\leq c,$

we obtain $c>a\omega_{2}^{\alpha_{2}}q_{1}+b\omega_{2}^{\alpha_{1}}q_{2},$
which together with (24) leads to $h_{1}(\omega_{2})>0.$

(ii) If $aq_{2}+bq_{1}\leq 1$, then $h_{2}(1)\geq 0.$ Thus $0<\omega_{2}\leq
1.$ Due to $c\geq bq_{2}$, we see that
$c>a\omega_{2}^{\alpha_{2}}q_{1}+b\omega_{2}^{\alpha_{1}}q_{2}$. The proof is
completed. ∎

##  4 Estimates for the functions $\mathcal{R}^{\lambda}$ and
$\mathcal{S}^{\beta}$

This section is devoted to the study of some important estimates of the
functions $\mathcal{R}^{\lambda}$ and $\mathcal{S}^{\beta}$ on $(0,\infty)$.
We recall their definitions from (9), viz.

$\displaystyle\mathcal{R}^{\lambda}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{s^{l(\alpha)-\lambda}}{sQ(s)}\right\\}(t),\quad\lambda\in\left\\{0,\alpha_{1},\alpha_{2}\right\\},$
$\displaystyle\mathcal{S}^{\beta}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{s^{l(\alpha)-\beta}}{Q(s)}\right\\}(t),\quad\beta\in\left\\{\alpha_{1},\alpha_{2},l(\alpha)\right\\},$

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

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\\}$:

$\displaystyle\mathcal{R}^{\lambda}(t)$ $\displaystyle=O(t^{-\nu})\quad\mbox{
as }t\to\infty,$ (31) $\displaystyle\mathcal{S}^{\beta}(t)$
$\displaystyle=O(t^{-\nu-1})\quad\mbox{ as }t\to\infty,$ (32)
$\displaystyle\mathcal{S}^{\beta}(t)$ $\displaystyle=O(t^{\nu-1})\quad\mbox{
as }t\to 0.$ (33)

Moreover,

$\int_{0}^{\infty}|\mathcal{S}^{\beta}(t)|dt<\infty.$ (34)

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
$\mathcal{S}^{\beta}$ and a continuous function.

######  Theorem 4.2.

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

$\displaystyle F_{g}^{\beta}(t):=\mathcal{S}^{\beta}\ast
g(t)=\int_{0}^{t}\mathcal{S}^{\beta}(t-s)g(s)ds$ (35)

and

$\displaystyle\bar{F}_{g}^{\beta}(t):=|\mathcal{S}^{\beta}|\ast|g|(t)=\int_{0}^{t}\mathcal{|}S^{\beta}(t-s)|\cdot|g(s)|ds.$
(36)

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

  * (i)

If $g$ is bounded then $\bar{F}_{g}^{\beta}$ and $F_{g}^{\beta}$ are also
bounded.

  * (ii)

If $\lim_{t\to\infty}g(t)=0$ then
$\lim_{t\to\infty}\bar{F}_{g}^{\beta}(t)=\lim_{t\to\infty}F_{g}^{\beta}(t)=0$.

  * (iii)

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

$\bar{F}_{g}^{\beta}(t)=O(t^{-\mu})\mbox{ and
}F_{g}^{\beta}(t)=O(t^{-\mu})\mbox{ for }t\to\infty$

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)|\leq\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)\geq 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

$\bar{F}_{g}^{\beta}(t)\leq\varepsilon\quad\mbox{ for all }t>\widetilde{T}.$
(37)

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

Our first observation is then that, from (32) and (33), we know that there
exists some constant $C>0$ such that

$|S^{\beta}(t)|\leq\begin{cases}Ct^{-\nu-1}&\mbox{ for }t\geq 1,\\\
Ct^{\nu-1}&\mbox{ for }t\leq 1.\end{cases}$ (38)

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/(3C)$ for all
$t>\hat{T}$. Using these values $\hat{T}$ and $C$, we then define

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

For $t>\widetilde{T}\geq\hat{T}+1$, we can then write

$\displaystyle\bar{F}_{g}^{\beta}(t)$
$\displaystyle=\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+\\!\int_{t-1}^{t}\\!|\mathcal{S}^{\beta}(t-s)|\cdot|g(s)|ds$
$\displaystyle=F_{1}(t)+F_{2}(t)+F_{3}(t).$ (39)

Our goal now is to show that, under these assumptions,
$F_{j}(t)\leq\varepsilon/3$ for $j=1,2,3$, which implies (37) and thus
suffices to prove part (ii) of the Theorem. In this context, we see that

$\displaystyle F_{1}(t)$
$\displaystyle=\int_{0}^{\hat{T}}|\mathcal{S}^{\beta}(t-s)|\cdot|g(s)|ds\leq\|g\|_{\infty}C\int_{0}^{\hat{T}}(t-s)^{-\nu-1}ds$
$\displaystyle\leq\|g\|_{\infty}C\hat{T}(t-\hat{T})^{-\nu-1}<\frac{\varepsilon}{3}.$

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

$\displaystyle F_{2}(t)$
$\displaystyle=\int_{\hat{T}}^{t-1}|\mathcal{S}^{\beta}(t-s)|\cdot|g(s)|ds\leq\frac{\nu\varepsilon}{3C}C\int_{\hat{T}}^{t-1}(t-s)^{-\nu-1}ds$
$\displaystyle=\frac{\nu\varepsilon}{3}\frac{1}{\nu}\left(1-(t-\hat{T})^{-\nu}\right)<\frac{\varepsilon}{3}$

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

Finally,

$\displaystyle F_{3}(t)$
$\displaystyle=\int_{t-1}^{t}|\mathcal{S}^{\beta}(t-s)|\cdot|g(s)|ds\leq\frac{\nu\varepsilon}{3C}C\int_{t-1}^{t}(t-s)^{\nu-1}ds=\frac{\nu\varepsilon}{3}\frac{1}{\nu}=\frac{\varepsilon}{3}$

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

For the proof of (iii), we note that (38) 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\geq 2$. Then we write

$\displaystyle\bar{F}_{g}^{\beta}(t)$
$\displaystyle=\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$
$\displaystyle\phantom{=}{}\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$
$\displaystyle=F_{4}(t)+F_{5}(t)+F_{6}(t)+F_{7}(t),$ (40)

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,

$|g(t)|\leq C^{\prime}t^{-\eta}\quad\forall t\geq 1$ (41)

with some $C^{\prime}>0$, so that the upper branch of (38) implies

$\displaystyle 0\leq F_{4}(t)$ $\displaystyle\leq
C\|g\|_{\infty}\int_{0}^{1}(t-s)^{-\nu-1}ds=\frac{C\|g\|_{\infty}}{\nu}\left((t-1)^{-\nu}-t^{-\nu}\right)$
$\displaystyle<\frac{C\|g\|_{\infty}}{\nu}(t-1)^{-\nu}=O(t^{-\nu})=O(t^{-\mu})$

and

$\displaystyle 0\leq F_{5}(t)$ $\displaystyle\leq
CC^{\prime}\int_{1}^{t/2}(t-s)^{-\nu-1}s^{-\eta}ds\leq
CC^{\prime}\int_{1}^{t/2}(t-s)^{-\nu-1}ds$
$\displaystyle=\frac{CC^{\prime}}{\nu}\left((t-1)^{-\nu}-(t/2)^{-\nu}\right)<\frac{CC^{\prime}}{\nu}(t-1)^{-\nu}=O(t^{-\nu})=O(t^{-\mu})$

as well as

$\displaystyle 0\leq F_{6}(t)$ $\displaystyle\leq
CC^{\prime}\int_{t/2}^{t-1}(t-s)^{-\nu-1}s^{-\eta}ds\leq
CC^{\prime}\left(\frac{t}{2}\right)^{-\eta}\int_{t/2}^{t-1}(t-s)^{-\nu-1}ds$
$\displaystyle=2^{\eta}\frac{CC^{\prime}}{\nu}t^{-\eta}\left(1-(t/2)^{-\nu}\right)<2^{\eta}\frac{CC^{\prime}}{\nu}t^{-\eta}=O(t^{-\eta})=O(t^{-\mu}).$

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

$\displaystyle 0\leq F_{7}(t)$ $\displaystyle\leq
CC^{\prime}\int_{t-1}^{t}(t-s)^{\nu-1}s^{-\eta}ds\leq
CC^{\prime}(t-1)^{-\eta}\int_{t-1}^{t}(t-s)^{\nu-1}ds$
$\displaystyle=\frac{CC^{\prime}}{\nu}(t-1)^{-\eta}=O(t^{-\eta})=O(t^{-\mu}),$

thus completing the proof. ∎

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

$\bar{F}^{\beta}_{g_{1}}(t)=O(t^{-\nu})\quad\mbox{ as }t\to\infty\qquad\mbox{
for }g_{1}(t)=\min\\{1,t^{-\nu}\\}.$ (42)

Moreover, assuming $\nu<1$ and setting

$g_{2}(t)=t^{-\nu}-g_{1}(t)=\begin{cases}t^{-\nu}-1&\mbox{ for }t\in[0,1],\\\
0&\mbox{ for }t>1,\end{cases}$

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\geq 2$ then we have

$\displaystyle\bar{F}^{\beta}_{g_{2}}(t)$
$\displaystyle=\int_{0}^{1}|\mathcal{S}^{\beta}(t-s)|(s^{-\nu}-1)ds\leq
C\int_{0}^{1}(t-s)^{-\nu-1}s^{-\nu}ds=Ct^{-2\nu}B_{1/t}(1-\nu,-\nu)$
$\displaystyle=\frac{C}{1-\nu}t^{-\nu-1}{}{}_{2}F_{1}(1-\nu,1+\nu;2-\nu;t^{-1})\leq
C^{\prime}t^{-\nu-1}$ (43)

with some $C^{\prime}>0$.

  * •

If $t\in[1,2]$ then

$\displaystyle\bar{F}^{\beta}_{g_{2}}(t)$
$\displaystyle=\int_{0}^{1}|\mathcal{S}^{\beta}(t-s)|(s^{-\nu}-1)ds\leq
C\int_{0}^{1}(t-s)^{\nu-1}s^{-\nu}ds=CB_{1/t}(1-\nu,\nu)$
$\displaystyle=\frac{C}{1-\nu}t^{\nu-1}{}{}_{2}F_{1}(1-\nu,1-\nu;2-\nu;t^{-1})\leq
C^{\prime\prime}t^{-\nu}$ (44)

with some $C^{\prime\prime}>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. (42), (43) and (44) in the following
way:

######  Remark 4.3.

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

$\displaystyle
t^{\nu}\int_{0}^{t}|\mathcal{S}^{\beta}(t-s)|\frac{1}{s^{\nu}}ds\leq C.$ (45)

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

###  5.1 Asymptotic behaviour of solutions to fractional linear planar systems

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

$\displaystyle^{C}D^{\alpha}_{0^{+}}x(t)$ $\displaystyle=Ax(t)+f(t),\;t>0,$
(46a) $\displaystyle x(0)$ $\displaystyle=x^{0}\in\mathbb{R}^{2},$ (46b)

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

Suppose that all zeros of the characteristic function
$Q=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 following statements
hold.

  * (i)

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

  * (ii)

If $\lim_{t\to\infty}f(t)=0$ then the solution to (46) tends to $0$ when
$t\to\infty$ for any $x^{0}\in\mathbb{R}^{2}$.

  * (iii)

If $\|f(t)\|=O(t^{-\eta})$ as $t\to\infty$ with some $\eta>0$ then every
solution $x$ of (46a) behaves as $\|x(t)\|=O(t^{-\mu})$ for $t\to\infty$ where
$\mu=\min\left\\{\alpha_{1},\alpha_{2},\eta\right\\}.$

###### Proof.

The proof is straightforward by combining Lemma 2.3, Lemma 4.1 and Theorem
4.2. ∎

Based on Theorem 5.1 and Lemmas 3.3, 3.4, 3.5, 3.6 and 3.7, we obtain the
following corollary.

######  Corollary 5.2.

Let

$q_{1}=\frac{\sin\frac{\alpha_{1}\pi}{2}}{\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}}\quad\mbox{
and }\quad
q_{2}=\frac{\sin\frac{\alpha_{2}\pi}{2}}{\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}}.$

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

  * (i)

$a_{11},a_{22}\leq 0$ and $\det A>0$.

  * (ii)

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

$(a_{22}q_{1})^{\alpha_{1}/\alpha_{2}}a_{22}q_{2}<\det A.$

  * (iii)

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

$(a_{11}q_{2})^{\alpha_{2}/\alpha_{1}}a_{11}q_{1}<\det A.$

  * (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}\leq\det
A$;

    * (iv)2

$a_{11}q_{2}+a_{22}q_{1}\leq 1$ and $a_{11}q_{1}+a_{22}q_{2}<\det A.$

  * (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}\leq\det A$;

    * (v)2

$a_{11}q_{2}+a_{22}q_{1}\leq 1$ and $a_{22}q_{2}\leq\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

$\displaystyle^{C}D^{\alpha}_{0^{+}}x(t)$ $\displaystyle=Ax(t)+f(x(t)),\quad
t>0,$ (47a) $\displaystyle x(0)$
$\displaystyle=x^{0}\in\Omega\subset\mathbb{R}^{2},$ (47b)

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,
$\Omega$ is an open subset of $\mathbb{R}^{2}$ containing the origin and
$f:\Omega\rightarrow\mathbb{R}^{2}$ is locally Lipschitz continuous at the
origin such that $f(0)=0$ and $\lim_{r\to 0}l_{f}(r)=0$ with

$l_{f}(r):=\sup_{x,y\in B(0,r),\;x\neq y}\frac{\|f(x)-f(y)\|}{\|x-y\|}.$

######  Definition 5.3.

The trivial solution of (47a) 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 (47) exists globally on the interval $[0,\infty)$ and

$\displaystyle\max\\{\sup_{t\in[0,1]}\|\varphi(t,x^{0})\|,\sup_{t\geq
1}t^{\gamma}\|\varphi(t,x^{0})\|\\}\leq m.$

Our aim is to prove the following theorem.

######  Theorem 5.4.

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 (47a) is Mittag-Leffler stable. More precisely, there
exist constants $\delta,\varepsilon>0$ such that for any $\|x^{0}\|<\delta$,
the unique solution $\varphi(\cdot,x^{0})$ of the initial value problem (47)
exists globally on $[0,\infty)$ and $\sup_{t\geq
1}t^{\nu}\|\varphi(t,x^{0})\|\leq\varepsilon$ with
$\nu=\min\\{\alpha_{1},\alpha_{2}\\}.$

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

$q_{1}=\frac{\sin\frac{\alpha_{1}\pi}{2}}{\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}}\quad\mbox{
and }\quad
q_{2}=\frac{\sin\frac{\alpha_{2}\pi}{2}}{\sin\frac{(\alpha_{1}+\alpha_{2})\pi}{2}}.$

The statement of Theorem 5.4 is true if one of the following conditions is
satisfied.

  * (i)

$a_{11},a_{22}\leq 0$ and $\det A>0$.

  * (ii)

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

$(a_{22}q_{1})^{\alpha_{1}/\alpha_{2}}a_{22}q_{2}<\det A.$

  * (iii)

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

$(a_{11}q_{2})^{\alpha_{2}/\alpha_{1}}a_{11}q_{1}<\det A.$

  * (iv)

$a_{11},a_{22},detA>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}\leq\det
A$;

    * (iv)2

$a_{11}q_{2}+a_{22}q_{1}\leq 1$ and $a_{11}q_{1}+a_{22}q_{2}<\det A.$

  * (v)

$a_{11}<0$, $a_{22},detA>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}\leq\det A$;

    * (v)2

$a_{11}q_{2}+a_{22}q_{1}\leq 1$ and $a_{22}q_{2}\leq\det A$.

###### Proof of Theorem 5.4.

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

$\displaystyle^{C}D^{\alpha}_{0^{+}}x(t)$
$\displaystyle=Ax(t)+\hat{f}(x(t)),t>0,$ (48a) $\displaystyle x(0)$
$\displaystyle=x^{0}.$ (48b)

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}))^{\rm
T}$ on $[0,\infty)$ staifies the relationships

$\displaystyle\hat{\varphi}_{1}(\cdot,x^{0})$
$\displaystyle=\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}$
$\displaystyle\phantom{=}\quad+\left((\mathcal{S}^{\alpha_{1}}-a_{22}\mathcal{S}^{l(\alpha)})\ast\hat{f}_{1}(\hat{\varphi}(\cdot,x^{0}))\right)(t)+a_{12}\left(\mathcal{S}^{l(\alpha)}\ast\hat{f}_{2}(\hat{\varphi}(\cdot,x^{0}))\right)(t),$
(49a) $\displaystyle\hat{\varphi}_{2}(t,x^{0})$
$\displaystyle=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}$
$\displaystyle+a_{21}\left(\mathcal{S}^{l(\alpha)}\ast\hat{f}_{1}(\hat{\varphi}(\cdot,x^{0}))\right)(t)+\left((\mathcal{S}^{\alpha_{2}}-a_{11}\mathcal{S}^{l(\alpha)})\ast\hat{f}_{2}(\hat{\varphi}(\cdot,x^{0}))\right)(t).$
(49b)

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 (48) is contained in the space
$C_{\infty}([0,\infty);\mathbb{R}^{2})$ which is equipped the norm

$\|\xi\|_{w}:=\max\\{\sup_{t\in[0,1]}\|\xi(t)\|,\sup_{t\geq
1}t^{\nu}\|\xi(t)\|\\}.$

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

Based on the representation in (49), we establish a Lyapunov–Perron type
operator $\mathcal{T}_{x^{0}}$ on the space $C_{w}([0,\infty);\mathbb{R}^{2})$
as follows. For any $\xi\in C_{w}([0,\infty);\mathbb{R}^{2})$, let

$\displaystyle(\mathcal{T}_{x_{0}}\xi)_{1}(t)$
$\displaystyle:=\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}$
$\displaystyle\phantom{=}\quad+\left((\mathcal{S}^{\alpha_{1}}-a_{22}\mathcal{S}^{l(\alpha)})\ast\hat{f}_{1}(\xi(\cdot))\right)(t)+a_{12}\left(\mathcal{S}^{l(\alpha)}\ast\hat{f}_{2}(\xi(\cdot))\right)(t),$
$\displaystyle(\mathcal{T}_{x_{0}}\xi)_{2}(t)$
$\displaystyle:=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}$
$\displaystyle\phantom{=}\quad+a_{21}\left(\mathcal{S}^{l(\alpha)}\ast\hat{f}_{1}(\xi(\cdot))\right)(t)+\left((\mathcal{S}^{\alpha_{2}}-a_{11}\mathcal{S}^{l(\alpha)})\ast\hat{f}_{2}(\xi(\cdot))\right)(t).$

On the interval $[0,1]$, we have

$\displaystyle|(\mathcal{T}_{x^{0}}\xi)_{1}(t)|$
$\displaystyle\leq(|\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}|$
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(\|\xi\|_{\infty})\|\xi\|_{w}\int_{0}^{t}\left(|\mathcal{S}^{\alpha_{1}}(s)|+(|a_{22}|+|a_{12}|)|\mathcal{S}^{l(\alpha)}(s)|\right)ds$
$\displaystyle\leq\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)$
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(\|\xi\|_{\infty})\|\xi\|_{w}C(1+|a_{22}|+|a_{12}|)\int_{0}^{1}\frac{1}{s^{1-\nu}}ds$
$\displaystyle\leq\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)$
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(\|\xi\|_{\infty})\|\xi\|_{w}\frac{C(1+|a_{22}|+|a_{12}|)}{\nu},$
(50) $\displaystyle|(\mathcal{T}_{x^{0}}\xi)_{2}(t)|$
$\displaystyle\leq|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}|$
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(\|\xi\|_{\infty})\|\xi\|_{w}\int_{0}^{t}\left(|\mathcal{S}^{\alpha_{2}}(s)|+(|a_{12}|+|a_{11}|)|\mathcal{S}^{l(\alpha)}(s)|\right)ds$
$\displaystyle\leq\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)$
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(||\xi||_{\infty})||\xi||_{w}C(1+|a_{21}|+|a_{11}|)\int_{0}^{1}\frac{1}{s^{1-\nu}}ds$
$\displaystyle\leq\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)$
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(\|\xi\|_{\infty})\|\xi\|_{w}\frac{C(1+|a_{21}|+|a_{11}|)}{\nu}.$
(51)

For $t\in[1,\infty)$, then

$\displaystyle t^{\nu}|(\mathcal{T}_{x_{0}}\xi)_{1}(t)|$
$\displaystyle\leq\sup_{t\geq
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)||x^{0}_{2}|\Big{)}$
$\displaystyle+l_{\hat{f}}(||\xi||_{\infty})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}s^{\nu}|\xi(s)|ds$
$\displaystyle\leq\sup_{t\geq
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)||x^{0}_{2}|\Big{)}$
$\displaystyle+l_{\hat{f}}(||\xi||_{\infty})\|\xi\|_{w}\sup_{t\geq
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,$
(52) $\displaystyle t^{\nu}|(\mathcal{T}_{x_{0}}\xi)_{2}(t)|$
$\displaystyle\leq\sup_{t\geq
1}\Big{(}|a_{21}|t^{\nu}|\mathcal{R}^{\alpha_{2}}(t)||x^{0}_{1}|+(t^{\nu}|\mathcal{R}^{0}(t)|||+|a_{11}|t^{\nu}|\mathcal{R}^{\alpha_{1}}(t)|)|x_{2}^{0}|\Big{)}$
$\displaystyle+l_{\hat{f}}(||\xi||_{\infty})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}s^{\nu}\|\xi(s)\|ds$
$\displaystyle\leq\sup_{t\geq
1}\Big{(}|a_{21}|t^{\nu}|\mathcal{R}^{\alpha_{2}}(t)||x^{0}_{1}|+(t^{\nu}|\mathcal{R}^{0}(t)|||+|a_{11}|t^{\nu}|\mathcal{R}^{\alpha_{1}}(t)|)|x_{2}^{0}|\Big{)}$
$\displaystyle+l_{\hat{f}}(||\xi||_{\infty})\|\xi\|_{w}\sup_{t\geq
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.$
(53)

From (52) and (53), we obtain the estimates

$\displaystyle\|(\mathcal{T}_{x^{0}}\xi)_{1}\|_{w,1}$
$\displaystyle\leq(\|\mathcal{R}^{0}\|_{w,1}+|a_{22}|\cdot\|\mathcal{R}^{\alpha_{2}}\|_{w,1}+|a_{12}|\cdot\|\mathcal{R}^{\alpha_{1}}\|_{w,1})\|x^{0}\|$
(54)
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(||\xi||_{\infty})\|\xi\|_{w}\left(\frac{C}{\nu}(1+|a_{12}|+|a_{22}|)+M_{\alpha_{1}}+(|a_{22}|+|a_{12}|)M_{l(\alpha)}\right)$

and

$\displaystyle\|(\mathcal{T}_{x^{0}}\xi)_{2}\|_{w,1}$
$\displaystyle\leq(\|\mathcal{R}^{0}\|_{w,1}+|a_{21}|\cdot\|\mathcal{R}^{\alpha_{2}}\|_{w,1}+|a_{11}|\cdot\|\mathcal{R}^{\alpha_{1}}\|_{w,1})\|x^{0}\|$
(55)
$\displaystyle\phantom{=}\quad+l_{\hat{f}}(||\xi||_{\infty})\|\xi\|_{w}\left(\frac{C}{\nu}(1+|a_{21}|+|a_{11}|)+M_{\alpha_{2}}+(|a_{21}|+|a_{11}|)M_{l(\alpha)}\right)$

where $\|\xi\|_{w,1}:=\max\\{\sup_{t\in[0,1]}|\xi(t)|,\sup_{t\geq
1}t^{\nu}|\xi(t)|\\}$ for any $\xi\in C_{\infty}([0,\infty);\mathbb{R})$,
$M_{\beta}=\sup_{t\geq
1}t^{\nu}\int_{0}^{t}|\mathcal{S}^{\beta}(t-s)|s^{-\nu}ds$, for
$\beta\in\\{\alpha_{1},\alpha_{2},l(\alpha)\\}$, and By (54) and (55), we see
that

$\displaystyle\|(\mathcal{T}_{x^{0}}\xi)\|_{w}$
$\displaystyle\leq(2\|\mathcal{R}^{0}\|_{w,1}+(|a_{11}|+|a_{12}|)\|\mathcal{R}^{\alpha_{1}}\|_{w,1}+(|a_{21}|+|a_{22}|\cdot\|\mathcal{R}^{\alpha_{2}}\|_{w,1})\|x^{0}\|$
$\displaystyle+l_{\hat{f}}(||\xi||_{\infty})\|\xi\|_{w}\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).$

On the other hand, by virtue of the assumption that $\lim_{r\to 0}l(r)=0$, we
can choose $\varepsilon\in(0,\varepsilon_{0})$ so that

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

Take
$\delta=\varepsilon(1-r)/\left(2\|\mathcal{R}^{0}\|_{w,1}+(|a_{11}|+|a_{12}|)\|\mathcal{R}^{\alpha_{1}}\|_{w,1}+(|a_{21}|+|a_{22}|)\|\mathcal{R}^{\alpha_{2}}\|_{w,1}\right)$,
then for any initial condition $x^{0}\in B(0,\delta)$, we have

$\|\mathcal{T}_{x^{0}}\xi\|_{w}\leq\varepsilon,\;\forall\xi\in
B_{C_{w}}(0,\varepsilon),$

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)$,

$\displaystyle\|\mathcal{T}_{x^{0}}\xi-\mathcal{T}_{x^{0}}\hat{\xi}\|_{w}$
$\displaystyle\leq\left[\frac{2C}{\nu}+M_{\alpha_{1}}+M_{\alpha_{2}}+\sum_{i,j=1}^{2}|a_{ij}|(M_{l(\alpha)}+\frac{C}{\nu})\right]l_{\hat{f}}(\varepsilon)\|\xi-\hat{\xi}\|_{w}$
(56) $\displaystyle=r_{0}\|\xi-\hat{\xi}\|_{w}.$ (57)

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 (48) 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 (47).
This completes the proof. ∎

##  6 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

$f_{i}(t)=\begin{cases}1&\text{ if }0\leq t<1,\\\ \frac{1}{t^{2i}}&\text{ if
}t\geq 1\end{cases}\quad(i=1,2).$

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

######  Example 6.1.

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

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (58)

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 (58) tends to the origin as
$t\to\infty$ 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\to\infty$ and $j=1,2$, thus demonstrating
that the decay behaviour of $x_{j}(t)$ is indeed $O(t^{-1/3})$.

Figure 1: Solution to the differential equation (58) from Example 6.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)$.

######  Example 6.2.

Consider the two-component incommensurate fractional-order nonlinear system

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (59)

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

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.

Figure 2: Solution to the differential equation (59) from Example 6.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)$.

Figure 3: Solution to the differential equation (59) from Example 6.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)$.

######  Example 6.3.

Consider the fractional linear system

$\displaystyle\begin{cases}{}^{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{cases}\quad t>0.$ (60)

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\to\infty$ 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.

Figure 4: Solution to the differential equation (60) from Example 6.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)$.

######  Example 6.4.

Consider the system

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (61)

Based on Lemma 3.5 and Theorem 5.4, we see that the trivial solution of (61)
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.

Figure 5: Solution to the differential equation (61) from Example 6.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)$.

######  Example 6.5.

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

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (62)

The system (62) 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\to\infty$
as $O(t^{-0.3})$. Once again, our numerical results, shown in Figure 6,
support this statement.

Figure 6: Solution to the differential equation (62) 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

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (63)

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.

Figure 7: Solution to the differential equation (63) from Example 6.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)$.

######  Example 6.7.

Consider the two-component incommensurate fractional-order linear system

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (64)

The system (64) 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.

Figure 8: Solution to the differential equation (64) from Example 6.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)$.

######  Example 6.8.

Consider two-component incommensurate fractional-order nonlinear system

$\displaystyle\begin{cases}{}^{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{cases}\quad
t>0.$ (65)

Its characteristic function $Q(s)=s^{0.9}+2s^{0.5}-s^{0.3}+2$. According to
Lemma 3.7(ii) and Theorem 5.4, the trivial solution of (65) is Mittag-Leffler
stable as illustrated graphically in Figure 9.

Figure 9: Solution to the differential equation (65) from Example 6.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)$.

## Appendix: Proof of Lemma 4.1

###### Proof of Lemma 4.1.

Due to the fact that there are no zeros of the characteristic function $Q$ in
the closed right half of the complex plane, from Lemma 3.1(iii), we can find
$\delta>0$ (which is small enough) such that all zeros of $Q$ are not in the
domain $|\arg(s)|\leq\frac{\pi}{2}+\delta.$ Let $R>0$ be a large enough
constant such that

$\displaystyle|Q(s)|\geq\frac{1}{2}|s|^{\alpha_{1}+\alpha_{2}}\quad\mbox{
whenever }|s|\geq R.$ (66)

For $\mu>0$ and $\theta\in(0,\pi)$, we establish an oriented contour
$\gamma(\mu,\theta)$ formed by three segments:

  * •

$\left\\{s\in\mathbb{C}:|s|\geq\mu,\arg s=\theta\right\\}$,

  * •

$\left\\{s\in\mathbb{C}:|s|=\mu,|\arg s|\leq\theta\right\\}$,

  * •

$\left\\{s\in\mathbb{C}:|s|\geq\mu,\arg s=-\theta\right\\}$.

(i) Because all zeros of $Q$ (if they exist) lie on the left of the contour
$\gamma(R,\frac{\pi}{2}+\delta),$ using the same argument as in [21, Lemma
4.1], we obtain the representation

$\displaystyle\mathcal{R}^{\lambda}(t)=\frac{1}{2\pi
i}\int_{\gamma(R,\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds,\;t>0,\lambda\in\left\\{0,\alpha_{1},\alpha_{2}\right\\}.$
(67)

Choose $\varepsilon>0$ such that $Q$ has no zero in the ball
$\left\\{s\in\mathbb{C}:|s|\leq\varepsilon\right\\}.$ From (67), we have

$\displaystyle\mathcal{R}^{\lambda}(t)$ $\displaystyle=\frac{1}{2\pi
i}\int_{\Lambda_{t}^{\prime}}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\gamma(\frac{\varepsilon}{t},\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds,t\geq
1$ $\displaystyle=I_{1}(t)+I_{2}(t),$ (68)

where $\Lambda_{t}^{\prime}$ is the clockwise oriented contour bounding the
domain

$\Omega_{t}:=\left\\{s\in\mathbb{C}:\frac{\varepsilon}{t}<|s|<R,|\arg
s|<\frac{\pi}{2}+\delta\right\\},$

see Figure 10. Notice that $(s^{l(\alpha)-\lambda-1}e^{st})/Q(s)$ is analytic
on $\Omega_{t}\cup\Lambda_{t}^{\prime}$ for all $t\geq 1$. Thus, by applying
Cauchy’s theorem, we obtain

$I_{1}(t)=0\mbox{ for all }t\geq 1.$

Figure 10: The contours and sets used in the proof of Lemma 4.1: The radii of
the green and magenta circular arcs are $\varepsilon/t$ and $R$, respectively.
The contour $\gamma(R,\pi/2+\delta)$ comprises the upper blue ray, the magenta
circular arc, and the lower blue ray and is traversed from top to bottom;
$\Omega_{t}$ is the open set bounded by the magenta and green boundary lines
(so these magenta and green lines together form the contour
$\Lambda_{t}^{\prime}$). $\Lambda_{1}$ comprises the upper blue ray and the
upper green line; $\Lambda_{2}$ denotes the union of the lower blue ray and
the lower green line, and $\Lambda_{3}$ is the green circular arc.

Therefore, for each $t\geq 1$, we see that

$\displaystyle\mathcal{R}^{\lambda}(t)$ $\displaystyle=\frac{1}{2\pi
i}\int_{\gamma(\frac{\varepsilon}{t},\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds$
$\displaystyle=\frac{1}{2\pi
i}\int_{\Lambda_{1}}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\Lambda_{2}}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\Lambda_{3}}\frac{s^{l(\alpha)-\lambda-1}e^{st}}{Q(s)}ds$
$\displaystyle=I_{3}(t)+I_{4}(t)+I_{5}(t)$ (69)

with

$\displaystyle\Lambda_{1}$
$\displaystyle:=\left\\{s\in\mathbb{C}:|s|\geq\frac{\varepsilon}{t},\arg
s=\frac{\pi}{2}+\delta\right\\},$ $\displaystyle\Lambda_{2}$
$\displaystyle:=\left\\{s\in\mathbb{C}:|s|\geq\frac{\varepsilon}{t},\arg
s=-(\frac{\pi}{2}+\delta)\right\\},$ $\displaystyle\Lambda_{3}$
$\displaystyle:=\left\\{s\in\mathbb{C}:|s|=\frac{\varepsilon}{t},|\arg
s|\leq\frac{\pi}{2}+\delta\right\\}.$

Put

$\displaystyle\eta:=\underset{s\in\gamma(\varepsilon,\frac{\pi}{2}+\delta)\cup
B(0,\varepsilon)}{\inf}|Q(s)|.$ (70)

For $s\in\Lambda_{1},$
$s=re^{i(\frac{\pi}{2}+\delta)}=r(-\sin\delta+i\cos\delta)$ with
$r\geq\frac{\varepsilon}{t}$, and therefore

$\displaystyle I_{3}(t)=\frac{1}{2\pi
i}\int_{\varepsilon/t}^{\infty}\frac{r^{l(\alpha)-\lambda-1}e^{i(l(\alpha)-\lambda-1)(\frac{\pi}{2}+\delta)}e^{rt(-\sin\delta+i\cos\delta)}(-\sin\delta+i\cos\delta)}{\widetilde{Q}(r)}dr.$
(71)

Here, $\widetilde{Q}(r)=Q(re^{i(\frac{\pi}{2}+\delta)})$. From (70), we have
the estimate $|\widetilde{Q}(r)|\geq\eta$ for all
$r\geq\frac{\varepsilon}{t}$. This implies that

$\displaystyle|I_{3}(t)|\leq\frac{1}{2\pi\eta}\int_{\varepsilon/t}^{\infty}r^{l(\alpha)-\lambda-1}e^{-rt\sin\delta}dr.$
(72)

By the change of variable $r=u/(t\sin\delta)$,

$\displaystyle\int_{\varepsilon/t}^{\infty}r^{l(\alpha)-\lambda-1}e^{-rt\sin\delta}dr$
$\displaystyle\leq\frac{1}{(t\sin\delta)^{l(\alpha)-\lambda}}\int_{0}^{\infty}u^{l(\alpha)-\lambda-1}e^{-u}du$
$\displaystyle=\frac{1}{(t\sin\delta)^{l(\alpha)-\lambda}}\Gamma(l(\alpha)-\lambda)$
$\displaystyle\leq\frac{C_{1,1}}{t^{l(\alpha)-\lambda}}.$ (73)

Hence,

$\displaystyle|I_{3}(t)|\leq\frac{C_{1,1}}{t^{l(\alpha)-\lambda}}.$ (74)

Similarly, there is a $C_{1,2}>0$ such that

$\displaystyle|I_{4}(t)|\leq\frac{C_{1,2}}{t^{l(\alpha)-\lambda}}.$ (75)

For $s\in\Lambda_{3},$ $s=(\varepsilon/t)e^{i\varphi}$ with
$|\varphi|\leq\frac{\pi}{2}+\delta$, and so

$\displaystyle I_{5}(t)$ $\displaystyle=-\frac{1}{2\pi
i}\int_{-(\frac{\pi}{2}+\delta)}^{\frac{\pi}{2}+\delta}\frac{\left(\varepsilon/t\right)^{l(\alpha)-\lambda-1}e^{i\varphi(l(\alpha)-\lambda-1)}e^{\varepsilon(\cos\varphi+i\sin\varphi)}i(\varepsilon/t)e^{i\varphi}}{\hat{Q}(\varphi)}d\varphi$
(76)

where $\hat{Q}(\varphi)=Q(\frac{\varepsilon}{t}e^{i\varphi})$. From (70), we
know that $|\hat{Q}(\varphi)|\geq\eta$ for all
$\varphi\in[-(\frac{\pi}{2}+\delta),\frac{\pi}{2}+\delta]$. Thus

$\displaystyle|I_{5}(t)|$
$\displaystyle\leq\frac{1}{2\pi\eta}\left(\frac{\varepsilon}{t}\right)^{l(\alpha)-\lambda}\int_{-(\frac{\pi}{2}+\delta)}^{\frac{\pi}{2}+\delta}e^{\varepsilon\cos\varphi}d\varphi\leq\frac{C_{1,3}}{t^{l(\alpha)-\lambda}}.$
(77)

From (Proof of Lemma 4.1.), (74), (75) and (77), we obtain

$\displaystyle|\mathcal{R}^{\lambda}(t)|$
$\displaystyle\leq\frac{C}{t^{l(\alpha)-\lambda}}\leq\frac{C}{t^{\nu}}$ (78)

for all $t\geq 1$ and all
$\lambda\in\left\\{0,\alpha_{1},\alpha_{2}\right\\}$, with
$C:=C_{1,1}+C_{1,2}+C_{1,3}$.

(ii) For the proof of the seocnd statement, we first look at the case
$\beta\in\left\\{\alpha_{1},\alpha_{2}\right\\}$. Here, we apply the arguments
as in the proof of the part (i) above to obtain

$\displaystyle\mathcal{S}^{\beta}(t)$ $\displaystyle=\frac{1}{2\pi
i}\int_{\gamma(\frac{\varepsilon}{t},\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds$
$\displaystyle=\frac{1}{2\pi
i}\int_{\Lambda_{1}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\Lambda_{2}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\Lambda_{3}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds$

with each $t\geq 1$. In the same way as above, we can find a constant
$C_{2,1}$ so that the estimate

$\displaystyle|\mathcal{S}^{\beta}(t)|\leq\frac{C_{2,1}}{t^{l(\alpha)-\beta+1}}$
(79)

holds for all $t\geq 1$ and all
$\beta\in\left\\{\alpha_{1},\alpha_{2}\right\\}$. Clearly,
$t^{l(\alpha)-\beta+1}\geq t^{\nu+1}$ for all $t\geq 1$ and all
$\beta\in\left\\{\alpha_{1},\alpha_{2}\right\\}$. Thus,

$\displaystyle|\mathcal{S}^{\beta}(t)|\leq\frac{C_{2,1}}{t^{\nu+1}},\;\beta\in\left\\{\alpha_{1},\alpha_{2}\right\\},\;t\geq
1.$ (80)

Next, we consider the remaining case $\beta=l(\alpha)$. For $t\geq 1,$ we see

$\displaystyle\mathcal{S}^{l(\alpha)}(t)$ $\displaystyle=\frac{1}{2\pi
i}\int_{\gamma(\frac{\varepsilon}{t},\frac{\pi}{2}+\delta)}\frac{e^{st}}{Q(s)}ds$
$\displaystyle=\frac{1}{2\pi
i}\int_{\gamma(\frac{\varepsilon}{t},\frac{\pi}{2}+\delta)}\frac{1}{\det
A}e^{st}ds-\frac{1}{2\pi
i}\int_{\gamma(\frac{\varepsilon}{t},\frac{\pi}{2}+\delta)}\frac{(s^{\alpha_{1}+\alpha_{2}}-a_{11}s^{\alpha_{2}}-a_{22}s^{\alpha_{1}})e^{st}}{(\det
A)Q(s)}ds$ $\displaystyle=I_{6}(t)+I_{7}(t)$ (81)

By using the same estimates as in the proof of the part (i) above, there exist
constants $C_{2,2},C_{2,3}$ and $C_{2,4}$ such that

$\displaystyle|I_{7}(t)|\leq\frac{C_{2,2}}{t^{l(\alpha)+1}}+\frac{C_{2,3}}{t^{\alpha_{1}+1}}+\frac{C_{2,4}}{t^{\alpha_{2}+1}}\leq\frac{\sum_{2\leq
i\leq 4}C_{2,i}}{t^{\nu+1}},\;\forall t\geq 1.$ (82)

On the other hand, by the change of variable $s=\frac{u^{1/\mu}}{t}$ with some
$\mu\in(0,1)$, we find

$\displaystyle I_{6}(t)$ $\displaystyle=\frac{1}{2\pi\mu i}\frac{1}{t\det
A}\int_{\gamma(\varepsilon^{\mu},\mu(\frac{\pi}{2}+\delta))}e^{u^{1/\mu}}u^{(1-\mu)/\mu}du=\frac{1}{t\det
A}\frac{1}{\Gamma(0)}=0,$ (83)

the last equality being deduced from [15, eq. (1.52)]. Combining (Proof of
Lemma 4.1.) and (82), for each $t\geq 1$, we conclude

$\displaystyle|\mathcal{S}^{l(\alpha)}(t)|\leq\frac{C}{t^{\nu+1}}.$ (84)

(iii) For each $t\in(0,1)$ and
$\beta\in\left\\{\alpha_{1},\alpha_{2},l(\alpha)\right\\}$, we have

$\displaystyle\mathcal{S}^{\beta}(t)$ $\displaystyle=\frac{1}{2\pi
i}\int_{\gamma(R,\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds$
$\displaystyle=\frac{1}{2\pi
i}\int_{\Psi_{t}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\gamma(\frac{R}{t},\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds$
$\displaystyle=I_{8}(t)+I_{9}(t),$ (85)

where $\Psi_{t}$ is the boundary of the domain
$U_{t}:=\left\\{s\in\mathbb{C}:R<|s|<R/t,|\arg
s|<\frac{\pi}{2}+\delta\right\\}$. Since $s^{l(\alpha)-\beta}e^{st}/Q(s)$ is
analytic on $U_{t}\cup\Psi_{t}$ for $t\in(0,1)$, by applying Cauchy’s theorem,
we obtain $I_{8}(t)=0$ for all $t\in(0,1)$. Thus,

$\displaystyle\mathcal{S}^{\beta}(t)$ $\displaystyle=\frac{1}{2\pi
i}\int_{\gamma(\frac{R}{t},\frac{\pi}{2}+\delta)}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds$
$\displaystyle=\frac{1}{2\pi
i}\int_{\Psi_{1}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\Psi_{2}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds+\frac{1}{2\pi
i}\int_{\Psi_{3}}\frac{s^{l(\alpha)-\beta}e^{st}}{Q(s)}ds$
$\displaystyle=I_{10}(t)+I_{11}(t)+I_{12}(t)$ (86)

where

$\displaystyle\Psi_{1}$
$\displaystyle:=\left\\{s\in\mathbb{C}:|s|\geq\frac{R}{t},\arg
s=\frac{\pi}{2}+\delta\right\\},$ $\displaystyle\Psi_{2}$
$\displaystyle:=\left\\{s\in\mathbb{C}:|s|\geq\frac{R}{t},\arg
s=-(\frac{\pi}{2}+\delta)\right\\},$ $\displaystyle\Psi_{3}$
$\displaystyle:=\left\\{s\in\mathbb{C}:|s|=\frac{R}{t},|\arg
s|\leq\frac{\pi}{2}+\delta\right\\}.$

For $s\in\Psi_{1},$
$s=re^{i(\frac{\pi}{2}+\delta)}=r(-\sin\delta+i\cos\delta)$ with $r\geq R/t$,
and so

$\displaystyle I_{10}(t)=\frac{1}{2\pi
i}\int_{R/t}^{\infty}\frac{r^{l(\alpha)-\beta}e^{i(l(\alpha)-\beta)(\frac{\pi}{2}+\delta)}e^{rt(-\sin\delta+i\cos\delta)}}{\widetilde{Q}(r)}(-\sin\delta+i\cos\delta)dr$
(87)

where $\widetilde{Q}(r)=Q(re^{i(\frac{\pi}{2}+\delta)})$. From (66), we have
the estimate

$\displaystyle|\widetilde{Q}(r)|=|Q(re^{i(\frac{\pi}{2}+\delta)})|\geq\frac{1}{2}|re^{i(\frac{\pi}{2}+\delta)}|^{\alpha_{1}+\alpha_{2}}=\frac{1}{2}r^{\alpha_{1}+\alpha_{2}},\;\forall
r\geq\frac{R}{t}.$ (88)

This implies

$\displaystyle|I_{10}(t)|\leq\frac{1}{\pi}\int_{R/t}^{\infty}\frac{r^{l(\alpha)-\beta}e^{-rt\sin\delta}}{r^{l(\alpha)}}dr\leq\frac{1}{\pi
t\sin\delta}\int_{R/t}^{\infty}\frac{1}{r^{\beta+1}}dr=\frac{C_{3,1}}{t^{1-\beta}}.$
(89)

The second inequality here is obtained by applying the relation $e^{-x}\leq
1/x$ for $x>0$. Similarly,

$\displaystyle|I_{11}(t)|$ $\displaystyle\leq\frac{C_{3,2}}{t^{1-\beta}}.$
(90)

For $s\in\Psi_{3}$, $s=(R/t)e^{i\varphi}$ with
$|\varphi|\leq\frac{\pi}{2}+\delta$, thus

$\displaystyle I_{12}(t)=\frac{1}{2\pi
i}\int_{-(\frac{\pi}{2}+\delta)}^{\frac{\pi}{2}+\delta}\left(\frac{R}{t}\right)^{l(\alpha)-\beta+1}\frac{e^{i(l(\alpha)-\beta)\varphi}e^{R(\cos\varphi+i\sin\varphi)}}{\hat{Q}(\varphi)}ie^{i\varphi}d\varphi,$
(91)

where $\hat{Q}(\varphi)=Q(\frac{R}{t}e^{i\varphi})$. From (66), we have

$\displaystyle|\hat{Q}(\varphi)|=\left|Q\left(\frac{R}{t}e^{i\varphi}\right)\right|\geq\frac{1}{2}\left|\frac{R}{t}e^{i\varphi}\right|^{\alpha_{1}+\alpha_{2}}=\frac{1}{2}\left(\frac{R}{t}\right)^{\alpha_{1}+\alpha_{2}},\;\forall\varphi\in[-(\frac{\pi}{2}+\delta),\frac{\pi}{2}+\delta],$
(92)

and thus

$\displaystyle|I_{12}(t)|$
$\displaystyle\leq\frac{1}{2\pi}\int_{-(\frac{\pi}{2}+\delta)}^{\frac{\pi}{2}+\delta}\left(\frac{R}{t}\right)^{l(\alpha)-\beta+1}\frac{|e^{i(l(\alpha)-\beta)\varphi}|\cdot|e^{R(\cos\varphi+i\sin\varphi)}|}{|\hat{Q}(\varphi)|}|ie^{i\varphi}|d\varphi$
$\displaystyle\leq\frac{R^{1-\beta}}{t^{1-\beta}}\frac{1}{\pi}\int_{-(\frac{\pi}{2}+\delta)}^{\frac{\pi}{2}+\delta}e^{R\cos\varphi}d\varphi\leq\frac{C_{3,3}}{t^{1-\beta}}.$
(93)

From (Proof of Lemma 4.1.), (89), (90) and (Proof of Lemma 4.1.), we obtain

$\displaystyle|\mathcal{S}^{\beta}(t)|\leq\frac{C}{t^{1-\beta}},\forall
t\in(0,1),\beta\in\left\\{\alpha_{1},\alpha_{2},l(\alpha)\right\\}.$ (94)

Finally, (34) is an immediate consequence of (32) and (33). ∎

## References

  * [1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions. 10th printing with corrections, Washington: National Bureau of Standards, 1972.
  * [2] D. Bǎleanu and A.M. Lopes, Handbook of Fractional Calculus with Applications: Applications in Engineering, Life and Social Sciences, Part A. Berlin, Boston: De Gruyter, 2019.
  * [3] D. Bǎleanu and A.M. Lopes, Handbook of Fractional Calculus with Applications: Applications in Engineering, Life and Social Sciences, Part B. Berlin, Boston: De Gruyter, 2019.
  * [4] C. Bonnet and J.R. Partington, Coprime factorizations and stability of fractional differential systems. Systems and Control Letters, 41 (2000), no. 3, pp. 167–174.
  * [5] O. Brandibur and E. Kaslik, Stability of two-component incommensurate fractional-order systems and applications to the investigation of a FitzHugh-Nagumo neuronal model. Math. Methods in the Appl. Sci., 41 (2018), no. 17, pp. 7182–7194.
  * [6] O. Brandibur and E. Kaslik, Exact stability and instability regions for two-dimensional linear autonomous multi-order systems of fractional-order differential equations. Fractional Calculus and Applied Analysis, 24 (2021), no. 1, pp. 225–253.
  * [7] N.D. Cong, T.S. Doan, S. Siegmund and H.T. Tuan, Linearized asymptotic stability for fractional differential equations. Electron. J. Qualitative Theory Differ. Equations 39 (2016), 1–13.
  * [8] N.D. Cong, H.T. Tuan and H. Trinh, On asymptotic properties of solutions to fractional differential equations. J. Math. Anal. Appl. 484 (2020), Article ID 123759.
  * [9] K. Diethelm,  The Analysis of Fractional Differential Equations. An application-oriented exposition using differential operators of Caputo type.  Lecture Notes in Mathematics 2004.  Berlin: Springer, 2010.
  * [10] K. Diethelm, S. Siegmund and H.T. Tuan, Asymptotic behaviour of solutions of linear multi-order fractional differential systems. Fractional Calculus and Applied Analysis, 20 (2017), no. 5, pp. 1165–1195.
  * [11] R. Garrappa, Numerical solution of fractional differential equations: A survey and a software tutorial. Mathematics 6 (2018), 16.
  * [12] E. Ivanova, X. Moreau and R. Malti, Stability and resonance conditions of second-order fractional systems. Journal of Vibration and Control, 24 (2016), no. 4, pp. 659–672.
  * [13] D. Matignon,  Stability results for fractional differential equations with applications to control processing.  Computational Eng. in Sys. Appl. 2 (1996), 963–968.
  * [14] I. Petráš, Handbook of Fractional Calculus with Applications: Applications in Control. Berlin, Boston: De Gruyter, 2019.
  * [15] I. Podlubny.  Fractional Differential Equations. An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of their Solution and some of their Applications. Mathematics in Science and Engineering, 198.  San Diego, CA: Academic Press, 1999.
  * [16] J. Sabatier, C. Farges and J.-C. Trigeassou, A stability test for non-commensurate fractional order systems. Systems and Control Letters, 62 (2013), no. 9, pp. 739–746.
  * [17] J. Shen and J. Lam, Stability and performance analysis for positive fractional-order systems with time varying delays. IEEE Trans. Automat. Control., 61 (2016), no. 9, pp. 2676–2681.
  * [18] V.E. Tarasov, Handbook of Fractional Calculus with Applications: Applications in Physics, Part A. Berlin, Boston: De Gruyter, 2019.
  * [19] V.E. Tarasov, Handbook of Fractional Calculus with Applications: Applications in Physics, Part B. Berlin, Boston: De Gruyter, 2019.
  * [20] J. Trigeassou, A. Benchellal, N. Maamri and T. Poinot, A frequency approach to the stability of fractional differential equations. Transactions on Systems, Signals and Devices, 4 (2009), no. 1, pp. 1–25.
  * [21] H.T. Tuan and S. Siegmund, Stability of scalar nonlinear fractional differential equations with linearly dominated delay. Fractional Calculus and Applied Analysis, 23 (2020), no. 1, pp. 250–267.
  * [22] H.T. Tuan, H. Trinh, A qualitative theory of time delay nonlinear fractional-order systems. SIAM Journal on Control and Optimization, 3 (2020), pp. 1491–1518.
  * [23] H.T. Tuan and H. Trinh, Global attractivity and asymptotic stability of mixed-order fractional systems. IET Control Theory and Applications, 14 (2020), pp. 1240–1245.
  * [24] H.T. Tuan, H. Trinh, J. Lam, Necessary and sufficient conditions of the positivity and stability to mixed fractional-order systems. International Journal of Robust and Nonlinear Control, 31 (2021), no. 1, pp. 37–50.
  * [25] G. Vainikko, Which functions are fractionally differentiable? Z. Anal. Anwend., 35 (2016), no. 4, pp. 465–487.
  * [26] S. Zhang, L. Liu, D. Xue and Y.Q. Chen, Stability and resonance analysis of a general non-commensurate elementary fractional-order system. Fractional Calculus and Applied Analysis, 23 (2020), no. 1, pp. 183–210.