On Solutions of the Initial Value Problem for the Three-Term Fractional Di ﬀ erential Equation with Caputo Derivatives

: In this paper, two forms of an exact solution and an analytical–numerical solution of the three-term fractional di ﬀ erential equation with the Caputo derivatives are presented. The Prabhakar function and an asymptotic expansion are utilized to present the double series solution. Using properties of the Pochhammer symbol, a solution is obtained in the form of an inﬁnite series of generalized hypergeometric functions. As an alternative for the series solutions of the fractional commensurate equation, a solution received by an analytical–numerical method based on the Laplace transform technique is proposed. This solution is obtained in the form of a ﬁnite sum of the Mittag-Le ﬄ er type functions. Numerical examples and a discussion are presented.


Introduction
The calculus of non-integer order and fractional differential equations has been applied in the mathematical modeling of physical phenomena in the different fields of science and engineering [1][2][3][4][5][6][7]. An example of the use of fractional derivatives in mathematical modeling is the fractional differential equation describing heat conduction in a medium. This equation was derived by using a non-local dependence between the heat flux vector and the temperature gradient [4]. As a result, the time-fractional partial differential equation of the heat conduction was obtained. The separation of variables in this equation leads to a time-fractional three-term ordinary differential equation [8,9].
The exact solutions of initialvalue problems of the fractional order are often expressed by the special functions [10,11]. In particular, there are the Mittag-Leffler function and the Prabhakar function, as well as the generalized Wright and generalized hypergeometric functions. The properties of the Mittag-Leffler and the Prabhakar functions are the subject of [12][13][14]. The generalized Wright function and hypergeometric functions are discussed in [15][16][17][18]. In the context of numerical calculations, it should be pointed out that the Prabhakar function and the generalized Wright function are not implemented in the software packages Mathematica and Maple, while the Mittag-Leffler function and the generalized hypergeometric function are predefined in these packages. Therefore, if a solution of the problem is expressed by the Prabhakar or generalized Wright function, then the usage of a relationship between this function and the generalized hypergeometric function allows one to convert the solution to the form with predefined hypergeometric function. In this paper, a relationship is used for presenting the solution of the three-term fractional equation in terms of the generalized hypergeometric functions.
Linear initial value problems are often solved by using the Laplace transform method. A stage of the usage of this method is determining the inverse Laplace transform. If an analytical inverse transform in a closed form is not possible to find, then the methods for numerical inversion can be used [19][20][21]. However, the inversion of the Laplace transform is an ill-conditioned problem [21,22]. This means that arbitrarily small changes in the transform can generate arbitrarily large changes in the value of the inverse Laplace transform. If the Laplace transform is given in the form of a rational function then, in order to find the inverse transform, an analytical-numerical method can be applied. In this case, the singularity points of the Laplace transform as roots of a polynomial can be numerically determined and used to create the inverse transform. A similar approach can be applied for the Laplace transform, which is obtained as a solution in the transform domain of a commensurate fractional differential equation [23]. The inverse Laplace transform is obtained in terms of the Mittag-Leffler and the Prabhakar functions.
The ill conditioning of the numerical inversion of the Laplace transforms is a disadvantage of the Laplace transform technique used to solve differential equations. In turn, a disadvantage of the series solutions is the slow convergence of many functional series or the loss of the convergence of these series for some arguments. Due to these disadvantages, a control of the numerical results of the solutions of the differential equations is needed. For the purpose of such control, various mathematical methods derived from different theoretical bases may be applied simultaneously.
In this paper, the two forms of an exact solution of a three-term fractional differential equation with the Caputo derivatives is presented. The solution of the differential equation, for which the difference of fractional derivative orders is a rational number, was derived in the form of a series of generalized hypergeometric functions. A solution of the commensurate fractional equation in the form of a finite sum of the Mittag-Leffler and the Prabhakar functions using an analytical-numerical method is presented. A numerical comparison and discussion of these methods and their effectiveness are carried out.

Exact Solution of the Three-Term Fractional Equation
We consider the differential equation in the form a D α t y(t) + A a D β t y(t) + B y(t) = f (t) (1) where A and B are real constants and a D α t and a D β t are the fractional Caputo derivatives of order α and β, respectively. The Caputo derivative is defined by where Γ is the gamma function. We continue for the rest of the paper with a = 0 and 0 < β < α ≤ 2. The fractional differential Equation (1) of order α ∈ (0, 1] is complemented by the initial condition. If 1 < α ≤ 2, then in addition to the condition (3a), the initial condition for the derivative is also required [3].
∂ y ∂ t t=0 = y 1 (3b) The general solution of the three-term fractional equation using the generalized Wright function is presented in [3]. The solution of the initial value problem (1), (3a,b) can be written in the form and whereas p ψ q is the generalized Wright function, which is defined as [13]: It can be noted that G 0 α, β is Green's function for the differential Equation (1) with zero initial conditions.
The function M j (α, β, γ, t) can be expressed by the Prabhakar function E γ α,β , which is defined by The Prabhakar function (8) . Using Equations (6)-(9), we obtain that For certain values of orders α and β, the function M j (α, β, γ, t) can be represented in terms of the generalized hypergeometric function p F q , which is defined by On the basis of Equations (8)-(11), we obtain for α − β = 1: For α − β = p/q, where p and q are positive integer numbers, using Equations (8) and (10), we can write: Another form of the expression (13) can be obtained by utilizing the following properties of the Pochhammer symbol: where a ∈ R, i, n, q ∈ N. Using these symmetric (14a) and non-symmetric (14b) properties in Equation (13), one obtains where Taking into account Equation (11), we find that U i,j (t) is the generalized hypergeometric function As a result, the solution of the initial value problem (1), (3) was obtained in the form (4) with function G γ α, β (t) given by Equation (5) and functions M j (α, β, γ, t) defined by (15), whereas U i,j (t) are the generalized hypergeometric functions (17).
As mentioned earlier, the generalized hypergeometric functions are implemented in the software packages Mathematica and Maple. Therefore, the presentation of a solution of an initial value problem by generalized hypergeometric functions is a benefit from the viewpoint of numerical calculations with the help of these software packages.

A Solution of the Three-Term Fractional Equation Using the Laplace Transform Technique
We consider the fractional order initial value problem (1), (3) with the derivative orders satisfying the inequalities: 0 < β < α ≤ 2. In order to solve the problem, we use the Laplace transform method.
where s is a complex parameter. The Laplace transform f (s) exists, if f (t) is an exponentially bounded function [24], i.e., if there exist real constants M > 0, t 0 > 0 and c such that For exponentially bounded functions, we utilize the property of the Laplace transform of the Caputo derivative [11]: Applying the Laplace transformation to Equation (1), using initial conditions Equations (3a) and (3b) and Equation (19) we obtain the Laplace transform of the solution y(s) in the form (20) From Equation (19), it follows that the term s α−2 q on the right-hand side of Equation (20) occurs only if 1 < α ≤ 2. Similarly, the term A s β−2 q occurs only if 1 < β ≤ 2. The inverse Laplace transform, y(t) = L −1 [y(s)], can be written as where g γ α, β (t) = L −1 g γ α, β (s) and The inverse transform of Equation (22) can be determined by an expansion of this function in a series. Using the expansion proposed by Podlubny [10], the function G γ α, β (s) can be expressed in the following series form Utilizing the Laplace transform pair, the inverse Laplace transform of the function (23) is obtained as it can be noted that the same form has the function G γ α, β (t), which was introduced in the previous section (Equations (5) and (10)) by using the relationship between the generalized Wright function and the Prabhakar function. Hence, we have G γ α, β (t) = g γ α, β (t) for all t for which the series (25) is convergent.
The solution of the fractional differential Equation (1) is presented in terms of the infinite series (5) of the generalized Wright functions (6) or the Prabhakar functions (10) or the generalized hypergeometric functions (15). We now derive a solution of this initial value problem in the form of a finite sum of the Mittag-Leffler and the Prabhakar functions. We consider the fractional differential Equation (1) with commensurate derivative orders α and β. In this case, there exists such a real number ν and two integer numbers m and n that α = m ν and β = n ν. Firstly, we define the auxiliary function where k = γ ν , m = α ν , n = β ν and · is the floor function. The denominator of this function ϕ m, n (z) = z m + A z n + B is the characteristic polynomial of the fractional commensurate Equation (1). Using Equations (22) and (26), we can write Applying the procedure of partial fraction decomposition, we rewrite function H k m,n (z) for 0 ≤ k < m in the form of a sum Using the formula (24), we find the inverse Laplace transform for 0 ≤ k < m as If k < 0, then the function (26) can be rewritten in the form In this case, z 0 = 0 is the root of the denominator of the function (31) and µ 0 = −k is the multiplicity of this root. Taking this into account in Equations (28) and (29), we write l = 0 instead of l = 1 in these equations. As a result, we obtain the function G γ α, β (t) as Finally, the solution of the initial problem (1), (3) is given by equation (21), whereas the function G

Numerical Analysis
The function (4) is a solution of the initial value problem (1), (3). This solution is expressed by the function G γ α, β (t), which is defined by Equation (5). Taking into account Equations (8) and (10) in Equation (5), we write the function G The series (33) can be used for the numerical calculation of the function G γ α, β (t) for a small argument t. For large arguments, the asymptotic expansion of the Prabhakar function can be utilized, which is presented in [13]. Using this asymptotic expansion in Equation (10), we obtain The results of the numerical calculations of the function G 0 α, β (t) by using three representations of this function are presented below. The first representation is the series (5) with coefficients (15) expressed by the generalized hypergeometric functions (17). The second method for the calculation of this function relies on using the double series (33) and (34), which are applicable for small and large arguments, respectively. The third representation of the function G 0 α, β (t) was obtained by using the analytical-numerical method in the form of the finite sum (30). The function G 0 α, β (t) is Green's function for the fractional differential Equation (1). For α = 2, this equation is known as the Bagley-Torvik equation [25].
We consider the differential Equation (1) assuming: α = 2, β = 3/2, A = 2, B = 1/2 and f (t) = sin t. Parameters occurring in Equation (26) are: m = 4, n = 3, k = 0 and ν = 1/2. In order to calculate the values of the function G 0 2, 3/2 (t) by using Equation (30), the roots z l of the characteristic polynomial and the corresponding coefficients C r l must be determined. The characteristic polynomial ϕ 4, 3 (z) = z 4 + 2 z 3 + 1/2 has two single real roots and a pair of conjugate complex roots. The roots z l and the corresponding coefficients C 1 l are Using Equation (30), we get the function G 0 2, 3/2 (t) in the form Numerical values of the function G 0 2;3/2 (t) computed using the infinite series (5) with coefficients (15) expressed by the generalized hypergeometric functions (17), the double series (33), (34) and the finite sum (35) are listed in Table 1. To receive proper accuracy of the values of this function represented by the series of the generalized hypergeometric functions (15), computations with a high precision should be made. As follows from Table 1, the values of the functions (33) and (34) in a certain interval are approximately equal, however, the first series can be only applied for small arguments and the second for large arguments. For the calculation of the function with the Formula (35), firstly the roots of the polynomial equation and the coefficients of the partial fraction decomposition must be determined. If the computations are made with the help of the software packages Mathematica or Maple, then the roots and the coefficients can be appointed with arbitrary precision. The Mittag-Leffler function occurring in the Formula (35) is predefined in both software packages and the arbitrary precision evaluation for the complex arguments is provided. The presented results show that the correctness of the received numerical values should be confirmed by using different forms of the analytical solution. Table 1. Values of the function G 0 2;3/2 (t) computed using the infinite series and the analytical-numerical method: (a) the series of generalized hypergeometric functions (Equations (5), (15) and (17)  Assuming f (t) = sin t in Equation (1), a particular solution of this equation satisfying zero initial conditions can be written in the convolution form For the calculation of this integral, the Laplace transform property can also be used. The Laplace transform y(s) of this function can be written as Introducing z = s 1/2 , we obtain an auxiliary function in which the denominator is ψ(z) = z 4 + 2z 3 + 2 z 4 + 1 . It follows that the four roots determined earlier should be complemented by the roots of the equation z 4 + 1 = 0. Then calculating the coefficients C 1 l , we can write the solution as a sum In the second example, we consider differential Equation (1), assuming that A = where w is a constant, w > 1. Moreover, we assume that the function f (t) is given in the form In this case, the exact solution of Equation (1) with zero initial conditions is the function Considering the differential Equation (1)  , where y(t) is obtained by using (42) and y exact (t) is the value of the function (40), are presented in Table 2 for w = 1.5; 2.0; 2.5. The small relative errors confirm the usefulness of the analytical-numerical method for also solving the commensurate fractional differential equation for the high degree of the characteristic polynomial. The fractional differential Equation (1) for α = 1 and 0 < β < 1 is known as the Basset equation [11]. We consider this equation, assuming A = 5, B = 2 and f (t) in the following form The exact solution of this equation for arbitrary β ∈ (0, 1), satisfying zero initial conditions is the function y exact (t) = sin λ t. The solution derived on the basis of Equations (4) and (30) for a rational order β = n/m has the form where m and n are positive integers, relatively prime numbers and m > n. The results of numerical calculations obtained by using the solution (44) for λ = 0.2 were compared with the exact solution. The absolute values of relative errors δ(t) for different values of the variable t and β = 3/4; 4/5; 9/10; 19/20 are given in Table 3. The accuracy of the results depends on errors of the numerical determining of roots of the characteristic polynomial, errors of calculations of the coefficients and errors of numerical calculations of the integrals occurring in the solution (44).

Conclusions
The exact solution of the three-term fractional differential equation has been presented in the form of the infinite series of the generalized Wright functions and the Prabhakar functions. The properties of the Pochhammer symbol allow one to express the solution of the problem in terms of the generalized hypergeometric function, which is predefined in the Mathematica and Maple software packages. It should be noted that using the power series representation of the Prabhakar function in numerical calculations is acceptable for small arguments, while for large arguments, asymptotic expansions should be utilized. An alternative for the series solution is the solution obtained by using the analytical-numerical method based on the Laplace transform technique. This solution is expressed in the form of a finite sum, using the Mittag-Leffer type functions. Numerical values obtained using the presented procedure and the values calculated using the series solution for selected parameters of the fractional equation have been compared. The accordance of the results was obtained in the case when the calculations of the partial sums of the infinite series were done with high precision. Numerical examples show that the method can be used for high degrees of the characteristic polynomials obtained in solving the fractional commensurate equations. Although the presented approach concerns the three-term time-fractional equations, the method can be used, after modification, in solving multi-term fractional commensurate differential equations.