A Chebyshev spectral method based on operational matrix for fractional differential equations involving non-singular Mittag-Leffler kernel

In this paper, we solve a system of fractional differential equations within a fractional derivative involving the Mittag-Leffler kernel by using the spectral methods. We apply the Chebyshev polynomials as a base and obtain the necessary operational matrix of fractional integral using the Clenshaw–Curtis formula. By applying the operational matrix, we obtain a system of linear algebraic equations. The approximate solution is computed by solving this system. The regularity of the solution investigated and a convergence analysis is provided. Numerical examples are provided to show the effectiveness and efficiency of the method.


Introduction
Despite the long history of fractional calculus in the field of mathematics, a large amount of real world applications of this field has appeared mainly during the last decades. This type of calculus has become so wide that almost no branch of science and engineering cannot be found without fractional calculus and a lot of books have been written in these regards (see for example Refs. [1][2][3][4] and the references therein). Increasing the use of fractional calculations has increased the variety of questions and resulted in various basic definitions for fractional integral and derivative. We recall that the Riemann-Liouville definition entails physically unacceptable initial conditions [1]; conversely for the Liouville-Caputo fractional derivative, the initial conditions are expressed in terms of integer-order derivatives having direct physical significance [1,5]. A few years ago Caputo and Fabrizio [6] have opened the following subject of debate within the mathematical community: is it possible to describe all nonlocal phenomena within the same basic kernels, namely the power kernel involved within the definition of Riemann-Liouville derivative and some other few basic fractional derivatives. If we analyze, step by step, the way Caputo has introduced his classical fractional derivative [5], we will realize that during the last step he generalized the classical integral to the fractional Riemann-Liouville integral (see for more details [5]). After that, about almost 50 years later, he kindly asked the mathematical community how the Gamma function appears into the description of real phenomena, and why only some existing fractional operators are required by experiments [6]. Immediately, Nieto and Losada found, by using the Laplace transform, the associated integral of the so-called Caputo-Fabrizio fractional derivative [7].
The dynamics of many applied physical or biological problem can be modeled by a system of fractional differential equations (FDEs) (for example see [24] for a relaxation system). A system of Mittag-Leffler non-singular FDEs can be described by where A is a constant matrix of dimension ν × ν, ν ∈ N is the dimension of the system, f : R → R ν is a known vector-valued function, ABC 0 D α t y(t) is a fractional derivative involving Mittag-Leffler functions (also known as AB type [10]) and y : R → R ν is the unknown function. Recently, it was observed that the system (1) is more successful for modeling of suspension concentration distribution in turbulent flows than other models [25].
The conditions for the existence and uniqueness of the solution to exponential nonsingular system can be found in [26]. The consistency condition Ay 0 + f(0) = 0, is one of them. It seems that this condition is also important for system (1) with Mittag-Leffler non-singular kernels and as mentioned in [27], we should consider the initial condition carefully. This imposes some restriction on system (1). However, due to the important dynamics of the solutions of system (1), it is significant to solve the system (1) analytically or numerically [28]. Because of the novelty and the newness of this topic there are a few articles on this subject. We found only, a linear piecewise polynomial base method for solving this system numerically [29]. However, solving this system with Chebyshev base polynomials has not been studied yet.
The spectral methods using Chebyshev polynomials are well known for differential and partial differential equations [30][31][32][33]. For smooth problems in simple geometries, they offer exponential rates of convergence or spectral accuracy. An important advantage of these methods over finite-difference methods is that computing the coefficient of the approximation, completely determines the solution at any point of the desired interval. Therefore, numerical solution of the system (1) using operational matrix spectral methods based on Chebyshev polynomials is very important.
The discrete orthogonality properties of the Chebyshev polynomials are the advantages over other orthogonal polynomials like Legendre polynomials. Also, the zeros of the Chebyshev polynomials are known analytically. These properties lead to the Clenshaw-Curtis formula which makes integration easy. We use this formula to obtain the operational matrix of the fractional integration.
The aim of this paper is to obtain an efficient numerical method to solve the system (1) using the operational matrix based on Chebyshev polynomials. For this purpose, we obtain the operational matrix approximation for fractional integral operator. We transform the system (1) to a system of weak singular integral equation and then using an operational matrix we obtain a system of linear algebraic equations. Solving the obtained algebraic system, we get the numerical approximation. We investigate the existence and convergence of the numerical solution. To this end, we also study the regularity of the exact solutions.
The structure of this paper is as follows. In Sect. 2, we review new definition of the fractional calculus and related results and the Chebyshev polynomials. In Sect. 3, we review the approximations of multi-variable functions in terms of the shifted Chebyshev polynomials and we obtain the operational matrix approximation for fractional integral operator. In Sect. 4, we propose a spectral method based on the operational matrix for solving system of Mittag-Leffler, non-singular FDEs. In Sect. 5, we obtain the regularity of the solutions. In Sect. 6, we study a convergence analysis for proposed method without discretization. In Sect. 7, we obtain the convergence results for discretized version. Finally, in Sect. 8, we provide some numerical examples to show the efficiency of the introduced method, and we present a comparison between the solution of the AB type FDEs and the Liouville-Caputo type FDEs.

Definitions and preliminaries
In this section, we first recall some basic definitions related and results to the Mittag-Leffler function [34]. Then we recall some basic definitions and results related to the new non-singular fractional derivative and integral formulas [10].

The Mittag-Leffler function
The Mittag-Leffler function is the cornerstone of fractional calculus. Several books and excellent papers [34][35][36][37] describe the importance of these types of operators. The concept of Mittag-Leffler calculus was introduced in [10] and the integral associated to the nonsingular fractional operator with Mittag-Leffler kernel was found by using the Laplace transform [38].
Throughout the paper, the symbol E α shows the one parameter Mittag-Leffler function [39] defined by The two-parameter Mittag-Leffler function is defined as Here, the notation denotes the gamma function. An interesting book containing the history, applications and the effect of the gamma functions on the progress in mathematics and the progress in describing the real phenomena can be found in [40].

The non-singular fractional derivative and integral involving Mittag-Leffler kernel
We use a Sobolev space defined by to define the fractional derivative as follows.
The associated fractional integral is also defined by [10] AB The fractional integral of (tt 0 ) β (β > -1) for α > 0 is The Newton-Leibniz formula for this fractional derivative and integral is obtained in [38,41].
From Theorem 2.1, the fractional derivative of a monomial t β (β > 0) is

Chebyshev polynomials
Here, we review some basic definitions and results related to the Chebyshev polynomials [30,42].
The Chebyshev polynomials are orthogonal with respect to the weight function w(x) = 1 √ 1-x 2 and the corresponding inner product is The well-known recursive formula with T 0 (x) = 1 and T 1 (x) = x is important for numerical computing of these polynomials, whereas we may use to compute Chebyshev polynomials in analysis. Since the range of interest of the problem (1) is [0, T], we can define the shifted Chebyshev polynomials T * n (x) by [30], Sect. 1.3) we could compute the shifted Chebyshev polynomials by The discrete orthogonality of Chebyshev polynomials leads to the Clenshaw-Curtis formula: where x k for k = 1, . . . , N + 1 are zeros of T N+1 (x). Therefore, we have Also, the norm of T * n (x), will be of importance later.

Function approximation
A function f (t) defined over the interval [0, T], may be expanded as where p N : , is an orthogonal projection, π N is the space of polynomials with degree not exceeding N , C and are the matrices of size (N + 1) × 1 and The following error estimate for the Dini-Lipschitz continuous function f provides the convergence of approximation by Chebyshev polynomials.
where ω is the modulus of continuity. Then gp n g ∞ → 0 as n → ∞.
where P α,M = (p n,r ) is the operational matrix of dimension N × N and its elements can be computed using for r = 1, . . . , N , and for n = 1, . . . , N and r = 0, . . . , N .
Proof Taking the fractional integral on both sides of (11), we get x n-k+α for n > 0, and For n = 0, it can easily be checked that Now, applying (15) to the f (x) = x n-k+α , we obtain By substituting the coefficients of T * r (x) from (19) into (17) and (18) we obtain the desired result.
, the maximum error of Clenshaw-Curtis formula is less than 4 fp N f ∞ [43]. Hence, the Clenshaw-Curtis formula for x n-k+α in the proof of 3.2 is convergent and we conclude that where P α = lim M→∞ P α,M .

Constructing the method
Taking fractional integration from both sides of system (1) and using (5), the system (1) can be written in the following form: Let E = I -1-α B(α) A. Then, using the following lemma, one can guarantee the invertibility of E.
A < α, and I denotes the identity matrix. Then the matrix Proof The proof is a direct result of the geometric series theorem. Now, multiplying by E -1 both sides of (21), we obtain the second kind of weakly singular integral equations of the form In order to obtain a numerical method, we suppose to be the approximate solution, f(t) = F T (t) and Y 0 = [y 0 , 0, . . . , 0] T . Substituting them into (22) and using the operational matrix of Theorem 3.2, we obtain where Solving the linear system (24), we obtain Y T , and finally we obtain the approximate solution using (23).
To solve (24), we can use the vectorization operators to obtain system of linear algebraic equations in the standard form. We denote by vec the vectorization of a matrix I ν×ν is the identity matrix and the notation ⊗ is for the Kronecker product. By the vec notation, the system (24) can be transformed to the standard form which it can be solved by mathematic software like MATLAB.

Regularity of the solution
For the simplicity of notation and analysis, we may write the system (22) as follows: where The system (25)  for i = 0, . . . , m and a real number c ∈ R such that g = t λ g 0 (t) + c and It is straightforward to show that the space C m,λ (0, T] equipped with the norm is a Banach space. We note that, for 0 < λ 1 < λ 2 ≤ 1, we have Remark 5.1 We note that, for f ∈ C m,λ (0, T] there exists a positive constant c > 0 such that and for f ∈ C 0,λ (0, T] the norms f ∞ and f 0,λ are equivalent. Proof By integral substitution (τ = tz), we obtain where Here, we should be concerned that the systems we investigated are of dimension ν ≥ 1 and we use the norm Proof By using Lemma 5.2,f(t) ∈ (C m,α (0, T]) ν . Consider a Picard iteration corresponding to the system (25) with y 0 =f . The first iteration can be written in the form , and the second iteration can be written in the form Using the variable transformation s = τ + (tτ )z, we obtain where Q 2 (t, τ ; α) : Proceeding by this procedure and by an argument similar to [44], Chap. 6 (note that we have used 1α instead of α), one can show that and Q n (t, s; α) := αE -1 A B(α) (α) 1 0 (1z) α-1 z (n-1)(α)-1 Q n-1 s + (ts)z, s; α dz.
Therefore, (I -T ) -1 is a bounded operator in (C m,α (0, T]) ν . Other parts of the proof are straightforward.

Convergence analysis
The orthogonal projection p N : (C[0, T]) ν → (π N ) ν (m, N ∈ N ) can be defined by where p N is defined by (13). The introduced method can be written in the form where N f(t)) and y N ∈ (π N ) ν . It is well known that the operator I α is compact on C[0, T] (see [45]) and hence is Proof Since t α satisfies the Dini-Lipschitz condition, f satisfies the Dini-Lipschitz condition and, by Theorem 3.1, f -p n f ∞ → 0, as N → ∞. The latter can be obtained by equivalency of the norms. Setting X = (C 0,α [0, T]) ν , T = T and using Theorem 6.2, we have as N → ∞. Here, the notation L((C 0,α [0, T]) ν ) shows the space of linear operators on (C 0,α [0, T]) ν , and the operator norm is induced norm. The operator T is compact and due to the fact that compact linear operators are bounded (see [45]), the operator T is also bounded. In order to obtain the convergence of the proposed method, we can use the following lemma to show that the operator (I -p N T ) -1 exists and is bounded for all sufficiently large N . Then, for all sufficiently large N , say N > N 0 , the operator (Ip N T ) -1 exists as a bounded operator from X to X. Moreover, it is uniformly bounded: Taking into account we can write Now, taking the norm from both sides of Eq. (31), we obtain Finally, we note that where c is a constant number, and we can state the following theorem.
Theorem 6.5 Assume that f ∈ (C 0,α (0, T]) ν , for α ∈ (0, 1). Then, for sufficiently large N , the approximated solution of system (1), say y N , obtained by (23) and (24) exists and converges to the exact solution y. Furthermore, we have where c is a constant number.
Also, the result of Theorem 6.5 is true with the norm · ∞ . This holds by considering the equivalency of the norms · 0,α,ν and · ∞ .

Existence and convergence results for discretized version
Often, we cannot compute the infinite series of P α , and instead we use P α,M M ∈ N. We call the obtained method the discretized version, because we discretized the corresponding integral by the Clenshaw-Curtis formula. For y N (t) = Y T (t), we can define T M y N by and we have The numerical approximation by the discretized version can now be obtained by solving where y N,M (t) = Y T M (t). According to Remark 3.3, we have as M → ∞ on C 0,λ with λ > 0. Hence, by Theorem 6.2, as M → ∞ on Banach space C 0,λ . In order to show that I -p N T M is invertible, we note that Regarding the arguments of previous section, (I -p N T ) is invertible for sufficiently large N . Also, (I + (I -p N T ) -1 (p N T y N -p N T M )) is invertible by geometric series theorem for sufficiently large M, since (I - Thus, we conclude that, for all sufficiently large M and N , say M > M 0 and N > N 0 , the operator I -p N T M is invertible. This fact guarantees the existence of the numerical solution and we can write In order to provide the convergence analysis, we note that Therefore, it remains to provide the convergence for y N,M -y N ∞ . Since

Numerical examples
To show the effectiveness and efficiency of the method some examples are presented for illustration. In the rest of the paper, we assume that B(α) = 1, and we obtain the maximum error by         Tables 2 and 3 show the maximum error for α = 0.5 and α = 0.9 and various N . As these tables show, with increasing N the maximum error decreases rapidly and the proposed method converges to the exact solution. By Theorem 6.5, y -y N ∞ is less than c( y -p N y ∞ + f -p N f ∞ ), Therefore, we set c = e 3   We can obtain the exact solution using the Laplace transform as follows: Table 4 shows the maximum error for α = 0.9 and various N . Figures 3 and 4 show the numerical solution for various α on [0, 15]. We observe that, as α approaches 0, the solution   Figure 5 shows the error behavior of this example, for α = 0.5. Theoretically, we expect yy N ∞ is less than c( yp N y ∞ + fp N f ∞ ) up to a constant multiplier c. We observe a similar behavior for both components of the error with c = e 0.5 . This is in complete agrement with the error analysis.

Application
Since dynamical systems with ordinary or fractional differential equations are abundant in natural phenomenon, we expect that, like other type of FDEs, the AB type FDEs will also be successful in modeling of natural phenomenon. This was confirmed with the comparison in the previous section. Beside the fact that the AB type derivative has been recently introduced, many applications in this field can be found in the literature we talked about in the introduction section. One of many examples for modeling of natural phenomenon using AB type FDEs is the amount of drug lidocaine in the bloodstream and body tissue [47]. The human disease of ventricular arrythmia or irregular heartbeat is treated clinically using the drug lidocaine. Let X(t) be the amount of lidocaine in the bloodstream and Y (t) be the amount of lidocaine in body tissue. Then the dynamics of the drug therapy This system is equivalent to the system of Example 8.3, with f = [0, 0] T , and was solved by ordinary and Liouville-Caputo fractional derivatives in the literature. The solution of this system with AB type fractional derivative is illustrated in Figs. 6 and 7 for various value of α.

A comparison of the proposed method with other methods
Due to the novelty of the subject, there are few numerical methods available in the literature for solving the system (1). We found only a type of predictor-corrector method intro-  [29] for solving such systems numerically. We use the numerical values reported in that paper to compare our proposed method with that one. Consider Example 4.1 of [29]. It corresponds to A = 0, and f (t) = 0, with the notations of this paper. As we see in Table 5, the result of our method, even with small N = 1, 2, 3, is much better than the result of [29].

Conclusion
We used the Clenshaw-Curtis quadrature to obtain an operational matrix of a fractional integral based on Chebyshev polynomials. Then, by taking the fractional integral from both sides of the system of FDEs involving the Mittag-Leffler kernel, we obtained a system of second kind weakly singular equations involving the fractional integral of the nonhomogeneous part. This system was transformed to a system of linear algebraic equations, and using vectorization operator we obtained the system of standard linear algebraic equations. Solving it we obtained an approximate solution of the system of AB type FDEs. Table 5 A comparison between the proposed method of this paper and the method proposed in [29], (Method 1), with N = 1, 2, 3 for α = 0.5 Numerical examples showed that the proposed method effectively and accurately solves the system of FDEs. The successfulness of the new definition of the fractional derivative and integral, involving the Mittag-Leffler kernel, makes it important to improve numerical methods for nonlinear system of FDEs and to implement numerical methods in other bases. Therefore, future studies with this new type of calculus are to find more experimental data and to compare the related results to the Liouville-Caputo and Riemann-Liouville derivatives. Also, the stability of the related fractional differential equations needs more investigation.