A New Second Order Approximation for Fractional Derivatives with Applications

We propose a generalized theory to construct higher order Grünwald type approximations for fractional derivatives. We use this generalization to simplify the proofs of orders for existing approximation forms for the fractional derivative. We also construct a set of higher order Grünwald type approximations for fractional derivatives in terms of a general real sequence and its generating function. From this, a second order approximation with shift is shown to be useful in approximating steady state problems and time dependent fractional diffusion problems. Stability and convergence for a Crank-Nicolson type scheme for this second order approximation are analyzed and are supported by numerical results.


Introduction
ractional calculus has a history that goes back to L'Hospital, Leibniz and Euler [1,2].A historical account of early works on fractional calculus can be found, for example, in [3].Fractional integral and fractional derivative are extensions of the integer order integrals and derivatives to a real or complex order.Various definitions of fractional derivatives have been proposed in the past, among which the Riemann-Liouville, Grünwald-Letnikov and Caputo derivative are common and established.Each definition characterizes certain properties of the integer order derivatives.Recently, fractional calculus found its way into the application domain in science and engineering.The field of application includes, but is not limited to, oscillation phenomena [4], visco-elasticity [5], control theory [6] and transport problems [7].Fractional derivatives are also found to be suitable to describe anomalous transport in an external field derived from the continuous time random walk [8], resulting in a fractional diffusion equation.The fractional diffusion equation involves fractional derivative either in time, in space or in both variables.
Fractional derivative is approximated by the Grünwald approximation obtained from the equivalent Grünwald-Letnikov formula of the Riemann-Liouville fractional derivative.Numerical experience and theoretical justifications have shown that application of this approximation as it is in the space fractional diffusion equation results in unstable solutions when explicit, implicit and even when the Crank-Nicolson (CN) type schemes are used [9].

F
The latter two schemes are popular for their unconditional stability for classical diffusion equations.This peculiar phenomenon for the implicit and CN type schemes is corrected and the stability is restored when a shifted form of the Grünwald approximation is used [9,10].
The Grünwald approximation is known to be of first order with the space discretization size h in the shifted and non-shifted form and is, therefore, useful only in first order schemes such as explicit Euler (forward) and implicit Euler (backward) for the fractional diffusion equation.
Since the CN approximation scheme is of second order in time step  , Tadjeran et al. [11] used extrapolation improvement for the space discretization to obtain a second order accuracy.Subsequently, second order approximations for the space fractional derivatives were obtained through some manipulations on the Grünwald approximation.Nasir et al. [12] obtained a second order accuracy through a non-integer shift in the Grünwald approximation, displaying super convergence.Convex combinations of various shifts of the shifted Grünwald approximation were used to obtain higher order approximations in Chinese schools [13,14,15,16,17], some of which are unconditionally stable for the space fractional diffusion equation with CN type schemes.Zhao and Deng [18] extended the concept of super convergence to derive a series of higher order approximations.Earlier, Lubich [19] obtained some higher order approximations for the fractional derivative without shift for orders up to 6. Numerical experiments show that these approximations are also unstable for the fractional diffusion equation when the methods mentioned above are used.Shifted forms of these higher order approximations diminish the order to one, making them unusable as Chen and Deng [20,21] noted.
In this paper, we construct a new second order Grünwald type approximation which can be used with shifts without reducing its order.We apply this second order approximation in steady state problems and time dependent fractional diffusion problems.A CN type scheme is devised for this approximation and a justification for stability and convergence is given.
The rest of the paper is organized as follows.In Section 2, definitions and notations are introduced.In Section 3, the main results of generalization are presented.In Section 4, approximations of order one and two with shifts are constructed.In Section 5 the constructed second order approximation is applied to devise numerical schemes for steady state and time dependant space fractional diffusion equations.Stability and convergence for the scheme of fractional diffusion equations are analysed in Section 6. Supporting numerical results are presented in Section 7 and conclusions are drawn in Section 8.

Denote by
, the space of Lebesgue integrable functions: and assume that it is sufficiently differentiable so that the following definitions hold.The left and right Riemann-Liouville fractional derivatives of order R   are defined respectively as , ) ( where The corresponding left and right Grünwald-Letnikov (GL) fractional derivatives are given respectively by ), ( 1) It is known that the Riemann-Liouville and Grünwald-Letnikov definitions are equivalent [22].Hence, from now onwards, we denote the left and right fractional derivatives as The Fourier transform (FT) and inverse FT of an integrable function , the FTs of the left and right fractional derivatives are given by For a fixed h , the Grünwald approximations for the fractional derivatives in (1) and ( 2) are obtained by simply dropping the limit in the GL definitions (3) and ( 4) as  respectively, where [ b a , it is zero-extended outside the interval to adopt the definitions of fractional derivatives and their approximations.The sums are restricted to a finite up to N which grows to infinity as 0  h .Often, N is chosen to be The Grünwald approximations are of first order accuracy and display unstable solutions in the approximation of fractional diffusion equation by implicit and CN type schemes [9].As a remedy, a shifted form of Grünwald formula with shift where the upper limits of the summations have been adjusted to cover the shift r .
Meerchaert et al. [9] showed that for a shift 1 = r , the shifted approximations are also first order approximations with unconditional stability restored in implicit Euler and CN type schemes for space fractional diffusion equations.
As for higher order approximations, Nasir et al. [12] derived a second order approximation with a non-integer shift, Tian et al. [13] used convex combinations of different shifted Grünwald forms to obtain two second order approximations: Hao et al. [14] obtained an quasi-compact order 4 approximation.Higher order approximations were also obtained earlier by Lubich [19], establishing a connection with the characteristic polynomials of multistep methods for ordinary differential equations.Specifically, if ) ( ), ( z z   are the characteristic polynomial [23] of a multistep method of order p , then

All the above approximations are based on the shifted Grünwald approximations
gives the weights for the Grünwald type approximation of the same order for the fractional derivative of order  .From the backward multistep methods, Lubich [19] derived higher order approximations of up to order six in the form

Generalization of Grünwald approximation
In this section, we generalize the concept of shifted Grunwald approximation to an arbitrary sequence of weights and analyse its properties.This generalization is then used to construct higher order approximations for fractional derivatives.
For a sufficiently smooth function ) (x f , denote the left Grünwald type operator with shift r and weights We consider the generating function 2. The property (10) is the consistency condition for the Grünwald type approximation operator  r h   , to the fractional differential operator.This consistency condition tells that the order of approximation is at least one.
3. For convenience, we do not distinguish the sequence of weights r k w , and its generating function ) (z W in the above definitions, and in the discussions and results that follow.
The following theorem establishes an equivalent characterization of the generator ) (z W for an approximation of fractional differential operator with order 1  p and shift r . Moreover, if 9), with the help of linearity, we obtain where we have used Now, (11) The equivalent condition (11) for order p approximation in Theorem 1 holds for the right fractional derivative as well.In fact, the condition on the generating function which is equivalent to (11).
One of the consequences of Theorem 1 is the following consistency condition.Using Theorem 1, one can check algebraically the following propositions for the generating function of approximation operators known previously in [12,13,19,20].Proposition 1 The approximation given by ( 7) is of second order.with shifts p and q respectively, check the Taylor series of

Proof. The coefficients
reducing the orders to 1. □

Construction of higher order approximation
We construct higher order approximation generating functions for fractional derivatives with shifts by the use of Theorem 1.The importance of Theorem 1 is that the construction process is entirely confined to algebraic manipulation with the aid of Taylor series expansion.
In this paper, we choose the form First and second order shifted Grünwald approximations have been constructed by the above algebraic method and we have the following.

Theorem 2
The first and second order shifted Grünwald approximations with shift r are given respectively by The generating function 5) and ( 6) for the left and right fractional derivatives respectively, and is independent of the shift r .) ( 2, z W r gives the second order Grünwald type approximations given by (9)   The Grünwald weights r k w , can be computed by a recurrence formula [24,25].

Applications to fractional differential equations
In this section, we apply the fourth order Grünwald type approximation derived in the previous section to fractional differential equations.

Steady state problems
In order to keep the discrete values where i u ˆ are the solutions of (19) and hence the approximation of the exact solution . Then, the matrix formulation of ( 19) is given by , = 2,1  A reduced at both ends as above.

Approximation of fractional diffusion equation
We consider the numerical approximation of the space fractional diffusion equation defined in the domain (20) with the initial and boundary conditions  .We use the following notations for conciseness: . We present the CN type scheme with the order 2 approximation in space using 2,1  .


, where Thus, the Crank-Nicolson type scheme in matrix form reads 1, ,0 where Let  B ˆ be the reduced matrix from  B and 1/2 ˆ m F be the reduced vector from 1/2  m F as was in Section 5.1.
After imposing the boundary conditions, equation (23) reduces to the ready-to-solve form I is the unit matrix of appropriate size.

Stability and convergence
We establish the stability of the Crank-Nicolson type scheme (23).We closely follow the analysis in [14] and present some required results.and the computed convergence orders are listed in Table 2.This test confirms the theoretical justification of second order for the approximation ) ( 2,1 z W .We consider the following test example for the fractional diffusion problem (20).

Proof. For any
and right fractional derivatives respectively, to cover the sum up to the boundary of these domain intervals, where ] [ y is the integer part of y .The left and right fractional derivatives, in this case, are denoted by r k w , which will play a central role in constructing approximations.Remark 1 1.Analogous definitions hold for the right fractional derivative with the same weights r k w , .We denote the right Grünwald

Corollary 1
If the generating function ) (z W gives a consistent approximation of the left and right fractional differential operators, then By Remark 1 and 2, the order p is at least one.When 0  z , the condition (11) becomes

Proposition 2
The approximation given by (8) is of second order accuracy.
an order p approximation, we set the first p coefficients ) (r a l to satisfy(11) in Theorem 1.That is, we impose conditions ( ), = pr u x p   for left and right fractional derivatives respectively.
power of a polynomial of degree 2 , the upper limit of the summation in(17) goes up to only

.
conditions, we obtain the ready-to-solve second order scheme as

F
U , respectively by deleting their first and last boundary entries.and last column vectors of the matrix 2,1 T into a uniform partition of size M with subintervals of length M T/ =  .These two partitions form a uniform partition of the 2

B
reduced again as before, and

.
be the space of grid functions in the computational domain in the space interval ] the discrete inner product and its corresponding norm respectively as .This matrix of size 1 is negative definite.For any positive integer N and for any non-


By virtue of Lemma 2, it is enough to prove that the coefficients

and 2 Lemma 4 2 , 1 
are all non-positive for 6  m .From Corollary 1, we have the consistency condition 0 inequality follows from the consistency condition.□ The approximation operator ' is negative definite.

Theorem 4 .. 5 .
From the above estimates, we have the following stability result.□The CN type difference scheme(23) of order 2 is unconditionally stable for For the convergence of the approximate solution from the CN type scheme, we have the following.Theorem The approximate solutions of the CN type scheme(23) with the given initial and boundary conditions are convergent for We test the approximation scheme devised in Subsection 5.1 using the steady state test problem

Table 1
are of order p accurate without shift.Moreover, if a non-zero shift r is introduced to the approximation, the orders reduce to one.

Table 1 .
with coefficientsWe denote the Grünwald approximation of order p by ,

Table 2 .
Second order approximation with shift