Solution of Fractional Differential Equation Systems and Computation of Matrix Mittag–Leffler Functions

In this paper, solutions for systems of linear fractional differential equations are considered. For the commensurate order case, solutions in terms of matrix Mittag–Leffler functions were derived by the Picard iterative process. For the incommensurate order case, the system was converted to a commensurate order case by newly introducing unknown functions. Computation of matrix Mittag–Leffler functions was considered using the methods of the Jordan canonical matrix and minimal polynomial or eigenpolynomial, respectively. Finally, numerical examples were solved using the proposed methods.


Introduction
In fractional calculus, different possibilities of defining integral and derivative of non-integral order, which generalize the classical integration operator J f (t) = t 0 f (s)ds and differentiation operator D f (t) = d dt f (t), are investigated.In this respect, the Grünwald-Letnikov fractional integral and derivative, the Riemann-Liouville fractional integral and derivative and the Liouville-Caputo fractional derivative are commonly used [1][2][3][4][5][6][7].
For a function f (t) defined on (0, +∞), the Riemann-Liouville integral of order β is given as and J β t f (t) = f (t) for β = 0, where Γ(•) is the Euler's gamma function.Equation (1) can be expressed as a convolution as The Riemann-Liouville derivative of order α is defined in the form while the Liouville-Caputo fractional derivative has the definition The existence, uniqueness, and stability of solutions for fractional differential equations were studied [3,[5][6][7][15][16][17].Analytical and numerical methods are presented in [4,5,11,13,[18][19][20][21].Solutions of some linear fractional differential equations may be expressed in terms of the Mittag-Leffler functions, which are defined as [5,22] If both of the parameters are 1, the exponential function is reduced: E 1,1 (z) = e z .In [23,24], the Mittag-Leffler functions were applied to generalize fractional calculus.In [25], multi-parameter Mittag-Leffler functions were studied and applied to investigate different models of fractional calculus.In [26,27], a family of generalized fractional derivative operators and the relationships with the Mittag-Leffler type functions were considered.
Since the convergence radius of the power series on the right-hand side of Equation ( 5) is infinity, the matrix Mittag-Leffler function is well-defined for any n × n matrix A, where A 0 = I, the n × n unit matrix.In [28], the system of fractional differential equations was applied to analyze lateral motion of an elastic column.In [29], finite difference scheme for time-fractional advection-diffusion-reaction equation leads to a semilinear system of fractional differential equations.Existence, uniqueness, and stability of solutions for system of linear fractional differential equations are discussed in [30,31].Proposed analytical methods for the solutions include the Laplace transform [32,33], the Adomian decomposition method [34], the exponential integrators [29], and the multi-variable Mittag-Leffler functions [35].
The solution of linear fractional system can be expressed in terms of the matrix Mittag-Leffler function.For the computation of matrix Mittag-Leffler functions, the Krylov subspace methods were employed in [36], and convergence results are presented in [37].The method based on the Jordan canonical form representation is described in [38], and the Bagley-Torvik equation was solved with the help of MATLAB.
Evaluation of the matrix Mittag-Leffler function requires evaluation of the original scalar Mittag-Leffler function and its derivatives, and the latter was investigated in [39], where a straightforward method is the series expansion of the derivatives of the Mittag-Leffler functions, i.e., the three-parameter Mittag-Leffler functions [40].We remark that Garrappa [40] presented the optimal parabolic contour algorithm of the inversion of the Laplace transform for the computation of the two-/three-parameter Mittag-Leffler functions.
We note that the generalized exponential time differencing methods for fractional differential equations [41] and the solutions of multi-term fractional differential equations [42] also made use of the matrix Mittag-Leffler functions.
In this paper, we consider the system of linear fractional differential equations in the Liouville-Caputo sense.In the next section, we derive solutions for a system of commensurate orders using the Picard iterative process in terms of matrix Mittag-Leffler functions.In Section 3, we consider the system of incommensurate orders and convert it to a commensurate order case by introducing new unknown functions.In Section 4, we discuss two computational methods of matrix Mittag-Leffler functions, i.e., the Jordan canonical matrix and the minimal polynomial or eigenpolynomial.In Section 5, we solve numerical examples using the proposed methods.Section 6 concludes.

Solutions of System of Fractional Differential Equations of Commensurate Orders
Consider the initial value problem for system of fractional differential equations in the Liouville-Caputo sense where 0 T is piecewise continuous and integrable on the considered interval, the initial condition x(0) = (x 1 (0), x 2 (0), . . ., x n (0)) T is specified, and the Liouville-Caputo fractional derivative D α t x(t) operates componentwise, i.e., D α t x(t) = (D α t x 1 (t), D α t x 2 (t), . . ., D α t x n (t)) T .We use the Picard iterative process to derive the series solution for the initial value problem.Operating Equation (7) with the integral operator J α t and using the formula We denote by ϕ ϕ ϕ k (t) the kth approximate solution, where the zeroth approximate solution is taken as From the recurrent formula, we calculate successively . . ., Taking notice of the formula for the fractional integral we derive the series expression for the solution by taking the limit k → ∞ for ϕ ϕ ϕ k (t), In terms of matrix Mittag-Leffler functions, the solution has the form where and are the transient solution and the steady-state solution, respectively.We remark that x t (t) in Equation ( 12) is consistent with Daftardar-Gejji and Babakhani's result [30].In [29], an analogy of Equation ( 11) was given for a semilinear fractional system.What follows will be developed based on the result (11).
If α = 1 in Equation ( 11), we obtain the results for the integer-order case By operating term-by-term, we can check that the two matrix Mittag-Leffler functions in Equation ( 11) satisfy the relations If the Liouville-Caputo fractional derivative is applied in Equation ( 15), the relation becomes where I is the n × n unit matrix.Note that the Mittag-Leffler function satisfies the equation Therefore, we can express the convergent improper integral in Equation ( 13) to a normal definite integral if the derivative f (t) is continuous on the interval [0, +∞), For numerical computation, Equation ( 20) is more practical than Equation (13) and will be used in Example 3.

The System of Fractional Differential Equations of Incommensurate Orders
The system of fractional differential equations of incommensurate orders has the form where a ik are constants, not all zeroes, f i (t) are prescribed functions, D α i t are the Liouville-Caputo fractional derivative operators, the orders α i are rational numbers satisfying 0 < α i ≤ 1, and x i (t) are functions to be solved with the specified initial values x i (0), i = 1, 2, . . ., n.
In order to derive a system with the same order for each equation, i.e., commensurate orders, we first give the composition rule for the Liouville-Caputo fractional derivative: If 0 < α, β < 1 and In fact, this rule can be verified using the Laplace transform as We express the orders as , where k i and m i are irreducible positive integers for i = 1, 2, . . ., n. Denote by m the least common multiple of the denominators of α i 's, i.e., m = l.c.m.{m 1 , m 2 , . . ., m n }.
Then, for each m i , there exists an integer l i such that m i = m/l i .Denote by d the greatest common divisor of the numerators of α i 's, i.e., Then, for each k i , there exists an integer g i such that k i = g i d.Now let β = d/m.Then we have α i = g i d m/l i = p i β, where p i = l i g i for i = 1, 2, . . ., n.In order to derive a system of commensurate orders in Section 2, we decompose each derivative operators D α i t into a composition of p i fractional derivative operators D β t .For this aim, we introduce and the new variables satisfy the following equation system: . . .
This is a system of fractional differential equations of commensurate orders, subject to the initial value conditions where the accessorial zero initial values are required from the composition rule (21).

Computation of Matrix Mittag-Leffler Functions
With the foregoing discussion, the solution of system of fractional differential equations involves the matrix Mittag-Leffler functions E α,1 (At α ) and E α,α (At α ).We consider the computation of the matrix Mittag-Leffler functions E α,ρ (At α ) in the following two methods, where the matrix Mittag-Leffler functions are represented by matrix polynomials with the coefficients in terms of scalar Mittag-Leffler functions E α,ρ (λ t α ).

Method by the Jordan Canonical Matrix
Let the Jordan canonical form of n × n matrix A be J, J = diag(J 1 , J 2 , • • • , J s ), and A = PJP −1 , where the Jordan blocks have the form Thus, from the matrix theory, we have Further using the formula we set f (z) = E α,ρ (zt α ) to obtain where

Method by the Minimal Polynomial or Eigenpolynomial
The minimal polynomial of an n × n matrix A over a field F is a single-variable polynomial m A (λ), in which the coefficient of the highest-degree term is required to be 1, of the least degree such that m A (A) = O.That is, it is an annihilation polynomial of matrix A with the least degree and leading coefficient 1.
From the matrix theory, for any n × n matrix A, the minimal polynomial m A (λ) exists uniquely.The minimal polynomial and the eigenpolynomial of a matrix A have the same roots.The degree of the minimal polynomial of the matrix A is no more than n by the Cayley-Hamilton theorem.
Let the minimal polynomial of nth-order matrix A be where λ 1 , λ 2 , . . ., λ s are all of the different eigenvalues of A. By the theory of matrix analysis [43], E α,ρ (At α ) can be expressed uniquely as a matrix polynomial of degree m − 1, the well-known Lagrange-Sylvester interpolation polynomial.
For the present case, we let the polynomial Then ϕ(A; t) = E α,ρ (At α ) holds if and only if the polynomial ϕ(λ; t) and the Mittag-Leffler function E α,ρ (λt α ) are consistent in the spectrum of matrix A, that is The coefficients a i (t), i = 0, 1, . . ., m − 1, in Equation ( 26) are determined from Equation ( 27).Thus we obtain the matrix Mittag-Leffler function E α,ρ (At α ) by a finite sum as If it is difficult to determine the minimal polynomial m A (λ), we can instead use the eigenpolynomial, We determine the polynomial of degree n − 1, such that ψ(A; t) = E α,ρ (At α ).This comes true from the system of equations

Numerical Examples
In the following four examples, we express matrix Mittag-Leffler functions E α,ρ (At α ) as matrix polynomials of A, where the coefficients are given in terms of scalar Mittag-Leffler functions E α,ρ (λ t α ), which may be calculated using Podlubny's MATLAB code [44] or Garrappa's MATLAB code [45].For Figures 1-3 in Examples 2-4, we checked that the curves plotted by the two MATLAB codes are almost consistent.
We note that Garrappa's code can calculate the three-parameter Mittag-Leffler functions, hence the derivatives of the two-parameter Mittag-Leffler functions [39].Moreover, due to implementing the optimal parabolic contour algorithm [40], this routine is more efficient in some cases.

Example 1. Consider the matrix Mittag
By using the Jordan canonical form A = PJP −1 , where we obtain Alternatively, by using the eigenpolynomial of matrix A, det(λI − A) = (λ − 1)(λ − 2) 2 , we let the polynomial a 0 (t) + a 1 (t)λ + a 2 (t)λ 2 and the function E α,ρ (λt α ) be identical in the spectrum of matrix A, i.e., Solving the system leads to Thus, the matrix Mittag-Leffler function also has the form We checked that the results of Equations ( 32) and (33) are consistent.We note that the minimal polynomial and the eigenpolynomial are identical for this example.

Example 2. Consider the matrix Mittag
and the solution of the initial value problem is determined as follows: The eigenvalues of the matrix A are λ 1 = −1 − i and λ 2 = −1 + i.By using the Jordan canonical form A = PJP −1 , where we obtain Alternatively, letting the polynomial a 0 (t) + a 1 (t)λ and the function E α,ρ (λt α ) be identical in the spectrum of matrix A, we obtain Solving the system leads to Thus, the matrix Mittag-Leffler function also has the form We checked that the results of Equations ( 35) and ( 36) are consistent.The solution of the initial value problem of Equation ( 34) parameterized by the order α is If α = 1, the solution becomes the classical result In Figure 1, we show the curves of the solutions x 1 (t) and x 2 (t) versus t for different values of α.
Example 3. Consider the solution of the nonhomogeneous initial value problem where f(t) = [e −t , sin t] T and A is the same second-order matrix as in Example 2.
If α = 1, the exact solution is From Equation ( 20), we have the solution of the nonhomogeneous initial value problem of Equation ( 39): where f(0) = [1, 0] T and f (t) = [−e −t , cos t] T .Here, we use Equation ( 20) instead of Equation ( 13) for the sake of numerical computation of the convolution.From Equation ( 36), the matrix Mittag-Leffler function has the expression In order to numerically calculate the convolution in Equation ( 41), we denote t n = nh, n = 0, 1, . . ., N. Thus, we have x(t 0 ) = x(0) = [0, 0] T .For 1 ≤ n ≤ N, we apply the composite trapezoidal rule to the convolution integral in Equation (41): where In Figure 2, we plot the curves of the solutions x 1 (t) and x 2 (t) versus t on the interval 0 ≤ t ≤ 6 for different values of α, where we take h = 0.01.Example 4. Consider the system of fractional differential equations subject to the initial condition x 1 (0) = 2, x 2 (0) = 3.

Conclusions
The solutions of system of linear fractional differential equations were considered.For the case of commensurate orders, we derived the solution in terms of matrix Mittag-Leffler functions using the Picard iterative process.For the case of incommensurate orders, through introducing new unknown functions, we constructed a system of commensurate orders.We then investigated how a matrix Mittag-Leffler function can be calculated.We discussed two methods, i.e., the Jordan canonical matrix and the minimal polynomial or eigenpolynomial, to express the matrix Mittag-Leffler function E α,ρ (At α ) into a matrix polynomial, where the coefficients were given in terms of scalar Mittag-Leffler functions E α,ρ (λt α ), which may be calculated using the MATLAB code [44,45].Finally, we solved numerical examples by using the proposed methods.Solutions were plotted using MATLAB.