Numerical Solutions for Multi-Term Fractional Order Differential Equations with Fractional Taylor Operational Matrix of Fractional Integration

In this article, we propose a numerical method based on the fractional Taylor vector for solving multi-term fractional differential equations. The main idea of this method is to reduce the given problems to a set of algebraic equations by utilizing the fractional Taylor operational matrix of fractional integration. This system of equations can be solved efficiently. Some numerical examples are given to demonstrate the accuracy and applicability. The results show that the presented method is efficient and applicable.


Introduction
Fractional calculus is an emerging field of mathematics, which is a generalisation of differentiation and integration to non-integer orders. The history of fractional calculus is almost as long as the history of classical calculus, beginning with some speculations of Leibniz (1695, 1697) and Euler (1730). However, fractional calculus and fractional differential equations (FDEs) are increasingly becoming popular in recent years. The progressively developing history of this old and yet novel topic can be found in [1][2][3][4][5]. In fact, fractional calculus provides the mathematical modeling of some important phenomena like social and natural in a more powerful way than the classical calculus. During the last few decades, many applications were reported in many branches of science and engineering such as chaotic systems [6,7], fluid mechanics [8], viscoelasticity [9], optimal control problems [10,11], chemical kinetics [12,13], electrochemistry [14], biology [15], physics [16], bioengineering [17], finance [18], social sciences [19], economics [20,21], optics [22], chemical reactions [23], rheology [24], and so on. Due to the importance of FDEs, the solutions of them are attracting widespread interest. On the other hand, analytical solutions are not always possible for solving them. Therefore, numerical techniques becomes more important for solving such equations.
Multi-term fractional differential equations are one of the most important type of FDEs, which is a system of mixed fractional and ordinary differential equations and involving more than one fractional differential operators. Nowadays, they are widely appearing for modelling of many important processes, especially for multirate systems. Their numerical solution is then a strong subject that deserves high attention. In this paper, motivated by the results reported in [40,41] for solving a smaller class of problems where the highest order of derivative is an integer and involving at most one noninteger order derivative, we go further and establish a method for numerical solutions for higher order and arbitrary multi-term fractional differential equations which have a general form where D α representing the Caputo fractional derivative of order α > 0 and we assume that 0 Multi-term fractional order differential equations also have useful properties and they can describe complex multi-rate physical processes in a various way and can be applied in many fields, see e.g., [2,4,26,42]. Basset [43] and Bagley-Torvik [44] equations can be given as important examples for smaller class of multi-term fractional differential equations. Existence, uniqueness and stability of solution for multi-term fractional differential equations are discussed in [45][46][47][48][49]. Because of difficulty of finding the exact solutions for such equations, many new numerical techniques have been developed to investigate the numerical solutions such as Adams method [50], Haar wavelet method [51], differential transform method [52], Adams-Bashforth-Moulton method [53], collocation method based on shifted Chebyshev polynomials of the first kind [54], Boubaker polynomials method [55], matrix Mittag-Leffler functions [56], differential transform method [57] and so on.
Our main purpose is to present an effective, reliable method to approximate initial value problem for the Equation (1). In order to reach this aim, we rewrite and focus the general type of Caputo multi-term fractional differential equation given in Equation (1) in the following linear form subject to the y (p) (0) = Y p , p = 0, 1, ..., n − 1 where n − 1 < α < n u i (i = 0, 1, ..., k) are known coefficients and (3) Here, we also state that the highest order α need not to be an integer. This equation is important in applications due to the fact it can treat the problems with fractional force, therefore it is suitable for being treated within fractional operators of Caputo type.
In this work, a numerical approach based on fractional Taylor vector is proposed to solve the initial value problem of general type of multi-term fractional differential equations which is given in Equations (2) and (3). The core idea of this method is to employ the operational matrix of fractional integration based on fractional Taylor vector to given problem and reduce it to a set of algebraic equations which can be efficiently solved.
The structure of the manuscript is organized as follows. In Section 2, we briefly introduce some preliminary ideas of fractional calculus and necessary definitions. In Section 3, an operational matrix of fractional integration based on fractional taylor vector is derived. In Section 4, we present the numerical algorithm to solve the given equation and a pseudo-code for matlab is also provided in Algorithm 1. In Section 5, the presented method is applied to six examples to demonstrate the efficiency. A final conclusion is presented in the last section.

Preliminary Knowledge
In this section, we recall some fundamental definitions and preliminary facts of fractional calculus.

The Fractional Integral and Derivative
Definition 1. The Riemann-Liouville fractional integral to order α of an integrable function y(t) is defined to be When applied to a power function, it yields the following result: The operator has a semigroup property, namely and it is linear, namely for any two functions y 1 ,y 2 and constants A 1 ,A 2 .

Definition 2.
The fractional derivative of y(t) of the order α in the Caputo sense is given as

1.
The Riemann-Liouville fractional integral and Caputo fractional derivative do not usually commute with each other. The following Newton-Leibniz identity gives an important relation between them: 2.
The Caputo fractional derivative also has the following substitution identity. If we write y 1 (q) = y(qR) and q = t/R, then where j − 1 < α ≤ j, j ∈ N

Fractional Taylor Basis Vector
We shall make use of the fractional Taylor vector, for m ∈ N and δ > 0 in the work of this paper.

Approximation of Function
Suppose that T mδ (t) ⊂ H, where H is the space of all square integrable functions on the interval [0, 1]. For any y ∈ H, since S = span 1, t δ , t 2δ , ..., t mδ is a finite dimensional vector space in H, then, y has a unique best approximation y * ∈ S, so that ∀ y ∈ S, y − y * ≤ y − y Therefore, the function y is approximated by fractional Taylor vector as following where T mδ (t) denote the fractional Taylor vector and are the unique coefficients.

The Numerical Algorithm
In this section, to solve the given multi-term fractional differential equation in Equations (2) and (3), we employ the fractional Taylor method. The algorithm of method is given below.
Firstly, by using the transformation q = t/R, we replace the variable t ∈ [0, R] with q ∈ [0, 1]. Now, by using Equation (8) in Equation (2), we get where f 1 (q) = f (qR) and y 1 (q) = y(qR). Similar to Equation (10) we approximate the y 1 (q) as (15) such that T mδ (q) = [1, q δ , q 2δ , ..., q mδ ] T is the fractional Taylor vector and the unique coefficients C T is given in Equation (11). Next, applying the Riemann-Liouville fractional integral on both side of (14), we get where Hence, by substituting initial conditions (3), we get j! . Now, by using the Equation (12), we approximate the fractional order integrals in Equation (17) and we have Finally, by taking the collocation points q j = j/m (j = 0, 1, ..., m) in Equation (18), we get m + 1 linear algebraic equations. This linear system can be solved for the unknown vector C T . Consequently, y 1 (q) can be approximated by Equation (15).

MATLAB Implementation of Method
The pseudocode given in Algorithm 1 below allows us to use proposed method in MATLAB for obtain a numerical solution of given problem [58].

Algorithm 1: Fractional Taylor Method
[A, b] = f ractionalTaylor(al pha, beta, Uk, f unc, t0, R, y0, m, delta) % Input % al pha is the highest order of fractional derivative of given equation % beta is the order of fractional derivatives other than alpha. beta must be a vector with decending ordered values % Uk is the vector of coefficients % f unc is defining the right hand side of given problem % t0 and R denotes the left and right endpoints % y0 is the initial conditions % m denotes the number of steps % delta is a real number greater than zero. We usually take delta = 1 or delta =fractional part of al pha where command f ractionalTaylor.m is defined by the Equation (18), gives us the linear system AC = B which is (m + 1) % algebraic equations with unknown coefficients C T % Next step is to use matlab function linsolve(A, b) to solve obtained algebraic equation for unknown coefficient vector C T with dimension (m + 1).
% C is an (m + 1) x 1 matrix which is the solution of linear system AC = B % Next step is substituting obtained coefficients to approxSoln() as input, where the command approxSoln() defined by Equation (15), we get the approximate solution of given problem [s, y] = approxSoln(C) % Input % C is the vector of coefficients obtained in previous step.
% Output % s is the nodes on [t0, R] in which the approximate solution calculated % y is the numerical solution evaluated in the points of s.

Illustrative Examples
To illustrate the applicability and effectiveness of the presented method, we give six examples in this section. In each example, we apply the fractional Taylor operational matrix method which is presented in previous section and the approximate results compared with analytical solutions.
Obtained results indicate that the proposed technique is very effective for multi-term fractional differential equations. In order to solve the numerical computations, MATLAB version R2015a has been used.

Example 1
Consider the following form of multi-order fractional differential equation [59] where the exact solution is We apply the given procedure which is implemented in previous section for solving the Equation (19) step by step.
Firstly, change variable t ∈ [0, R] to q ∈ [0, 1] by using q = t/R. Now, we use the Equation (8) and get where 0 ≤ q ≤ 1. Next, using Equation (7) we get Now, using Equation (21) and substituting initial conditions y(0) = V 0 , y (0) = V 1 into equation From Equation (12), we have Now, taking R = 1 in Equation (23) and putting the given values for V 0 , V 1 , u i , β i where i = 0, 1, 2, 3 into this equation, we get Finally, taking the collocation points q j = j/m (j = 0, 1, ..., m) generates a linear algebraic system of dimension m + 1 with unknown vector C T . In order to solve this system by using presented method and comparing the results, we choose δ = 1 and different values of m.
To show the efficiency, we compared the numerical results with the method given in [59]. Table 1, compares the obtained results for absolute error with m = 4, 6, 7. We observe from Table 1 that, the absolute errors for presented method are smaller and the numerical solution is more accurate for the same size of m. In Figures 1-3, we present the graphical representation of comparison between exact solution and the numerical solutions obtained by proposed method and the method of [59] for the problem (19) with m = 4, 6, 7 respectively. From these results, we can conclude that m = 4 and m = 6 give larger absolute error, while m = 7 gives smaller absolute error (10 −16 ) and more precise numerical solution. These comparisons also shows that the results obtained by proposed method is closer to the exact solution than the results of [59].   In Figure 4, we show the graphical representation of absolute errors obtained by using proposed method and the method of [59] with m, n = 6. From Figure 4, we can conclude that the absolute error obtained by our method is remaining smaller and stable while the absolute error of other method is increasing in the interval [0, 1].
In Figures 5 and 6, we give the graphical representation of absolute errors obtained by using proposed method with m = 4, 7 respectively.

Example 2
In this example, we consider the Equation (19) with α = 2, V 0 = V 1 = 0, the coefficients u 0 = u 2 = −1, u 1 = 0, u 3 = 2 and β 0 = 0, β 2 = 2 3 , β 3 = 5 3 and the function is The exact solution of this equation is y(t) = t 3 [59]. Applying the same procedure to given problem as presented in Example 1, we get the following equation C T T mδ = 2q 1/3 C T (G 1/3 * T mδ (q)) − q 4/3 C T (G 4/3 * T mδ (q)) − q 2 C T (G 2 * T mδ (q)) + I 2 f 1 (q) (25) As we stated in previous example, collocating this equation at the nodes q j = j/m (j = 0, 1, ..., m) generates a system of algebraic equations. In this example, to solve this sysem for C T , we choose δ = 1, 1.5 and different values of m. Table 2 shows the results for obtained absolute errors by using presented method with m = 2, 3. From these results, we can see that, there is satisfactory agreement between the exact solution and numerical solutions. The absolute error is achieved about 10 −15 . We also note that, the proposed method gives better results for m = 2 by taking δ = 1.5. In Figure 7a, we show the graphical representation of obtained numerical solution and the exact solution of the given problem. Figure 7b presents the obtained absolute error by using proposed method with m = 3.

Example 3
Consider the multi-term fractional order initial value problem [54] D (2.2) y(t) + 1.3D (1.5) y(t) + 2.6y(t) = sin(2t), (26) with initial conditions y(0) = y (0) = y (0) = 0, where the equation have the series solution given by [52] y s (t) = 28561 3600000 In order to solve this problem, we choose δ = 1, and m = 10. We give the comparison of series solution and the numerical solution obtained by presented method in Table 3. Table 4 compares the obtained absolute errors by using presented method with the results of [54]. From this compared results, it can be seen that the approximate solution is very close to series solution for a small number of m for the given method.
From the compared results of Table 4, we can conclude that the proposed method has better approach to series solution with a smaller m.  The graphical representation of comparison between series solution and numerical solutions obtained by presented method and the method of [54] in the interval [0, 1] is illustrated in Figure 8. In Figure 9, we show present graphical representation of absolute errors obtained by using proposed method and the method of [54] with m = 10. In Figure 10, we show the graphical representation for series solution and the numerical results of presented method for the interval [0, 10]. The results plotted in Figure 10 are in a very good and satisfactory agreement with the series solution given in [52] and the results of [60].
In order to apply the presented method to Equation (28) and compare the results with methods of [54,61,62], we solve this problem with α = 0.3, 0.5, 0.7, 1.25, 1.5, 1.85, and different values for δ and m. The obtained results are presented as below.
In Table 5, we list the results of obtained absolute errors for α = 0.3, 0.5, 0.7 by use of presented method. Also, the results for α = 1.25, 1.5, 1.85 are given in Table 6.   In Figure 11a,b, we present the graphical representation of obtained results for numerical and exact solution of the given problem and absolute error for α = 1.5 in the interval [0, 1].   Table 7 lists the obtained absolute errors for the given problem (28) at t = 1, 5, 10, 50 and α = 1.5 by use of presented method and some other methods in literature [54,61,62]. From this compared results, we can say that the numerical solution obtained by use of proposed method is in better agreement with the exact solution and obtained absolute error is smaller. In Figure 13, the behaviour of absolute error for α = 1.5 with m = 4 and δ = 1/2, 1 at t ∈ [0, 50] is presented. From this graph, it can be seen that we get better results by taking δ = 1/2 for this example and the numerical solution is very close to exact solution for a small number of m.

Example 5
In this example, we consider the following form of linear multi-term fractional differential equation with variable coefficients [65] with, y(0) = 2, y (0) = 0 where 0 < β 2 < 1, 1 < β 1 < 2 and whose the exact solution is given by y(t) = 2 − t 2 2 . We give the numerical solution for the given problem by proposed method for a = 1, b(t) = √ t, c(t) = t   Table 8, we give the results for maximum errors obtained by use of proposed method and comparison with the results of [65,66]. From this compared results, we can see that the numerical solution obtained by use of proposed method is closer to the exact solution.  Figure 14 presents the graphical representation for behaviour of numerical and exact solutions with m = 6. From this representation, we can see that the numerical solution is in a very good agreement with exact solution.

Example 6
For the last example, let us consider the below fractional differential equation [63] y(0) = 1 which arises, for example, in the study of generalized Basset force occuring when a spherical object sinks in a (relatively dense) incompressible viscous fluid; see [43,67]. By use of Laplace transformation of Caputo derivatives, we get the analytical solution as following where the Mittag-Leffler function E λ,µ (t) with parameters λ, µ > 0 is given as This Mittag-Leffler function and its variations are very significant in fractional calculus and fractional differential equations [68].
In order to solve given problem by use of proposed method and compare the results, we take t ∈ (0, 5] and use different values of δ and m. Table 9 lists the exact and obtained numerical solutions by use of presented method and method of [63] for the given problem for m = 5, 10, 15, 20. Comparison of this results shows that, even for small values of m, the numerical solution obtained by use of presented method is in a better agreement with exact solution. In Figures 15a, 16a and 17a, we present the graphical representation of comparison between exact solution and the numerical solutions obtained by using proposed method and the method of [63] with taking m = 5, 10, 20 respectively. Also in Figures 15b, 16b and 17b we show the behaviour of absolute errors obtained by proposed method and the method of [63] in the interval [0, 1] with m = 5, 10, 20.   From these graphical results represented in Figures 15-17, we can conclude that the absolute error obtained by our method is remaining smaller when compared the absolute error of method given in Reference [63].

Conclusions
In this work, an operational matrix based on the fractional Taylor vector is used to numerically solve the multi-term fractional differential equations by reducing them to a set of linear algebraic equations, which simplifies the problem. From comparison of the obtained results with exact solutions and also with results of other methods in the literature, we conclude that the proposed method provides the solution with high accuracy. The findings also show that, even for the small number of steps, we can get satisfactory results by using presented method. All computational results are obtained by using MATLAB.