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 ef-ﬁciency 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 is 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 have gained considerable importance due to their frequent appearance applications in fluid flow, polymer rheology, economics, biophysics, control theory, psychology and so on 10 [5, 6,7].
The main reason that fractional differential equations (FDEs) are being used to modeling 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 the previous time history [8]. In other words, fractional derivative provides a perfect tool when it is used 15 to describe the memory and hereditary properties of various materials and processes [9,10]. Some of the other main differences between fractional calculus and classical calculus are: (1) FDEs are, at least, as stable as their integer order counterpart [11,12]; (2) Using FDEs can help to reduce the errors arising from the neglected parameters in modeling real-life phenomena [13,14]; (3) In some situations, FDEs models seem 20 more consistent with the real phenomena than the integer-order models [15,16]; (4) Fractional order models are more general [17] and in the limit results obtained from FC coincide with those obtained from classical calculus [18] 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 the FDEs. It is well 25 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 make evaluation very expensive [19? , 20]. For this reason, necessarily, one may need an efficient approximate and numerical technique for the solution of FDEs [21]. 30 Odibat et al. constructed a numerical scheme for the numerical solution of FDEs based on the modified trapezoidal rule and the fractional Euler's method [22]. To obtain a numerical solution scheme for the fractional differential equations, authors of [23] divided the time interval into a set of small subintervals, and utilized quadratic interpolation polynomial between two successive intervals to approximate unknown 35 functions. Cao and Xu applied quadratic interpolation polynomial to construct a high order scheme based on the so-called block-by-block approach, for the fractional ordinary differential equations [24]. 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 approximated the 40 Hadamard finite-part integral [25]. Yan et al. designed a high order numerical scheme for solving a linear fractional differential equation by approximating the Hadamard finite-part integral with the quadratic interpolation polynomials [26]. 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 45 for solving a nonlinear FDEs is also obtained in this paper. Pal et al. designed an extrapolation algorithm for solving linear FDEs based on the direct discretization of the fractional differential operator [27].
In this paper, we will introduce two new numerical algorithms for solving the nonlinear FDEs, which are expressed in terms of Caputo type fractional derivatives. In 50 these algorithms properties of the Caputo derivative are used to reduce the FDE into a Volterra type integral equation of the second kind. We then use the Lagrange interpolation polynomials of degree three and four to approximate the integral and the proposed numerical algorithms has the truncation error O(h 4 ) and O(h 5 ) for all α > 0. The stability of the numerical method is proved based on the properties of the weights in the 55 numerical algorithm under the assumption that the time T > 0 is sufficiently small.
Such properties are used in the first time to prove the stability of the numerical methods for solving fractional differential equations. To our best knowledge, there is no numerical algorithm for solving nonlinear fractional differential equation with the convergence order greater than 4 in the literature. We also introduce a new way to analyze 60 the stability of the numerical methods for solving fractional differential equations.
The outline of the paper is as follows. Numerical algorithms are presented in Section 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 in Section 4. Linear stability analysis of the proposed schemes is given in Section 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 Section 6. Finally, Some concluding remarks are given in Section 7.

70
Consider the nonlinear fractional differential equation, with α > 0, where C 0 D α t denotes the Caputo fractional derivative and f (t, u) satisfies the Lipschitz condition with respect to the second variable, i.e., there exists a constant L > 0 such It is well known that the initial value problem (1) is equivalent to the Volterra integral equation where h(t) = (1) at the point t j is denoted by y j . For 75 notational convenience, let F(τ) = f (τ, y(τ)) and F j = f (t j , y 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ê 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 (6), thus need to be solved simultaneously. An explicit solution of these three equations is given in Appendix A section.

80
To construct the scheme for the next steps, I k+1 , k ≥ 3 is descritized as follows in which like as (4), for the first three integrals ( j = 0, 1, 2, 3),F is the piecewise Lagrange cubic interpolation polynomial of F(τ) associated with the nodes t 0 , t 1 , t 2 and t 3 . For the reminder integrals ( j = 3, 4, . . . , k + 1),F j+1 is chosen to be the piecewise Lagrange cubic interpolation polynomial of F(τ) associated with the nodes t j−2 , t j−1 , t j and t j+1 . In this way, for k ≥ 3 we havê The following special cases should be excluded: For this reason, after some explicit calculations y k+1 for k ≥ 3 can be approximated as follows: To summarize, we obtain the following novel scheme: where d k+1 j are defined as above.

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 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 (14), thus need to be solved simultaneously. An explicit solution of these four equations is given in Appendix B section.

85
To design the schema for the next steps, I k+1 , k ≥ 4 is descritized as follows in which like as (12), for the first four integrals ( j = 0, 1, 2, 3),F is the piecewise Lagrange interpolation polynomial of degree four associated with the nodes t 0 , t 1 , t 2 , , t 3 and t 4 . For the reminder integrals ( j = 4, 5, . . . , k + 1),F j+1 is the piecewise Lagrange interpolation polynomial of degree four associated with the nodes t j−3 , t j−2 , t j−1 , t j and t j+1 . In this way, for k ≥ 4 we have the following weights The following special cases should be excluded: Therefore y k+1 for k ≥ 4 can be approximated as follows: Thus a new numerical algorithm II is described by (14) and (19) with the weights b k+1 j defined as above.

Error analysis
For the numerical algorithm I the truncation error at the step k + 1 is defined by whereỹ k+1 is an approximation to y(t k+1 ), evaluated by using the algorithm I (11) with exact previous solutions, i.e. for k ≥ 3, For the numerical algorithm II (19), the definition of truncation error is the same as (20), whereỹ k+1 for k ≥ 4 is as follows: Proof. We have, by Eqs.
Theorem 2. Let r k+1 (h) be the truncation error defined in (20). If F(τ) ∈ C 5 [0, T ] for some suitable chosen T , then for the numerical algorithm II (14) and (19) there exists a positive constant C > 0, independent of h, such that Proof. The details of the proof is similar to that of Theorem 1 so are neglected. We have, by Eqs.

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 cause small errors in the numerical solution [28,29]. Suppose that y k+1 andỹ k+1 are numerical solutions in (11), and the initial conditions are given by y then it concluded that the scheme (11) is stable [30]. It is similar to define the numerical stability for the numerical algorithm II (14) and (19). 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.

95
Lemma 1. For the weights of the novel scheme (11) we have 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 exist C α , which is dissimilar values at each cases, such that the following inequality is holds.
For j = 4, 5, . . . , k − 2 we have Hence, above equation has the simplify form, Combining all above results, by choosing sufficiently large C α , and also sufficiently small T one can get (24) to complete the proof of the Lemma.

Lemma 2.
For the weights of the novel scheme (19) we have where C α only depends on α.
Proof. The idea of the proof is similar to that of Lemma 1, so is omitted.
Proof. Suppose that y k+1 andỹ k+1 are numerical solutions in (11), 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 true. We have, using the Lipschitz condition assumption (2), By Lemma 1 one can get which implies that Now for sufficiently small T , one can complete the proof by using the mathematical induction (27) and by choosing constant C α,T sufficiently large.
Proof. The proof is similar to the proof of Theorem 3.

Linear stability analysis
Consider the following test problem to investigate stability region of the presented methods: The new method (11) gives the following iteration formula for solving (28): Denoting z = λ h α , we get Let y j = ξ j , then by assuming ξ = e iθ with 0 ≤ θ ≤ 2π we get the following stability region for the scheme (11) The stability region of the algorithm II (14) and (19) can be achieved in a quite similar way. The stability region of the numerical algorithm I is obtained in Figs. (1) and (2)

Numerical results
To check the numerical errors between the exact and the numerical solution, nu- Example 1. Consider the following fractional differential equation where The exact solution is y(t) = t 4 . At the time t = 1 for different step sizes h and different The exact solution to this initial value problem is y(t) = t 2 − t. Example 3. Consider the following fractional differential equation The exact solution is y(t) = t 4 − 1 2 t 3 . Table 7 shows the absolute errors of the presented schemes and the method reported in Ref. [27] at the time t = 1. From this Table it is observed that the error of presented method is decreased significantly. we shall try to follow this idea to construct higher order schemes for solving nonlinear FDEs.

Appendix A
The idea of solving y 1 , y 2 and y 3 form (6) is as follows. For simplicity, we assume that f (t, y) = µy + g(t) for understanding the idea of the numerical method. We have the following linear system of equations, from (6), , in which Substituting (35) and (38) to (37) leads to in which Now, firstly one can calculate y 3 from given initial conditions and known function g(t).
Then y 2 and y 1 can be calculated by (38) and (35), respectively.