# Asymptotic behavior of solutions of linear multi-order fractional differential systems

Kai Diethelm, Stefan Siegmund and H.T. Tuan

## Abstract

In this paper, we investigate some aspects of the qualitative theory for multi-order fractional differential equation systems. First, we obtain a fundamental result on the existence and uniqueness for multi-order fractional differential equation systems. Next, a representation of solutions of homogeneous linear multi-order fractional differential equation systems in series form is provided. Finally, we give characteristics regarding the asymptotic behavior of solutions to some classes of linear multi-order fractional differential equation systems.

### 1 Introduction

In recent years, fractional calculus has received increasing attention due to its applications in a variety of disciplines such as mechanics, physics, chemistry, biology, electrical engineering, control theory, material science, mathematical psychology. For more details, we refer the reader to the monographs [2, 7, 15, 19, 20].

A particularly interesting aspect in this connection that does not pertain to classical mathematical models using integer-order differential operators has recently been discussed in the context of a number of applications in the life sciences [2, 8, 9]: It appears that certain real world problems can be described by a system of fractional differential equations where each equation may have an order that differs from the orders of the other equations of the system. We shall call such systems multi-order fractional differential systems.

Among the published papers, it seems that the authors mainly concentrated on approximating solutions of multi-order fractional differential equations, see e.g. [1, 10, 11, 12, 14, 17, 18, 21, 22, 25, 26]. The investigation of the analytical properties of such systems is often restricted to the case where the orders of the differential operators are rational [5, 6, 7, 16]. For the general case, rigorous mathematical studies of even the most fundamental questions in this context do not seem to be readily available.

Therefore, in this paper we consider d-dimensional linear multi-order fractional differential equation systems

Dαixi(t)=j=1daijxj(t)+gi(t),i=1,2,,d,(1.1)

with orders αi ∈ (0, 1], coefficients aij ∈ ℂ, gi : [0, ∞) → ℂ continuous, i, j = 1, …, d, and the Caputo differential operator of order α > 0

Dαy(t):=JααDαy(t)

which is defined for Cα functions y : [0, T] → ℂd, T > 0, with the classical derivative D and the Riemann-Liouville operator

Jβy(t):=1Γ(β)0t(ts)β1y(s)ds

for β > 0 and J0y(t) := y(t) (see e.g. [7]). Note that Dαy can also be defined for not necessarily differentiable functions, e.g. if α ∈ (0, 1), for continuous functions y for which limt→0tα(v(t) − v(0)) exists, is finite, and

limθ1supt[0,T]θtt(ts)(α+1)(v(t)v(s))ds=0,

cf. [23, Theorem 5.2].

For convenience, we use the notation

D(α1,,αd)x(t):=Dα1Dαdx(t)=Dα1x1(t)Dαdxd(t).(1.2)

Using the notation A = (aij)i,j=1,…,d∈ℂd×d, g = (g1, …, gd) and x(t) = (x1(t), …, xd(t)), the system (1.1) can then be rewritten as

Dα1Dαdx(t)=Ax(t)+g(t).(1.3)

Of central importance are the two-parameter Mittag-Leffler functions Eα,β : ℂ → ℂ, α > 0, β ≥ 0, with

Eα,β(z):=j=0zjΓ(αj+β)(1.4)

and the one-parameter Mittag-Leffler functions Eα : ℂ → ℂ, α > 0, defined by Eα := Eα,1 (see e.g. [13]).

The structure of the paper is as follows. In Section 2, we first introduce a result on the existence and uniqueness of solutions to multi-order fractional differential equations. Then, we give a representation of solutions to homogeneous linear multi-order fractional differential equations in series form. Section 3 is devoted to the study of the asymptotic behavior of solutions of linear multi-order fractional differential equations. More precisely, we obtain some criterion on the asymptotic behavior of solutions to these equations. Some auxiliary results concerning the Mittag-Leffler functions and the asymptotic behavior of solutions of scalar linear fractional differential equations are shown in Appendix A.

### 2 Fundamental theory of multi-order fractional differential equations

In this section we provide some fundamental results regarding multi-order fractional differential equations. Specifically, we shall prove a Picard-Lindelöf type existence and uniqueness result in Subsection 2.1, and Subsection 2.2 will be devoted to a description of the structure of the associated solutions in the linear case.

#### 2.1 Existence and uniqueness of solutions to a class of multi-order fractional differential equations

Let T > 0. In this subsection we consider the existence and uniqueness of solutions to the multi-order fractional differential equation

Dαx(t)=f(t,x(t)),t(0,T],(2.1)

where α := (α1, …, αd) ∈ (0, 1]d and f = (f1, …, fd) : [0, T] × ℂd → ℂd is continuous. With similar arguments as in [7, Chapter 6] or [15, § 3.5], one can show that for x0=(x10,,xd0) ∈ ℂd and a continuous function x : [0, T] → ℂd for which Dαx is defined (cf. [23, Theorem 5.2]), the following two statements are equivalent:

1. x satisfies the d-dimensional differential equation system (2.1) together with the initial condition x(0) = x0,

2. x satisfies the Volterra integral equation

x(t)=x0+Jα[f(,x())](t)t[0,T],(2.2)

where Jα[f(⋅, x(⋅))] (t):= (Jα1[f1(⋅, x(⋅))](t),…,Jαd [fd(⋅, x(⋅))] (t)).

Following the usual convention, we define solutions of (2.1) by considering (2.2) for continuous functions.

#### Definition 2.1

A continuous function x : [0, T] → ℂd is called a solution to the differential equation (2.1) with initial condition x(0) = x0 if it satisfies the integral equation (2.2).

#### Remark 2.2

Because we assume the function f to be continuous, we can see that, for every solution x of (2.1) in the sense of Definition 2.1, the function f(⋅, x(⋅)) is continuous, too. Therefore, in view of the fact that the solution x satisfies the integral equation (2.2), it follows for i ∈ {1, 2, …, d} that the i-th component of x can be written as the sum of a constant and the Riemann-Liouville integral of order αi of a continuous function. Using the arguments of [7, proof of Theorem 3.7], we can then conclude that xi fulfils the conditions of [23, Theorem 5.2] and thus that Dαixi exists and is continuous. Therefore, under the continuity assumption on f, a solution to (2.1) in the sense of Definition 2.1 is automatically a strong solution to the differential equation in the classical sense.

Our basic assumption on the given function f will be that all its components fi : [0, T] × ℂd→ ℂ are continuous and satisfy a Lipschitz condition with respect to the second variable, i.e.

|fi(t,y)fi(t,z)|Lyz,t[0,T],y,zCd(2.3)

with some constant L > 0, where ∥⋅∥ is the max norm on ℂd, i.e., ∥y∥ = max {|y1|, …, |yd|} for all y = (y1, …, yd) ∈ ℂd.

We are now in a position to formulate a result on unique existence of solutions of initial value problems.

#### Theorem 2.3

Consider the equation(2.1). Assume that the function f is continuous and satisfies the Lipschitz condition(2.3). Then, for any x0=(x10,,xd0) ∈ ℂd, the differential equation (2.1) has a unique solution in the set C([0, T]; ℂd) that satisfies the initial condition x(0) = x0.

#### Proof

Let λ > 0 be a constant such that

max1idLλαi<1.

On the space C([0, T];ℂd) we define a new norm ∥⋅∥λ as

φλ:=sup0tTφ(t)exp(λt).

Using standard arguments, it is easy to see that (C([0, T];ℂd), ∥⋅∥λ) is a Banach space. For any x0=(x10,,xd0) ∈ ℂd, we define an operator 𝓣x0 : C([0, T]; ℂd) → C([0, T];ℂd) by

Tx0φ(t):=((Tx0φ)1(t),,(Tx0φ)d(t)),

where for 1 ≤ id, t ∈ [0, T] and φC([0, T];ℂd),

(Tx0φ)i(t):=xi0+1Γ(αi)0t(tτ)αi1fi(τ,φ(τ))dτ.

We see that for every φ, φ^C([0, T];ℂd), every t ∈ [0, T] and all 1 ≤ id,

|(Tx0φ)i(t)(Tx0φ^)i(t)|exp(λt)LΓ(αi)exp(λt)0t(tτ)αi1exp(λτ)φ(τ)φ^(τ)exp(λτ)dτLΓ(αi)exp(λt)0t(tτ)αi1exp(λτ)dτsup0θTφ(θ)φ^(θ)exp(λθ)LΓ(αi)0tuαi1exp(λu)duφφ^λ=LΓ(αi)λαi0λtvαi1exp(v)dvφφ^λLΓ(αi)λαi0vαi1exp(v)dvφφ^λ=Lλαiφφ^λ.(2.4)

It is clear that the operator 𝓣x0 maps the space (C([0, T];ℂd), ∥⋅∥λ) to itself; moreover, from (2.4) we obtain the estimate

Tx0φTx0φ^λLλαiφφ^λ,φ,φ^C([0,T];Cd)

which, by definition of λ, shows that this operator is a contractive mapping on this space. Due to the fact that (C([0, T];ℂd), ∥⋅∥λ) is a Banach space, by Banach’s fixed point theorem, there exists a unique fixed point φ of 𝓣x0 in this space. This fixed point is the unique solution of the Volterra equation (2.2) and hence, as stated above, also the unique solution to the initial value problem consisting of the differential equation (2.1) and the initial condition x(0) = x0 in C([0, T]; ℂd). The proof is complete.□

In Section 3 we shall look at the behavior of solutions to multi-order systems as the independent variable goes to infinity. For this purpose, it is important to have an existence and uniqueness result that is not restricted to functions defined on bounded intervals. Fortunately, the following result immediately follows from Theorem 2.3.

#### Corollary 2.4

Let f : [0, ∞) × ℂd → ℂdbe continuous and satisfy a Lipschitz condition with respect to the second variable. Moreover, let α ∈ (0, 1]dand x0 ∈ ℂd. Then, the initial value problem

Dαx(t)=f(t,x(t))(t>0),x(0)=x0,

has a unique solution in C([0, ∞); ℂd).

#### 2.2 A representation of solutions to homogeneous linear multi-order fractional differential equations

In this subsection we concentrate on a particularly important and fundamental special case of the class of differential equations discussed in Subsection 2.1, namely we shall look at the solutions to homogeneous linear equations with constant coefficients, i.e. to differential equations of the form

Dα1Dαdx(t)=Ax(t),(2.5)

which is the special case of (1.3) where g(t) = 0 for all t.

Our basic result in this section, Theorem 2.6, provides some information about the structure of the solutions to the system (2.5) in the case of an arbitrary matrix A ∈ ℂd×d and an arbitrary vector (α1, …, αd) ∈ (0, 1]d.

In order to motivate our results, we start with the case d = 2. In this case, the system (2.5) has the form

Dα1x1(t)=a11x1(t)+a12x2(t),(2.6a)
Dα2x2(t)=a21x1(t)+a22x2(t).(2.6b)

First of all, Corollary 2.4 asserts that, for any initial condition (x1(0), x2(0)) = (x10,x20) ∈ ℂ2, this system has a unique continuous solution x = (x1, x2) on [0, ∞). Moreover, for equations of this structure, the fractional version of the variation-of-constants method [7, Theorem 7.2 and Remark 7.1] provides the relations

x1(t)=x10Eα1(a11tα1)+a120t(ts)α11Eα1,α1a11(ts)α1x2(s)ds,(2.7a)
x2(t)=x20Eα2(a22tα2)+a210t(ts)α21Eα2,α2a22(ts)α2x1(s)ds,(2.7b)

for all t ≥ 0. This representation indicates that we should seek the solution components in the class of generalized power series of the form

x1(t)=x10+k=1=0b1ktkα1+α2,(2.8a)
x2(t)=x20+k=0=1b2ktkα1+α2.(2.8b)

Assuming a suitable convergence behavior of these series, we may differentiate in a termwise manner and obtain

Dα1x1(t)=k=1=0b1kΓ(kα1+α2+1)Γ((k1)α1+α2+1)t(k1)α1+α2=k=0=0b1,k+1,Γ((k+1)α1+α2+1)Γ(kα1+α2+1)tkα1+α2,Dα2x2(t)=k=0=1b2kΓ(kα1+α2+1)Γ(kα1+(1)α2+1)tkα1+(1)α2=k=0=0b2,k,+1Γ(kα1+(+1)α2+1)Γ(kα1+α2+1)tkα1+α2.

Plugging these representations into the differential equation system (2.5), we find

a11x10+a11k=1=0b1ktkα1+α2+a12x20+a12k=0=1b2ktkα1+α2=k=0=0b1,k+1,Γ((k+1)α1+α2+1)Γ(kα1+α2+1)tkα1+α2,a21x10+a21k=1=0b1ktkα1+α2+a22x20+a22k=0=1b2ktkα1+α2=k=0=0b2,k,+1Γ(kα1+(+1)α2+1)Γ(kα1+α2+1)tkα1+α2.

A comparison of coefficients of t1+α2 then yields the equations

b110=1Γ(α1+1)(a11x10+a12x20),(2.9a)
b201=1Γ(α2+1)(a21x10+a22x20),(2.9b)
b1,k+1,0=Γ(kα1+1)Γ((k+1)α1+1)a11b1k0(k=1,2,),(2.9c)
b1,1,=Γ(α2+1)Γ(α1+α2+1)a12b20(=1,2,),(2.9d)
b1,k+1,=Γ(kα1+α2+1)Γ((k+1)α1+α2+1)(a11b1k+a12b2k)(k,=1,2,),(2.9e)
b2,0,+1=Γ(α2+1)Γ((+1)α2+1)a22b20(=1,2,),(2.9f)
b2,k,1=Γ(kα1+1)Γ(kα1+α2+1)a21b1k0(k=1,2,),(2.9g)
b2,k,+1=Γ(kα1+α2+1)Γ(kα1+(+1)α2+1)(a21b1k+a22b2k)(k,=1,2,).(2.9h)

Formally introducing the quantities

b10=0for =1,2,andb2k0=0for k=1,2,,b100=x10andb200=x20,(2.10a)

we see that the system (2.9) can be simplified to

b1,k+1,=Γ(kα1+α2+1)Γ((k+1)α1+α2+1)(a11b1k+a12b2k)(k,=0,1,2,),(2.10b)
b2,k,+1=Γ(kα1+α2+1)Γ(kα1+(+1)α2+1)(a21b1k+a22b2k)(k,=0,1,2,).(2.10c)

A brief inspection of these formulas reveals that, given the initial values x10 and x20, they can indeed be used to compute all coefficients that appear in the representation (2.8) in a recursive manner. Specifically, the coefficients b1kℓ and b2kℓ for k + = μ can be computed via eqs. (2.10b) and (2.10c), respectively, and this computation only requires the knowledge of b1kℓ and b2kℓ with k + = μ − 1. Thus one can first compute all b1kℓ and b2kℓ with k + = 1, then with k + = 2, etc.

A closer look at the recurrence relations (2.10) allows us to prove that the series from (2.8) converge for all t ≥ 0. To this end we first state a preliminary result.

#### Lemma 2.5

Let the values b1kℓand b2kℓbe defined as in (2.10) with arbitraryx10,x20 ∈ ℂ. Then, for j ∈ {1, 2} the series

sj(z):=k=0=0|bjk|zk+

is convergent for all z ∈ ℂ.

Actually, it is immediately clear that the desired convergence property is a consequence of this lemma, since the series

k=0=0|bjk||t|kα1+α2

is, on the one hand, a majorant for xj(t) and is, on the other hand, convergent for all t > 0 according to

k=0=0|bjk||t|kα1+α2k=0=0|bjk||t|(k+)max{α1,α2}=sj(tmax{α1,α2})for t1,k=0=0|bjk||t|(k+)min{α1,α2}=sj(tmin{α1,α2})for t<1.

#### Proof of lemma 2.5

Since the series in question does not have any negative summands, we may rearrange the terms according to powers of z; this yields

sj(z)=k=0μ=0k|bj,μ,kμ|zk.

It is therefore evident that, in order to investigate the convergence radius of this series, we need to estimate expressions of the form

βjk:=μ=0k|bj,μ,kμ|.

In fact, we shall demonstrate that for sufficiently large k

0β1k+β2kc1c2kΓ(kα+1),(2.11)

where c1 and c2 are certain positive constants and

α:=min{α1,α2}.

Equation (2.11) tells us that the classical power series for the Mittag-Leffler function Eα* — that is well known to be convergent on the entire complex plane — evaluated at c2 |z| is a majorant for the series s1 and s2 that we are interested in, and hence the series expansions for s1(z) and s2(z) also converge for all z as required.

Thus, it only remains to prove (2.11). The left inequality is clear by definition. To prove the right inequality, we employ the relations (2.10a), (2.10b) and (2.10c) and see, using the notation a : = maxi,j ∈ {1, 2} |aij|, that we have for k ≥ 2 the following chain of inequalities:

μ=0k|b1,μ,kμ|+|b2,μ,kμ|μ=1k|b1,μ,kμ|+μ=0k1|b2,μ,kμ|a¯μ=1kΓ((μ1)α1+(kμ)α2+1)Γ(μα1+(kμ)α2+1)(|b1,μ1,kμ|+|b2,μ1,kμ|)+a¯μ=0k1Γ(μα1+(kμ1)α2+1)Γ(μα1+(kμ)α2+1)(|b1,μ,kμ1|+|b2,μ,kμ1|)=a¯μ=1kΓ((μ1)α1+(kμ)α2+1)Γ(μα1+(kμ)α2+1)(|b1,μ1,kμ|+|b2,μ1,kμ|)+a¯μ=1kΓ((μ1)α1+(kμ)α2+1)Γ((μ1)α1+(kμ+1)α2+1)(|b1,μ1,kμ|+|b2,μ1,kμ|)=a¯μ=1kwk,μ(α1,α2)(|b1,μ1,kμ|+|b2,μ1,kμ|)

with

wk,μ(α1,α2)=Γ((μ1)α1+(kμ)α2+1)Γ(μα1+(kμ)α2+1)+Γ((μ1)α1+(kμ)α2+1)Γ((μ1)α1+(kμ+1)α2+1).

Both fractions on the right-hand side have the same numerator but their denominators differ by α2α1; the well known monotonicity of the Gamma function thus allows us to conclude that, for sufficiently large k, we have

wk,μ(α1,α2)2Γ((μ1)α1+(kμ)α2+1)Γ((μ1)α1+(kμ)α2+α+1)=2Γ(u+μ(α1α2))Γ(u+μ(α1α2)+α)(2.12)

with u := −α1 + k α2 + 1. For γ > 0 and z → ∞, Stirling’s formula yields the asymptotic relation Γ(z) / Γ(z + γ) = zγ (1 + o(1)) which is monotonically decreasing in z. Hence, for sufficiently large k, the quotient on the right-hand side of (2.12) is monotonically decreasing with respect to μ for α1α2 and monotonically increasing with respect to μ if α1 < α2. Therefore, the maximum of this expression over all admissible values of μ is attained at μ = 1 if α1α2 and at μ = k if α1 < α2. These observations may be summarized in the form

wk,μ(α1,α2)2Γ(k12(α1+α2|α1α2|)+1)Γ(k12(α1+α2|α1α2|)+α+1)=2Γ((k1)α+1)Γ(kα+1),

and this implies

β1k+β2k2a¯Γ((k1)α+1)Γ(kα+1)μ=1k(|b1,μ1,kμ|+|b2,μ1,kμ|)=2a¯Γ((k1)α+1)Γ(kα+1)μ=0k1(|b1,μ,k1μ|+|b2,μ,k1μ|)=2a¯Γ((k1)α+1)Γ(kα+1)(β1,k1+β2,k1),

if k is large enough. Thus, for a sufficiently large and fixed constant N and arbitrary k, by induction, we deduce the estimate

β1,N+k+β2,N+k(2a¯)kΓ(Nα+1)Γ((N+k)α+1)(β1,N+β2,N),

which shows (2.11) and completes the proof of Lemma 2.5. □

The same ideas and methods can be applied if the dimension of the fractional differential equation system is greater than 2. This leads to the following result:

#### Theorem 2.6

Letα = (α1, …, αd) ∈ (0, 1]dand A ∈ ℂd × d. Then, for eachx0 ∈ ℂd, the initial value problem

Dαx(t)=Ax(t),x(0)=x0,(2.13)

has a uniquely determined solution inC([0,∞); ℂd). The components of this solution can be expressed in the form

xj(t)=k=01,2,,j1,j+1,,d=1bk,1,2,,j1,j+1,,dtkαj+μ=1,μjdμαμ,(2.14)

and the series in eq. (2.14)converges for all t ≥ 0.

### 3 Asymptotic behavior of solutions of multi-order fractional differential equations

Having established these foundations, we can now come to the core of this paper, namely the discussion of the asymptotic behavior of solutions of linear multi-order fractional differential systems.

#### 3.1 Systems with (block) triangular coefficient matrices

Assume that αk ∈ (0,1] for 1 ≤ kd. In the case that the coefficient matrix A of the system (1.1) has a triangular structure, we provide a detailed investigation of the asymptotic behaviour of the system’s solutions. More precisely, we obtain a necessary and sufficient condition such that all solutions of the homogeneous system associated to (1.1) tend to zero at infinity, and we derive sufficient conditions for all solutions of the full inhomogeneous system (1.1) to have this property. In this context we stress that the αk may be completely arbitrary numbers from the interval (0,1]; in particular it is allowed that αk = αk for some kk′.

Thus, let us now first consider the system

Dαixi(t)=j=idaijxj(t),1id,(3.1a)

i.e. the case of a homogeneous system with an upper triangular matrix A, together with the initial condition

xi(0)=xi0,1id.(3.1b)

In order to exclude the pathological and practically irrelevant case where the right-hand sides of certain equations from the system (3.1a) do not depend on their respective unknown functions, we shall explicitly assume throughout this subsection that aii ≠ 0 for all i = 1,2,…, d. In other words, we assume the matrix A to be not only upper triangular but also nonsingular.

The case where A is of lower triangular form can be handled in a completely analog manner; we shall not treat this case explicitly. The associated inhomogeneous system will be discussed later; cf. Corollary 3.3.

In the simplest nontrivial case d = 2, the system (3.1a) has the form

Dα1x1(t)=a11x1(t)+a12x2(t),Dα2x2(t)=a11x1(t)+a22x2(t),

and it is a relatively simple matter to explicitly compute its solution. Specifically, in view of the triangular structure of the coefficient matrix, one can solve the second equation of the system directly and obtain the well known result [7, Theorem 4.3]

x2(t)=x20Eα2(a22tα2).(3.2a)

This result can be plugged into the system’s first equation which then takes the form

Dα1x1(t)=a11x1(t)+a12x20Eα2(a22tα2).

For equations of this structure, the fractional version of the variation-of-constants method [7, Theorem 7.2 and Remark 7.1] provides the solution

x1(t)=x10Eα1(a11tα1)=+a12x200t(ts)α11Eα1,α1a11(ts)α1Eα2(a22sα2)ds.(3.2b)

From the representation (3.2) it is evident that the solution vector (x1, x2) is an element of the function space C[0,∞). Moreover, the power series representations of the Mittag-Leffler functions imply that the component x2(t) can be written as a power series in tα2, and therefore its asymptotic behavior as t → 0+ is of the form

x2(t)=x20+c2a22Γ(α2+1)tα2+O(t2α2),

whereas the behavior of x1(t) in this respect can be described by

x1(t)=x10+c1a11Γ(α1+1)tα1+o(tα1)

with some constant c1 ∈ ℂ. The arguments employed in Subsection 2.2 can be used to derive more details.

These considerations can directly be generalized to homogeneous upper triangular systems of arbitrary dimension d. In this case we obtain the set of equations

xi(t)=xi0Eαi(aiitαi)+j=i+1daij0t(ts)αi1Eαi,αiaii(ts)αixj(s)ds(3.3)

for i = d, d−1, …, 1 which can be recursively evaluated to explicitly compute the solutions.

Some known results about the asymptotic behavior of the Mittag-Leffler functions admit to draw the conclusions required in the asymptotic behavior analysis. The main result in this context is the following theorem. The proof of its statements requires a number of auxiliary results that can be considered as minor extensions of already known theorems and lemmas. Since these extensions may be of a certain degree of independent interest, we have explicitly formulated and collected them, together with complete proofs, in Appendix A.

#### Theorem 3.1

1. Every solution of the system(3.1a)converges to zero at infinity if and only if

|arg(akk)|>αkπ2k{1,,d}.(3.4)
2. If there existsk ∈ {1, …, d} such that |arg{(akk)| < αkπ / 2 then there exists some x0such that the solution to the system(3.1a)that satisfies the initial condition x(0) = x0is unbounded.

#### Proof

For the proof of part (i), we will first show that the condition (3.4) is sufficient to assert that every solution of (3.1a) converges to zero at infinity. Indeed, for any initial value x0=(x10,,xd0)Cd, we denote the solution of (3.1a) starting from x0 by φ (⋅, x0) = (φ1(⋅, x0), …, φd(⋅, x0)). Our proof will use mathematical induction over the index j of the components of the solution vector in a backward direction. Thus, for our induction basis we consider j = d. Since the d-th equation of the system (3.1a) reads

it follows from Lemma A4(i) that the condition |arg add| > αdπ / 2 is sufficient to assert that φd (t, x0) → 0 as t → ∞ for all x0. For the induction step, we assume that we have already shown that the components d, d−1, …, j+1 of the solution tend to 0 as t → ∞ for any choice of the initial values. Then we need to prove that this is also true for the j-th component. To this end we recall that the j-th component of the differential equation system (3.1a) reads

Dαjxj(t)=ajjxj(t)+k=j+1dajkφk(t,x0).

All the terms in the sum are already known and, because of the induction hypothesis, they are continuous and tend to zero as t → ∞. Thus we may apply Lemma A4(i) and immediately deduce that xj has this property as well.

To conclude the proof of part (i) we now have to demonstrate that (3.4) is also necessary for all solutions of (3.1a) tend to zero as t→∞. To this end we assume that (3.4) does not hold. Then there exists an index k0 ∈ {1, 2, …, d} which satisfies

|arg(aii)|>αiπ2 for k0+1id and |arg(ak0,k0)|αk0π2,

i.e. k0 is the largest index for which (3.4) is violated. Consider the equation

Dαk0xk0(t)=ak0k0xk0(t)+f(t) with f(t):=i=k0+1dak0ixi(t).(3.5)

Since (3.4) is true for all i > k0, the arguments used above imply that f is continuous and tends to zero at infinity. As in the considerations above, we may use the fractional variation-of-constants method [7, Theorem 7.2 and Remark 7.1] to see that the set of all solutions to (3.5) consists of the functions

φk0(t,x0)=xk00Eαk0(ak0,k0tαk0)+h(t),(3.6)

where xk00 runs through the entire complex plane and where

h(t):=0t(tτ)αk01Eαk0,αk0(ak0,k0(tτ)αk0)f(τ)dτ.

The well known asymptotic behavior of the Mittag-Leffler functions [13, Proposition 3.6 and Theorem 4.3] then implies that Eαk0(ak0,k0tαk0) does not converge to 0 as t → ∞ because of our assumption on the relation of αk0 and |arg ak0, k0|. Now assume that there exists some xk00 ∈ ℂ such that φk0(t, x0) → 0 as t → ∞. Then, it follows that for every x~0Cd with xk0=x~k0 for k = k0+1, …, d and xk00x~k00, we have

φk0(t,x~0)=x~k00Eαk0(ak0,k0tαk0)+h(t)=(x~k00xk00)Eαk0(ak0,k0tαk0)+φk0(t,x0).

For t → ∞, the last summand on the right-hand side of this equality tends to zero but the other summand does not, and hence we conclude that φ(t, x~0) does not tend to zero as t → ∞ which yields our required contradiction.

For the proof of part (ii), we — much as above — know that there exists an index k0 ∈ {1, 2, …, d} which satisfies

|arg(aii)|αiπ2 for k0+1id and |arg(ak0k0)|<αk0π2.

We may then proceed in the same way as in the second part of the proof of (i). However, now we know that |Eαk0 (ak0k0tαk0)| → ∞ as t → ∞, and therefore we may even conclude that there exists some xk00 ∈ ℂ such that φk0(t, x0) is unbounded. □

#### Remark 3.2

The same arguments can be used if the coefficient matrix A of the system has a block-upper triangular structure and the differentiation matrix on the left-hand side of the differential equation has a block structure with identical block sizes where each block consists of differential operators of the same order, i.e. if the differential equation has the form

D1Dnx(t)=A11A12A1nA22A2nAnnx(t)(3.7)

where, using the notation Iμ for the μ-dimensional unit matrix,

Dj=DαjIdj,

Ajk ∈ ℂdj × dk and x = (x1, …, xd) with d=j=1ndj: In this case,

• all solutions of the system (3.7) converge to zero as t → ∞ if and only if, for all j = 1,2,…,n, all eigenvalues λjk, k = 1,2,…, dj, of the matrix Ajj satisfy |arg λjk| > αjπ /2, and

• whenever there exist some j ∈ {1, 2, …, n} and k ∈ {1, 2, …, dj} with |arg λjk| < αjπ /2, there exists an initial value whose corresponding solution is unbounded.

A close inspection of the proof of Theorem 3.1 reveals that the statement of its part (i) can easily be extended to cover a class of inhomogeneous problems:

#### Corollary 3.3

Consider the differential equation system

Dαixi(t)=j=idaijxj(t)+gi(t),1id,(3.8)

where, for all i = 1,2, …, d, the functionsgi : [0, ∞) → ℂ are continuous and satisfy

limtgi(t)=0.

Every solution of the inhomogeneous system(3.8)converges to zero at infinity if and only if all solutions of the associated homogeneous system(3.1a)tend to zero as t→∞, i.e. if and only if condition(3.4)is satisfied.

#### Proof

Assume first that every solution of (3.8) tends to zero as t → ∞. In order to prove that every solution of the corresponding homogeneous system (3.1a) converges to zero, we choose an arbitrary x0 ∈ ℂd. It is then sufficient to show that the solution of (3.1a) that starts at x0 converges to zero as t → ∞. To this end, we take the solutions φ(⋅,x0) and φ(⋅,0) of (3.8) that start at x0 and at 0, respectively. By assumption, both these functions tend to 0 as t → ∞. Thus, φ(t, x0) − φ(t, 0) tends to 0 as t → ∞ as well. But clearly, this difference is identical to the solution of the homogeneous system (3.1a) that starts at x0.

Regarding the proof of the other direction of the equivalence, we assume that the condition (3.4) is satisfied. Under this hypothesis, we may proceed as in the first part of the proof of Theorem 3.1(i). Using the argumentation via Lemma A4(i) employed in the induction step there, we can derive that φd(t, x0) → 0 as t → ∞ for any x0 ∈ ℂd. Then we can proceed inductively as in the first part of the proof of Theorem 3.1(i) and demonstrate that the other components of φ(⋅, x0) vanish near infinity as well. The proof is complete. □

#### Remark 3.4

Clearly, the same arguments can be used to extend the statement of Remark 3.2 regarding block triangular systems to the inhomogeneous case as well.

#### 3.2 Systems with general coefficient matrices

With respect to the stability theory for such systems of equations with general (not necessarily triangular or block triangular) coefficient matrices, we are not yet in a position to provide a comprehensive theory. We can, however, develop an approach that works under certain restrictions on the orders of the differential operators involved. Specifically we shall assume that αj ∈ (0, 1] for all j and that there exists some α* ∈ (0, 1] and some ρj ∈ ℚ such that αj = ρjα*.

In this case, there exist positive integers pj and qj (j = 1, 2, …, d) such that, for all j, gcd(pj, qj) = 1 and ρj = pj / qj. Then we define q to be the least common multiple of the qj. This allows us to deduce that for every j there exists some positive integer rj such that αj = α*rj / q (clearly, rj = pjq / qj). According to [7, Theorem 8.1], we can then rewrite the j-th equation of the original system (1.1) as an equivalent system of rj differential equations of order α* / q. Thus, the entire system (1.1) can be expressed as a system of d* = j=1drj equations of order α* / q. This new system has the form

Dα/qx(t)=Ax(t)+g(t)(3.9a)

where the matrix A* has the block structure

with matrices Ajk ∈ ℂrj×rk given by

Ajj=010000100001ajj000 for j=1,2,,d,(3.9c)

and

Ajk=0000000ajk00 for j,k=1,2,,d and jk,(3.9d)

and with the vector g* being defined by

g(t)=(0,,0r11 times,g1(t),0,,0r21 times,g2(t),0,,0rd1 times,gd(t)).(3.9e)

While the dimension d* of this new system is potentially very much larger than the dimension d of the original system, thus substantially increasing the complexity, we obtain a significant advantage because all equations of the system now have the same order, so that we may invoke the well known classical theory to investigate the asymptotic behavior of solutions of the system. Specifically, in view of this construction, we can immediately deduce from [7, Theorem 8.1]:

#### Theorem 3.5

Let the function g : [0, ∞) → ℂdbe continuous and satisfy g(t) → 0 for t → ∞. Moreover, assume that αj ∈ (0, 1] for all j and that there exist some α* ∈ (0, 1] and some ρj ∈ ℚ such that αj = ρjα*for all j. Then, all solutions of the original differential equation system(1.1)converge to zero at infinity if the eigenvaluesλjof the associated system’s coefficient matrix A*defined in eqs. (3.9b), (3.9c)and(3.9d)satisfy | arg λj | > πα* / (2q) for all j, where q is the least common multiple of the denominators of the ρj.

#### Proof

From [7, Theorem 8.1], we see that the systems (1.1) and (3.9a) are equivalent. Hence, we only concentrate on the system (3.9a). By changing variable x* = Ty, where T is the non-singular matrix which transforms A* into a Jordan normal form B, the system (3.9a) becomes

Dα/qy(t)=By(t)+g^(t),(3.10)

where B = T−1AT = diag(B1, …, Bj, …, Bs) where Bj is the Jordan block corresponding the eigenvalue λj of the matrix A* and g^ = T−1g*. Note that limt→∞g^(t) = 0. Now, using the same arguments as in the proof of Theorem 3.1 and Corollary 3.3, we see that every solution of the system (3.9a) tends to zero if and only if the eigenvalues λj of the associated system’s coefficient matrix A* satisfy | arg λj| > πα* / (2q) for all j. The proof is complete.□

Unfortunately, this criterion is based on the new system’s coefficient matrix A*, and thus it only indirectly makes use of the coefficients of the original matrix A. It would be useful to have a formulation that allows to directly draw such a conclusion from the original matrix without having to explicitly form the much larger new matrix and to compute its eigenvalues. However, the following example indicates that we can probably not expect to find a simple criterion that permits to immediately decide the question for the solution asymptotics for a given differential equation system.

#### Example 3.6

Consider the system

(D1/2x1(t)D1/4x2(t))=Ax(t), where A=(a11a12a21a22)=(0.0000110.00220.1).(3.11)

Following the development above, we may choose α* = 1 and q = 4 in this example, and thus this two-dimensional system can be rewritten as a three-dimensional system of order α*/q = 1/4 in the form

D1/4x(t)=Ax(t) with A=(0100.00001010.002200.1).(3.12)

The components x1andx3 of the solution to this new system are then identical to the two components x1 and x2, respectively, of the original system’s solution. The eigenvalues of A* are λ1 = −0.103917 and λ2/3 = 0.101958 ± 0.10385i so that arg λ1=πand|argλ2|=|argλ3|=0.79459>π/8 = πα* /(2q). Therefore, Theorem 3.5 asserts that all solutions of the system given in eq. (3.11) tend to zero at infinity.

However, this observation does not appear to be immediately deducible from the original matrix A. By a simple calculation, we see that the eigenvalues of this matrix are λ1 = 0.0673111 and λ2 = 0.0326989 and thus arg λ1 = arg λ2 = 0 — a property that one would normally associate with a system for which, in particular, unbounded solutions must be expected.

Similarly, the diagonal entries of A are real and positive as well, so their arguments are zero too. Thus, an argumentation based on the diagonal entries and not the eigenvalues like the one that we had shown to be valid for triangular systems in Subsection 3.1 is not directly applicable to the case of a general (non-triangular) coefficient matrix either.

This seemingly negative observation is not the final word though. Using different techniques we may actually derive a strategy that allows to investigate the stability question in a satisfactory manner at least for the case of a homogeneous system. Specifically, from the proof of Theorem 2.6 we see that all solutions of the homogeneous multi-order system (2.5) are exponentially bounded. (This essentially follows from the generalized power series representation of the solution components and the estimate (2.11) for the coefficients of these series.) Hence, we may take the Laplace transform on both sides of this system. This leads to

sαiXi(s)sαi1xi(0)=j=1daijXj(s),i=1,,d,(3.13)

where Xi(s) is the Laplace transform of the i-th component xi(t) of the solution x(t). The system (3.13) can be rewritten in the form

Δ(s)X1(s)X2(s)Xd(s)=b1(s)b2(s)bd(s)(3.14a)

where

bi(s)=sαi1xi(0),i=1,,d,

and

Using a standard result from the Laplace transform based stability theory [4, Theorem 1], we immediately obtain the following criterion on the asymptotic behavior of the system (2.5):

#### Theorem 3.7

Consider the homogeneous multi-order system(2.5)and let the function Δ be defined as in(3.14b). If all the roots of the characteristic equation det Δ(s) = 0 have negative real parts, then all solutions of the system(2.5)converge to zero at infinity.

#### Remark 3.8

In the triangular case considered in Subsection 3.1, we were able to extend our results derived for homogeneous equations also to the inhomgeneous case, cf. Corollary 3.3. This was possible mainly because the triangular structure allowed us to handle the individual equations of the given system in a step-by-step manner one at a time which made it possible to employ the variation-of-constants formula that is available for scalar equations or single-order systems. In the general case considered here, a suitable generalization of the variation-of-constants formula to the setting of multi-order systems is not readily available and does not appear to be straightforward to derive. The authors plan to address this question in a future work.

Dedicated to Professor Virginia Kiryakova on the occasion of her 65th birthday and the 20th anniversary of FCAA

### Appendix A. Auxiliary results

In this appendix we collect some auxiliary results that we used in the proofs of our theorems above. For the formulation of these auxiliary results we shall use the notation

Λαs:=λC{0}:|arg(λ)|>απ2

and

Λαu:=λC{0}:|arg(λ)|<απ2,

where the superscripts “s” and “u” can be interpreted as “stable region” and “unstable region”, respectively. We note that the lemmas below can be interpreted as generalizations of some results provided in [3] where similar statements have been derived under more restrictive assumptions on the parameter λ.

### Lemma A.1

Letλbe an arbitrary complex number andα ∈ (0, 1]. There exists a positive real number m(α, λ) such that for everyt > 0 the following estimates hold:

1. IfλΛαu, then

|Eα(λtα)1αexp(λ1/αt)|m(α,λ)min{tα,1},|tα1Eα,α(λtα)1αλ1/α1exp(λ1/αt)|m(α,λ)min{t1α,t1+α}.
2. If λΛαs, then

tα1Eα,α(λtα)m(α,λ)min{t1α,t1+α}.

### Proof

In the case α = 1 the results are trivially true because then Eα(z) = Eα, α(z) = exp(z). We therefore only have to deal with the case 0 < α < 1 explicitly.

Let us start with the case 0 < t ≤ 1. In this case, the minimum in the first claim of (i) has the value 1. Thus, this claim is an immediate consequence of the fact that the expression on its left-hand side is a continuous function of t ∈ [0, 1]. Similarly, we can see — in view of the continuity of the Mittag-Leffler functions and the exponentials on [0, 1] — that the expressions on the left-hand sides of the two other claims can be bounded by O(t− 1 + α) = O(min {t− 1 − α, t− 1 + α}).

The statements for t > 1 (where the minima are always attained by the first expression in the braces) immediately follow from well-known results about the asymptotic behavior of Mittag-Leffler functions; specifically, we have (cf., e.g., [Proposition 3.6 and Theorem 4.3] or [20, Theorems 1.3 and 1.4]) that

Eα,β(z)=1αz(1β)/αexp(z1/α)k=1pzkΓ(βαk)+O(|z|p1) for zΛαu(A.1)

and

Eα,β(z)=k=1pzkΓ(βαk)+O(|z|p1) for zΛαs,(A.2)

hold for arbitrary p ∈ ℕ and |z| → ∞. Upon choosing t > 0 and z := λtα, we then observe that the relation zΛαs holds if and only if λΛαs, and an analog equivalence exists for Λαu. Using this approach, the first statement of (i) follows from eq. (A.1) with p = 1. Similarly, the second statement of (i) and the statement of (ii) follow from eqs. (A.1) and (A.2), respectively, upon setting p = 2 and noticing that the summands for k = 1 vanish because they contain a factor 1/ Γ(αα) = 1/Γ(0) = 0.□

### Lemma A.2

Letλ ∈ ℂ ∖ {0} andα ∈ (0, 1]. There exists a positive constantK(α, λ) such that for allt ≥ 1 the following estimates hold:

1. IfλΛαu, then

tλ1/α1Eα(λtα)exp(λ1/ατ)dτK(α,λ),0t(tτ)α1Eα,α(λ(tτ)α)λ1/α1Eα(λtα)exp(λ1/ατ)dτK(α,λ).
2. IfλΛαs, then

0t(tτ)α1Eα,α(λ(tτ)α)dτK(α,λ).

### Proof

Once again the statements are trivially true for α = 1. The proof of the remaining cases is very similar to the proof of [3, Lemma 5].

For the first claim of part (i), the first statement of Lemma A.1(i) allows us to proceed as follows:

tλ1/α1Eα(λtα)exp(λ1/ατ)dτ|λ|1/α1t1αexp(λ1/αt)+m(α,λ)tαexp(λ1/ατ)dτ=|λ|1/α11αtexp(λ1/α(tτ))dτ+m(α,λ)tαtexp(λ1/ατ)dτ.

For the evaluation of these integrals we recall that λΛαu, and hence |arg λ1/α| < π/2 which implies that ℜ λ1/α > 0. Making use of this inequality in combination with the identity |exp(λ1/αz)| = exp(ℜ λ1/αz) for z ∈ ℝ, we conclude

texp(λ1/α(tτ))dτ=0exp(λ1/αu)du=1λ1/α

and

texp(λ1/ατ)dτ=texp(λ1/ατ)dτ<0exp(λ1/ατ)dτ=1λ1/α.

These estimates conclude this part of the proof.

The proof of the second claim of part (i) uses the second statement of Lemma A.1(i). Specifically, that result allows us to write

0t(tτ)α1Eα,α(λ(tτ)α)λ1/α1Eα(λtα)exp(λ1/ατ)dτ0t1αλ1/α1exp(λ1/α(tτ))λ1/α1Eα(λtα)exp(λ1/ατ)dτ+m(α,λ)0tmin{(tτ)1α,(tτ)1+α}dτ.(A.3)

Since we have assumed that t ≥ 1, we may bound the last integral as follows:

0tmin{(tτ)1α,(tτ)1+α}dτ=0tmin{τ1α,τ1+α}dτ=01τ1+αdτ+1tτ1αdτ=1α+1αtα1=2α1αtα<2α.(A.4)

Moreover, for the first integral on the right-hand side of eq. (A.3) we may invoke the first statement of Lemma A.1(i) and conclude that

0t1αλ1/α1exp(λ1/α(tτ))λ1/α1Eα(λtα)exp(λ1/ατ)dτ=|λ|1/α11αexp(λ1/αt)Eα(λtα)0t|exp(λ1/ατ)|dτ|λ|1/α1m(α,λ)tα0t|exp(λ1/ατ)|dτ=|λ|1/α1m(α,λ)tα0texp(λ1/ατ)dτ=|λ|1/α1m(α,λ)tα1λ1/α1exp(λ1/αt).

As above, our assumption that λΛαu implies that ℜ λ1/α > 0, and hence this last expression is uniformly bounded for all t ≥ 1. This completes the proof of the second statement of part (i).

Finally, for part (ii), Lemma A.1(ii) and the fact that t ≥ 1 allow us to estimate as follows:

0t(tτ)α1Eα,α(λ(tτ)α)dτ=0tτα1Eα,α(λτα)dτm(α,λ)0tmin{τ1α,τ1+α}dτ<m(α,λ)2α,

where the last estimate uses the result (A.4). Thus the desired result follows.□

### Lemma A.3

For any continuous and bounded function f: [0, ∞) → ℂ, α ∈ (0, 1] andλΛαu, we have

limt0t(tτ)α1Eα,α(λ(tτ)α)Eα(λtα)f(τ)dτ=λ1/α10exp(λ1/ατ)f(τ)dτ.(A.5)

### Proof

Again, the case α = 1 is trivial.

For 0 < α < 1, we first remark that the expression on the left-hand side of eq. (A.5) is well defined: The denominator is non-zero because, as shown by Wiman [24, pp.225–226], the Mittag-Leffler function Eα does not have any zeros in Λαu. Thus, since λΛαu implies that tλΛαu for all t > 0, we conclude that Eα(λtα) ≠ 0 for all t > 0.

Next we note that the integral on the right-hand side of eq. (A.5) exists because f is assumed to be continuous (which asserts the existence of the integral over any compact subinterval [0, T] with arbitrary T > 0) and bounded which admits us to bound the absolute value of the integrand by

|exp(λ1/ατ)f(τ)|exp(λ1/ατ)supt0|f(t)|.

As we already noted in earlier proofs, ℜλ1/α > 0, and hence this bound provides a convergent majorant for the integral over [0, ∞), thus asserting the existence and finiteness of the improper integral on the right-hand side of eq. (A.5).

Then, the first statement of Lemma A.1(i) implies that |Eα(λtα)| exhibits an unbounded growth as t → ∞ and hence that

limt1αexp(λ1/αt)Eα(λtα)=1.

It thus follows that

limt0t(tτ)α1Eα,α(λ(tτ)α)Eα(λtα)f(τ)dτ=limtα0t(tτ)α1Eα,α(λ(tτ)α)exp(λ1/αt)f(τ)dτ,

if one of the limits exists (which immediately implies the existence of the other one).

For t > 1 we see that

αt1t(tτ)α1Eα,α(λ(tτ)α)f(τ)dτsupu0|f(u)|sup0u1|Eα,α(λuα)|α01uα1du=supu0|f(u)|sup0u1|Eα,α(λuα)|.

Evidently, the upper bound depends on f, α and λ but not on t. It therefore follows, once again using the unbounded growth of |exp(λ1/αt)| for t → ∞, that

limtαt1t(tτ)α1Eα,α(λ(tτ)α)exp(λ1/αt)f(τ)dτ=0.

In order to complete the proof of Lemma A.3, it therefore suffices to show that

limtα0t1(tτ)α1Eα,α(λ(tτ)α)exp(λ1/αt)f(τ)dτ=λ1/α10exp(λ1/ατ)f(τ)dτ.(A.6)

To this end, we recall that the second statement of Lemma A.1(i) implies

0t1α(tτ)α1Eα,α(λ(tτ)α)λ1/α1exp(λ1/α(tτ))exp(λ1/αt)f(τ)dτsupu0|f(u)|0t1αm(α,λ)(tτ)1αexp(λ1/αt)dτsupu0|f(u)|αm(α,λ)|exp(λ1/αt)|1tτ1αdτsupu0|f(u)|m(α,λ)|exp(λ1/αt)|

for t > 1; in particular we once again see that the upper bound converges to zero as t → ∞, and therefore (A.6) follows as desired.□

Using Lemma A.1, Lemma A.2 and Lemma A.3, we obtain the asymptotic behavior of solutions to scalar linear fractional differential equations as described in the following result.

### Lemma A.4

Letα ∈ (0, 1], and let f : [0, ∞) → ℂ be a continuous function with the property limt → ∞ |f(t)| = 0. Consider the differential equation

Dαx(t)=λx(t)+f(t),t>0.(A.7)

The following statements hold:

1. If |arg(λ)| > απ/2 then all solutions of(A.7)tend to zero ast → ∞.

2. If |arg(λ)| < απ/2 then eq. (A.7)has a unique bounded solution. Moreover, this solution tends to zero ast → ∞.

### Proof

In either case, we start from the variation of constants formula [7, Theorem 7.2 and Remark 7.1] which tells us that the solution φ(⋅, x0) of (A.7) that satisfies the condition φ(0, x0) = x0 is given by

φ(t,x0)=x0Eα(λtα)+0t(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ.(A.8)

In order to prove part (i), let ε > 0 be arbitrarily small. We can find a constant T > 0 such that |f(t)| < ε for all tT. For t > T + 1 and x0 ∈ ℂ, we split up the integral on the right-hand side of eq. (A.8) according to

φ(t,x0)=x0Eα(λtα)+0T(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ=x0Eα(λtα)+Tt1(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ=x0Eα(λtα)+t1t(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ.

By virtue of Lemma A.1 (ii), we have

limtx0Eα(λtα)=0.(A.9)

On the other hand, by a simple computation, we obtain

|Tt1(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ|ε1tT|τα1Eα,α(λτα)|dτεm(α,λ)α,(A.10)

due to Lemma A.1(ii), and

|t1t(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ|ε01|τα1Eα,α(λτα)|dτεEα,α+1(|λ|)(A.11)

(see [20, eq. (1.99)]). Furthermore,

|0T(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ|supt0|f(t)|tTt|τα1Eα,α(λτα)|dτm(α,λ)supt0|f(t)|α(tT)α,(A.12)

due to Lemma A.1(ii). Since ε is arbitrarily small, from eqs. (A.9), ( A.10), ( A.11) and ( A.12), we get

limt|φ(t,x0)|=0,

and the proof of part (i) is complete.

For the proof of (ii), we note that Lemma A.3 admits us to precisely describe the asymptotic behavior of the integral on the right-hand side of eq. (A.8), namely

0t(tτ)α1Eα,α(λ(tτ)α)f(τ)dτ=Eα(λtα)λ1/α10exp(λ1/ατ)f(τ)dτ(1+o(1)).,

Thus, by (A.8), any solution to the differential equation behaves as

φ(t,x0)=Eα(λtα)[x0+λ1/α10exp(λ1/ατ)f(τ)dτ(1+o(1))](A.13)

for t → ∞. Since |arg λ| < απ/2, we know that Eα (λtα) is unbounded as t → ∞. Thus, a necessary condition for the entire expression on the right-hand side of (A.13) to be bounded is that the term in brackets converges to zero as t → ∞. Clearly, this is the case if and only if

x0=x¯0:=λ1/α10exp(λ1/ατ)f(τ)dτ.

Thus, the differential equation ( A.7) has at most one bounded solution, and it remains to prove that this solution has the property φ(t, x0) → 0 as t → ∞ (which, in particular, implies that the solution is bounded and hence that a bounded solution exists).

To this end, let ε > 0 be an arbitrary positive real number. Then there exists a positive constant T > 0 such that

|f(t)|ε for all tT.(A.14)

For any tT + 1, we put

H1(t)=Eα(λtα)λ1/α1texp(λ1/ατ)f(τ)dτ,H2(t)=0T[(tτ)α1Eα,α(λ(tτ)α)λ1/α1exp(λ1/ατ)Eα(λtα)]f(τ)dτ,H3(t)=Tt[(tτ)α1Eα,α(λ(tτ)α)λ1/α1exp(λ1/ατ)Eα(λtα)]f(τ)dτ.

It is then clear from eq. (A.8) and the definition of x0 that

φ(t,x¯0)=H1(t)+H2(t)+H3(t).

By virtue of (A.14) and the first statement of Lemma A.2(i), we have

|H1(t)|εK(α,λ).(A.15)

Using both statements of Lemma A.1(i), we obtain, since tT ≥ 1,

|H2(t)|supt0|f(t)|××0T[|λ|1/α1|exp(λ1/α(tτ))αEα(λtα)exp(λ1/ατ)|+m(α,λ)(tτ)1+α]dτsupt0|f(t)|[|λ|1/α10T|exp(λ1/ατ)||1αexp(λ1/αt)Eα(λtα)|dτsupt0|f(t)|[+m(α,λ)0Tdτ(tτ)1+α]m(α,λ)supt0|f(t)|××[|λ|1/α1tα0T|exp(λ1/ατ)|dτ+(tT)αtαα].(A.16)

Since λΛαu, we conclude once again that

0T|exp(λ1/ατ)|dτ=0Texp(λ1/ατ)dτ=1λ1/α[1exp(λ1/αT)]1λ1/α,

and thus we see from eq. (A.16) that

H2(t)0 as t.(A.17)

Furthermore, by (A.14) and the second statement of Lemma A.2(i), we have

|H3(t)|εK(α,λ).(A.18)

From (A.15), (A.17), (A.18) and the fact that ε > 0 can be made arbitrarily small, we conclude

limtφ(t,x¯0)=0.

The proof is complete.□

## Acknowledgement

The work of H.T. Tuan is supported by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under Grant No. 101.03-2017.01.

#### References

[1] M.H. Atabakzadeh, M.H. Akrami, G.H. Erjaee, Chebyshev operational matrix method for solving multi-order fractional ordinary differential equations. Appl. Math. Model. 37 (2013), 8903–8911.10.1016/j.apm.2013.04.019Search in Google Scholar

[2] D. Baleanu, K. Diethelm, E. Scalas, J.J. Trujillo, Fractional Calculus: Models and Numerical Methods. 2nd edition, World Scientific, Singapore, 2016.10.1142/10044Search in Google Scholar

[3] N.D. Cong, T.S. Doan, S. Siegmund, H.T. Tuan, On stable manifolds for planar fractional differential equations. Appl. Math. Comput. 226 (2014), 157–168.10.1016/j.amc.2013.10.010Search in Google Scholar

[4] W. Deng, C. Li, J. Lu, Stability analysis of linear fractional differential system with multiple time delays. Nonlinear Dyn. 48 (2007), 409–416.10.1007/s11071-006-9094-0Search in Google Scholar

[5] W. Deng, C. Li, Q. Guo, Analysis of fractional differential equations with multi-orders. Fractals15 (2007), 1–10.10.1142/S0218348X07003472Search in Google Scholar

[6] K. Diethelm, Multi-term fractional differential equations, multi-order fractional differential systems and their numerical solution. J. Européen des Systèmes Automatisés42 (2008), 665–676.10.3166/jesa.42.665-676Search in Google Scholar

[7] K. Diethelm, The Analysis of Fractional Differential Equations. Springer-Verlag, Berlin, 2010.10.1007/978-3-642-14574-2Search in Google Scholar

[8] K. Diethelm, A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlin. Dynamics71 (2013), 613–619.10.1007/s11071-012-0475-2Search in Google Scholar

[9] K. Diethelm, Properties of the solutions to “fractionalized’' ODE systems, with applications to processes arising in the life sciences. In: D.T. Spasic, N. Grahovac, M. Zigic, M. Rapaic, T.M. Atanackovic (Eds.): Proc. Internat. Conf. Fractional Differentiation and its Applications 2016, Vol. 1. Faculty of Technical Sciences, Novi Sad (2016), 32–44.Search in Google Scholar

[10] K. Diethelm, N.J. Ford, Multi-order fractional differential equations and their numerical solution. Appl. Math. and Comput. 154 (2004), 621–640.10.1016/S0096-3003(03)00739-2Search in Google Scholar

[11] A. El-Mesiry, A. El-Sayed, H. El-Saka, Numerical methods for multi-term fractional (arbitrary) orders differential equations. Comput. and Appl. Math. 160 (2005), 683–699.10.1016/j.amc.2003.11.026Search in Google Scholar

[12] V. Gejji, H. Jafari, Solving a multi-order fractional differential equation using Adomian decomposition. Appl. Math. and Comput. 189 (2007), 541–548.10.1016/j.amc.2006.11.129Search in Google Scholar

[13] R. Gorenflo, A.A. Kilbas, F. Mainardi, S.V. Rogosin, Mittag-Leffler Functions, Related Topics and Applications. Springer-Verlag, Berlin, 2014.10.1007/978-3-662-43930-2Search in Google Scholar

[14] E. Hesameddini, E. Asadollahifard, Numerical solution of multi-order fractional differential equations via the sinc collocation method. Iranian J. of Numer. Anal. and Optimiz. 5 (2015), 37–48.Search in Google Scholar

[15] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations. North-Holland Mathematics Studies, Vol. 204. Elsevier, Amsterdam, 2006.10.1016/S0304-0208(06)80001-0Search in Google Scholar

[16] C. Li, F. Zhang, J. Kurths, F. Zeng, Equivalent system for a multiple-rational-order fractional differential system. Phil. Trans. R. Soc. A371 (2013), 371: 20120156.10.1098/rsta.2012.0156Search in Google Scholar PubMed

[17] Y. Li, Solving a nonlinear fractional differential equation using Chebyshev wavelets. Commun. in Nonlin. Sci. and Numer. Simul. 15 (2010), 2284–2292.10.1016/j.cnsns.2009.09.020Search in Google Scholar

[18] Y. Liu, Existence of solutions of IVPs for differential systems on half line with sequential fractional derivative operators. African Diaspora J. of Math. (New Ser.)18 (2015), 27–54.Search in Google Scholar

[19] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations. John Wiley & Sons, New York, 1993.Search in Google Scholar

[20] 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. Academic Press, Inc., San Diego, 1999.Search in Google Scholar

[21] T. Sardar, S.S. Ray, R.K. Bera, B.B. Biswas, The analytical approximate solution of the multi-term fractionally damped van der Pol equation. Physica Scripta80, No 2 (2009), # 025003.10.1088/0031-8949/80/02/025003Search in Google Scholar

[22] N.H. Sweilam, M.M. Khader, R.F. Al-Bar, Numerical studies for a multi-order fractional differential equation. Phys. Letters A371 (2007), 26–33.10.1016/j.physleta.2007.06.016Search in Google Scholar

[23] G. Vainikko, Which functions are fractionally differentiable? Zeitschrift für Anal. und ihre Anwendungen35 (2016), 465–487.10.4171/ZAA/1574Search in Google Scholar

[24] A. Wiman, Über die Nullstellen der Funktionen Eα(x). Acta Math. 29 (1905), 217–234.10.1007/BF02403204Search in Google Scholar

[25] S. Yang, A. Xiao, H. Su, Convergence of the variational iteration method for solving multi-order fractional differential equations. Computers and Math. with Appl. 60 (2010), 2871–2879.10.1016/j.camwa.2010.09.044Search in Google Scholar

[26] Y. Yang, Solving a nonlinear multi-order fractional differential equation using Legendre pseudo-spectral method. Appl. Math. 4 (2013), 113–118.10.4236/am.2013.41020Search in Google Scholar