Stability of fractional-order systems with Prabhakar derivatives

Abstract

Fractional derivatives of Prabhakar type are capturing an increasing interest since their ability to describe anomalous relaxation phenomena (in dielectrics and other fields) showing a simultaneous nonlocal and nonlinear behaviour. In this paper we study the asymptotic stability of systems of differential equations with the Prabhakar derivative, providing an exact characterization of the corresponding stability region. Asymptotic expansions (for small and large arguments) of the solution of linear differential equations of Prabhakar type and a numerical method for nonlinear systems are derived. Numerical experiments are hence presented to validate theoretical findings.

1 Introduction

The Prabhakar function is named after the Indian mathematician Tilak Raj Prabhakar who introduced in 1971 a generalization to three parameters of the Mittag–Leffler function [18] and studied a convolution integral operator with this function as kernel [36].

After their introduction, Prabhakar’s function and integral have been overlooked for a long time until, in the first years of the twenty-first century, the connections with the Havriliak–Negami (HN) dielectric model [21] have been put in light. The HN model was introduced to incorporate the asymmetry and broadness observed in the dielectric dispersion of some polymers and today it is recognized as manifestation of the simultaneous nonlocality and nonlinearity [33, 38] in the response of complex and heterogeneous systems. For these reasons operators based on the Prabhakar function are employed to describe in the time-domain sophisticated relaxation models in several areas (see e.g. [2, 4, 5, 12,13,14,15,16, 19, 22, 30, 37, 41]).

In 2002 the Prabhakar integral was studied in the context of weakly-singular Volterra integral equations and an interpretation in the framework of fractional calculus was provided [23], thus leading 2 years later to the proposition of a left-inverse operator of the Prabhakar fractional integral [24]. A regularization of this inverse, known as the fractional Prabhakar derivative, was introduced in [8], and 1 year later all these preliminary ideas were incorporated in a more general theory [10], successively deepened in [9, 11]. We refer to the recent survey paper [17] for a comprehensive history and collection of background material and applications of the Prabhakar fractional calculus.

Theoretical aspects of the Prabhakar derivative have been studied in a fair number of works. However, there still persist some not completely clear aspects which must be deepened in order to profitably employ the fractional Prabhakar derivative in the analysis and simulation of linear and nonlinear systems.

This paper focuses on the asymptotic stability of fractional-order systems with Prabhakar derivatives. Due to the nonlinear dependence of this derivative on a certain number of parameters, this is a difficult and highly complex task, and there are only a couple of previously published papers which have tackled this issue [1, 6], obtaining some sufficient conditions for the asymptotic stability of linear systems with Prabhakar derivative. Consequently, our aim is to clarify several aspects presented in [1, 6] and to give a rigorous and complete characterization of the stability region of fractional-order systems with Prabhakar derivatives, essentially obtaining a generalization of the well-known Matignon theorem [32] for standard fractional calculus.

Our main result is formulated as a necessary and sufficient condition for the asymptotic stability of a linear autonomous systems with Prabhakar derivatives and an application to the study of nonlinear systems is also provided.

This paper is organized in the following way. Section 2 is devoted to present a short review of the basic material on the Prabhakar function and on the fractional Prabhakar calculus. Section 3 describes the main results concerning stability properties of systems of differential equations with the fractional Prabhakar derivative. A characterization of the corresponding stability region by means of the root locus method is presented in Sect. 4. In Sect. 5 we derive the asymptotic expansion of the solution of linear differential equations with the Prabhakar derivative together with a numerical method for solving nonlinear problems. Finally, some numerical experiments are presented in Sect. 6 with the aim of validating the theoretical results.

2 Preliminaries on Prabhakar function and Prabhakar calculus

Given three real parameters $\alpha$, $\beta$ and $\gamma$, the Prabhakar function is defined by the series representation

$$\begin{aligned} E_{\alpha ,\beta }^{\gamma } (z) = \frac{1}{\varGamma (\gamma )} \sum _{k=0}^{\infty } \frac{\varGamma (\gamma +k) z^{k}}{k! \varGamma (\alpha k + \beta )} , \quad \alpha > 0 , \quad z \in {\mathbb {C}}, \end{aligned}$$

where, as usual, $\varGamma (x) = \int _0^\infty t^{x-1} {\mathrm {e}}^{-t} {\mathrm {d}}t$ is the Euler-Gamma function. This is an entire function of order $\rho =1/\alpha$ and type $\sigma =1$.

More generally, the Prabhakar function is defined for complex parameters, provided that $\mathfrak {R}(\alpha )>0$; in this paper we prefer, however, to focus just on real parameters in view of their wider range of applications.

It is immediate to see that when $\gamma =1$ the function $E_{\alpha ,\beta }^{\gamma } (z)$ reduces to the standard two parameter ML function $E_{\alpha ,\beta }(z)$, when $\beta =\gamma =1$ to the one-parameter ML function $E_{\alpha }(z)$ and when $\alpha =\beta =\gamma =1$ the correspondence with the exponential function ${\mathrm {e}}^{z}$ is obtained. Whenever $\gamma =-j$, with $j \in {\mathbb {N}}$, it is easy to verify that the Prabhakar function is the j-th degree polynomial

$$\begin{aligned} E_{\alpha ,\beta }^{-j}(z) = \sum _{k=0}^{j} (-1)^k\left( {\begin{array}{c}j\\ k\end{array}}\right) \frac{z^{k}}{\varGamma (\alpha k + \beta )} . \end{aligned}$$

Although an analytical representation of the Laplace transform (LT) of $E_{\alpha ,\beta }^{\gamma }(z)$ is not known, it is possible to evaluate the LT of the generalization

$$\begin{aligned} e_{\alpha ,\beta }^{\gamma }(t; \omega ) = t^{\beta -1} E_{\alpha ,\beta }^{\gamma }(t^{\alpha } \omega ) , \quad t > 0, \quad \omega \in {\mathbb {C}}, \end{aligned}$$

which, for $\mathfrak {R}(s)> 0$ and $|s| > |\omega |^{\frac{1}{\alpha }}$, is

$$\begin{aligned} {{\mathcal {E}}}_{\alpha ,\beta }^{\gamma }(s; \omega ) {:}{=}{{\mathcal {L}}} \Bigl ( e_{\alpha ,\beta }^{\gamma }(t; \omega ) \, ; \, s \Bigr ) = \frac{s^{\alpha \gamma - \beta }}{(s^{\alpha } - \omega )^{\gamma }}. \end{aligned}$$

Since the Prabhakar function, and in particular its generalization $e_{\alpha ,\beta }^{\gamma }(t; \omega )$, is employed for the description of relaxation phenomena, it is of importance to identify the range of parameters for which it turns out to be completely monotonic (CM). We recall that a function $f:(0,+\infty ) \rightarrow {\mathbb {R}}$ is CM if it has derivatives of any order $k \in {\mathbb {N}}$ and $(-1)^k f^{(k)}(t) \ge 0$ on $(0,+\infty )$. The CM properties of the Prabhakar function have been studied in [31, 39] and it is possible to prove that $e_{\alpha ,\beta }^{\gamma }(t; \omega )$ is CM if

$$\begin{aligned} \omega<0,\quad 0< \alpha \le 1, \quad 0 < \alpha \gamma \le \beta \le 1. \end{aligned}$$

(1)

The asymptotic behaviour of the Prabhakar function for large arguments in the whole complex plane has been studied in [9, 34, 35]. In particular, for $0<\alpha \le 1$ it is

$$\begin{aligned} E_{\alpha ,\beta }^{\gamma }(z) \sim \left\{ \begin{array}{ll} {{\mathcal {F}}}_{\alpha ,\beta }^{\gamma }(z) + {{\mathcal {A}}}_{\alpha ,\beta }^{\gamma }(z {\mathrm {e}}^{\mp \pi {\mathrm {i}}}) &{} |\arg z |< \frac{\alpha \pi }{2} \\ {{\mathcal {A}}}_{\alpha ,\beta }^{\gamma }(z {\mathrm {e}}^{\mp \pi {\mathrm {i}}}) + {{\mathcal {F}}}_{\alpha ,\beta }^{\gamma }(z) &{} \frac{\alpha \pi }{2}< |\arg z |< \alpha \pi \\ {{\mathcal {A}}}_{\alpha ,\beta }^{\gamma }(z {\mathrm {e}}^{\mp \pi {\mathrm {i}}}) &{} \alpha \pi < |\arg z | \le \pi \\ \end{array}\right. \end{aligned}$$

as $|z| \rightarrow \infty$ and where the sign in ${\mathrm {e}}^{\mp \pi {\mathrm {i}}}$ must be understood as negative for z in the upper complex half-plane and positive otherwise. We have adopted the convention proposed in [17] by which in each sum it is first presented the dominant term. The exponential and algebraic expansions ${{\mathcal {F}}}_{\alpha ,\beta }^{\gamma }(z)$ and ${{\mathcal {A}}}_{\alpha ,\beta }^{\gamma }(z {\mathrm {e}}^{\mp \pi {\mathrm {i}}})$ are, respectively,

$$\begin{aligned} {{\mathcal {F}}}_{\alpha ,\beta }^{\gamma }(z) = \frac{1}{\varGamma (\gamma )} {\mathrm {e}}^{z^{1/\alpha }} z^{\frac{\gamma -\beta }{\alpha }} \frac{1}{\alpha ^{\gamma }} \sum _{k=0}^{\infty } c_{k} z^{-\frac{k}{\alpha }} \end{aligned}$$

and

$$\begin{aligned} {{\mathcal {A}}}_{\alpha ,\beta }^{\gamma }(z) = \frac{z^{-\gamma }}{\varGamma (\gamma )} \sum _{k=0}^{\infty } \frac{(-1)^{k} \varGamma (k+\gamma )}{k! \varGamma (\beta -\alpha (k+\gamma ))} z^{-k} , \end{aligned}$$

where $c_k$ are the coefficients in the inverse factorial expansion of

$$\begin{aligned} F_{\alpha ,\beta }^{\gamma }(s) = \frac{\varGamma (\gamma +s)\varGamma (\alpha s + 1-\gamma +\beta )}{\varGamma (s+1)\varGamma (\alpha s + \beta )} , \end{aligned}$$

(2)

as $|s| \rightarrow \infty$, with $|\arg (s)| \le \pi - \epsilon$ for any arbitrarily small $\epsilon >0$. The first few entries of coefficients $c_k$ are explicitly provided in [35], but they can be evaluated by an algorithm described in [34] and further explained in [9].

For $\alpha ,\beta ,\gamma > 0$ the Prabhakar fractional integral of a function $f \in L_1[0,T]$ is the convolution of f with the Prabhakar kernel $e_{\alpha ,\beta }^{\gamma }(t; \omega )$, namely

$$\begin{aligned} {}_0{{\mathcal {J}}}_{\alpha ,\beta ,\omega }^{\gamma } f(t) = \int _0^t {\mathrm{e}}_{\alpha ,\beta }^{\gamma }(t-\tau ;\omega ) f(\tau ) {\mathrm {d}}\tau . \end{aligned}$$

(3)

Its inverse operator regularized in Caputo’s sense provides, in the case $0<\beta < 1$ and for functions $f \in AC[0,T]$, the fractional Prabhakar derivative

$$\begin{aligned} {}^{{\tiny {\text {C}}}}_0\! {\mathcal {D}}^\gamma _{\alpha ,\beta ,\omega } f(t)=\int _0^t \mathrm{e}^{-\gamma }_{\alpha ,1-\beta }(t-\tau ;\omega )f'(\tau )\mathrm{d}\tau \end{aligned}$$

(4)

(we refer to [15] for a discussion of the special case $\beta =1$).

3 Asymptotic stability of linear systems of Prabhakar-type FDEs

Due to the main interest in practical applications, we will assume throughout this paper that the parameters $\alpha$, $\beta$ and $\gamma$ fulfil the condition (1) under which the Prabhakar kernel is CM.

Consider the following linear system of Prabhakar-type fractional-order differential equations:

$$\begin{aligned} {}^{{\tiny {\text {C}}}}_0\! {\mathcal {D}}^\gamma _{\alpha ,\beta ,\omega } y(t)=Ay(t), \end{aligned}$$

(5)

coupled with the initial condition $y(0)=y_0$, and where ${}^{{\tiny {\text {C}}}}_0\! {\mathcal {D}}^\gamma _{\alpha ,\beta ,\lambda }$ is the Prabhakar differential operator (regularized in the Caputo sense) defined according to (4).

System (5) is equivalent to the following system of weakly singular Volterra integral equations of convolution type (see, for example [15, 25]):

$$\begin{aligned} y(t) = y_0 + A \int _0^t \mathrm{e}_{\alpha ,\beta }^{\gamma }(t-\tau ;\omega ) y(\tau ) {\mathrm {d}}\tau . \end{aligned}$$

(6)

For the theory of linear Volterra integral equations, including the case when the convolution kernel is completely monotonic, we refer to [3, 20, 26, 40].

From the LT of the Prabhakar kernel we observe that the characteristic equation associated with system (5) is

$$\begin{aligned} \det \left( s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma I-A\right) =0 , \end{aligned}$$

(7)

where according to [7], the principal values (first branches) of the complex power functions are taken into account.

It is easy to see that s is a root of the characteristic equation (7) if and only if there exists an eigenvalue $\lambda$ of the matrix A such that

$$\begin{aligned} s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma =\lambda . \end{aligned}$$

(8)

We obtain the following characterization of the asymptotic stability of system (5), in terms of the roots of its characteristic equation.

Proposition 1

The linear system (5) is asymptotically stable if and only if

$$\begin{aligned} \sigma (A)\subset S_{\alpha ,\beta ,\omega }^\gamma \end{aligned}$$

where $\sigma (A)$ denotes the spectrum of the matrix A and

$$\begin{aligned}S_{\alpha ,\beta ,\omega }^\gamma =\{\lambda \in {\mathbb {C}}~:~ s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma \ne \lambda ,~\forall ~\mathfrak {R}(s)\ge 0\}.\end{aligned}$$

In what follows, we will give a complete characterization of the stability region $S_{\alpha ,\beta ,\omega }^\gamma$.

4 Stability region by the root locus method

The boundary of the stability region $S_{\alpha ,\beta ,\omega }^\gamma$ will be determined using the root locus method. We first give the following preliminary results

Lemma 1

The function $\varLambda :\left[ 0,\frac{\alpha \pi }{2}\right) \rightarrow {\mathbb {C}}$ given by:

$$\begin{aligned} \varLambda (\theta )=|\omega |^{\frac{\beta }{\alpha }}~\dfrac{(\sin \theta )^{\frac{\beta }{\alpha }-\gamma }\left( \sin \frac{\alpha \pi }{2}\right) ^\gamma }{\left( \sin \left( \frac{\alpha \pi }{2}-\theta \right) \right) ^\frac{\beta }{\alpha }}~ {\mathrm {e}}^{{\mathrm {i}}\left[ \gamma \theta +(\beta -\alpha \gamma )\frac{\pi }{2}\right] } \end{aligned}$$

is a $C^\infty$ function which satisfies

$$\begin{aligned}&\lim _{\theta \rightarrow 0} \varLambda (\theta )={\left\{ \begin{array}{ll} 0 ,&{}\quad {\mathrm{if}}\; \beta -\alpha \gamma >0,\\ |\omega |^\gamma ,&{}\quad {\mathrm{if}}\; \beta -\alpha \gamma =0. \end{array}\right. }\\&\lim _{\theta \rightarrow \frac{\alpha \pi }{2}}|\varLambda (\theta )|=\infty \quad {\mathrm{and}}\quad \lim _{\theta \rightarrow \frac{\alpha \pi }{2}}\mathrm{Arg}(\varLambda (\theta ))=\frac{\beta \pi }{2}. \end{aligned}$$

Moreover:

$$\begin{aligned} 0\le \mathrm{Arg}(\varLambda (\theta ))\le \frac{\beta \pi }{2}, \quad \forall ~ \theta \in \left[ 0,\frac{\alpha \pi }{2}\right) . \end{aligned}$$

The image of the function $\varLambda$ in the complex plane, i.e. the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ defined by the parametric equation

$$\begin{aligned} \varPsi _{\alpha ,\beta ,\omega }^\gamma : \quad \lambda =\varLambda (\theta ) ,\quad \theta \in \left[ 0,\frac{\alpha \pi }{2}\right) , \end{aligned}$$

is a simple curve, included in the first quadrant of the complex plane.

Proof

The first part of the proof is trivial. Moreover, as inequalities (1) hold and $0\le \theta <\frac{\alpha \pi }{2}$, we have:

$$\begin{aligned} 0\le \frac{ (\beta -\alpha \gamma )\pi }{2}<\mathrm{Arg}(\varLambda (\theta ))<\frac{\beta \pi }{2}\le \frac{\pi }{2}. \end{aligned}$$

Therefore, the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ is included in the first quadrant of the complex plane.

Assuming by contradiction that $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ is not simple, there exist $\theta$ and $\theta '$ such that $0\le \theta<\theta '<\frac{\alpha \pi }{2}$ and $\varLambda (\theta )=\varLambda (\theta ')$. Therefore, $|\varLambda (\theta )|=|\varLambda (\theta ')|$, or equivalently:

$$\begin{aligned} \frac{\sin (\frac{\alpha \pi }{2}-\theta )}{\sin (\frac{\alpha \pi }{2}-\theta ')}=\left( \frac{\sin \theta }{\sin \theta '}\right) ^{1-\frac{\alpha \gamma }{\beta }}. \end{aligned}$$

As $0\le \theta<\theta '<\frac{\alpha \pi }{2}$, the left-hand side of this equality is larger than 1, while the right-hand side is subunitary, which is absurd. Hence, $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ is a simple curve. $\square$

Clearly, when $\beta -\alpha \gamma =0$, the parametric equation of the curve given by Lemma 1 simplifies to

$$\begin{aligned} \varPsi _{\alpha ,\alpha \gamma ,\omega }^\gamma : \quad \lambda =\left( \dfrac{|\omega |\sin \frac{\alpha \pi }{2}}{\sin \left( \frac{\alpha \pi }{2}-\theta \right) }\right) ^\gamma \! {\mathrm {e}}^{{\mathrm {i}}\gamma \theta }, \quad \theta \in \left[ 0,\frac{\alpha \pi }{2}\right) . \end{aligned}$$

In what follows, ${\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$ denotes the complex conjugate of the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ defined in Lemma 1, i.e.

$$\begin{aligned} {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma =\{\lambda \in {\mathbb {C}}~:~{\overline{\lambda }}\in \varPsi _{\alpha ,\beta ,\omega }^\gamma \}. \end{aligned}$$

We obtain the following result, characterizing the root locus of the characteristic equation (8):

Proposition 2

The characteristic equation (8) has pure imaginary roots if and only if $\lambda \in \varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$.

Proof

Assuming that Eq. (8) has a root $s=i\mu$, with $\mu \ge 0$, let us consider $\rho >0$ and $\theta \in (-\pi ,\pi ]$ such that

$$\begin{aligned} ({\mathrm {i}}\mu )^\alpha -\omega =\rho {\mathrm {e}}^{{\mathrm {i}}\theta }. \end{aligned}$$

Hence, Eq. (8) has a pure imaginary root if and only if there exist $\mu \ge 0$ and $\rho >0$ and $\theta \in (-\pi ,\pi ]$ such that

$$\begin{aligned} {\left\{ \begin{array}{ll} ({\mathrm {i}}\mu )^{\beta -\alpha \gamma }(\rho {\mathrm {e}}^{{\mathrm {i}}\theta })^\gamma =\lambda \\ ({\mathrm {i}}\mu )^\alpha -\omega =\rho {\mathrm {e}}^{{\mathrm {i}}\theta } \end{array}\right. } \end{aligned}$$

(9)

Taking the real and imaginary parts in the second equation of system (9), it follows that

$$\begin{aligned} {\left\{ \begin{array}{ll} \mu ^\alpha \cos \frac{\alpha \pi }{2}=\rho \cos \theta +\omega \\ \mu ^\alpha \sin \frac{\alpha \pi }{2}=\rho \sin \theta \end{array}\right. } \end{aligned}$$

and hence:

$$\begin{aligned} \mu ^\alpha = \dfrac{|\omega |\sin \theta }{\sin \left( \frac{\alpha \pi }{2}-\theta \right) }\quad \text {and}\quad \rho =\dfrac{|\omega |\sin \frac{\alpha \pi }{2}}{\sin \left( \frac{\alpha \pi }{2}-\theta \right) }. \end{aligned}$$

(10)

It is obvious that since $\mu \ge 0$ and $\rho >0$, the following inequalities must be satisfied:

$$\begin{aligned} \sin \theta \ge 0\quad \text {and}\quad \sin \left( \frac{\alpha \pi }{2}-\theta \right) >0. \end{aligned}$$

which is equivalent to $\theta \in \left[ 0,\frac{\alpha \pi }{2}\right)$.

Replacing $\mu$ and $\rho$ given by (10) into the first equation of system (9), we deduce that $\lambda \in \varPsi _{\alpha ,\beta ,\omega }^\gamma$. In a similar way, assuming that Eq. (8) has a root $s=-{\mathrm {i}}\mu$, with $\mu \ge 0$, it follows that $\lambda \in {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$. $\square$

Let us denote by $N(\alpha ,\beta ,\gamma ,\omega ,\lambda )$ the number of unstable roots ($\mathfrak {R}(s)\ge 0$) of the characteristic equation (8), including their multiplicities. The following lemma shows that the function $N(\alpha ,\beta ,\gamma ,\omega ,\lambda )$ is well defined. Moreover, some important properties are also established, which are needed for the proof of the main results.

Lemma 2

Let $\lambda \in {\mathbb {C}}$ and $\alpha ,\beta ,\gamma ,\omega$ satisfy inequalities (1). The following statements hold:

• (i) The characteristic function equation (8) has at most a finite number of roots satisfying $\mathfrak {R}(s)\ge 0$.

• (ii) The function $\lambda \mapsto N(\alpha ,\beta ,\gamma ,\omega ,\lambda )$ is continuous at each $\lambda \notin \varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$, and hence, $N(\alpha ,\beta ,\gamma ,\omega ,\lambda )$ is constant on each connected component of the set ${\mathbb {C}}{\setminus } (\varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma )$.

Proof

Let us denote

$$\begin{aligned} \varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda )=s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma -\lambda . \end{aligned}$$

We will first show that the set of unstable roots of Eq. (8) is bounded. Indeed, if s is a root of (8) such that $\mathfrak {R}(s)\ge 0$, as $\alpha \in (0,1]$, it follows that $\mathfrak {R}(s^\alpha )\ge 0$. Moreover, as $\omega <0$, we have:

$$\begin{aligned} |s^\alpha -\omega |=\sqrt{|s|^{2\alpha }+\omega ^2-2\omega \mathfrak {R}(s^\alpha )}\ge |s|^{\alpha }. \end{aligned}$$

Therefore:

$$\begin{aligned} |\lambda |=|s|^{\beta -\alpha \gamma }|s^\alpha -\omega |^\gamma \ge |s|^{\beta -\alpha \gamma }|s|^{\alpha \gamma }=|s|^\beta , \end{aligned}$$

and therefore, $|s|\le |\lambda |^{\frac{1}{\beta }}$.

Proof of statement (i) Let us first consider $\lambda \ne 0$. Assuming that the characteristic equation (8) has an infinite number of unstable roots, the Bolzano–Weierstrass theorem implies that there exists a convergent sequence of unstable roots $(s_n)$ with the limit $s_0\ne 0$, such that $\mathfrak {R}(s_0)\ge 0$. Since the function $\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda )$ is analytic in ${\mathbb {C}}{\setminus }{\mathbb {R}}_-$, the principle of permanence implies that it is identically zero, which is absurd. Hence, the function $N(\alpha ,\beta ,\gamma ,\omega ,\lambda )$ is finite and well defined.

If $\lambda =0$, the number of unstable roots of (5) is finite because the equation $s^\alpha -\omega =0$ has a finite number of unstable roots. This can be shown in a similar way as above, by a simple application of the principle of permanence.

Proof of statement (ii) Let $\lambda _0\in {\mathbb {C}}{\setminus }(\varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma )$ and $r>0$ such that the open neighborhood $B_r(\lambda _0)=\{\lambda \in {\mathbb {C}}~:~|\lambda -\lambda _0|<r\}$ is included in the set ${\mathbb {C}}{\setminus } (\varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma )$.

For any $\lambda \in B_r(\lambda _0)$ we have that $|\lambda |<r+|\lambda _0|$, and hence, based on the first part of the proof, any unstable root of $\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda )$ satisfies:

$$\begin{aligned} |s|< (r+|\lambda _0|)^{\frac{1}{\beta }}. \end{aligned}$$

Let us denote by (c) the simple closed curve, oriented counterclockwise, bounding the open half-disk

$$\begin{aligned} D=\{s\in {\mathbb {C}}~:~\mathfrak {R}(s)> 0,~ 0<|s|<(r+|\lambda _0|)^{\frac{1}{\beta }}\}. \end{aligned}$$

By the above construction and Proposition 2, it is clear that for any $\lambda \in B_r(\lambda _0)$, all unstable roots of the characteristic function $\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda )$ belong to D.

As $\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda _0)\ne 0$ for any $s\in (c)$, it is easy to see that

$$\begin{aligned} m_0=\min \limits _{s\in (c)}|\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda _0)|>0. \end{aligned}$$

Considering $r'=\min \{m_0,r\}$ it follows that for any $s\in (c)$ and for any $\lambda \in B_{r'}(\lambda _0)\subset B_r(\lambda _0)$, we have:

$$\begin{aligned} |&\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda )-\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda _0)|\\&=|\lambda -\lambda _0|<r'\le m_0\le |\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda _0)|. \end{aligned}$$

By Rouché’s theorem, it follows that $\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda _0)$ and $\varDelta (s;\alpha ,\beta ,\gamma ,\omega ,\lambda )$ have the same number of zeros in the half-disk D, and hence

$$\begin{aligned} N(\alpha ,\beta ,\gamma ,\omega ,\lambda )=N(\alpha ,\beta ,\gamma ,\omega ,\lambda _0)\quad ,~\forall ~ \lambda \in B_{r'}(\lambda _0). \end{aligned}$$

Hence, the function $\lambda \mapsto N(\alpha ,\beta ,\gamma ,\omega ,\lambda )$ is continuous on ${\mathbb {C}}{\setminus } (\varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma )$, and from the fact that it is integer-valued, we deduce that it is constant on each connected component of ${\mathbb {C}}{\setminus } (\varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma )$. $\square$

We now give the main result which characterizes the stability region $S_{\alpha ,\beta ,\omega }^\gamma$ of system (5).

Theorem 1

The stability region $S_{\alpha ,\beta ,\omega }^\gamma$ of system (5) is the region of the complex plane which includes ${\mathbb {C}}_-=\{\lambda \in {\mathbb {C}}~:\mathfrak {R}(\lambda )<0\}$ and is bounded by $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ and its complex conjugate ${\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$.

Proof

Lemma 1 implies that $\varPsi _{\alpha ,\beta ,\omega }^\gamma \cup {\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$ partition the complex plane into two disjoint regions, which will be denoted by $D_-$ and $D_+$. As $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ and ${\overline{\varPsi }}_{\alpha ,\beta ,\omega }^\gamma$ are included, respectively, in the first and fourth quadrants of complex plane, one of these regions includes ${\mathbb {C}}_-$ (we will further assume that ${\mathbb {C}}_-\subset D_-$). Moreover, based on Lemma 2, these regions have the property that, for every $\lambda$ within a given region, the number of unstable roots of the characteristic equation (8) is constant.

In what follows, we will show that if $\lambda \in (-\infty ,0)$, the characteristic equation (8) does not have any roots with positive real part. Indeed, let us assume by contradiction that there exists $s\in {\mathbb {C}}$, $\mathfrak {R}(s)\ge 0$ such that

$$\begin{aligned} s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma =\lambda . \end{aligned}$$

As both s and ${\overline{s}}$ are roots of the above equation, we may further assume that $\mathrm{Arg}(s)\in \left[ 0,\frac{\pi }{2}\right]$.

On the one hand, we have:

$$\begin{aligned} \mathrm{Arg}\left[ s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma \right] =\mathrm{Arg}(\lambda )=\pi . \end{aligned}$$

(11)

On the other hand, as $\beta -\alpha \gamma \in [0,1)$, we have:

$$\begin{aligned} \mathrm{Arg}(s^{\beta -\alpha \gamma })&=(\beta -\alpha \gamma )\mathrm{Arg}(s)\\&\quad +2\pi \left\lfloor \frac{\pi -(\beta -\alpha \gamma )\mathrm{Arg}(s)}{2\pi } \right\rfloor \\&=(\beta -\alpha \gamma )\mathrm{Arg}(s)\in \left[ 0,\frac{\pi }{2}\right] . \end{aligned}$$

Moreover, as $\omega <0$ and $\alpha \in (0,1]$, we deduce that $\mathrm{Arg}(s^\alpha )\in \left[ 0,\frac{\pi }{2}\right]$ and:

$$\begin{aligned} 0<\mathrm{Arg}(s^\alpha -\omega )<\mathrm{Arg}(s^\alpha )=\alpha \mathrm{Arg}(s)\le \frac{\pi }{2} \end{aligned}$$

and hence, as $0<\alpha \gamma \le 1$, it follows that:

$$\begin{aligned} 0<\gamma \mathrm{Arg}(s^\alpha -\omega )<\alpha \gamma \mathrm{Arg}(s)\le \frac{\pi }{2}. \end{aligned}$$

Therefore:

$$\begin{aligned} \mathrm{Arg}\left( (s^\alpha -\omega )^\gamma \right)&=\gamma \mathrm{Arg}(s^\alpha -\omega )\\&\quad +2\pi \left\lfloor \frac{\pi -\gamma \mathrm{Arg}(s^\alpha -\omega )}{2\pi } \right\rfloor \\&=\gamma \mathrm{Arg}(s^\alpha -\omega )\in \left[ 0,\frac{\pi }{2}\right] . \end{aligned}$$

Finally, combining the previous results and taking into account that $\beta \le 1$, we get:

$$\begin{aligned} 0\le \mathrm{Arg}&\left[ s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma \right] \\&=\mathrm{Arg}(s^{\beta -\alpha \gamma })+ \mathrm{Arg}\left( (s^\alpha -\omega )^\gamma \right) \\&=(\beta -\alpha \gamma )\mathrm{Arg}(s)+\gamma \mathrm{Arg}(s^\alpha -\omega )\\&< (\beta -\alpha \gamma )\mathrm{Arg}(s)+\alpha \gamma \mathrm{Arg}(s)\\&=\beta \mathrm{Arg}(s)\\&\le \frac{\pi }{2}, \end{aligned}$$

which is in contradiction with (11). Therefore, all the roots of the characteristic equation (8) are in the left half plane, whenever $\lambda \in (-\infty ,0)$. Hence, based on Lemma 2, we obtain that $N(\alpha ,\beta ,\gamma ,\omega ,\lambda )=0$ for any $\lambda \in D_-$, implying that $D_-\subset S_{\alpha ,\beta ,\omega }^\gamma$.

Proof of the transversality condition.

Let us denote by $s(\lambda )$ the unique root of the characteristic equation (8) such that $\mathfrak {R}\left[ s(\lambda ^\star )\right] =0$ when $\lambda ^\star \in \varPsi _{\alpha ,\beta ,\omega }^\gamma$. Let $\theta \in \left( 0,\frac{\alpha \pi }{2}\right)$ such that $\lambda ^\star =\varLambda (\theta )$, and hence, based on the proof of Proposition 2, we have:

$$\begin{aligned} s(\lambda ^\star )={\mathrm {i}}\mu = {\mathrm {i}}\left( \dfrac{|\omega |\sin \theta }{\sin \left( \frac{\alpha \pi }{2}-\theta \right) }\right) ^{\frac{1}{\alpha }}. \end{aligned}$$

Let us denote $F(s)=s^{\beta -\alpha \gamma }(s^\alpha -\omega )^\gamma$. By the implicit function theorem, as $F(s(\lambda ))=\lambda$, we have:

$$\begin{aligned} \frac{\partial s}{\partial \mathfrak {R}(\lambda )}=\frac{1}{F'(s)}\quad \text {and}\quad \frac{\partial s}{\partial \mathfrak {I}(\lambda )}=\frac{{\mathrm {i}}}{F'(s)} \end{aligned}$$

and hence:

$$\begin{aligned} \frac{\partial \mathfrak {R}(s)}{\partial \mathfrak {R}(\lambda )}=\frac{\mathfrak {R}(F'(s))}{|F'(s)|^2}\quad \text {and}\quad \frac{\partial \mathfrak {R}(s)}{\partial \mathfrak {I}(\lambda )}=\frac{\mathfrak {I}(F'(s))}{|F'(s)|^2}, \end{aligned}$$

which leads to:

$$\begin{aligned} \nabla \mathfrak {R}(s):=\frac{\partial \mathfrak {R}(s)}{\partial \mathfrak {R}(\lambda )}+ {\mathrm {i}}\frac{\partial \mathfrak {R}(s)}{\partial \mathfrak {I}(\lambda )}=\frac{F'(s)}{|F'(s)|^2}. \end{aligned}$$

A simple computation shows that:

$$\begin{aligned} F'(s)=s^{-1}\left[ \beta +\alpha \gamma \omega (s^\alpha -\omega )^{-1}\right] F(s), \end{aligned}$$

which gives:

$$\begin{aligned} \begin{aligned} F'(s)\vert _{\lambda =\lambda ^\star }&=F'({\mathrm {i}}\mu ) \\&=-{\mathrm {i}}\mu ^{-1}\left[ \beta +\alpha \gamma \omega ((i\mu )^\alpha -\omega )^{-1}\right] \lambda ^\star . \ \end{aligned} \end{aligned}$$

Finally, taking into account that $\lambda ^\star =\varLambda (\theta )$ and expressing $({\mathrm {i}}\mu )^\alpha -\omega$ in terms of $\theta$, it follows that:

$$\begin{aligned} F'(s)\vert _{\lambda =\lambda ^\star }=-{\mathrm {i}}\mu ^{-1}\left[ \beta -\alpha \gamma \dfrac{\sin \left( \frac{\alpha \pi }{2}-\theta \right) }{\sin \frac{\alpha \pi }{2}}{\mathrm {e}}^{-{\mathrm {i}}\theta }\right] \varLambda (\theta ). \end{aligned}$$

On the other hand, a straightforward computation gives:

$$\begin{aligned} \varLambda '(\theta )= & {} \dfrac{\sin \frac{\alpha \pi }{2}}{\alpha \sin \theta \sin \left( \frac{\alpha \pi }{2}-\theta \right) } \\&\quad \left[ \beta - \alpha \gamma \dfrac{\sin \left( \frac{\alpha \pi }{2}-\theta \right) }{\sin \frac{\alpha \pi }{2}}{\mathrm {e}}^{-{\mathrm {i}}\theta } \right] \varLambda (\theta ) \end{aligned}$$

and therefore:

$$\begin{aligned} F'(s)\vert _{\lambda =\lambda ^\star }=-{\mathrm {i}}\mu ^{-1}\dfrac{\alpha \sin \theta \sin \left( \frac{\alpha \pi }{2}-\theta \right) }{\sin \frac{\alpha \pi }{2}}\varLambda '(\theta ). \end{aligned}$$

Hence, exploiting the ${\mathbb {R}}^2$ vector space structure which underlies ${\mathbb {C}}$ and considering that the parametrization of the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ is fixed in the direction of increasing $\theta$, it follows that the gradient vector $\nabla \mathfrak {R}(s)(\lambda ^\star )$ is in fact a right-pointing normal vector to the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$, pointing towards the region $D_+$. We deduce that as the parameter $\lambda$ crosses the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ from the region $D_-$ into the region $D_+$, $\mathfrak {R}(s(\lambda ))$ becomes positive, which ensures that the transversality condition holds.

Moreover, this shows that if $\lambda \in D_+$, the characteristic equation (8) has at least one root with positive real part, and hence, we finally obtain that $S_{\alpha ,\beta ,\omega }^\gamma =D_-$, which completes the proof. $\square$

It is important to emphasize that Theorem 1 gives a complete characterization of the stability region $S_{\alpha ,\beta ,\omega }^\gamma$ of system (5). Some examples are shown in Figs. 1 and 2.

Remark 1

Based on Lemma 1, it is clear that when $\gamma \rightarrow 0$ (i.e. when the Prabhakar derivative in (5) reduces to the standard Caputo derivative of order $\beta$), the curve $\varPsi ^\gamma _{\alpha ,\beta ,\omega }$ approaches the half line $\mathrm{Arg}(\lambda )=\frac{\beta \pi }{2}$ of the complex plane, and hence the stability region is indeed

$$\begin{aligned} S_{\alpha ,\beta ,\omega }^0=\{\lambda \in {\mathbb {C}}~:~|\mathrm{Arg}(\lambda )|>\frac{\beta \pi }{2}\}, \end{aligned}$$

which is in accordance with Matignon’s theorem [32]. Hence, Theorem 1 is a generalization of Matignon’s theorem for the case of systems of fractional differential equations with Prabhakar derivatives.

Remark 2

The transversality condition which was verified in the proof of Theorem 1 ensures that if $\lambda \in {\mathbb {C}}{\setminus }\{0\}$ is a simple eigenvalue of the matrix A of the linear system (5), when the parameters $(\alpha ,\beta ,\gamma ,\omega )$ of the Prabhakar derivative are varied and $\lambda$ crosses from $S_{\alpha ,\beta ,\omega }^\gamma$ to its open complementary $\text {int}({\mathbb {C}}{\setminus } S_{\alpha ,\beta ,\omega }^\gamma )$, exactly one root of the characteristic equation (7) crosses the imaginary axis from the left half-plane to the right half-plane of ${\mathbb {C}}$. More generally, based on a similar argument, we can express the number of unstable roots ($\mathfrak {R}(s)\ge 0$) of the characteristic equation (7) as follows:

$$\begin{aligned} N(\alpha ,\beta ,\gamma ,\omega )= & {} \sum _{\lambda \in \sigma (A)} N(\alpha ,\beta ,\gamma ,\omega ,\lambda )\\= & {} \sum _{\lambda \in {\mathbb {C}}{\setminus } S_{\alpha ,\beta ,\omega }^\gamma } m(\lambda ), \end{aligned}$$

where $m(\lambda )$ denotes the algebraic multiplicity of the eigenvalue $\lambda$.

Fig. 1

Stability region $S_{\alpha ,\beta ,\omega }^\gamma$ for fixed values of $\alpha , \beta , \omega$ and increasing values of $\gamma$. Last figure is for the special case $\beta -\alpha \gamma =0$

Fig. 2

Stability region $S_{\alpha ,\beta ,\omega }^\gamma$ for fixed values of $\alpha , \beta , \gamma$ and decreasing values of $\omega$

5 Solution of linear and nonlinear FDEs of Prabhakar type

Analytical solutions of (5) are not available in a simple and closed form. With the purpose of verifying the theoretical findings on stability regions, we derive here asymptotic representations of the exact solution for small and large arguments, together with a numerical method for solving more general nonlinear problems with the fractional Prabhakar derivative.

For convenience we derive asymptotic expansion just for the scalar case $y(t):[0,T]\rightarrow {\mathbb {C}}$ and $A\in {\mathbb {C}}$; the generalization to the vector case is, however, straightforward.

5.1 Asymptotic expansion for small arguments

By means of the LT we can rewrite (5) in the LT domain as

$$\begin{aligned} s^{\beta -\alpha \gamma -1}\bigl (s^{\alpha }-\omega \bigr )^{\gamma } \bigl (s {\hat{y}}(s) - y_0 \bigr ) = A {\hat{y}}(s) , \end{aligned}$$

with ${\hat{y}}(s)$ the LT of y(t). Therefore the solution in the LT domain is ${\hat{y}}(s) = \mathcal{H}(s) y_0$, where

$$\begin{aligned} \mathcal{H}(s) = \frac{s^{\beta -\alpha \gamma -1}(s^{\alpha }-\omega )^{\gamma } }{s^{\beta -\alpha \gamma }(s^{\alpha }-\omega )^{\gamma } - A}. \end{aligned}$$

(12)

Observe now that

$$\begin{aligned} \mathcal{H}(s) = \displaystyle \frac{s^{-1} \Bigl (1-\displaystyle \frac{\omega }{s^{\alpha }}\Bigr )^{\gamma } }{\Bigl (1-\displaystyle \frac{\omega }{s^{\alpha }}\Bigr )^{\gamma } - \frac{A}{s^{\beta }} } = s^{-1} \left( 1- \frac{A}{s^{\beta }\Bigl (1-\frac{\omega }{s^{\alpha }}\Bigr )^{\gamma }}\right) ^{-1} \end{aligned}$$

and, for sufficiently large |s|, we can expand

$$\begin{aligned} \mathcal{H}(s) = s^{-1}\sum _{j=0}^{\infty } \frac{A^j}{s^{j \beta }\Bigl (1-\frac{\omega }{s^{\alpha }}\Bigr )^{j \gamma }} = \sum _{j=0}^{\infty } \frac{s^{\alpha \gamma j - j \beta - 1} A^j}{(s^{\alpha }-\omega )^{j\gamma }} . \end{aligned}$$

Therefore, after inverting back each LT of Prabhakar functions in the series we are able to obtain

$$\begin{aligned} y(t) = \sum _{j=0}^{\infty } A^j t^{j \beta } E_{\alpha , j \beta +1}^{j \gamma } (\omega t^{\alpha }) y_0 \end{aligned}$$

which holds as $t \rightarrow 0$.

5.2 Asymptotic expansion for large arguments

To derive an asymptotic expansion of the solution y(t) of (5) as $t\rightarrow \infty$ we consider again the function ${{\mathcal {H}}}(s)$.

Since we are now interested in study the solution ${\hat{y}} = {{\mathcal {H}}}(s)y_0$ in the LT domain as $|s|\rightarrow 0$, we have to take into account the singularity of ${{\mathcal {H}}}(s)$.

Due to the transversality condition stated by Theorem 1, ${{\mathcal {H}}}(s)$ has just one singularity, say ${\bar{s}}$, which can be eliminated in the formula for the inversion of the LT by the residue subtraction

$$\begin{aligned} h(t) = \mathop {\mathrm{Res}} \bigl ( {\mathrm {e}}^{st} \mathcal{H}(s) , {\bar{s}} \bigr ) + {\hat{h}}(t) , \end{aligned}$$

where

$$\begin{aligned} {\hat{h}}(t) = \frac{1}{2\pi i} \int _{{\mathcal {C}}} {\mathrm {e}}^{st} \mathcal{H}(s) {\mathrm {d}}s \end{aligned}$$

and ${\mathcal {C}}$ is any contour in the complex plane leaving ${\bar{s}}$ at its right and not crossing the branch-cut placed on the negative real semi-axis.

An analytical expression of ${\bar{s}}$ seems not available and therefore ${\bar{s}}$ must be evaluated numerically after solving the equation $s^{\beta -\alpha \gamma }(s^{\alpha }-\omega )^{\gamma } = A$ or, equivalently,

$$\begin{aligned} s^{\mu } - \omega s^{\mu -\alpha } - B = 0 , \quad \mu = \frac{\beta }{\gamma }, \quad B = A^{\frac{1}{\gamma }} . \end{aligned}$$

The corresponding residue can be instead evaluated by simple derivations. Indeed, since we assume $0<\alpha <1$ and real $\omega <0$, it is ${\bar{s}}^{\alpha }-\omega \not =0$; therefore, after standard derivations one obtains

$$\begin{aligned} \mathop {\mathrm{Res}} \bigl ( {\mathrm {e}}^{st} \mathcal{H}(s) , {\bar{s}} \bigr ) = C^{\gamma }_{\alpha ,\beta ,\omega }({\bar{s}}) {\mathrm {e}}^{{\bar{s}}t} , \end{aligned}$$

where $C^{\gamma }_{\alpha ,\beta ,\omega }(s)$ is constant with respect to t and

$$\begin{aligned} C^{\gamma }_{\alpha ,\beta ,\omega }(s) = \frac{ (s^{\alpha } - \omega )}{\beta s^{\alpha } - (\beta -\alpha \gamma ) \omega }. \end{aligned}$$

To evaluate ${\hat{h}}(t)$ we first observe from (12) that

$$\begin{aligned} \mathcal{H}(s) = - s^{\beta -\alpha \gamma -1} \frac{(s^{\alpha }-\omega )^{\gamma }}{A} \left( 1 - s^{\beta -\alpha \gamma } \frac{(s^{\alpha }-\omega )^{\gamma }}{A} \right) ^{-1} , \end{aligned}$$

and hence, for sufficiently small |s|, it is possible to consider the expansion

$$\begin{aligned} \begin{aligned} \mathcal{H}(s)&=- s^{\beta -\alpha \gamma -1} \frac{(s^{\alpha }-\omega )^{\gamma }}{A} \sum _{k=0}^{\infty } s^{k\beta -k\alpha \gamma }\frac{(s^{\alpha }-\omega )^{k \gamma }}{A^k}\\&= - \sum _{k=0}^{\infty } s^{(k+1)\beta -(k+1)\alpha \gamma -1 }\frac{(s^{\alpha }-\omega )^{(k+1) \gamma }}{A^{(k+1)}} \\&= - \sum _{k=1}^{\infty } s^{k\beta -k\alpha \gamma -1}\frac{(s^{\alpha }-\omega )^{k \gamma }}{A^{k}} \end{aligned} \end{aligned}$$

We can therefore transform back $\mathcal{H}(s)$ from the LT domain to the time domain to obtain ${\hat{h}}(t)$ and hence the expansion of the solution y(t) as $t \rightarrow \infty$

$$\begin{aligned} y(t) = \left[ C^{\gamma }_{\alpha ,\beta ,\omega }({\bar{s}}) {\mathrm {e}}^{{\bar{s}}t} - \sum _{k=1}^{\infty } \frac{t^{-k\beta }}{A^k} E_{\alpha ,1-k\beta }^{-k\gamma }(t^{\alpha } \omega ) \right] y_0 . \end{aligned}$$

This formula can be exploited, in connection with the asymptotic representation of the Prabhakar function introduced in Sect. 2 to provide a representation of the solution y(t) of (5).

5.3 Numerical solution

To devise an effective method for solving not only the linear system (5) but, more generally, any nonlinear system such as

$$\begin{aligned} \left\{ \begin{array}{l} {}^{{\tiny {\text {C}}}}_0\! {\mathcal {D}}^\gamma _{\alpha ,\beta ,\omega } y(t) = f(t,y(t)) \\ y(0) = y_0 \ \end{array}\right. , \end{aligned}$$

(13)

we use as starting point the standard trapezoidal rule

$$\begin{aligned} y_{n+1} - y_{n} = \frac{h}{2} \bigl (f(t_n,y_n) + f(t_{n+1},y_{n+1})\bigr ) \end{aligned}$$

for ordinary differential equations and its generalization to our problem is made in the framework devised by Lubich [27,28,29]. The choice of the trapezoidal rule as starting point for devising a numerical method for the solution of (13) is motivated by its excellent stability properties. Given the generating function of the trapezoidal rule

$$\begin{aligned} \delta (\xi ) = \frac{2(1-\xi )}{1+\xi } , \end{aligned}$$

a corresponding trapezoidal convolution quadrature rule for (13) evaluates the approximation $y_n$ of $y(t_n)$ by the formula

$$\begin{aligned} y_n = y_0 + h^{\beta }\sum _{j=0}^{s} w_{n,j} f(t_j,y_j) + h^{\beta } \sum _{j=0}^{n} c_{n-j} f(t_j,y_j) . \end{aligned}$$

Convolution weights $c_n$ are the coefficients in the asymptotic expansion of

$$\begin{aligned} \sum _{n=0}^{\infty } c_n \xi ^n = \frac{1}{h^{\beta }} G(\xi ) , \quad G(\xi ) = {{\mathcal {E}}}_{\alpha ,\beta } \Bigl ( \frac{\delta (\xi )}{h} ; \omega \Bigr ) , \end{aligned}$$

with ${{\mathcal {E}}}_{\alpha ,\beta } (s ; \omega )$ the LT of ${{\mathcal {E}}}_{\alpha ,\beta } (s ; \omega )$, and can be evaluated with high accuracy by a quadrature rule applied to the Cauchy integral

$$\begin{aligned} c_n = \frac{1}{2 \pi {\mathrm {i}}} \int _{{\mathcal {C}}} \xi ^{-n-1} G(\xi ) {\mathrm {d}}\xi \end{aligned}$$

with ${{\mathcal {C}}}$ a suitably selected closed contour encircling the origin but not any singularity of $G(\xi )$. Starting weights $w_{n,j}, j=0,1,\dots ,s$ are instead introduced to deal with the lack of smoothness at 0 of the solution and evaluated after imposing that exact solutions are obtained when $f(t,y(t))=t^{\nu }$, with $\nu$ multiple of $\beta$ less than 1. We refer again to [27,28,29] for a more detailed description.

6 Numerical experiments

6.1 Asymptotic stability

To verify the theoretical findings on the asymptotic stability we consider here, for the selection of the parameters $\alpha =0.8$, $\beta =0.9$, $\gamma =0.8$ and $\omega =-1.0$, the solution of (5) in the scalar case, for three distinct values of the coefficient $A\in {\mathbb {C}}$.

The three values of the coefficient A are selected, respectively, just inside, on the border and just outside the stability region $S_{\alpha ,\beta ,\omega }^\gamma$ determined by Theorem 1 and depicted in the left plot of Fig. 3. Since the three values of A are almost indistinguishable in the small box in the first quadrant, an enlarged view of this box is provided in the right plot of the same Fig. 3.

To observe the asymptotic behaviour of the solution of (5) we have considered both the asymptotic expansion (for large arguments) and the numerical method devised in the previous section. For large t the two approaches provide overlapping results, thus showing their reliability. We therefore report, in the following plots, only the outcomes from the numerical method which hold for small and large t.

The first plot illustrates the solution of (5) with the coefficient $A_1 = 0.866 + 1.171{\mathrm {i}}$ inside the stability region. As expected from the theory, the solution illustrated in Fig. 4 shows a stable behaviour decaying to zero.

Fig. 3

Stability region $S_{\alpha ,\beta ,\omega }^\gamma$ for $\alpha =0.8$, $\beta =0.9$, $\gamma =0.8$ and $\omega =-1$ (left plot) and zoom of the box with the three values $A_1$, $A_2$, $A_3$ near the border of the stability region (right plot)

The second plot shows the solution of the same problem when the coefficient $A_2 = 0.901 + 1.161{\mathrm {i}}$ is instead used. Since $A_2$ is on the border of the stability region, we expect that, after a transient phase, the solution presents sustained oscillations which neither decay nor amplify. This behaviour is indeed confirmed by the numerical experiment reported in Fig. 5.

Finally the third experiment concerns the coefficient $A_3 = 0.936 + 1.151{\mathrm {i}}$ outside the stability region. In accordance with theoretical expectations, the plot in Fig. 6 shows an unstable solution with oscillations of growing amplitude as t increases.

6.2 A nonlinear example

It is of interest to provide an example of an application to nonlinear systems of the theoretical results on the asymptotic stability of linear systems of differential equations with the Prabhakar derivative.

Fig. 4

Solution of the linear scalar equation (5) with $A_1=0.866 + 1.171{\mathrm {i}}$ inside the stability region

Fig. 5

Solution of the linear scalar equation (5) with $A_2=0.901 + 1.161{\mathrm {i}}$ on the border of the stability region

To this purpose we consider a dynamical system of Brusselator type

$$\begin{aligned} {\left\{ \begin{array}{ll} {}^{{\tiny {\text {C}}}}_0\! {\mathcal {D}}^\gamma _{\alpha ,\beta ,\omega } x(t)=1 - (b+1)x(t) + a x(t)^2 y(t) \\ {}^{{\tiny {\text {C}}}}_0\! {\mathcal {D}}^\gamma _{\alpha ,\beta ,\omega }y(t)=bx(t)-a x^2(t)y(t) \end{array}\right. } , \end{aligned}$$

(14)

describing an autocatalytic and oscillating chemical reaction, with the integer-order derivative replaced by the fractional Prabhakar derivative.

It is easy to compute that (1, b/a) is the equilibrium of the system and the eigenvalues of the Jacobian evaluated at the equilibrium point are

$$\begin{aligned} \lambda _{1,2} = \frac{b-a-1 \pm \sqrt{(b-a-1)^2-4a}}{2} . \end{aligned}$$

Let us take into consideration the values of the coefficients $a=10$ and $b=14$ for which the corresponding eigenvalues are $\lambda _{1,2} = 1.5 \pm 2.7839{\mathrm {i}}$.

Fig. 6

Solution of the linear scalar equation (5) with $A_3 = 0.936 + 1.151{\mathrm {i}}$ outside the stability region

Depending on the choice of parameters $\alpha$, $\beta$, $\gamma$ and $\omega$ of the Prabhakar derivative, the eigenvalues $\lambda _1$ and $\lambda _2$ can lay inside or outside the corresponding stability region $S_{\alpha ,\beta ,\omega }^\gamma$. In Fig. 7 we show the location of $\lambda _1$ and $\lambda _2$ with respect to the stability region of the Prabhakar derivative when $\alpha =0.9$, $\beta =0.95$ and $\gamma =0.8$ and $\omega =-4.0$ (left plot) or $\omega =-0.5$ (right plot).

The solution of the system (14) when $\omega =-4.0$, namely when $\lambda _{1,2}$ lie inside the stability region, is shown in Fig. 8. The two components x(t) and y(t) of the solution approach the equilibrium state (the dotted line), although in a quite slow way, in accordance with the behaviour expected from the asymptotic stability theory.

Fig. 7

Stability regions for $\alpha =0.9$, $\beta =0.95$, $\gamma =0.8$ and $\omega =-4.0$ (left plot) or $\omega =-0.5$ (right plot) and location of the eigenvalues $\lambda _{1,2}$ of the linearized Brusselator system

Fig. 8

Solution of the Brusselator system (14) for $\lambda _1$ and $\lambda _2$ in the stability region

When $\omega =-0.5$, and hence $\lambda _{1,2}$ are outside the stability region, the equilibrium point is instead unstable and, indeed, the solution of (14) oscillates without ever reaching the equilibrium point as shown in Fig. 9.

Fig. 9

Solution of the Brusselator system (14) for $\lambda _1$ and $\lambda _2$ outside the stability region

In fact, fixing the parameters $\alpha =0.9$, $\beta =0.95$ and $\gamma =0.8$ and numerically solving equation $\varLambda (\theta )=\lambda _1$, where $\varLambda (\theta )$ defines the parametric equation of the curve $\varPsi _{\alpha ,\beta ,\omega }^\gamma$ as given in Lemma 1, we obtain the critical value $\omega ^\star =-1.58444$, in correspondence of which the eigenvalues $\lambda _{1,2}$ belong to the boundary of the stability region $S_{\alpha ,\beta ,\omega ^\star }^\gamma$. We may consider that in this case, the critical value $\omega ^\star =-1.58444$ of the parameter $\omega$ corresponds to a Hopf bifurcation in the Brusselator-type system (14), resulting in the loss of asymptotic stability of the equilibrium for $\omega >\omega ^\star$ and the appearance of an attracting quasi-periodic orbit, as shown in Fig. 9. However, we emphasize that to the best of our knowledge, at present, the bifurcation theory of fractional-order differential equations with Prabhakar derivatives has not been investigated.

7 Concluding remarks

In this paper we have studied asymptotic stability properties of systems of fractional differential equations with the Prabhakar derivative. A complete characterization of the stability region was obtained, in terms of the eigenvalues of the system’s matrix, thus generalizing classical results for the stability of fractional-order systems.

We have also obtained the asymptotic representations of the solution of linear fractional Prabhakar differential equations for small and large arguments, and we have presented a numerical method for (linear and nonlinear) differential equations of fractional Prabhakar type.

Some numerical experiments using the asymptotic expansion and the numerical method have been presented in order to validate the theoretical findings. An application to the study of a nonlinear system has also been discussed.