Numerical Schemes for Fractional Ordinary Differential Equations

Fractional calculus, which has almost the same history as classic calculus, did not attract enough attention for a long time. However, in recent decades, fractional calculus and fractional differential equations become more and more popular because of its powerful potential applications. A large number of new differential equations (models) that involve fractional calculus are developed. These models have been applied successfully, e.g., in mechanics (theory of viscoelasticity), biology (modelling of polymers and proteins), chemistry (modelling the anomalous diffusion behavior of Brownian particles), electrical engineering (electromagnetic waves) etc (Bouchaud & Georges, 1990; Hilfer, 2000; Kirchner et al., 2000; Metzler & Klafter, 2000; Zaslavsky, 2002; Mainrdi, 2008). Meanwhile, some rich fractional dynamical behavior which reflect the inherent nature of realistic physical systems are observed. In short, fractional calculus and fractional differential equations have played more andmore important role in almost all the scientific fields. One of the most important fractional models is the following initial value problem


Introduction
Fractional calculus, which has almost the same history as classic calculus, did not attract enough attention for a long time. However, in recent decades, fractional calculus and fractional differential equations become more and more popular because of its powerful potential applications. A large number of new differential equations (models) that involve fractional calculus are developed. These models have been applied successfully, e.g., in mechanics (theory of viscoelasticity), biology (modelling of polymers and proteins), chemistry (modelling the anomalous diffusion behavior of Brownian particles), electrical engineering (electromagnetic waves) etc (Bouchaud & Georges, 1990;Hilfer, 2000;Kirchner et al., 2000;Metzler & Klafter, 2000;Zaslavsky, 2002;Mainrdi, 2008). Meanwhile, some rich fractional dynamical behavior which reflect the inherent nature of realistic physical systems are observed. In short, fractional calculus and fractional differential equations have played more and more important role in almost all the scientific fields. One of the most important fractional models is the following initial value problem D α * y(t)= f (t, y(t)), y (k) (x 0 )=y (k) 0 , k = 0, 1, ··· , ⌈α⌉−1, where α ∈ (0, ∞), y (k) 0 can be any real numbers, and D α * denotes the fractional derivative in the Caputo sense, defined by D α * y(t)=J n−α D n y(t)= here n := ⌈α⌉ is the smallest integer not less than α, y (n) (x) is the classical nth-order derivative of y(x) and for μ > 0,J μ is the μ-order Riemann-Liouville integral operator expressed as follows The use of Caputo derivative in above equation is partly because of the convenience to specify the initial conditions. Since the initial conditions are expressed in terms of values of the unknown function and its integer-order derivatives which have clear physical meaning (Podlubny, 1999;. However, from pure mathematical viewpoint, the Riemann-Liouville derivative is more welcome and many earlier research papers use it instead of Caputo derivative (Podlubny, 1999). In general, specifying some additional conditions is necessary to make sure the discussed equation has a unique solution. These additional conditions, in many situations, describe some properties of the solution at the initial time (Heymans & Podlubny, 2005), the fractional derivative does not have convenient used physical meaning (there are already some progress in the geometric and the physical interpretation of the fractional calculus (Podlubny, 2002) and the physical interpretation of the initial conditions in terms of the Riemann-Liouville fractional derivatives of the unknown function has also been discussed in (Podlubny, 2002)).
Just like the classic calculus and differential equations, the theories of fractional differentials, integrals and differential equations have been developing. With the development of the theories of fractional calculus, many research monographs are published, e.g., (Oldham & Spanier, 1974;Podlubny, 1999;Samko et al., 1993). In the literatures, several analytical methodologies, such as, Laplace transform, Mellin transform, Fourier transform, are restored to obtain the analytical solutions of the fractional equations by many authors (Metzler & Klafter, 2000;Podlubny, 1999;Samko et al., 1993;Zaslavsky, 2002;Mainrdi, 2008), however, similar to treating classical differential equations, they can mainly deal with linear fractional differential equation with constant coefficients. Usually, for nonlinear systems these techniques do not work. So in many cases the more reasonable option is to find its numerical solution. As is well known, the difficulty of solving fractional differential equations is essentially because fractional calculus are non-local operators. This non-local property means that the next state of a system not only depends on its current state but also on its historical states starting from the initial time. This property is closer to reality and is the main reason why fractional calculus has become more and more useful and popular. In other words, this non-local property is good for modeling reality, but a challenge for numerical computations. Much effort has been devoted during the recent years to the numerical investigations of fractional calculus and fractional dynamics of (1) (Lubich, 1985;1986;Podlubny, 1999). More recently, Diethelm et al successfully presented the numerical approximation of (1) using Adams-type predictor-corrector approach (Diethelm & Ford, 2002a) and the detailed error analysis of this method was given in (Diethelm et al., 2004). The convergent order of Diethelm's predictor-corrector approach was proved to be min(2, 1 + α) (Diethelm et al., 2004). Because of the non-local property of the fractional derivatives, the arithmetic complexity of their algorithm with step size h is O(h −2 ), whereas a comparable algorithm for a classical initial value problem only gives rise to O(h −1 ). To improve the accuracy and reduce the arithmetic complexity, some techniques such as the Richardson extrapolation, short memory principle and corresponding mixed numerical schemes are developed. In (Deng, 2007a), we present an improved version of the predictor-corrector algorithm with the accuracy increased to min(2, 1 + 2α) and half of the computational cost is reduced comparing to the original one in (Diethelm et al., 2004). Furthermore, we apprehend the short memory principle from a new viewpoint (Deng, 2007b); after using the nested meshes presented in (Ford & Simpson, 2001) and combining the short memory principle and the predictor-corrector approach, we minimize the computational complexity to O(h −1 log(h −1 )) at preserving the order of accuracy 2.
This chapter briefly reviews the recent development of the predictor-corrector approach for fractional dynamic systems. The plan of this chapter is as follows. In Section 2, we briefly discuss the short memory principle and the nested meshes. In Section 3, the predictor-corrector schemes and its improved versions are presented, meanwhile the convergent order and arithmetic complexity are also proposed. In Section 4, we provide two numerical examples to illustrate the performance of our numerical schemes. Conclusions are given in the last section.

Short memory principle
We can see that the fractional derivative (2) is an operator depending on the past states of the process y(t) (see Fig 1). However, for t ≫ 0 and 0 < α < 1, the behavior of y(t) far * () Dyt a 0 t the "" past of () yt } Fig. 1. The fractional derivative operating on the "past" of y(t) (the red part).
from upper terminal suggests that the original fractional derivative can possibly be replaced by the fractional derivative with a moving lower terminal (Podlubny, 1999). This means that the history of the process y(t) can be approached only over a fixed period of recent history. Inspired by this kind of idea, Podlubny in (Podlubny, 1999) introduces a fixed integral length memory principle for Riemann-Liouville derivative (see Fig 2). He shows that the truncation The "short-memory" principle with "fixed memory length" L(< t), where 0 D α t denote the Riemann-Liouville derivative operator.
error gives E < ML −α /Γ(1 − α) with a fixed integral length L. To accelerate the computation without significant loss of accuracy, in (Ford & Simpson, 2001) Ford and Simpson present a short-memory principle for Caputo derivative. In (Deng, 2007b) we apprehend the short memory principle from a new viewpoint, and then it is closer to reality and extends the effective range of short memory principle from α ∈ (0, 1) to α ∈ (0, 2). In view of the scaling property of fractional integral, the nested meshes are possibly used. And the nested meshes can be produced by splitting the integral interval [0, where T = ωh, h ∈ R + , m, ω, p ∈ N and p m T≤t n < p m+1 T . Denote M h = {hn, n ∈ N} and l 1 , l 2 ∈ N, l 1 > l 2 , then M l 2 h ⊃ M l 1 h .

Predictor-corrector schemes
The problem of determining y(t) by means of the information of initial value y 0 , is called initial value problem. As the classic theory of ordinary differential equation, if a function y(t) satisfies the initial value problem (1) pointwisely, then y(t) is called the solution of the fractional differential equation (1). For the existence and uniqueness of the solutions of fractional differential equations, one can see (Podlubny, 1999). But the explicit formula of the solution y(t) can't usually be obtained in spite of we can prove the existence of solution. The most 357 Numerical Schemes for Fractional Ordinary Differential Equations www.intechopen.com reasonable way is to use numerical methods, and the obtained solution is called the numerical approximation of the solution of the differential equations (1) and denoted by y h in the following sections.

Numerical schemes
In this section, we show the predictor-corrector schemes of (1). Using the Laplace transform formula for the Caputo fractional derivative (Podlubny, 1999) From (1), we have or where F(s, Y(s)) = L( f (t, y(t))). Applying the inverse Laplace transform gives where the fact are used. The approximation is based on the equivalent form of the Volterra integral equation (7). A fractional Adams-predictor-corrector approach was firstly developed by Diethelm et al (Diethelm et al., 2002b) to numerically solve the problem (7). Using the standard quadrature techniques for the integral in (7), denote g(τ)= f (τ, y(τ)), the integral is replaced by the trapezoidal quadrature formula at point t n+1 where g n+1 is the piecewise linear interpolation of g with nodes and knots chosen at t j , j = 0, 1, 2, ··· , n + 1. After some elementary calculations, the right hand side of (8) gives where the uniform mesh is used and h is the stepsize. And if we use the product rectangle rule, the right hand of (8) can be written as where Then the predictor and corrector formulae for solving (7) are given, respectively, by and The approximation accuracy of the scheme (11)-(12) is O(h min{2,1+α} ) (Diethelm et al., 2004). Now we make some improvements for the scheme (11)-(12). We modify the approximation of (8) as (Deng, 2007a) where g n is the piecewise linear interpolation for g with nodes and knots chosen at t j , j = 0, 1, 2, ··· , n. Then using the standard quadrature technique, the right hand of (13) can be recast as Hence, this algorithm for the predictor step can be improved as (Deng, 2007b) The new predictor-corrector approach (15) and (12) has the numerical accuracy O(h min{2,1+2α} ) (the detailed analysis is left in Section 3.2). Obviously half of the computational cost can be reduced, for 0 < α 1, if we reformulate (15) and (12) as and The arithmetic complexity of the above two predictor-corrector schemes ((11)-(12), and (16)-(17) For further reducing the computational cost, we understand the short memory principle from a new viewpoint, and then go to design the predictor-corrector scheme (Deng, 2007b). For α ∈ (0, 1) and α ∈ (1, ∞), we rewrite (7) as, respectively, and By observation of (18) and (19), we can see that the non-local property of D α * induces the term In fact, if α ∈ (0, 2) the integration kernel of (20) fades faster when the time history becomes longer, more concretely, Obviously, we can see that the integration (20)'s kernel (t n+1 − τ) α−1 − (t n − τ) α−1 decays (algebraically) by the order 2 − α when α ∈ (0, 2), but in (Ford & Simpson, 2001) the integral kernel (t n+1 − τ) α−1 decays with the order 1 − α when α ∈ (0, 1). This is the main reason why we can extend the range of the short memory principle of fractional differential equations from α ∈ (0, 1) to α ∈ (0, 2). For the nested mashes defined by (3), we can take the step length h in the integral [t n − pT , t n ] and in the subsequent is ignored, it does not destroy the computational accuracy in general. We will pay our attention to the numerical approximation by using short memory principle in the integration (20) in the following part.
The predictor and corrector formulae based on the analytical formula (19) is fully described by (28) and (29) with the weighted a j,n defined by (24). It is apparent that the same term ∑ n j=0 a j,n f (t j , y h (t j )) exists in both predictor (28) and corrector formulae (29), we just need to compute one times at each predictor-corrector iteration step.
In order to reduce the computational cost, in conjunction with the nested meshes memory principle (3), decompose the integral and still use the product trapezoidal quadrature formula at each subinterval but with different step lengths.
Employing above analysis we obtain the following predictor-corrector algorithm and y h (t n+1 )=y (1) In this case, the nested meshes predictor-corrector algorithm has the computational cost of O(h −1 log(h −1 )) for α ∈ (1, 2).

Error analysis and convergent order
In this section, we make the local truncation error and convergent order analysis for the improved predictor-corrector approaches (16)- (17), (26)-(29), and (32)-(33). First we present several lemmas which will be used later.
With the similar method, we can obtain the following lemma.
where C α only depends on α.
The error of the product rectangle rule is given as Lemma 3. Suppose that ∂ f (τ, y(τ))/∂t ∈ C[0, t), for some suitable t, then we have where C α only depends on α. Proof.
Lemma 4. Suppose that ∂ 2 f (τ, y(τ))/∂ 2 τ ∈ C[0, t), for some suitable t, then we have we have (41)) Here ξ, η ∈ [t n , t n+1 ] and in the above equalities the second integral mean value theorem and the properties of second divided differences are used.
Proof. The idea of this lemma's proof is similar to the above lemmas, namely Lemma 6. (Diethelm et al., 2004) Assume that the solution of the initial value problem satisfies with some γ 1 , γ 2 ≥ 0 and δ 1 , δ 2 ≥ 0. Then, for some suitably chosen T ≥ 0, we have with q = min{δ 1 + α, δ 2 } and N = ⌊T/h⌋. Theorem 1. For the fractional initial problem (1), if D α * y(t) ∈ C 2 [0, T] for some suitable T, then the convergent order of our algorithm with the predictor and corrector formulae (16) and (17) gives and applying lemmas 2, 4, 6, we have the above result.

Numerical Schemes for Fractional Ordinary Differential Equations 13
Proof. This proof will be used based on mathematical induction. In view of the given initial condition, the induction basis j = 0 is presupposed, it has convergent order 2. Now, assuming that the convergent order is 2 for j = 0, 1, ··· , k, k ≤ n, we have the local truncation error a j,n f (t j , y(t j )) a j,n f (t j , y h (t j )) .

Then we have
where z * ∈ (t n , t n+1 ), Lemmas 4 and 5 in the above proof are used, and we also utilize the result |y(t n+1 ) − y P h (t n+1 )| = O(h min{α+1,3} ) which can be proved by using Lemmas 4 and 5 and the similar idea to above proof, its sketch proof is given as We have proved that the local truncation error of our algorithm (26)- (27) and (28)-(29) is O(h 3 ) when α > 1, so the convergent order is 2.
Lemma 7. (Ford Simpson [6,Theorem 1]) The nested mesh scheme preserves the order of the underlying quadrature rule on which it is based.
Because of Theorem 2, Lemma 7, and the analysis in above section, we have Theorem 3. When α > 1,if∂ 2 f (τ, y(τ))/∂ 2 τ ∈ C[0, t) for some suitable t, then the local truncation error of our algorithm with the predictor and corrector formulae (32)-(33) (α ∈ (1, 2))i sO(h 3 ) and the convergent order is 2, i.e., max j=0,1,··· ,n+1 The comparison of the local truncation error, convergent order and arithmetic complexity for different predictor-corrector schemes are presented in the following table. So far, our convergence results are obtained under the smoothness assumptions of D α * y(t) and f , which obviously depend on the smoothness properties of the solution y(t) (Deng, 2010).

Numerical results
In this section we present two numerical examples to illustrate the performance of our proposed predictor-corrector schemes. We show the order of convergence in the the absolute Schemes Local truncation error Convergent order Arithmetic complexity where error(h) is the absolute error |y(t) − y h (t)| with the step h.
Example 1. Consider the following fractional differential equation with the initial conditions y(0)=0, α ∈ (0, 1), (50) or y(0)=0, y ′ (0)=0, α ∈ (1, 2), (51) note that the exact solution to this problem is The numerical results for the predictor-corrector scheme (16)-(17) at time t = 1 with different steps and different α are reported in Tables 2 and 3. Table 1 shows the numerical errors at time t = 1 between the exact solution and the numerical solution of the predictor-corrector schemes (16)-(17) for (49) with different α ∈ (0, 1) for various step sizes. It is seen that the rate of convergence of the numerical results is of order O(τ min{2,1+2α} ) for the predictor-corrector schemes (16)- (17). Tables 2 and 3 show the ratio of the error reduction as the grid is refined. The last row states the theoretical orders of convergence (abbreviated as TOC), which are the results we theoretically prove in the above theorems. From Tables 2 and 3, it can be noted that the numerical results are in excellent agreement with the theoretical ones given in Theorem 1.
Example 2. Now, consider the following fractional differential equation (Diethelm et al., 2004) with the initial conditions y(0)=0, α ∈ (0, 1), (54) or y(0)=0, y ′ (0)=−1, α ∈ (1, 2), (55) the exact solution to this problem is y(t)=t 2 − t.  Table 3. Absolute errors and convergence orders for predictor-corrector schemes (16)-(17) at time t = 1 with different 1 < α < 2. Table 4 shows the numerical errors at time t = 1 between the exact solution and the predictor-corrector schemes (16)-(17) for (53) with different α ∈ (0, 1) for various step sizes. Again we numerically verify that the rate of convergence for predictor-corrector schemes (16)- (17) is O(τ min{2,1+2α} ). From  Table 5. Error behavior versus the variation of p and T (the definition of p and T are given in (3)) at time t = 50 with exact (analytical) value y(50)=2450, fractional order α = 1.5, step length h = 1/80. and the errors of the improved predictor-corrector scheme (16)- (17) are improved significantly compared with the Tables 3 and 4 given in (Diethelm et al., 2004). Table 5 shows the computing values, the absolute numerical errors, and the relative numerical errors by using the scheme (32) and (33) to compute Example 2 for different values of p and T defined in (3). According to the numerical results we can see that the computing errors are generally acceptable in engineering when the computational cost is greatly minimized, especially the computing error is not sensitive to the value of p. On the other hand, this numerical example also illuminates that the algorithm is numerically stable.

Conclusions
We briefly review the numerical techniques for efficiently solving fractional ordinary differential equations. The possible improvements for the fractional predictor-corrector algorithm are presented. Even though the algorithm is designed for scalar fractional ordinary differential equation, it can be easily extended to the systems. In fact, based on this algorithm, we have simulated the dynamics of fractional systems (Deng, 2007c;d). In addition, the fractional predictor-corrector algorithm combining with the scheme of generating stochastic variables also works well for the stochastic fractional ordinary differential equations (Deng & Barkai, 2009). This book demonstrates applications and case studies performed by experts for professionals and students in the field of technology, engineering, materials, decision making management and other industries in which mathematical modelling plays a role. Each chapter discusses an example and these are ranging from wellknown standards to novelty applications. Models are developed and analysed in details, authors carefully consider the procedure for constructing a mathematical replacement of phenomenon under consideration. For most of the cases this leads to the partial differential equations, for the solution of which numerical methods are necessary to use. The term Model is mainly understood as an ensemble of equations which describe the variables and interrelations of a physical system or process. Developments in computer technology and related software have provided numerous tools of increasing power for specialists in mathematical modelling. One finds a variety of these used to obtain the numerical results of the book.