High order algorithms for numerical solution of fractional differential equations

In this paper, two novel high order numerical algorithms are proposed for solving fractional differential equations where the fractional derivative is considered in the Caputo sense. The total domain is discretized into a set of small subdomains and then the unknown functions are approximated using the piecewise Lagrange interpolation polynomial of degree three and degree four. The detailed error analysis is presented, and it is analytically proven that the proposed algorithms are of orders 4 and 5. The stability of the algorithms is rigorously established and the stability region is also achieved. Numerical examples are provided to check the theoretical results and illustrate the efficiency and applicability of the novel algorithms.


Introduction
The subject of fractional calculus (theory of integration and differentiation of arbitrary order) can be considered as an old yet novel topic. It has been an ongoing topic for more than 300 years; however, since the 1970s, it has been gaining increasing attention [1]. Firstly, there were almost no practical applications of fractional calculus (FC), and it was considered by many as an abstract area containing only mathematical manipulations of little or no use [2][3][4]. Recently, FC has been widely used in various applications in almost every field of science, engineering, and mathematics, and it has gained considerable importance due to its frequent applications in fluid flow, polymer rheology, economics, biophysics, control theory, psychology, and so on [5,6].
The main reason that fractional differential equations (FDEs) are being used to model real phenomena is that they are nonlocal in nature, that is, a realistic model of a physical phenomenon depends not only on the time instant but also on the previous time history [7]. In other words, fractional derivative (FD) provides a perfect tool when it is used to describe the memory and hereditary properties of various materials and processes [8,9]. Some of the other main differences between fractional calculus and classical calculus are: (i) FDEs are, at least, as stable as their integer order counterparts [10,11]; (ii) using FDEs can help to reduce the errors arising from the neglected parameters in modeling reallife phenomena [12,13]; (iii) in some situations, the FDE models seem more consistent with the real phenomena than the integer order models [14,15]; (iv) modeling of real-life processes with FD is more reliable and accurate than the classical derivative [16,17]; (v) fractional order models are more general [18] and in the limit results obtained from FC coincide with those obtained from classical calculus [19], and so on.
The wide applicability of FC in the field of science and engineering motivates researchers to try to find out the analytical or numerical solutions for FDEs. It is well known that the analytical and closed solutions of FDEs cannot generally be obtained and if luckily obtained always contain some infinite series (such as Mittag-Leffler function), which makes evaluation very expensive [20,21]. For this reason, necessarily, one may need an efficient approximate and numerical technique for the solution of FDEs [22].
Odibat et al. constructed a numerical scheme for the numerical solution of FDEs based on the modified trapezoidal rule and the fractional Euler method [23]. To present a numerical solution scheme for the fractional differential equations, authors of [24] divided total time into a set of small intervals and utilized quadratic interpolation polynomial between two successive intervals to approximate unknown functions. Cao and Xu applied quadratic interpolation polynomial to construct a high order scheme based on the socalled block-by-block approach for the fractional ordinary differential equations [25]. The convergence order of this scheme is 3 + α for 0 < α ≤ 1 and 4 for α > 1. Diethelm proposed an implicit numerical algorithm for solving FDEs by using piecewise linear interpolation polynomials to approximate the Hadamard finite-part integral [26]. Yan et al. designed a high order numerical scheme for solving a linear fractional differential equation defined in terms of the Riemann-Liouville fractional derivative [27]. This method is based on a direct discretisation of the fractional differential operator, and the order of convergence of the method is O(h 3-α ). A high order fractional Adams-type method for solving nonlinear FDEs is also obtained in this paper. Asl et al. [28] applied the piecewise cubic interpolation polynomial to construct a numerical scheme for solving linear FDEs. The convergence order of this scheme is O(h 4-α ) for 0 < α ≤ 1.
In recent years a lot of attention has been paid to finding effective numerical methods to simulate the fractional partial differential equations, such as meshless methods [29][30][31] and homotopy analysis method [32]. Dehghan et al. presented a fourth order numerical method for solving a fractional partial integro-differential equation [33]. Authors of [34] proposed a high accuracy numerical method to simulate a nonlinear fractional partial integro-differential equation.
The present study proposes to develop two numerical algorithms for the numerical approximation of a class of FDEs which are expressed in terms of Caputo-type fractional derivatives. In these algorithms the properties of the Caputo derivative are used to reduce the FDE into a Volterra-type integral equation of the second kind. The outline of the paper is as follows. Numerical algorithms are presented in Sect. 2 by using the piecewise Lagrange interpolation polynomial of degree three and degree four. Section 3 deals with the error analysis of the presented algorithms, and stability analysis of these algorithms is given on Sect. 4. Linear stability analysis of the presented schemes is given in Sect. 5 to achieve stability region of these methods. To demonstrate the effectiveness and high accuracy of the proposed methods, some numerical examples are provided in Sect. 6. The experimental orders of convergence are also investigated in this section. Finally, some concluding remarks are given in Sect. 7.

Numerical algorithms
The present paper further discusses numerical algorithms for the initial value problem for a FDE in the form where C 0 D α t denotes the Caputo fractional derivative(FD). The reason for choosing Caputo derivative in (1) is partly due to the fact that the Caputo derivative of a constant is zero.
More importantly, the Caputo FD determines an ordinary derivative followed by a fractional integral to obtain the preferred order of FDs. That means the Caputo FD allows us to couple FDEs with traditional (local) boundary and initial conditions. These conditions have clear physical interpretation and are familiar to us [8,35,36]. Motivated by these advantages, the Caputo FD has been applied to model various physical and real-life problems such as epidemic model for HIV transmission [8,9], symbiosis system [16], multi-pulse splitting process [17], reaction-diffusion problems [37,38], and tumor-host models [39].
Presented numerical algorithms are based on the analytical property that the initial value problem (1) is equivalent to the Volterra integral equation [40,41] where h(t) = α -1 j=0 t j j! y (j) (0), in the sense that a continuous function solves (2) if and only if it solves (1). The piecewise Lagrange interpolation polynomials of degree three and degree four are used to approximate the integral in (2). For an integer N and the given time T, the interval [0, T] is divided into t j = jh, j = 0, . . . , N , where h = T/N is the step length. The numerical solution of Eq. (1) at the point t j is denoted by y j . For notational convenience, let F(τ ) = f (τ , y(τ )) and F j = f (t j , y j ) = βy j + g(t j ).

Numerical algorithm I
We start with computing the value of y(t) at t 1 , t 2 , and t 3 , simultaneously. Consider the following integral for the first three steps (k = 0, 1, 2): whereF(τ ) is chosen to be the piecewise Lagrange cubic interpolation polynomial of F(τ ) associated with the nodes t 0 , t 1 , t 2 , and t 3 . In this way, we havê To explain further, for y 1 , we have For this reason, after some elementary calculations, y k+1 for the first three steps k = 0, 1, 2 can be approximated as follows: As it is mentioned above, the first three step solutions y 1 , y 2 , and y 3 are coupled in (5), thus need to be solved simultaneously. This procedure leads to a linear system of equations, which can be solved explicitly by means of back-substitution.
The following special cases should be excluded: In this way, after some explicit calculations, y k+1 for k ≥ 3 can be approximated as follows: Summing up the above arguments, we have the following novel scheme:

Numerical algorithm II
Consider the following integral for the first four steps (k = 0, 1, 2, 3): whereF(τ ) is the piecewise Lagrange interpolation polynomial of degree four associated with the nodes t 0 , t 1 , t 2 " t 3 , and t 4 . Therefore, one can achieve the following weights: j . Hence y k+1 for the first four steps k = 0, 1, 2, 3 can be determined as follows: It is obvious that the first four step solutions y 1 , y 2 , y 3 , and y 4 are coupled in (13), thus need to be solved simultaneously. This procedure leads to a linear system of equations, which can be solved explicitly by means of back-substitution.

Error analysis
For numerical algorithm I, the truncation error at step k + 1 is defined by [25] r k+1 (h) := y(t k+1 ) -ỹ k+1 , whereỹ k+1 is an approximation to y(t k+1 ), evaluated by using algorithm I (10) with exact previous solutions, i.e., for k ≥ 3, For numerical algorithm II (18), the definition of truncation error is the same as (19), wherẽ y k+1 for k ≥ 4 is as follows: Theorem 1 Let r k+1 (h) be the truncation error defined in (19). If F(τ ) ∈ C 4 [0, T] for some suitable chosen T , then for numerical algorithm I (10) there exists a positive constant C > 0, independent of h, such that Proof By the definition of truncation error (19) we have We compute y(t k+1 ) andỹ k+1 from (2) and (20), respectively, which yields We have, by (6), whereF(τ ) andF j+1 are defined by (6). Thus, according to the property of cubic interpolation polynomials, Here,τ j ∈ [t j , t j+1 ] and the second integral mean value theorem is used. By the definition of t j = jh and after evaluating the integrals in the above inequality, we get The series in the above expression can be easily summed up, which yields Theorem 2 Let r k+1 (h) be the truncation error defined in (19). If F(τ ) ∈ C 5 [0, T] for some suitable chosen T , then for numerical algorithm II (13) and (18) there exists a positive constant C > 0, independent of h, such that Proof The proof is similar to the proof of Theorem 1. We omit the proof here.

Stability analysis
The stability of a numerical scheme mainly refers to that if there is a perturbation in the initial condition, then the small change causes small errors in the numerical solution [40,41]. Suppose that y k+1 andỹ k+1 are numerical solutions in (10), and the initial conditions are given by y (i) 0 andỹ (i) 0 respectively. If there exists a positive constant C independent of h such that then we conclude that scheme (10) is stable [42]. It is similar to defining the numerical stability for numerical algorithm II (13) and (18). Assume that F(τ ) is sufficiently smooth and C α > 0 is independent of all discretization parameters. Firstly, we introduce two lemmas which will be used in stability analysis. (10), we have k+1 j=0 d k+1

Lemma 1 For the weights of novel scheme
where C α only depends on α.
Proof For d k+1 0 , we have where j -1 ≤ ξ j ≤ j, j = 1, 2, 3. Therefore we have Using similar analysis it can be shown that for j = 1, 2, 3, k -1, k, k + 1 there exists C α , which has dissimilar values at each case, such that the following inequality holds: For j = 4, 5, . . . , k -2, we have , andτ 4 ∈ [t j+2 , t j+3 ]. Hence, the above equation has the simplified form Combining all the above results, by choosing sufficiently large C α and also sufficiently small T, one can get (23) to complete the proof of the lemma.

Lemma 2 For the weights of the novel scheme
where C α only depends on α.
Proof The idea of the proof is similar to that of Lemma 1, so it is omitted.
Proof Suppose that y k+1 andỹ k+1 are numerical solutions in (10) and the initial conditions are given by y (i) 0 andỹ (i) 0 respectively. We shall use mathematical induction. Assume that is true for (j = 0, 1, . . . , k). We must prove that this also holds for j = k + 1. Note that, by assumptions of the given initial conditions, the induction basis (j = 0) is presupposed.
Utilizing Lemma 1, one can get Now, for sufficiently small T, one can end the proof using mathematical induction (26) and by choosing constant C α,T sufficiently large.
Then the presented algorithm II is stable.
Proof The proof is similar to Theorem 3.

Linear stability analysis
Following the ideas of [43,44], consider the following test problem to investigate stability region of the presented methods: where λ ∈ C is a complex number. The new method (13) and (18) gives the following iteration formula for solving (27): Denoting z = λh α , we get Let y j = ξ j , then by assuming ξ = e iθ with 0 ≤ θ ≤ 2π we get the following stability region for scheme (10): The stability region of algorithm I (10) can be achieved in a quite similar way. The stability region of numerical algorithm I is obtained in Figs

Numerical results
To check the numerical errors between the exact and the numerical solution, numerical experiments are carried out in this section. The presented examples have exact solutions and also have been solved by other numerical methods from literature. This allows one to compare the numerical results obtained by the presented schemes with the analytical solutions or those obtained by other methods. We also aim to graphically investigate the convergence order of the presented algorithms using the idea of [27] and [28]. The experimental orders of convergence (EOC) of the algorithms are measured by where E(h) is the absolute error |y(t j )y j | at t = 1 for the step size h. By the results of Theorems 1 and 2, we have where (order) is 4 for the presented algorithm (10) and is 5 for the presented algorithm (13) and (18). Before proceeding, we first define some notations. Let Y = log 2 (|E(h)|) and X = log 2 (h), also we use algorithm I to denote the presented algorithm (10) and algorithm II to denote the presented algorithm (13) and (18). We use E I and E II to denote the absolute error of algorithm I and algorithm II respectively. To graphically examine the order of convergence, we shall plot the straight line y = (order) × x and also the function Y = Y (X). The parallelism of these two lines can confirm the results of Theorems 1 and 2. Numerical computations are carried out with the help of MATLAB R2014a on an Intel(R) Core(TM)2 Duo CPU (2.67 GHz, 4 GB RAM) Windows 7 system.
Example 1 Consider the following fractional differential equation: where The exact solution is y(t) = t ν . At the time t = 1, for different step sizes h and different α, the approximate solutions for equation (31) are obtained for two values of ν = 4 and ν = 3.45 by using the presented algorithms. First, we choose ν = 4, which for this case y(t) is a sufficiently smooth function. The absolute errors of the presented algorithms and of the method reported in Ref. [45] are shown in Table 1. This table shows that novel schemes are valid methods in solving a fractional differential equation. Although the convergence ratios in Table 1 are a bit slower, they are not too far from our analytical results. This can be mainly due to some round-off   Example 2 Consider the following fractional differential equation: The exact solution is y(t) = t 4 -1 2 t 3 . Table 3 shows the absolute errors of the presented schemes and the method reported in Ref. [46] at the time t = 1. From this table it is observed that the error of the presented method has decreased significantly. In Figs. 9-10,  Example 3 In this example we aim to examine the validity and efficiency of the presented methods when they are applied to solve a fractional order system of differential equations.
We solve system (33) using the presented algorithms numerically. Figure 11 shows the difference between the curves of exact and approximate solutions for α = 0.925 with h = The experimentally determined orders of convergence for algorithm I at t = 1 for (32) 1/80. It is observed that the numerical solutions obtained by algorithm I and algorithm II are in excellent agreement with the exact solution. Figure 12(a) shows the exact and approximate solutions of system (33) obtained by algorithm I for t = 50 with α = 0.95, h = 1/40. Figure 12(b) shows the exact and approximate solutions of system (33) obtained by algorithm II for t = 50 with α = 0.975, h = 1/40. From Fig. 12 it is observed that the numerical solutions are at very good agreement with the exact solution and the presented algorithms are still valid and capable of solving problems in a large scale computational domain.

Conclusion
This paper provides two high order numerical schemes with theoretically proved convergence order of 4 and 5 for solving FDEs. The properties of the Caputo derivative are used to reduce FDEs into a Volterra integral equation. After dividing the total domain into a set of grid points, the piecewise Lagrange interpolation polynomials of degree three and degree four are utilized to approximate unknown functions. The stability and error estimate of the methods are investigated. Moreover, graphical illustrations for the stability region of the schemes are derived. In the numerical results, the accuracy and convergence orders Figure 10 The experimentally determined orders of convergence for algorithm II at t = 1 for (32) of the developed algorithms have been studied. Numerical results confirm the theoretical results of the algorithms. Furthermore, we extended the presented algorithms to solve a 3D fractional order system of equations. It is shown that the algorithms are valid and efficient for solving a linear fractional system, and it is observed that the algorithms are capable of solving large domain problems. In the future, we shall try to follow this idea to construct higher order schemes for solving nonlinear FDEs.