A novel numerical approach for solving fractional order differential equations using hybrid functions

Abstract: This article presents a novel numerical method for seeking the numerical solutions of fractional order differential equations using hybrid functions consisting of block-pulse functions and Taylor polynomials. The fractional integrals operational matrix of the hybrid function is conducted through projecting the hybrid functions onto block-pulse functions. Then, the fractional order differential equations are converted to a set of algebraic equations via the derived operational matrix. Then, the numerical solutions are obtained via solving the algebraic equations. Moreover, we perform error analysis of the algorithm and gives the upper bound of absolute error. Finally, numerical examples are given to show the effectiveness of the proposed method.


Introduction
Fractional order differential equations (FODEs) are the extension of classical integer order differential equations (IODEs) using the concept of fractional calculus (FC) [1]. Compared to integer calculus (FC), FC is non-local and long-time memory, which enables FODEs be an alternative powerful mathematical tool for modeling the nature of many physical systems and real processes. For example, the relaxation modulus in viscoelastic material exhibits power-law behavior, which indicates that in the relaxation process the material has history memory. In this case, the relationship between stress and strain can be modeled using a FODE as σ(t) = Eτ α d α ε(t) dt α [2]. Also, the Warburg impedance has a fractional-power-law dependence on angular frequency, which is modeled as Z(s) = As −1/2 [3].
In the last decades, linear or nonlinear FODEs have emerged in many scientific and engineering fields such as neural networks [4,5], chaos [6], diffusion process [7], fluid flows in porous media [8] and optimal control [9]. Now, more and more applications of FODEs can be found in practice, there is a great need of finding the solution of FODEs. However, it is not an easy thing to obtain the exact or analytic solution for most FODEs, pursuing numerical solution of FODEs is an important thing. In the past, the finite difference method [10], predictor-corrector method [11], Adomian decomposition method [12], variational iteration method [13] and Homotopy analysis method [14] had been proposed to solve FODES.
Though those methods are efficient and can provide good approximation to the real solution, however, they possess high computational complexities, which hampers its real application. So, developing efficient and effective methods for solving FODEs becomes an urgent task. Recently, the operational matrices had been widely adopted by researchers to solve different kinds of FODEs [15]. The core idea of operational matrices based methods is to transform the FODEs to a set of algebraic equations. As a result, the problem is dramatically simplified. In the literature, the operational matrices of fractional operators of block pulse functions (BPFs) [16], Bernstein polynomials [17], Bernoulli polynomials [18], Taylor polynomials [19] had been developed and adopted by several scholars to solve FODEs numerically. Among the above different polynomials or functions, the BPFs' fractional integral operational matrix is upper triangular matrix, in which the k−th row can be obtained by shifting the (k − 1)−th row to the right. This an appealing property which enable one to solve the set of algebraic equations more easily and drastically reduce the computation burden involved in the computation of matrix inverse. However, BPFs are not enough smooth, more BPFs are required if high approximation accuracy is desired, and in turn, the dimension of operational matrix increases.
In recent years, the hybrid functions, which mainly combines BPFs or Haar function with other polynomials, had been widely adopted by researchers to solve different kind of FODEs. For example, in [20], the hybrid of BPFs and Bernoulli polynomials is adopted to find the numerical solution of nonlinear fractional integro-differential equations. In [21], the hybrid of Legendre polynomial and Haar function, called Legendre wavelet is used to solve distributed order FODEs. In [22], Chebyshev wavelets is adopted to solve nonlinear variable order FODEs. One of the advantage of hybrid functions is that they are piecewise polynomial instead of piecewise constant in each interval as BPFs. Thus, hybrid functions are smoother than BPFs, which leads to a more accurate approximation to a function than BPFs if an equal number of basis functions are used. Therefore, one can obtain more accurate solution for a FODEs using hybrid functions than BPFs.
Observing the aforementioned facts, this paper presents a new numerical method for solving FODEs based on hybrid of BPFs and Taylor polynomials (HBT). The fractional integrals operational matrix of HBT is derived through projecting the HBT functions onto a set of BPFs. Then, the FODEs to be solved is transformed into a set of algebraic equations using the derived operational matrix. Through solving the algebraic equations, one can obtain the numerical solution of FODEs.
The organization of this paper is arranged as follows. Some basics about fractional calculus are briefly reviewed in Section 2. In Sections 3, the basic formulation of HBT functions is given, and the calculation of fractional integral operational matrix of HBT is given in Section 4. The error analysis of the proposed method is given in Section 5. In Section 6, numerical experiments are presented to verify the effectiveness of the proposed method. In the last, the conclusive remarks are presented in Section 7.

Definition of fractional calculus
Unlike integer calculus, the definition of FC is not unique. There are several definitions for fractional integration and derivatives. Among these definitions, the widely used one are the Riemann-Liouville (R-L) definition and the Caputo definition.
In R-L definition, the fractional integral of function f (t) is given as [23] where α > 0 is the order, 0 and t are the low and upper limits of fractional integrals, Γ(·) denotes the Gamma function. The fractional derivative of a function f (t) in the sense of Caputo definition is (2.2) R-L fractional integral and Caputo fractional derivative has the following relationship and Since the derivation is the inverse operation of integration, the derivation of a function can easily be achieved from the result of integration.

BPFs and HBT functions
The Taylor polynomials defined on the interval [0, T ] is [24] T j (t) = t j , j = 0, 1, 2, · · · (3.1) A set of BPFs ψ i (t), i = 1, 2, · · · N is given as follows [25] The BPFs ψ i (t) are disjoint and orthogonal, that is to say  A function f (t) in C n+1 [0, T ) can be expanded into an N-term BPF series as T , the constant coefficients c i in Eq (3.5) are given by On the interval [0, T ), the hybrid of BPFs and Taylor polynomials are defined as [24] where i and j are the order of BPFs and the degree of Taylor polynomials respectively.

Function approximation
Let Since f * (t) ∈ X, there exist the unique coefficients c 10 , c 11 , · · · , c N(M−1) such that where T denotes transpose of vector or matrix, C and H(t) are N × M column vectors and To evaluate C, firstly let Using Eq (3.8) one gets Therefore, For HBT functions D has the following form and Hence, C T in Eq (3.8) is given by (3.14)

Fractional integral operational matrix of HBT
Here, the fractional integral operational matrix of HBT is derived. To simplify the derivation process, the HBT functions are first projected onto a set of BPFs, and the projection matrix is calculated. The HBT operational matrix is obtained by using matrix operations of the projection matrix and BPF operational matrix.
The R-L fractional integration operator of hybrid function vector H(t) can be written as where matrix P α is the fractional integral operational matrix of HBT. To simplify the derivation process, the hybrid functions are projected onto a set of BPFs ψ i (t), i = 1, 2, · · · , m, and m = N ×M. Specifically, the hybrid function vector defined in Eq (3.10) can be expressed as where Ψ(t) = [ψ 1 (t), ψ 2 (t), · · · , ψ m (t)] T and Φ is the projection matrix which transform the hybrid functions onto BPFs. In general, for arbitrary M and N, one has The expressions of a i j , b i j and d i j are given in Appendix A. As an example, let N = 2, M = 3 and T = 1, one has the function vector (4.4) can be expressed as In this case, The projection matrix is derived in detail in Appendix A and omitted here. It is easy to see from the above example that the obtained projection matrix is in block-diagonal form. This elegant property is useful for simplifying the computation of the operational matrices. From Eq (4.1) and Eq (4.2), one has According to Ref. [25], where is the BPFs operational matrix of fractional integration with Substituting Eq (4.7) into Eq (4.6), one can get Therefore, P α = ΦF α Φ −1 . (4.10)

Error analysis
Denote the HBT approximation to function f (t) at the level l as f where l is determined by N and M (l = l N,M ). We replace f (t) with f * (t) in Eq (2.1) and call the resulting integral HBT approximation of the α− order R-L fractional integral of f (t).
then the absolute error between I α f (t) and I α f * (t) is the truncation error is where ξ i lies between i−1 N and t. Since the best approximation is unique [26], for all t ∈ [ i−1 N T, i N T ), the absolute error between f (t) and f * (t) is given as . The absolute error between I α f (t) and I α f * (t) is according to the assumption | f (M) (t) |≤ K, the absolute error between I α f (t) and I α f * (t) can be estimated as where K is a finite positive value.
In order to verify the maximum absolute error arised by HBT approximation is smaller than the upper bound in theory given in Theorem 1, the function f (t) = t is selected as an example. The analytic expression of α order fractional integral of f (t) = t is given as Using Eq (4.9) with α = 0.5, T = 1 and N = 2, M = 3, the HBT estimation of ( 0 I α t )(t) is obtained as  Table 1 presents the absolute errors obtained by the approximation of BPFs and HBT functions. The piecewise-polynomial property of HBT functions lead to a more accurate approximation of the fractional integral even with the equal number of basis functions as BPFs. Therefore, HBT approximation of R-L fractional integral is effective.

Results
To show the effectiveness of our presented method, the fractional integral operational matrix of HBT is used to solve several FODEs. The solutions obtained by our proposed method are compared with their exact solutions, those by block-pulse operational matrix (BPOM) method and fractional Taylor operational matrix (FTOM) method [19].

Example 1
A fractional Riccati equation as [28,29] together with the initial condition, one has Since H(t) = ΦΨ(t), from Eq (6.3) one has , a 2 , · · · , a m ], (6.5) and using Eqs (3.3) and (3.4) one can obtain Substituting Eqs(4.2), (6.2), (6.4) and (6.6) into Eq (6.1), one has Equation (6.7) is a set of algebraic equations. In this paper, the Matlab function fsolve is used to solve Eq (6.7). The numerical solution, for N = 4, M = 6 is shown in Figure 1. When α = 1, the exact solution of this equation is  On the other hand, the numerical solutions obtained by the present method with N = 4, M = 6 for α = 0.7, 0.8, 0.9, 1 are shown in Figure 1. It is demonstrate that the approximate solution obtained by the present method is in good agreement with the exact solution when α = 1. Moreover, as indicated in [28], the solution continuously depends on the time-fractional derivative.

Example 2
A fractional differential equation as and subject to x(0) = x (0) = 0 is considered [30]. The exact solution of this equation is (6.9) and take the initial states into consideration, one has D α 2 x(t) = C T P 2−α 2 H(t), (6.10) Since H(t) = ΦΨ(t), from Eq (6.9) we have Expanding f (t) onto HBT functions, one has where f T is a known constant vector. Substituting Eq (6.9)-(6.11) and Eq (6.14) into Eq (6.8), together with H(t) = ΦΨ(t) , then   Table 3. The computational results illustrate that, compared with BPOM method and FTOM method, the present method can get more accurate approximate solutions although it has a large amount of calculation. Moreover, for the present method, the error is smaller and smaller when N and M increase. Therefore for higher accuracy of the approximation, lager N and M are recommended.

A. The projection matrix of HBT functions
In this section, the derivation process of the projection matrix from HBT functions to BPFs is given. Note that Eq (3.10) is the HBT function vector. First, elements h 1 j (t) of H(t) are considered and can be expanded into BPFs as follows: where j = 0, 1, . . . , (M − 1) and m = N × M. The coefficient a ( j+1)k can be computed according to Eq (3.6) Since h 1 j (t) is nonzero when 0 ≤ t < T N , and ψ k (t) is defined on the interval [ k−1 N M T, k N M T ), h 1 j (t)ψ k (t) is nonzero only when 1 ≤ k ≤ M. Observing this fact, Eq (A.2) can be written as otherwise.
After some manipulations, it is easy to obtain otherwise. Define