On the analysis of mixed-index time fractional differential equation systems

In this paper we study the class of mixed-index time fractional differential equations in which different components of the problem have different time fractional derivatives on the left hand side. We prove a theorem on the solution of the linear system of equations, which collapses to the well-known Mittag-Leffler solution in the case the indices are the same, and also generalises the solution of the so-called linear sequential class of time fractional problems. We also investigate the asymptotic stability properties of this class of problems using Laplace transforms and show how Laplace transforms can be used to write solutions as linear combinations of generalised Mittag-Leffler functions in some cases. Finally we illustrate our results with some numerical simulations.

Time fractional models typically have solutions with heavy tails as described by the Mittag-Leffler matrix function [13] that naturally occurs when solving time fractional linear systems. However such models are usually only described by a single fractional exponent, α, associated with the fractional derivative. The fractional exponent can allow the coupling of different processes that may be occurring in different spatial domains by using different fractional exponents for the different regimes. One natural application here would be the coupling of models describing anomalous diffusion of proteins on the plasma membrane of the cell with the behaviour of other proteins in the cytosol of the cell. Tian et al [14] addressed this problem by coupling a stochastic model (based on the Stochastic Simulation Algorithm [15]) for the plasma membrane with systems of ordinary differential equations describing reaction cascades within the cell. It may also be necessary to couple more than two models and so in this paper we introduce a formulation that focuses on coupling an arbitrary number of domains in which dynamical processes are occurring described by different anomalous diffusive processes. This leads us to consider the r index time fractional differential equation problem in Caputo form A ij y j + F i (y), y i (0) = z i , y i ∈ R mi , i = 1, · · · , r, (1) or in vector form D α t y = A y + F (y). Here the A ij are m i × m j matrices, while A is the associated block matrix of dimension r j=1 m j and α = (α 1 , · · · , α r ) has all components α i ∈ (0, 1]. We believe that a modelling approach based on this formulation has not been fully developed before. We note that scalar linear sequential fractional problems have been considered whose solution can be described by multi-indexed Mittag-Leffler functions [16], and there are a number of articles on the numerical solution of multi-term fractional differential equations [17,18,19], and while mixed index problems can, in some cases, be written in the form of linear sequential problems, namely R i=1 D βi t y = f (y), we claim that it is inappropriate to do so in many cases.
Therefore in this paper we develop a new theorem that gives the analytical solution of equations such as (1) that reduces to the Mittag-Leffler expansion in the case that all the indices are the same (section 3) and generalises the class of linear sequential problems. We then analyse the asymptotic stability properties of these mixed index problems using Laplace transform techniques (section 4), relating our results with known results that have been developed in control theory. In section 5 we show that, in the case that the α i are all rational, the solutions to the linear problem can be written as a linear combination of generalised Mittag-Leffler functions, again using ideas from control theory and transfer functions. In section 6 we present some numerical simulations illustrating the results in this paper and give some discussion on how these ideas can be used to solve semi-linear problems either by extending the methodology of exponential integrators to Mittag-Leffler functions, or by writing the solution as sums of certain Mittag-Leffler expansions.

Background
We consider the linear system given in (1) with r = 2. It will be convenient to let We will call such a system a time fractional index-2 system. Here the Caputo time fractional derivative with starting point at t = 0 is defined (see Podlubny [20], for example), as Furthermore, given a fixed mesh of size h then a first order approximation of the Caputo derivative [21] is given by If β = α then the solution to (1) is given by the Mittag-Leffler expansion where Γ(x) is the Gamma function. If the problem is completely decoupled, say A 2 = 0, then from (3) the solution to (1) and (2) satisfies In order to solve (4), this requires us to solve problems of the form Before making further headway, we need some additional background material.
Definition 1. Generalisations of the Mittag-Leffler functions are given by The following result will be important in section 5. .
The Caputo derivatives satisfy the following relationships.
Lemma 4. The solution of the scalar, linear, non-homogeneous problem is Proof: Using the integral form from Lemma 3, (7) can be rewritten as We now apply a Picard-style iteration of the form where y 0 (t) = y 0 , ∀t.
It can be shown that this iteration will converge to (8) -see [16].
Proof: Use Definition 1 and integrate the left hand side term by term.

Remark 1.
The function multiplying f (s) in the integrand of (8), namely can be viewed as a Green function. For example, when α = 1, G 1 (t − s) = e λ(t−s) .
The generalisation of the class of problems given by (7) to the systems case takes the form In the case that F (t) = 0 the solution of the linear homogeneous system is From this we can prove Theorem 6. The solution of (9) is given by Proof: We can use the idea of a Green function. But first of all it is trivial to see from Definition 1 that Now the solution of (9) can be written as where G α (t − s) is a matrix Green function, or alternatively Finally it is clear from (12) that is the Green function and the result is proved by using Lemma 5.
Note that the proofs of Lemma 4 and Theorem 6 can be found in Podlubny [20].

The solution of mixed index linear systems
The main focus of this paper is to consider generalisations of (9) where there are differing values of α on the left hand side. In its general form, we will let y = (y 1 , · · · , y r ) ∈ R m where y i ∈ R mi and m = r i=1 m i . We will also assume F (t) = (F 1 (t) , · · · , F r (t) ) and that A can be written in block form A = (A ij ) r i,j=1 , A ij ∈ R mi×mj . We will also let α = (α 1 , · · · , α r ) and consider a class of linear, non-homogeneous multi-indexed systems of FDEs of the form that we interpret as the system The index of the system is said to be r.
In the case that F = 0, then by letting we can rewrite (13) as M y = 0, where M is the block matrix, whose determinant must be zero, with Thus, in the case all m i = 1, so that the individual components are scalar and so m = r, (15) implies Det(M ) y r = 0.
For example, when r = 2 this becomes While, for r = 3 this gives after some simplification Clearly there is a general formula for arbitrary r in terms of the cofactors of A. In particular, it can be fitted into the framework of linear sequential FDEs [16,20,21,22,23]. These take the form However, this characterisation is not particularly simple, useful, or computationally expedient. Furthermore when the m i are not 1, so that the individual components are not scalar, then there is no simple representation such as (16) and new approaches are needed. Before we consider this new approach we note the converse, namely that (16) can always be written in the form of (13) for a suitable matrix A with a special structure. In particular we can write (16) in the form of (13) with p = r − 1 as an r dimensional, r index problem with α = (β 0 , β 1 , · · · , β p ), and For completeness we note in the case that d = 0 and f (t) = 0, an explicit solution to this problem was given in Podlubny [20]. This can be found by considering the transfer function (see section 4) given by By finding the poles of this function and converting back to the untransformed domain, Podlubny gives the solution as .
We now return to the index-2 problem (1) and (2). We first claim that the solution takes the matrix form where the α n,j , β n,j are appropriate matrices, of size m 1 × m and m 2 × m, respectively, that are to be determined. We now use the fact that Using (17) and (18) the left hand side of (1) is If we define α n,n+1 = 0, β n0 = 0, n = 1, 2, · · · (20) then the right hand side of (1) is Equating (19) and (21) we find along with (20) that for n = 0, 1, 2, · · · α 00 β 00 In order to get a succinct representation of the solution based on (17) and (22), it will be convenient to write so p n (t) ∈ R m(n+1)×m , and let p 0 (t) = I m .
Thus we can state the following theorem.
(i) In the case α = β, (24) reduces, as expected, to (ii) It will be convenient to define the matrix so that the solution (24) can be expressed as (iii) If the fractional index-2 system has initial condition y(t 0 ) = z then the solution is We note that in solving (9) an equivalent solution to (11) is where G α is the Green function satisfying This leads us to give a general result on the solution of the mixed index problem with a time-dependent forcing function D α,β t y(t) = Ay(t) + F (t), but first we need the following definition. Definition 2. Let y(t) = (y 1 (t), y 2 (t)) , then define I α,β t y(s)ds = I α t y 1 (s)ds, I β t y 2 (s)ds .
Theorem 8. The solution to the fractional index-2 problem is given by Proof: The result follows from the above discussion and noting that We now turn to analysing the asymptotic stability of linear fractional index-2 systems.

Asymptotic stability of multi-index systems
The first contribution to the asymptotic stability analysis of time fractional linear systems was by Matignon [24]. Given the linear system D α t y(t) = Ay(t) in Caputo form, then taking the Laplace transform and using the definition of the Caputo derivative gives Here X(s) is the Laplace transform of y(t). If we write w = s α , then the matrix s α I − A will be nonsingular if w is not an eigenvalue of A. In the w-domain this will happen if Re(σ(A)) ≤ 0, where σ(A) denotes the spectrum of A. In the s-domain this will happen if |Re(σ(A))| ≥ απ 2 . That is, the eigenvalues of A lie in the complex plane minus the sector subtended by angle απ symmetric about the positive real axis -see Figure 1.  In fact Laplace transforms are a very powerful technique for studying the asymptotic stability of mixed index fractional systems. Deng et al. [25] studied the stability of linear time fractional systems with delays using Laplace transforms. Given the delay system then the Laplace transforms results in ∆(s) X = b ∆(s) = Diag(s α1 , · · · , s αm ) − L (33) L ij = a ij e −sτij , i, j = 1, · · · , m.
Hence, Deng et al. proved: Theorem 9. If all the zeros of the characteristic polynomial of ∆(s) have negative real part then the zero solution of (32) is asymptotically stable.
Deng et al. also proved a very nice result in the case that all the indices α 1 , · · · , α m are rational.
Theorem 10. Consider (32) with no delays and all the α i ∈ (0, 1) and are rational. In particular let and let M be the lowest common multiple of all the denominators and set γ = 1 M . Then the problem will be asymptotically stable if all the roots, λ, of Remarks 3.
(i) If α i = α, i = 1, · · · , m then Theorem 10 reduces to the result of Matignon. The proof of Theorem 10 comes immediately from (33) where p(λ) is the characteristic polynomial of ∆(s).
(ii) A nice survey on the stability (both linear and nonlinear) of fractional differential equations is given in Li and Zhang [26], while Saberi Najafi et al. [27] has extended some of these stability results to distributed order fractional differential equations with respect to an order density function. Zhang et al [28] consider the stability of nonlinear fractional differential equations.
(iii) Radwan et al. [29] note that the stability analysis of mixed index problems reduces to the study of the roots of the characteristic equation In the case that the α i are arbitrary real numbers, the study of the roots of (34) is difficult. By letting s = e z , we can cast this in the framework of quasi (or exponential) polynomials (Rivero et al. [30]). The zeros of exponential polynomials have been studied by Ritt [31].
The general form of an exponential polynomial with constant coefficients is a j e αj z .
An analogue of the fact that a polynomial of degree k can have up to k roots is expressed by a Theorem due to Tamarkin, Pólya and Schwengler (see [31]).
Theorem 11. Let P be the smallest convex polygon containing the values α 1 , · · · , α k and let the sides of P be s 1 , · · · , s k . Then there exist k half strips with half rays parallel to the outer normal to b i that contain all the zeros of f . If |b i | is the length of b i , then the number of zeros in the i th half strip with modulus less than or equal to r is asymptotically r |bi| 2π .
If the α i are rational and with M the lowest common multiple of the denominators, this reduces to the polynomial This leads us to think about stability from a control theory point of view. Thus given the system where α n > · · · > α 0 , β M > · · · β 0 then the solution of (35) can be written in terms of the transfer function where s is the Laplace variable (see Rivero et al. [30], Petras [32]).
This leads to the following result, given in Cěrmák and Kisela [33].
Theorem 12. Equation (38) is asymptotically stable with α > β > 0 real and α β rational if We now follow this idea but for arbitrarily sized systems in our mixed index format, and this leads to slight modifications to (38). We first make a slight simplification and take m 1 = m 2 and we also assume that A 2 is nonsingular, then problem (1) leads to and substituting into the equation for y 1 gives This leads us to consider the roots of the characteristic function In the scalar case this gives an extension to (38) where the characteristic equation is Now reverting to Laplace transforms of (1) and (2) then This can be written in systems form as where or alternatively as This can now be considered as a generalised eigenvalue problem. From (41) we require D 1 − A to be nonsingular. That is Let us write v = (v 1 , v 2 ) and assume α ≥ β and that s β I − B 2 is nonsingular, so that from the previous analysis this means Hence

Thus (43) and
Det(( define the asymptotic stability boundary -see also (39). In order to make this more specific, let m 1 = m 2 = 1 and Note that σ(A) = {d ± √ ab}. Then (44) becomes Furthermore, let b = −a = θ, so that the eigenvalues of A are d ± iθ and (46) becomes If we now assume that which defines the asymptotic stability boundary (the imaginary axis) when α = β = 1, then (47) becomes Now since θ and d are real, the imaginary part of the right hand side of (48) must be zero, so that Equations (49) and (50) will define the asymptotic stability boundary with θ as a function of d. Rewriting (49) as d = r α+β sin α+β 2 π r α sin απ 2 + r β sin βπ 2 .
(51) and substituting (50) leads after simplification to Using the relationships gives Since cos α − β 2 π − cos α + β 2 π = 2 sin απ 2 sin βπ 2 and letting x = r α−β , then we can write (52) as Furthermore, we can write (51) as It is easily seen that as a function of x the minimum of (53) is when x = 1. (56) Remarks 4. We have the following results forθ in three particular cases: In the case α + β = 1 we see from (53) that Letting α = 1 2 + with > 0 small, then x = r 2 . This means that x 2 +1 2x , as a function of r, is very shallow apart from when r is near the origin or very large. Hence the asymptotic stability boundary will be almost constant over long periods of d when α and β are close together.
we can write (53) and (54) as where Due to the nonlinearities in (58) it is hard to determine an explicit simple relation between φ and d except if α = 2β. In this case we make use of the following Lemma.
Lemma 14. If x 2 − ax + b = 0 and x 2 − cx + d = 0 then there is a solution Proof: By subtraction of the two equations and substitution.
In the case of (57) and (58) then (59) becomes Note that Some manipulation from (60) leads to Now since α = 2β, this reduces to By takingθ = arctan(φ) this gives an explicit relationship betweenθ and d for the case α = 2β.

Remarks 5.
• β = 1 2 , α = 1 gives (62) It is clear from (61) that when d = 0 and d = ∞, then θ = π 2 and then the angle will make an excursion from π 2 down to a minimum value and back to π 2 as d increases. For example, in the case of β = 1 2 , α = 1 we can show from (62) that the minimum value of the angle is when Some of these aspects are shown in the Simulations and Results section.

Further analysis
Returning to (41) and taking m 1 = m 2 = 1 and then the Laplace transform in (42) is where Now if α and β are rational (α ≤ β) α = m n , β = p q , m ≤ n, p ≤ q, positive integers and with z = s 1 nq , then Hence (64) gives From Descartes rule of sign, then (65) will have at most 4 real zeros if mq + np is even, and at most 5 real zeros if mq + np is odd. Now factorise Det(z) = Π N j=1 (z − λ j ), N = mq + np, where there are at most 4 real zeros if N is even and at most 5 real zeros if N is odd. Then using (66) and (67) we can write where the A (i) j can be found by writing Using Lemma 2 withα = 1 nq ,β = 1 − α +α leads to the following result.
Theorem 15. The solution of the mixed index 2 problem with α = m n , β = p q , m ≤ n, p ≤ q all positive integers is, with N = mq + np, given by j ) , where the λ j are the zeros of (65) and the A j are the coefficients in the partial fraction expansion.

Remarks 6.
(i) In the case that α = β then (68) should collapse to the solution and this is not immediately clear. However, in this case, mq = np and so which is a quadratic function in z np while the equivalent p 1 and p 2 numerator functions are linear in z np . Thus in (68) N is replaced by 2, 1 nq is replaced by α, and 1 − α + 1 nq becomes 1. Thus (68) reduces to that then becomes (69). (ii) In the case that α is rational and β = Kα, K a positive integer, then If we factorise Det(s) = Π K+1 j=1 (s α − λ j ) and find A (1) j , A (2) j , j = 1, · · · , K + 1 by then we have the following Corollary.
Corollary 16. The solution of the mixed index 2 problem with α rational and β = Kα, K a positive integer, is given by where the vectors A j and "eigenvalues" λ j satisfy (71). As a particular example, take K = 2, α = p q , then the λ j and A j in Corollary 16 satisfy

In other words
Clearly in the case described by Corollary 16, writing the solution as a linear combination of generalised Mittag-Leffler functions makes the evaluation of the solution much more computationally efficient.  In Figures 2 and 3 we plot the asymptotic stability boundary of the two dimensional, index-two problem given by (1) where

Simulations and results
for the two cases considered in section 4, namely β = 1, α = 1 2 ( Figure 2) and β = 2 3 , α = 1 3 ( Figure  3). Since the eigenvalues of A are d ± iθ, we plot on the vertical axis the angleθ in radians, where   θ = 1 π arctan( θ λ ), as a function of d. In Figure 2 we see thatθ ∈ ( 1 4 , 1 2 ) corresponding to an angle lying between 45 • and 90 • , as expected from the theory. We also plot the angle, in green, corresponding to the midpoint between these two extremes, i.e. 3 8 π. We see that for the most part the asymptotic stability angle lies above this midpoint except for the values of d, as shown in the right hand figure.
In the case of Figure 3, we give a similar plot as Figure 2. We also plot in green the midpoint between the two lines subtended by angles 1 3 π and 1 6 π, namely 1 4 π. As with Figure 2 there is a small range of d for which the asymptotic stability angle drops beneath 1 4 π. Furthermore, it is clear from Remarks 4(ii) that as α and β approach one another, the asymptotic stability boundary will be almost constant over increasingly longer periods of d and will only asymptotically approach the angle π 2 for very small and very large values of d -see Remark part (ii).
In Figure 4 we confirm the asymptotic stability analysis showing sustained and decaying oscillations with α = 1 2 , β = 1 (top panel) and α = 1 3 , β = 2 3 (bottom panel). In all four cases, d = 1 while for the top panel we take θ = 2(1 + √ 33 − 1 + 0.3. In Figure 5 we present phase plots of y 1 versus y 2 for the two decaying oscillations cases. The figures confirm our theoretical results on the asymptotic stabiity boundary and also show the effects that the fractional indices have on the period of the solutions. As α approaches β we expect the oscillatory behaviour to disappear.
Finally, in Figure 6 we consider the problem in which case the eigenvalues of A are d ± θ. We take d = −1, θ = 1 2 and present the solutions for four pairs of indices, namely (α, β) = (0.85, 0.95), (0.5, 0.95), (0.2, 0.05), (0.15, 0.95). The simulations show that the components of the solution y 1 and y 2 seem to pick up "energy" from one another due to the coupling and that as the distance between α and β grows there is a greater separation between the two components. Finally, as α gets smaller, the solutions appear to "flat-line" more quickly.

Conclusions
In this paper we have studied mixed index fractional differential equations with coupling between the different components. We find an analytical expression for the solution of the linear system that generalises the Mittag-Leffler expansion of a matrix and the solution of linear sequential fractional differential equations. We can use this result to derive new numerical methods that generalise the concept of exponential methods used in the approximation of the Mittag-Leffler matrix function, see [34,35,36], for example, and exponential integrators [37], [38]. The second element would deal with developing numerical techniques for the integration component that incorporates the integral of a function times a Green function. We also use Laplace transform techniques to find the asymptotic stability domain in terms of the eigenvalues Figure 6: For A given by (73) with d = −1, θ = 1 2 so that the eigenvalues are − 3 2 , − 1 2 , showing the effect of variation of α with fixed β on the system dynamics.