*Review* **Stability Properties of Multi-Term Fractional-Differential Equations**

**Oana Brandibur 1,† and Éva Kaslik 1,2,\*,†**


**Abstract:** Necessary and sufficient stability and instability conditions are reviewed and extended for multi-term homogeneous linear fractional differential equations with Caputo derivatives and constant coefficients. A comprehensive review of the state of the art regarding the stability analysis of two-term and three-term fractional-order differential equations is provided, which is then extended to the case of four-term fractional-order differential equations. The stability and instability properties are characterized with respect to the coefficients of the multi-term fractional differential equations, leading to both fractional-order-dependent and fractional-order-independent characterizations. In the general case, fractional-order-independent stability and instability properties are described for fractional-order differential equations with an arbitrary number of fractional derivatives.

**Keywords:** multi-order fractional differential equation; stability; instability; Caputo derivative

### **1. Introduction**

Fractional-order differential equations have been extensively used in modelling realworld phenomena from a wide range of domains [1–5], considering the memory and hereditary properties that the fractional derivatives are provided with [1,5]. Two-term fractional differential equations describe the modelling of the motion of a rigid plate that is immersed in a viscous liquid, also known as the Bagley–Torvik equation [6], whereas the Basset equation [7] is also described using linear two-term fractional differential equation. Other examples worth mentioning are the equation of the inextendible pendulum [8] and the fractional harmonic oscillator [9,10], which are both described by multi-term fractionalorder differential equations with three fractional derivatives of Caputo type.

The most well-known types of fractional derivatives, which are typically not equivalent, are the Grünwald–Letnikov derivative, the Riemann–Liouville derivative, and the Caputo derivative. Given that it only requires initial conditions defined in terms of integerorder derivatives, which represent well-understood characteristics of physical phenomena, the Caputo derivative seems to be more applicable to real-world problems.

The existence and uniqueness of solutions for multi-term fractional differential equations have been explored the past years [11–13]. Stability theory has been broadly explored as well. Two-term fractional-order differential equations have been investigated in terms of their stability properties. Refs. [14–18] studied the boundary-value problems for multiterm fractional-differential equations. The asymptotic behaviour of solutions of linear multi-order fractional differential systems are determined in Ref. [19], whereas Ref. [20] summarize the recent stability results of fractional differential equations, most of them in comparison to their classical integer-order counterparts.

The stability of linear systems of fractional-order differential equations is also considered in Ref. [21], where the fractional orders are proportional and transformations form the physical plane to the fractional-domain is discussed. Moreover, recent results [22]

**Citation:** Brandibur, O.; Kaslik, E. Stability Properties of Multi-Term Fractional-Differential Equations. *Fractal Fract.* **2023**, *7*, 117. https:// doi.org/10.3390/fractalfract7020117

Academic Editors: Angelo B. Mingarelli, Leila Gholizadeh Zivlaei and Mohammad Dehghan

Received: 20 December 2022 Revised: 22 January 2023 Accepted: 24 January 2023 Published: 26 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

generalize the formulation, the existence of the solution, and provide stability properties for multi-term systems of fractional-order differential equations. Some stability theorems are also obtained in Ref. [23] for both systems of fractional differential equations with multi-order and for multi-term fractional differential equations.

The purpose of this paper is to investigate multi-term fractional-order equations in terms of their stability and instability properties, both with respect to the roots of the corresponding characteristic equations and with respect to the coefficients of the considered equations. Moreover, we obtain sufficient conditions for the asymptotic stability and instability of the fractional-order differential equation, independent of the choice of the considered fractional orders. Furthermore, we also determine the stability and instability conditions for multi-term fractional differential equations with four fractional-order derivatives of Caputo type, by extending the recent results obtained for multi-term fractionalorder equations with three derivatives [24], considering the analogy of the corresponding characteristic equations.

The structure of this paper is as follows. The preliminaries regarding basic definitions in fractional calculus are included in Section 2. The main results obtained for fractional-order differential equations with an arbitrary number of derivatives, as well as fractional-order-independent stability and instability results, are provided in Section 3. As a particular case, we recall the main stability and instability conditions for two-term fractional differential equations described in Section 4. These are followed by an overview of the main results recently obtained in the case of three-term fractional-order differential equations, which are listed in Section 5. The next section aims to present and prove the main novel results of this work, by extending the stability and instability conditions to the case of multi-term fractional differential equations with four derivatives of Caputo type, in Section 6. Both fractional-order dependent and independent stability and instability conditions are obtained. The conclusions are drawn in Section 7.

#### **2. Preliminaries**

**Definition 1.** *The Caputo fractional differential operator of order q* > 0 *is defined as*

$${}^{c}D^{q}x(t) := \begin{cases} \frac{1}{\Gamma(n-q)} \int\_{0}^{t} \frac{\mathfrak{x}^{(n)}(\tau)}{(t-\tau)^{q+1-n}} d\tau, & \text{if } n-1 < q < n\\ \frac{d^{n}\mathfrak{x}(t)}{dt^{n}}, & \text{if } q = n \end{cases}$$

*where n* <sup>∈</sup> <sup>N</sup> *and the Gamma function is* <sup>Γ</sup>(*z*) = <sup>∞</sup> 0 *yz*−1*e* <sup>−</sup>*ydy, for* (*z*) <sup>&</sup>gt; <sup>0</sup>*.*

The Laplace transform of the Caputo derivative or an arbitrary fractional order *q* is given as in Ref. [12]:

**Proposition 1.** *The Laplace transform for the fractional-order Caputo derivative of order q* ∈ (*n* − 1, *n*]*, n* ∈ N∗*, of a function x is:*

$$
\mathcal{L}(^cD^q \mathfrak{x})(s) = s^q X(s) - \sum\_{k=0}^{n-1} s^{q-k-1} \mathfrak{x}^{(k)}(0),
$$

*where X*(*s*) *represents the Laplace transform of the function x.*

Consider the following general linear multi-term fractional differential equation

$$\,^cD^q\mathbf{x}(t) + \sum\_{i=1}^{n-1} \boldsymbol{\kappa}\_i \,^cD^{q\_i}\mathbf{x}(t) + \beta \mathbf{x}(t) = \mathbf{0} \tag{1}$$

where *αi*, *β* are real coefficients, for any *i* ∈ {1, 2, ... , *n* − 1}, and *q*, *q*1, *q*2, ... , *qn*−<sup>1</sup> are the fractional orders of the Caputo derivatives, such that 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < ... < *qn*−<sup>1</sup> < *q* ≤ 2.

The linear homogeneous fractional-order differential Equation (1) can be seen as the linearization at the trivial equilibrium of a nonlinear autonomous fractional-order differential equation of the following form

$$\,^cD^qx(t) = f(x(t), ^cD^{q\_1}x(t), ^cD^{q\_2}x(t), \dots, ^cD^{q\_{n-1}}x(t))\tag{2}$$

where *<sup>f</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup> is a continuously differentiable function such that *<sup>f</sup>*(0, 0, . . . , 0) = 0. The initial condition associated to Equation (1) is either


In what follows, *ϕ*(*t*, *xi*) represents the unique solution of (1), which satisfies one of the initial considered above. One can prove the existence and uniqueness of the solution of the initial value problem associated to the equation (1) in a similar way to the case of fractional-order systems of differential equations [8,11]. As a matter of fact, we can express the solution of the linear constant coefficient fractional differential Equation (1) using the fractional Green's function as presented in Ref. [12], or by employing the fractional meta-trigonometric approach in the commensurate case [25].

Henceforward, we consider a non-exponential asymptotic stability concept, named Mittag–Leffler stability [26], in particular O(*t* <sup>−</sup>*q*)-asymptotic stability, as it reflects the algebraic decay of the solution. This is justified by the fact that fractional derivatives have a memory effect and hereditary properties and therefore the asymptotic stability of the trivial solution of equation (1) is not of exponential type [15,27].

#### **Definition 2.**


$$\left\| \,\varphi(t, \mathbf{x}\_0) \right\| = \mathcal{O}(t^{-q}) \quad \text{as } t \to \infty.$$

#### **3. General Fractional-Order-Independent Stability and Instability Results**

In what follows, we determine the characteristic equation associated to the multiterm fractional-order differential Equation (1), which we will analyse in order to obtain fractional-order-independent stability and instability results.

Applying the Laplace transform, Equation (1) becomes

$$\left(\mathbf{s}^{\eta} + \sum\_{i=1}^{n-1} \alpha\_i \mathbf{s}^{\eta\_i} + \boldsymbol{\beta}\right) \mathbf{X}(\mathbf{s}) = F(\mathbf{s})\_{\boldsymbol{\epsilon}}$$

where *X*(*s*) represents the Laplace transform of the function *x*, *sq*<sup>∗</sup> represents the first branch of the complex power function [28], with *q*<sup>∗</sup> ∈ {*q*, *q*1, *q*2, ... , *qn*−1} and *F*(*s*) incorporates the initial conditions.

Thus, we obtain the characteristic equation

$$s^q + \sum\_{i=1}^{n-1} \mathfrak{a}\_i \mathfrak{s}^{q\_i} + \beta = 0. \tag{3}$$

Our aim is to analyse the distribution of the roots of Equation (3). Therefore, we consider the complex valued function

$$
\Delta(\mathbf{s}; \mathfrak{a}\_1, \dots, \mathfrak{a}\_{n-1}, \mathfrak{b}, \mathfrak{q}\_1, \dots, \mathfrak{q}\_{n-1}, \mathfrak{q}) = \mathbf{s}^q + \sum\_{i=1}^{n-1} \mathfrak{a}\_i \mathbf{s}^{q\_i} + \beta.
$$

The next result gives us the necessary and sufficient conditions regarding the asymptotic stability and instability of the trivial solution of Equation (1) with respect to the sign of the real parts of the roots of its associated characteristic in Equation (3). A similar result has been formulated for the case of three-term linear fractional-order differential equations [24].

#### **Theorem 1.**


In what follows, let us consider the set

$$D = \{ (q\_1, q\_2, \dots, q\_{n-1}, q) \in \mathbb{R}^n \, : \, 0 < q\_1 < q\_2 < \dots < q\_{n-1} < q \le 2 \}.$$

The following lemma provides a sufficient condition for the instability of the Equation (1) with respect to its coefficients, regardless of the fractional orders (*q*1, *q*2, ... , *qn*−1, *q*) ∈ *D*.

**Lemma 1.** *If β* < 0 *or α*<sup>1</sup> + *α*<sup>2</sup> + ... + *αn*−<sup>1</sup> + *β* + 1 ≤ 0*, the trivial solution of Equation* (1) *is unstable, regardless of the fractional orders q*1*, q*2*,...,qn*−<sup>1</sup> *and q.*

**Proof.** If *α*<sup>1</sup> + *α*<sup>2</sup> + ... + *αn*−<sup>1</sup> + *β* + 1 ≤ 0 we have that

$$
\Delta(1; \mathfrak{a}\_1, \dots, \mathfrak{a}\_{n-1}, \mathfrak{b}, \mathfrak{q}\_1, \dots, \mathfrak{q}\_{n-1}, \mathfrak{q}) = 1 + \mathfrak{a}\_1 + \mathfrak{a}\_2 + \dots + \mathfrak{a}\_{n-1} + \mathfrak{b} \le 0.
$$

For *β* < 0, it is obvious that

$$
\Delta(0; \alpha\_1, \dots, \alpha\_{n-1}, \beta, q\_1, \dots, q\_{n-1}, q) = \beta < 0.
$$

On the other hand, we notice that

Δ(*s*; *α*1,..., *αn*−1, *β*, *q*1,..., *qn*−1, *q*) → ∞ for *s* → ∞.

Hence, for both cases the function *s* → Δ(*s*; *α*1, ... , *αn*−1, *β*, *q*1, ... , *qn*−1, *q*) has at least one strictly positive real root. Therefore, Equation (1) is unstable, regardless of the fractional orders *q*1, *q*2,. . . , *qn*−<sup>1</sup> and *q*.

In what follows, we only consider the case *β* > 0, as we have previously seen in Lemma 1 that if *β* < 0, the Equation (1) is unstable, regardless of the fractional orders (*q*1, *q*2, ... , *qn*−1, *q*) ∈ *D*. Therefore, we can state the following result, which provides a sufficient condition for the asymptotic stability of the trivial solution of Equation (1), regardless of its fractional orders.

**Lemma 2.** *If α<sup>i</sup>* > 0*, i* ∈ {1, 2, ... , *n* − 1} *and β* > 0*, the trivial solution of Equation* (1) *is asymptotically stable, for any fractional orders q*1*, q*2*,...,qn*−<sup>1</sup> *and q.*

**Proof.** Let *α<sup>i</sup>* > 0, *i* ∈ {1, 2, ... , *n* − 1} and *β* > 0. Assuming the contrary, that there exists *s*<sup>0</sup> = 0 a root of the characteristic function Δ(*s*; *α*1, ... , *αn*−1, *β*, *q*1, ... , *qn*−1, *q*) such that (*s*0) ≥ 0, we have

$$s\_0^{q} + \alpha\_1 s\_0^{q\_1} + \alpha\_2 s\_0^{q\_2} + \dots + \alpha\_{n-1} s\_0^{q\_{n-1}} + \beta = 0.1$$

Multiplying this previous relation by *s* − *q* 2 <sup>0</sup> , it follows that

$$s\_0^{\frac{q}{2}} + \alpha\_1 s\_0^{q\_1 - \frac{q}{2}} + \alpha\_2 s\_0^{q\_2 - \frac{q}{2}} + \dots + \alpha\_{n-1} s\_0^{q\_{n-1} - \frac{q}{2}} + \beta s\_0^{-\frac{q}{2}} = 0.$$

Since <sup>±</sup>*<sup>q</sup>* <sup>2</sup> <sup>∈</sup> [−1, 1] and *<sup>q</sup>*<sup>1</sup> <sup>−</sup> *<sup>q</sup>* <sup>2</sup> , *<sup>q</sup>*<sup>2</sup> <sup>−</sup> *<sup>q</sup>* <sup>2</sup> , ... , *qn*−<sup>1</sup> <sup>−</sup> *<sup>q</sup>* <sup>2</sup> ∈ (−1, 1), we have that the real parts of each term from the left-hand side of the equality from above are positive. Therefore,

$$\Re \left( s\_0^{\frac{q}{2}} + \alpha\_1 s\_0^{q\_1 - \frac{q}{2}} + \alpha\_2 s\_0^{q\_2 - \frac{q}{2}} + \dots + \alpha\_{n-1} s\_0^{q\_{n-1} - \frac{q}{2}} + \beta s\_0^{-\frac{q}{2}} \right) > 0,$$

which is a contradiction.

Thus, (*s*) < 0, for any *s* being a root of the characteristic function Δ, which shows that the trivial solution of Equation (1) is asymptotically stable, regardless of the fractional orders *q*1, *q*2,. . . , *qn*−<sup>1</sup> and *q*.

It is important to emphasize that Lemmas 1 and 2 generalize the recent results from Ref. [24], which strictly refer to three-term fractional differential equations. The techniques employed in the proofs presented above allow for a simple generalization to the multi-term case, leading to easily verifiable inequalities involving the constant coefficients, which guarantee the instability/stability of the trivial solution of Equation (1), for any choice of the fractional orders of the Caputo derivatives.

In the following sections, we explore fractional-order-dependent stability and instability conditions for particular cases of multi-term fractional-order differential equations. We will first review some results from the literature for the cases of two-term and three-term fractional-order differential equations, followed by an extension of the results to four-term fractional differential equations

#### **4. Two-Term Fractional-Order Differential Equations**

In this section, we review some stability properties of two-term fractional differential equations similar to those obtained by Cermák and Kisela [ ˘ 15], as a particular case of the multi-term fractional differential Equation (1).

Consider the following two-term fractional differential equation

$$
\hbar^c \mathbf{D}^q \mathbf{x}(t) + \mathbf{a}^c \mathbf{D}^{q\_1} \mathbf{x}(t) + \beta \mathbf{x}(t) = \mathbf{0},\tag{4}
$$

where the coefficients *α* and *β* are the real numbers and *q*<sup>1</sup> and *q* are the fractional orders of the Caputo derivatives, with 0 < *q*<sup>1</sup> < *q* ≤ 2.

The associated characteristic equation for the two-term fractional differential Equation (4) is

$$s^q + \mathfrak{a}s^{q\_1} + \mathcal{B} = 0.\tag{5}$$

In what follows, we give a characterisation of the stability properties of Equation (4), which may be obtained as a particular case of Theorem 2. For the proof of this theorem we refer to Ref. [24].

**Proposition 2.** *Let* 0 < *q*<sup>1</sup> < *q* ≤ 2*. For β* > 0*, let*

$$\alpha^\* = -\frac{\beta^{1-\frac{q\_1}{q}} \sin \frac{q\pi}{2}}{\left(\sin \frac{q\_1\pi}{2}\right)^{\frac{q\_1}{q}} \left(\sin \frac{(q-q\_1)\pi}{2}\right)^{1-\frac{q\_1}{q}}}.$$

#### *Then*


**Remark 1.** *The previous proposition gives us sufficient conditions of instability of the considered equation regardless of the fractional orders q*<sup>1</sup> *and q and sufficient conditions of both asymptotic stability and instability dependent of the fractional orders of Equation* (4)*. We notice that the results obtained are in accordance to the conditions of stability and instability determined by Cermak and Kisela in Ref. [15], where particular rational fractional orders were considered.*

In what follows, we will recall some well-known examples of two-term fractional-order differential equations, which were comprehensively described in Ref. [24].

**Example 1.** *The Basset equation is described by the following linear fractional-order differential equation*

$$
\dot{\mathbf{x}}(t) + \mathbf{a}^c \mathbf{D}^{q\_1} \mathbf{x}(t) + \beta \mathbf{x}(t) = \mathbf{0}, \quad \mathbf{0} < q\_1 < 1,\tag{6}
$$

*which originates in the study of the generalised Basset force ensuing when a spherical object is submerged in an incompressible viscous fluid [7].*

*It is clear that the characteristic equation associated to* (6) *is*

$$s + \alpha s^{q\_1} + \beta = 0.$$

*Based on the theoretical results previously mentioned we notice that:*


$$\alpha > \alpha^\*(\beta, q\_1) = -\beta^{1-q\_1} \left( \cot \frac{q\_1 \pi}{2} \right)^{q\_1} \sec \frac{q\_1 \pi}{2}.$$

**Example 2.** *The Bagley–Torvik equation is described by following fractional-order differential equation:*

$$\mathbf{x}^{\prime\prime}(t) + \mathbf{a}^{\varepsilon} \mathbf{D}^{q\_1} \mathbf{x}(t) + \beta \mathbf{x}(t) = \mathbf{0}, \quad \mathbf{0} < q\_1 < \mathbf{2} \tag{7}$$

*which arises when modelling the motion of a rigid plate that immerses in a viscous liquid [6]. The characteristic equation associated to Equation* (7) *is*

$$s^2 + \mathfrak{a}s^{q\_1} + \beta = 0.$$

*The critical value for the parameter α is α*∗(*β*, *q*1) = 0*, for any fractional order q*<sup>1</sup> *and any β* > 0*. The trivial solution of Equation* (7) *will then be asymptotically stable of order* O(*t* −{*q*1}) *if and only if α* > 0 *and β* > 0*, which is in accordance to the results obtained in Ref. [29].*

#### **5. Three-Term Fractional-Differential Equations**

The results presented in this section were rigorously investigated in Ref. [24] and will be enumerated in the following.

We consider the following multi-term fractional differential equation

$${}^{c}D^{q}\mathbf{x}(t) + a\_{1}^{c}D^{q\_{1}}\mathbf{x}(t) + a\_{2}^{c}D^{q\_{2}}\mathbf{x}(t) + \beta \mathbf{x}(t) = \mathbf{0} \tag{8}$$

where *α*1, *α*2, *β* are the real numbers and *q*, *q*1, *q*<sup>2</sup> are the fractional orders of the Caputo derivatives, with 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < *q* ≤ 2.

The characteristic equations becomes

$$
\Delta(s; \mathfrak{a}\_1, \mathfrak{a}\_2, \mathfrak{b}, q\_1, q\_2, q) := s^q + \mathfrak{a}\_1 s^{q\_1} + \mathfrak{a}\_2 s^{q\_2} + \mathcal{B} = 0. \tag{9}
$$

Let *<sup>D</sup>* <sup>=</sup> {(*q*1, *<sup>q</sup>*2, *<sup>q</sup>*) <sup>∈</sup> <sup>R</sup><sup>3</sup> : 0 <sup>&</sup>lt; *<sup>q</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>q</sup>*<sup>2</sup> <sup>&</sup>lt; *<sup>q</sup>* <sup>≤</sup> <sup>2</sup>} and let *<sup>β</sup>* <sup>&</sup>gt; 0.

**Lemma 3.** *Let* (*q*1, *q*2, *q*) ∈ *D and β* > 0 *be arbitrarily fixed. We consider the smooth parametric curve in the* (*α*1, *α*2)*-plane defined by*

$$\Gamma(\beta, q\_1, q\_2, q) \;:\quad \begin{cases} \alpha\_1 = \beta^{1 - \frac{q\_1}{q}} h(\omega, q\_1, q\_2, q) \\ \alpha\_2 = \beta^{1 - \frac{q\_2}{q}} h(\omega, q\_2, q\_1, q) \end{cases} \; \omega > 0,$$

*where h* : (0, ∞) × *D* → (0, ∞) *is given by:*

$$h(\omega, q\_1, q\_2, q) = \omega^{-\frac{q\_1}{q}} \left[ \omega \rho(q - q\_2, q\_2 - q\_1) - \rho(q\_2, q\_2 - q\_1) \right]$$

*with the function <sup>ρ</sup> defined as <sup>ρ</sup>*(*a*, *<sup>b</sup>*) = sin *<sup>a</sup><sup>π</sup>* 2 sin *<sup>b</sup><sup>π</sup>* 2 *, for any a* ∈ [0, 2], *b* ∈ [−1, 0) ∪ (0, 1].

*The following statements hold:*


**Theorem 2** (Fractional-order-dependent results)**.**

*Let β* > 0*,* 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < *q* ≤ 2 *be arbitrarily fixed. Consider the curve* Γ(*β*, *q*1, *q*2, *q*) *and the function φβ*,*q*1,*q*2,*<sup>q</sup>* : R → R *previously defined.*


$$
\alpha\_2 > \phi\_{\beta, q\_1, q\_2, q}(a\_1).
$$

*iii. If α*<sup>2</sup> < *φβ*,*q*1,*q*2,*q*(*α*1)*, the trivial solution of equation (8) is unstable.*

**Remark 2.** *If q*<sup>1</sup> = *q*<sup>2</sup> =: *q*∗*,* Γ(*β*, *q*1, *q*2, *q*) *becomes the line*

$$\alpha\_1 + \alpha + 2 = -\frac{\beta^{1 - \frac{q^\*}{q}} \sin \frac{q \pi}{2}}{\left(\sin \frac{q^\* \pi}{2}\right)^{\frac{q^\*}{q}} \left(\sin \frac{(q - q^\*) \pi}{2}\right)^{1 - \frac{q^\*}{q}}},$$

*which is in accordance with the results from Proposition 2.*

**Remark 3.** *We emphasize that in the particular case of multi-term fractional-order differential equations with three derivatives, the fractional-order independent asymptotic stability and instability conditions described by Lemmas 1 and 2 become necessary and sufficient conditions, compared to the general case of fractional-order differential equations of n fractional derivatives.*

**Remark 4.** *In Figure 1, we have plotted three parametric curves* Γ(*β*, *q*1, *q*2, *q*) *for β* = 3 *and different values of the fractional orders q*1, *q*2, *q such that* 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < *q* ≤ 2 *in the* (*α*1, *α*2) *plane, along with the fractional-order-independent stability and instability regions, respectively (plotted with red and blue). In fact, we have proved in Ref. [24] that the reunion of all the parametric*

*curves* Γ(*β*, *q*1, *q*2, *q*) *fill the space between the fractional-order-independent stability and instability regions and will not intersect neither of them.*

**Figure 1.** Parametric curves Γ(*β*, *q*1, *q*2, *q*) (black) for *β* = 3 and several values of the fractional orders *q*1, *q*<sup>2</sup> and *q*. The blue and red regions represent the *fractional-order-independent* instability and asymptotic stability regions, given by Lemma 1 and Lemma 2, respectively.

Similar to the case of two-term fractional-order differential equations described in Section 3, we will next enumerate some examples of practical applications described by three-term fractional-order differential equations, which were presented in detail in [24].

**Example 3.** *The equation of motion of the inextendible pendulum [8] is described by the following three term fractional-differential equation*

$$
\phi^{\prime\prime} + \mu \tau^{q\_1} \cdot ^c D^{q\_1} \phi + \nu \tau^{q\_2} \cdot ^c D^{q\_2} \phi + \frac{\mathcal{g}}{L} \sin \phi = 0. \tag{10}
$$

*The trivial equilibrium of the equation is asymptotically stable, for any fractional orders q*<sup>1</sup> *and q*2*. Moreover, Theorem 1 gives us that*

$$\phi(t) = \mathcal{O}(t^{-q'}) \quad \text{as } t \to \infty, \quad \text{where } q' = \min\{\{q\_1\}, \{q\_2\}\}.$$

**Example 4.** *The fractional harmonic oscillator is modelled by the following fractional differential equation:*

$$\mathbf{x}^{\prime\prime} + \mathbf{a}\_1^{\varepsilon} \mathbf{D}^{q\_1} \mathbf{x} + \mathbf{a}\_2 \mathbf{x}^{\prime} + \beta \mathbf{x} = \mathbf{0}, \quad \text{with } q\_1 \in (0, 1) \tag{11}$$

*where α*<sup>1</sup> *and α*<sup>2</sup> *represent the friction coefficient and viscous damping coefficient [9,10], and β* > 0*.*


#### **6. Four-Term Fractional-Differential Equations**

In this section, our main objective is to generalize the previously presented results of Theorem 2, by considering a linear fractional differential equations with four Caputo derivatives:

$$\mathbf{u}^{\varepsilon}\mathbf{D}^{q}\mathbf{x}(t) + a\_{1}^{\varepsilon}\mathbf{D}^{q\_{1}}\mathbf{x}(t) + a\_{2}^{\varepsilon}\mathbf{D}^{q\_{2}}\mathbf{x}(t) + a\_{3}^{\varepsilon}\mathbf{D}^{q\_{3}}\mathbf{x}(t) + \beta\mathbf{x}(t) = \mathbf{0},\tag{12}$$

where *α*1, *α*2, *α*3, *β* are the real numbers and *q*, *q*1, *q*2, *q*<sup>3</sup> are the fractional orders of the Caputo derivatives, with 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < *q*<sup>3</sup> < *q* ≤ 2.

By means of the Laplace transform technique, we obtain the characteristic equation

$$
\Delta(\mathbf{s}; \mathfrak{a}\_1, \mathfrak{a}\_2, \mathfrak{a}\_3, \mathfrak{b}, q\_1, q\_2, q\_3, q) = \mathbf{s}^q + \mathfrak{a}\_1 \mathbf{s}^{q\_1} + \mathfrak{a}\_2 \mathbf{s}^{q\_2} + \mathfrak{a}\_3 \mathbf{s}^{q\_3} + \mathcal{J} = 0. \tag{13}
$$

Next, we will assume that the fractional orders *q*1, *q*2, *q*3, *q* are arbitrarily fixed inside the domain

$$D = \{ (q\_1, q\_2, q\_3, q) \in \mathbb{R}^4 \, : \, 0 < q\_1 < q\_2 < q\_3 < q \le 2 \}.$$

Moreover, as *β* < 0 implies that the Equation (12) is unstable (according to Lemma 1), for any choice of the fractional orders *q*1, *q*2, *q*3, *q*, we will further assume that *β* > 0.

**Lemma 4.** *Let* (*q*1, *q*2, *q*3, *q*) ∈ *D and β* > 0 *be arbitrarily fixed. Consider the parametric surface in the* (*α*1, *α*2, *α*3)*-space defined by*

$$S(\beta, q\_1, q\_2, q\_3, q) \quad : \quad \begin{cases} \alpha\_1 = \beta^{1 - \frac{q\_1}{q}} h(\omega, \eta, q\_1, q\_2, q\_3, q) \\ \alpha\_2 = \beta^{1 - \frac{q\_2}{q}} h(\omega, \eta, q\_2, q\_1, q\_3, q) \\ \alpha\_3 = \beta^{1 - \frac{q\_3}{q}} \eta \omega^{-\frac{q\_3}{q}} \end{cases} \quad \omega > 0, \eta \in \mathbb{R}$$

*where h* : (0, ∞) × R × *D* → (0, ∞) *is given by:*

$$\ln h(\omega, \eta, q\_1, q\_2, q\_3, q) = \omega^{-\frac{q\_1}{q}} [\eta \rho(q\_3 - q\_2, q\_2 - q\_1) + \omega \rho(q - q\_2, q\_2 - q\_1) - \rho(q\_2, q\_2 - q\_1)]$$

*with the function ρ defined as*

$$\rho(a,b) = \frac{\sin\frac{a\pi}{2}}{\sin\frac{b\pi}{2}} \quad \text{, } \forall a \in [0,2], \ b \in [-1,0) \cup (0,1].$$

*The surface S*(*β*, *q*1, *q*2, *q*3, *q*) *lies outside the first octant (α*<sup>1</sup> > 0, *α*<sup>2</sup> > 0, *α*<sup>3</sup> > 0*) of the* (*α*1, *α*2, *α*3) *space.*

**Proof.** Assuming the contrary, that there exist *ω* > 0 and *η* ∈ R such that *α*<sup>1</sup> > 0, *α*<sup>2</sup> > 0 and *α*<sup>3</sup> > 0, it follows that *η* > 0 and

$$\begin{cases} & \eta \rho (q\_3 - q\_2, q\_2 - q\_1) + \omega \rho (q - q\_2, q\_2 - q\_1) - \rho (q\_2, q\_2 - q\_1) > 0 \\ & \eta \rho (q\_3 - q\_1, q\_1 - q\_2) + \omega \rho (q - q\_1, q\_1 - q\_2) - \rho (q\_1, q\_1 - q\_2) > 0 \end{cases}$$

which is equivalent to

$$\begin{cases} & \eta \cdot \sin(q\_3 - q\_2)\frac{\pi}{2} + \omega \cdot \sin(q - q\_2)\frac{\pi}{2} - \sin\frac{q\_2\pi}{2} > 0\\ & -\eta \cdot \sin(q\_3 - q\_1)\frac{\pi}{2} - \omega \cdot \sin(q - q\_1)\frac{\pi}{2} + \sin\frac{q\_1\pi}{2} > 0 \end{cases}$$

Eliminating *ω*, it leads to the following inequality:

$$\begin{aligned} &\eta \left( \sin(q\_3 - q\_2) \frac{\pi}{2} \sin(q - q\_1) \frac{\pi}{2} - \sin(q\_3 - q\_1) \frac{\pi}{2} \sin(q - q\_2) \frac{\pi}{2} \right) > 0\\ &> \sin \frac{q\_2 \pi}{2} \sin(q - q\_1) \frac{\pi}{2} - \sin \frac{q\_1 \pi}{2} \sin(q - q\_2) \frac{\pi}{2} .\end{aligned}$$

By applying elementary trigonometric identities, the previous relation is equivalent to the inequality:

$$
\eta \sin(q\_3 - q) \frac{\pi}{2} > \sin \frac{q \pi}{2},
$$

which is absurd, as *<sup>η</sup>* <sup>&</sup>gt; 0, sin(*q*<sup>3</sup> <sup>−</sup> *<sup>q</sup>*) *<sup>π</sup>* <sup>2</sup> <sup>&</sup>lt; 0 and sin *<sup>q</sup><sup>π</sup>* <sup>2</sup> > 0.

Thus, the surface *S*(*β*, *q*1, *q*2, *q*3, *q*) does not intersect the first octant of the (*α*1, *α*2, *α*3) space.

**Remark 5.** *If q*<sup>1</sup> = *q*<sup>2</sup> = *q*<sup>3</sup> =: *q*∗*, S*(*β*, *q*1, *q*2, *q*3, *q*) *represents the plane*

$$\alpha\_1 + \alpha\_2 + \alpha\_3 = -\frac{\beta^{1 - \frac{q^\*}{q}} \sin \frac{q \pi}{2}}{\left(\sin \frac{q^\* \pi}{2}\right)^{\frac{q^\*}{q}} \left(\sin \frac{(q - q^\*) \pi}{2}\right)^{1 - \frac{q^\*}{q}}}.$$

**Lemma 5.** *Let β* > 0*,* (*q*1, *q*2, *q*3, *q*) ∈ *D be arbitrarily fixed. Consider the surface S*(*β*, *q*1, *q*2, *q*3, *q*) *defined in Lemma 4. The characteristic Equation* (13) *has a pair of complex conjugated roots on the imaginary axis of the complex plane if and only if* (*α*1, *α*2, *α*3) ∈ *S*(*β*, *q*1, *q*2, *q*3, *q*)*.*

**Proof.** The characteristic Equation (13) has a pair of pure imaginary roots if and only if there exists *ω* > 0 such that

$$
\Delta(s; \alpha\_1, \alpha\_2, \alpha\_3, \beta, q\_1, q\_2, q\_3, q) = 0, \quad \text{where } s = i(\beta \omega)^{\frac{1}{q}},
$$

which is equivalent to

$$i^q \beta \omega + \mathfrak{a}\_1 i^{q\_1} \left(\beta \omega \right)^{\frac{q\_1}{q}} + \mathfrak{a}\_2 i^{q\_2} \left(\beta \omega \right)^{\frac{q\_2}{q}} + \mathfrak{a}\_3 i^{q\_3} \left(\beta \omega \right)^{\frac{q\_3}{q}} + \delta = 0.1$$

Dividing the previous relation by *β*, we have:

$$i^{\mathfrak{g}}\omega + \varkappa\_1 \beta^{\frac{q\_1}{q}-1} \omega^{\frac{q\_1}{q}} i^{q\_1} + \varkappa\_2 \beta^{\frac{q\_2}{q}-1} \omega^{\frac{q\_2}{q}} i^{q\_2} + \varkappa\_3 \beta^{\frac{q\_3}{q}-1} \omega^{\frac{q\_3}{q}} i^{q\_3} + 1 = 0.1$$

Taking the real and imaginary part from the relation from above, we obtain the following system:

$$\begin{cases} \omega \cos\frac{q\pi t}{2} + a\_1 \beta^{\frac{q\_1}{q}-1} \omega^{\frac{q\_1}{q}} \cos\frac{q\_1\pi t}{2} + a\_2 \beta^{\frac{q\_2}{q}-1} \omega^{\frac{q\_2}{q}} \cos\frac{q\_2\pi t}{2} + a\_3 \beta^{\frac{q\_3}{q}-1} \omega^{\frac{q\_3}{q}} \cos\frac{q\_3\pi t}{2} + 1 = 0 \\\\ \omega \sin\frac{q\pi t}{2} + a\_1 \beta^{\frac{q\_1}{q}-1} \omega^{\frac{q\_1}{q}} \sin\frac{q\_1\pi t}{2} + a\_2 \beta^{\frac{q\_2}{q}-1} \omega^{\frac{q\_2}{q}} \sin\frac{q\_2\pi t}{2} + a\_3 \beta^{\frac{q\_3}{q}-1} \omega^{\frac{q\_3}{q}} \sin\frac{q\_3\pi t}{2} = 0 \end{cases}$$

We solve the previous system for *α*<sup>1</sup> and *α*<sup>2</sup> and, denoting *η* = *α*3*ω q*3 *<sup>q</sup> δ q*3 *<sup>q</sup>* −1 , it follows that the characteristic Equation (13) has a pair of complex conjugated roots on the imaginary axis if and only if the triplet (*α*1, *α*2, *α*3) belongs to the parametric surface *S*(*β*, *q*1, *q*2, *q*3, *q*) defined in Lemma 4.

#### **Theorem 3** (Fractional-order-dependent results)**.**

*Let β* > 0*,* 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < *q*<sup>3</sup> < *q* ≤ 2 *be arbitrarily fixed. Consider the surface S*(*β*, *q*1, *q*2, *q*3, *q*) *previously defined. The trivial solution of Equation (12) is* O(*t* −*q* )*-asymptotically stable (where q* = min{{*q*}, {*q*1}, {*q*2}, {*q*3}}*) if and only if* (*α*1, *α*2, *α*3) *belongs to the connected component of the three-dimensional space bounded by the surface S*(*β*, *q*1, *q*2, *q*3, *q*)*, which contains the positive octant.*

**Proof.** If *α<sup>i</sup>* > 0, for *i* = 1, 3, it follows from Lemma 2 that the trivial solution of Equation (12) is asymptotically stable. As the parameters *αi*, *i* = 1, 3 are varied, the qualitative behaviour may change if and only if (*α*1, *α*2, *α*3) ∈ *S*(*β*, *q*1, *q*2, *q*3, *q*), which follows from Lemma 5. Therefore, it is clear that if (*α*1, *α*2, *α*3) belongs to the open connected region, which contains the positive octant, then the trivial solution of Equation (12) is asymptotically stable.

In what follows, we verify the transversality condition.

Let *β* > 0, (*q*1, *q*2, *q*3, *q*) ∈ *D* be arbitrarily fixed. Consider the surface *S*(*β*, *q*1, *q*2, *q*3, *q*) defined in Lemma 4. For simplicity, denote *Qj* <sup>=</sup> *qj q* , with *j* ∈ {1, 2, 3}.

We first compute the cross product of the vectors *∂α*<sup>1</sup> *∂ω* , *∂α*<sup>2</sup> *∂ω* , *∂α*<sup>3</sup> *∂ω* and *∂α*<sup>1</sup> *∂η* , *∂α*<sup>2</sup> *∂η* , *∂α*<sup>3</sup> *∂η* , which is the normal vector to the surface *S*(*β*, *q*1, *q*2, *q*3, *q*). *N* = *∂α*<sup>1</sup> *∂ω* , *∂α*<sup>2</sup> *∂ω* , *∂α*<sup>3</sup> *∂ω* × *∂α*<sup>1</sup> *∂η* , *∂α*<sup>2</sup> *∂η* , *∂α*<sup>3</sup> *∂η* = = *∂α*<sup>2</sup> *∂ω* · *∂α*<sup>3</sup> *∂η* <sup>−</sup> *∂α*<sup>2</sup> *∂η* · *∂α*<sup>3</sup> *∂ω* , *∂α*<sup>1</sup> *∂η* · *∂α*<sup>3</sup> *∂ω* <sup>−</sup> *∂α*<sup>1</sup> *∂ω* · *∂α*<sup>3</sup> *∂η* , *∂α*<sup>1</sup> *∂ω* · *∂α*<sup>2</sup> *∂η* − · *∂α*<sup>1</sup> *∂η* · *∂α*<sup>2</sup> *∂ω*

where

$$\begin{aligned} \frac{\partial a\_1}{\partial \omega} &= \beta^{1-Q\_1} \omega^{-Q\_1-1} \left[ -Q\_1 \eta \rho(q\_3 - q\_2, q\_2 - q\_1) - Q\_1 \omega \rho(q - q\_2, q\_2 - q\_1) + 1 \right. \\ &+ Q\_1 \rho(q\_2, q\_2 - q\_1) + \omega \rho(q - q\_2, q\_2 - q\_1)] \\ \frac{\partial a\_2}{\partial \omega} &= \beta^{1-Q\_2} \omega^{-Q\_2 - 1} \left[ -Q\_2 \eta \rho(q\_3 - q\_1, q\_1 - q\_2) - Q\_2 \omega \rho(q - q\_1, q\_1 - q\_2) + 1 \right. \\ &+ Q\_2 \rho(q\_1, q\_1 - q\_2) + \omega \rho(q - q\_1, q\_1 - q\_2)] \\ \frac{\partial a\_3}{\partial \omega} &= -Q\_3 \eta \beta^{1-Q\_3} \omega^{-Q\_3 - 1} \end{aligned}$$

$$\begin{aligned} \frac{\partial \alpha\_1}{\partial \eta} &= \beta^{1-Q\_1} \omega^{-Q\_1} \rho (q\_3 - q\_2, q\_2 - q\_1) \\ \frac{\partial \alpha\_2}{\partial \eta} &= \beta^{1-Q\_2} \omega^{-Q\_2} \rho (q\_3 - q\_1, q\_1 - q\_2) \\ \frac{\partial \alpha\_3}{\partial \eta} &= \beta^{1-Q\_3} \omega^{-Q\_3} \end{aligned}$$

Let *s*(*α*1, *α*2, *α*3, *β*, *q*1, *q*2, *q*3, *q*) denote the root of the characteristic function Δ(*s*; *α*1, *α*2, *α*3, *β*, *q*1, *q*2, *q*3, *q*) which satisfies *s*(*α*<sup>∗</sup> <sup>1</sup>, *α*<sup>∗</sup> <sup>2</sup>, *α*<sup>∗</sup> <sup>3</sup>, *β*, *q*1, *q*2, *q*3, *q*) = *i*(*βω*) 1 *<sup>q</sup>* , with (*α*∗ <sup>1</sup>, *α*<sup>∗</sup> <sup>2</sup>, *α*<sup>∗</sup> <sup>3</sup> ) ∈ *S*(*β*, *q*1, *q*2, *q*3, *q*).

Taking the derivative with respect to *α*<sup>1</sup> in the characteristic equation

$$s^{q} + \mathfrak{a}\_{1}s^{q\_{1}} + \mathfrak{a}\_{2}s^{q\_{2}} + \mathfrak{a}\_{3}s^{q\_{3}} + \beta = 0\_{r}$$

we obtain

$$q s^{q-1} \frac{\partial s}{\partial \alpha\_1} + s^{q\_1} + \alpha\_1 q\_1 s^{q\_1 - 1} \frac{\partial s}{\partial \alpha\_1} + \alpha\_2 q\_2 s^{q\_2 - 1} \frac{\partial s}{\partial \alpha\_1} + \alpha\_3 q\_3 s^{q\_3 - 1} \frac{\partial s}{\partial \alpha\_1} = 0, \quad \alpha\_1 = 0$$

which is equivalent to

$$\frac{\partial s}{\partial \alpha\_1} = \frac{-s^{q\_1}}{q s^{q\_1 - 1} + \alpha\_1 q\_1 s^{q\_1 - 1} + \alpha\_2 q\_2 s^{q\_2 - 1} + \alpha\_3 q\_3 s^{q\_3 - 1}}.$$

Taking the real part of the previous relation, it follows that:

$$\frac{\partial \mathcal{H}(s)}{\partial \mathcal{A}\_1} = \mathcal{H}\left(\frac{\partial s}{\partial \mathcal{a}\_1}\right) = \mathcal{H}\left(\frac{-s^{q\_1}}{q s^{q\_1 - 1} + \kappa\_1 q\_1 s^{q\_1 - 1} + \kappa\_2 q\_2 s^{q\_2 - 1} + \kappa\_3 q\_3 s^{q\_3 - 1}}\right).$$

.

We have:

*∂*(*s*) *∂α*<sup>1</sup> (*α*∗ <sup>1</sup> ,*α*<sup>∗</sup> <sup>2</sup> ,*α*<sup>∗</sup> 3 ) = ⎛ ⎝−(*i*(*βω*) 1 *<sup>q</sup>* )*q*<sup>1</sup> *P*(*i*(*βω*) 1 *q* ) ⎞ ⎠ = −(*βω*) *q*1 *<sup>q</sup>* ⎛ ⎝ *i <sup>q</sup>*1*P*(*i*(*βω*) 1 *q* ) |*P*(*i*(*βω*) 1 *<sup>q</sup>* )|<sup>2</sup> ⎞ ⎠ = <sup>=</sup> <sup>−</sup>(*βω*) *q*1 *q* |*P*(*i*(*βω*) 1 *<sup>q</sup>* )|<sup>2</sup> *i <sup>q</sup>*1*P*(*i*(*βω*) 1 *q* ) . Computing *i <sup>q</sup>*1*P*(*i*(*βω*) 1 *q* ) we obtain: *i <sup>q</sup>*1*P*(*i*(*βω*) 1 *q* ) = *i <sup>q</sup>*1*P*(*i*(*βω*) 1 *q* ) = = *i* <sup>−</sup>*q*<sup>1</sup> (*βω*) − 1 *q* ! *qiq*−1*βω* + *α*<sup>∗</sup> <sup>1</sup> *q*1*i <sup>q</sup>*1−1(*βω*)*Q*<sup>1</sup> + *α*<sup>∗</sup> <sup>2</sup> *q*2*i <sup>q</sup>*2−1(*βω*)*Q*<sup>2</sup> + *α*<sup>∗</sup> <sup>3</sup> *q*3*i <sup>q</sup>*3−1(*βω*)*Q*<sup>3</sup> " = *<sup>q</sup>β*1<sup>−</sup> <sup>1</sup> *<sup>q</sup> <sup>ω</sup>*<sup>−</sup> <sup>1</sup> *q i <sup>q</sup>*−*q*1−1*ω*+*Q*1*i* <sup>−</sup>1*ωQ*<sup>1</sup> *h*(*ω*, *η*, *q*1, *q*2, *q*3, *q*)+ + *Q*2*i <sup>q</sup>*2−*q*−1*ωQ*<sup>2</sup> *h*(*ω*, *η*, *q*2, *q*1, *q*3, *q*)+*Q*3*i <sup>q</sup>*3−*q*1−1*η* 

As the second term from the precious relation is purely imaginary and

$$\mathfrak{R}(i^{q-q\_1-1}) = \sin\frac{(q-q\_1)\pi}{2}, \quad \mathfrak{R}(i^{q\_2-q\_1-1}) = \sin\frac{(q\_2-q\_1)\pi}{2}, \quad \mathfrak{R}(i^{q\_3-q\_1-1}) = \sin\frac{(q\_3-q\_1)\pi}{2}, \quad \mathfrak{R}(i^{q\_4-q\_1}) = \sin\frac{(q\_4-q\_1)\pi}{2}, \quad \mathfrak{R}(i^{q\_5-q\_1}) = \sin\frac{(q\_5-q\_1)\pi}{2}$$
  $\text{ it follows that}$ 

$$\begin{split} \Re\left(\overline{i^{q\_1}}P(i(\mathcal{G}\omega)^{\frac{1}{q}})\right) &= q\beta^{1-\frac{1}{q}}\omega^{-\frac{1}{q}}\Big(\omega\sin(q-q\_1)\frac{\pi}{2} + Q\_2\sin(q\_2-q\_1)\frac{\pi}{2}\Big[\eta\rho(q\_3-q\_1,q\_1-q\_2) + \eta\rho(q\_3-q\_2)\Big] \\ &+ \omega\rho(q-q\_1,q\_1-q\_2) - \rho(q\_1,q\_1-q\_2)\Big] + Q\_3\eta\sin(q\_3-q\_1)\frac{\pi}{2} \end{split}$$

A simple computation leads to

$$\frac{\partial \mathfrak{R}(s)}{\partial \mathfrak{a}\_{1}}\Big|\_{\left(\mathfrak{a}\_{1}^{\tau},\mathfrak{a}\_{2}^{\tau},\mathfrak{a}\_{3}^{\tau}\right)} = \frac{q\mathcal{J}^{Q\_{1}+Q\_{2}+Q\_{3}-\frac{1}{q}-1}\omega^{Q\_{1}+Q\_{2}+Q\_{3}-\frac{1}{q}+1}}{|P(i(\mathfrak{k}\omega)^{\frac{1}{q}})|^{2}\sin(q\_{2}-q\_{1})\frac{\pi}{2}} \cdot \left(\frac{\partial \mathfrak{a}\_{2}}{\partial \omega}\cdot\frac{\partial \mathfrak{a}\_{3}}{\partial \eta} - \frac{\partial \mathfrak{a}\_{2}}{\partial \eta}\cdot\frac{\partial \mathfrak{a}\_{3}}{\partial \omega}\right).$$

Following similar steps, we can deduce that

$$\frac{\partial \Re(s)}{\partial \mathfrak{a}\_2}\Big|\_{(a\_1^\*,a\_2^\*,a\_3^\*)} = \frac{q\beta^{Q\_1+Q\_2+Q\_3-\frac{1}{q}-1}\omega^{Q\_1+Q\_2+Q\_3-\frac{1}{q}+1}}{|P(i(\beta\omega)^{\frac{1}{q}})|^2 \sin(q\_2-q\_1)\frac{\pi}{2}} \cdot \left(\frac{\partial a\_1}{\partial \eta}\cdot\frac{\partial a\_3}{\partial \omega} - \frac{\partial a\_1}{\partial \omega}\cdot\frac{\partial a\_3}{\partial \eta}\right)|$$

and

$$\frac{\partial \Re(s)}{\partial a\_3}\Big|\_{(a\_1^\*, a\_2^\*, a\_3^\*)} = \frac{q\beta^{Q\_1 + Q\_2 + Q\_3 - \frac{1}{q} - 1}\omega^{Q\_1 + Q\_2 + Q\_3 - \frac{1}{q} + 1}}{|P(i(\beta \omega)^{\frac{1}{q}})|^2 \sin(q\_2 - q\_1)\frac{\pi}{2}} \cdot \left(\frac{\partial a\_1}{\partial \omega} \cdot \frac{\partial a\_2}{\partial \eta} - \frac{\partial a\_1}{\partial \eta} \cdot \frac{\partial a\_2}{\partial \omega}\right).$$

Therefore, we obtain that

$$\begin{split} \nabla \Re(s)(a\_1^\*, a\_2^\*, a\_3^\*) &= \left( \frac{\partial \Re(s)}{\partial a\_1}, \frac{\partial \Re(s)}{\partial a\_2}, \frac{\partial \Re(s)}{\partial a\_3} \right) \bigg|\_{ (a\_1^\*, a\_2^\*, a\_3^\*) } = \\ &= \frac{q \beta^{Q\_1 + Q\_2 + Q\_3 - \frac{1}{7} - 1} \omega^{Q\_1 + Q\_2 + Q\_3 - \frac{1}{7} + 1}}{|P(i(\beta \omega)^{\frac{1}{7}})|^2 \sin(q\_2 - q\_1) \frac{\pi}{2}} \left( \frac{\partial a\_1}{\partial \omega}, \frac{\partial a\_2}{\partial \omega}, \frac{\partial a\_3}{\partial \omega} \right) \times \left( \frac{\partial a\_1}{\partial \eta}, \frac{\partial a\_2}{\partial \eta}, \frac{\partial a\_3}{\partial \eta} \right). \end{split}$$

Thus, the vector ∇(*s*)(*α*<sup>∗</sup> <sup>1</sup>, *α*<sup>∗</sup> <sup>2</sup>, *α*<sup>∗</sup> <sup>3</sup> ) is a normal vector to the surface *S*(*β*, *q*1, *q*2, *q*3, *q*). Moreover, the vector ∇(*s*)(*α*<sup>∗</sup> <sup>1</sup>, *α*<sup>∗</sup> <sup>2</sup>, *α*<sup>∗</sup> <sup>3</sup> ) points away from the open connected region which contains the positive octant and is bounded by the surface *S*(*β*, *q*1, *q*2, *q*3, *q*), as it is perpendicular to the tangent vector *∂α*<sup>1</sup> *∂η* , *∂α*<sup>2</sup> *∂η* , *∂α*<sup>3</sup> *∂η* , where *∂α<sup>i</sup> ∂η* <sup>&</sup>gt; 0, *<sup>i</sup>* <sup>=</sup> 1, 3.

**Remark 6.** *In Figure 2, we have plotted three parametric surfaces S*(*β*, *q*1, *q*2, *q*3, *q*) *for β* = 5 *and different values of the fractional orders q*1, *q*2, *q*3, *q such that* 0 < *q*<sup>1</sup> < *q*<sup>2</sup> < *q*<sup>3</sup> < *q* ≤ 2 *in the* (*α*1, *α*2, *α*3)*-space, along with the fractional-order-independent stability and instability regions, respectively (plotted red and blue). We theorize that the reunion of all the parametric surfaces S*(*β*, *q*1, *q*2, *q*3, *q*) *will fill the space between the fractional-order-independent stability and instability regions and will not intersect neither of them.*

**Figure 2.** Parametric surfaces *S*(*β*, *q*1, *q*2, *q*3, *q*) (yellow) for *β* = 5 and several values of the fractional orders: *q*<sup>1</sup> = 0.4, *q*<sup>2</sup> = 0.8, *q*<sup>3</sup> = 1.5, *q* = 1.8 (left), *q*<sup>1</sup> = 0.2, *q*<sup>2</sup> = 0.3, *q*<sup>3</sup> = 1.2, *q* = 1.5 (centre) and *q*<sup>1</sup> = 0.7, *q*<sup>2</sup> = 0.71, *q*<sup>3</sup> = 0.72, *q* = 1.5 (right). The blue and red regions represent the *fractionalorder-independent* instability and asymptotic stability regions, given by Lemma 1 and Lemma 2, respectively.

**Example 5.** *To exemplify the previous theoretical results, we consider a mathematical model of a damped unforced Duffing-type oscillator expressed as the following multi-term fractional differential equation:*

$$m\mathbf{x}''' + c\_1 \, ^\circ D^{q\_1} \mathbf{x} + c\_2 \, ^\circ D^{q\_2} \mathbf{x} + c\mathbf{x}' + k\_1 \mathbf{x} + k\_2 \mathbf{x}^3 = \mathbf{0}, \quad \text{with } q\_1, q\_2 \in (0, 1) \tag{14}$$

*where m* > 0 *is the mass, k*<sup>1</sup> *and k*<sup>2</sup> *are the linear and non-linear stiffness coefficients, respectively, while c*1*, c*<sup>2</sup> *and c are linear viscous damping coefficients [30]. The linearization of this equation at the trivial equilibrium represents a particular case of Equation* (12) *with q* = 2 *and q*<sup>3</sup> = 1*.*

*By Lemma 2 it follows that if c*<sup>1</sup> > 0*, c*<sup>2</sup> > 0*, c* > 0 *and k*<sup>1</sup> > 0*, then the trivial equilibrium of* (14) *is asymptotically stable, regardless of the choice of the fractional order q*1*. Nevertheless, if a negative damping coefficient c is chosen, the stability of the trivial solution depends on the fractional orders q*<sup>1</sup> *and q*2*.*

*For example, if m* = 1*, c*<sup>1</sup> = 1*, c*<sup>2</sup> = 0.1*, c*<sup>3</sup> = −1*, k*<sup>1</sup> = 4 *and k*<sup>2</sup> = 1*, for q*<sup>2</sup> = 0.95*, then the corresponding critical value of the fractional order q*<sup>1</sup> *is found numerically by solving the parametric equations of the surface defined in Lemma 4 for ω, η and q*1*, leading to q* <sup>1</sup> = 0.885052*. The trivial equilibrium of Equation* (14) *is asymptotically stable if and only if q*<sup>1</sup> > *q* <sup>1</sup>*. The numerical solutions [31,32] of Equation* (14) *with respect to several values of q*<sup>1</sup> *are plotted in Figure 3, showing that when q*<sup>1</sup> < *q* <sup>1</sup>*, the trivial solution of* (14) *becomes unstable and persistent oscillations appear.*

**Figure 3.** Numerical solutions of (14) for *m* = 1, *c*<sup>1</sup> = 1, *c*<sup>2</sup> = 0.1, *c*<sup>3</sup> = −1, *k*<sup>1</sup> = 4 and *k*<sup>2</sup> = 1, for *q*<sup>2</sup> = 0.95 and several values of the fractional order *q*1.

#### **7. Conclusions**

In this paper, on the one hand, we have focused our attention on determining general **fractional-order-independent** stability and instability conditions for multi-term fractionalorder differential equations of an arbitrary number of terms. In fact, the obtained conditions are formulated as easily verifiable inequalities involving the constant coefficients of the linear fractional differential equations.

On the other hand, we also presented a brief overview of **fractional-order-dependent** stability and instability conditions for two-term and three-term fractional-order differential equations, previously reported in the literature. As a generalisation, we provided an investigation of fractional-order dependent stability and instability conditions for four-term fractional-differential equations.

Employing similar techniques as in the proof of Theorem 3, these results can further be generalized to the case of fractional differential equations with multiple Caputo derivatives, leading to a characterisation of the hypersurface, which represents the boundary of the stability region in the hyperspace of coefficients. However, it becomes clear that fractional-order-dependent stability and instability conditions for multi-term fractional differential equations cannot be expressed in a simple way, and as expected, the complexity of the problem increases along with the number of Caputo derivatives included in the fractional differential equation. As a possible direction for future research, our aim is to find a way to overcome this limitation, possibly by finding a different parametrization of the hypersurface at the boundary of the stability region.

Moreover, at this stage, the results presented in this paper represent a timely review of the main contributions related to qualitatively characterizing the stability and instability properties of multi-term fractional-order differential equations with Caputo derivatives. Fractional-order dependent stability and instability properties for multi-order fractional differential equations having an arbitrary number of Caputo derivatives are still yet to be determined. Moreover, we foresee that this investigation shall become more arduous as we increase the number of fractional derivatives, or if we consider different types of fractional derivatives, such as the more general Prabhakar derivatives or derivatives with Sonine kernels.

**Author Contributions:** Conceptualization, O.B. and É.K.; methodology, O.B. and É.K.; software, O.B. and É.K.; validation, O.B. and É.K.; formal analysis, O.B. and É.K.; investigation, O.B. and É.K.; resources, O.B. and É.K.; data curation, O.B. and É.K.; writing—original draft preparation, O.B. and É.K.; writing—review and editing, O.B. and É.K.; visualization, O.B. and É.K.; supervision, O.B. and É.K.; project administration, O.B. and É.K.; funding acquisition, O.B. and É.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
