Abstract and Applied Analysis

Volume 2015 (2015), Article ID 952057, 9 pages

http://dx.doi.org/10.1155/2015/952057

## A New Grünwald-Letnikov Derivative Derived from a Second-Order Scheme

^{1}School of Computer Science and Applied Mathematics, University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa^{2}DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa

Received 25 May 2015; Accepted 16 August 2015

Academic Editor: Tiecheng Xia

Copyright © 2015 B. A. Jacobs. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A novel derivation of a second-order accurate Grünwald-Letnikov-type approximation to the fractional derivative of a function is presented. This scheme is shown to be second-order accurate under certain modifications to account for poor accuracy in approximating the asymptotic behavior near the lower limit of differentiation. Some example functions are chosen and numerical results are presented to illustrate the efficacy of this new method over some other popular choices for discretizing fractional derivatives.

#### 1. Introduction

Fractional differential equations and fractional order derivatives have become increasingly popular in recent years. The use of fractional order derivatives has become abundant in modeling techniques as well as computational methods for the numerical solution of these models. Zeng et al. [1] derive a generalized fractional diffusion equation by addressing deficiencies in previously proposed models for anomalous diffusion. An application of modelling plumes moving at a fractional rate was investigated by Benson et al. [2] which showed a fractional advection-dispersion model to be appropriate as the mean squared displacement followed . He [3] presented a more exact model for seepage flow by relaxing the continuity assumption and went on to solve the arising equation using the variational iteration method (VIM). Metzler et al. [4] provide a review of anomalous diffusion equations presented as well as the formulation of a general diffusion case overcoming previous inconsistencies. Moving past anomalous diffusion Gafiychuk et al. [5, 6] studied a fractional reaction-diffusion system modelling an activator-inhibitor dynamic. The application of fractional order models is various and with this increased complexity requires both analytic and numerical methods to find useful solutions.

Meerschaert and Tadjeran [7] construct an implicit scheme for two-sided space-fractional partial differential equations based on two first-order approximations to the left- and right-shifted Grünwald-Letnikov fractional derivative. Scherer et al. [8] also derive an implicit scheme for the fractional heat equation in the Caputo sense, making use of the standard Grünwald-Letnikov scheme as well as the relationship between the Riemann-Liouville derivative and the Caputo derivative. There has been work conducted on comparing numerical methods with analytical methods; El-Sayed et al. [9], for example, compare a modified numerical method, wherein the authors apply a transform to homogenize the initial conditions, with the Adomian decomposition method. Another approach taken by Karatay et al. [10] derives an unconditionally stable implicit scheme based on a first-order approximation to the Caputo derivative. An interesting study is presented by Su et al. [11] which combines the first-order shifted Grünwald-Letnikov scheme with characteristic methods and a departure from the Riemann-Liouville and Caputo derivatives leads Chen et al. [12] to derive a conditionally stable explicit scheme and unconditionally stable implicit scheme for the Caputo-Riesz fractional diffusion equation that are superlinearly convergent.

Tadjeran and Meerschaert [13] implement a second-order accurate Crank-Nicolson-type finite difference scheme for the space-fractional diffusion equation, making use of the right-shifted Grünwald-Letnikov fractional derivative. Diethelm et al. [14] propose a predictor-corrector-type approach for the numerical solution of both linear and nonlinear fractional differential equations including systems of equations. Mohebbi et al. [15] study a high-order accurate numerical method for obtaining the solution to the subdiffusion equation with a nonlinear source term. Cui [16] goes on to solve the two-dimensional time-fractional diffusion equation using a high-order compact alternating direction implicit method. Meerschaert et al. [17] solve the two-dimensional space-fractional dispersion equation using a right-shifted Grünwald-Letnikov discretization since the order of the fractional derivative is above 1 but less than 2. Lynch et al. [18] implement a finite difference approach for the space-fractional dispersion problem but implement discretization switching to maximize accuracy based on the order of derivative. Momani and Odibat [19] depart from the finite-difference approach and instead compare two analytic methods for the solution of linear fractional differential equations, namely, the variational iteration method and the Adomian decomposition method. These authors go on to explore the application of these methods to nonlinear partial differential equations of fractional order in [20].

The present work derives a difference approximation in a novel way that produces a form similar to the known result presented by Oldham and Spanier in [21]. This text cites an error in the approximation of the fractional derivative of a constant and linear function of order that is proportional to the number of terms in the series approximation. A novel technique is coupled with this approximation to attain order accuracy. Although the result obtained is known, to the best of the author’s knowledge this derivation is novel.

Another important technique used in the numerical application of fractional differential equations is the “short-memory” principle [22]. This technique truncates the approximation sum to avoid high computational costs. It is shown in Section 4 that the present method is highly accurate even with very few terms in the approximation sum.

The following section presents some preliminary results used in the remainder of the paper. Section 3 goes on to derive the new method of discretizing a fractional derivative with results on the performance of this technique reported in Section 4. We make some concluding remarks and discuss opportunities for further work in Section 5.

#### 2. Preliminaries

This section presents some preliminary results that are put to use in subsequent sections. However, in the interest of brevity, we omit some useful properties but direct the interested reader to the work by Li and Deng [23]. Often for numerical simulation of fractional order problems the Grünwald-Letnikov derivative is used as a means of discretizing a fractional order derivative. The Grünwald-Letnikov discretization has been used for numerical schemes for fractional partial differential equations [13, 24, 25] and is defined as follows.

*Definition 1. *The Grünwald-Letnikov derivative of order of a function is where and .

However the limit definition is only practically useful for finite-difference implementations; hence the following definition is used.

*Definition 2. *The Grünwald-Letnikov derivative of order of a function is [22] where .

Definition 2 is used as a means of calculating the exact fractional derivative of a function. The exact solution obtained here is then used as a measure for accuracy when compared with the discrete approximations. It is prudent to mention that if is sufficiently smooth, , then the Grünwald-Letnikov derivative is equivalent to the popular Riemann-Liouville defined below.

*Definition 3. *The Riemann-Liouville derivative of fractional order of a function is [22]

The right-shifted Grünwald-Letnikov approximation introduced by Meerschaert and Tadjeran [25], as used in [13], leads to stable numerical schemes and is defined below.

*Definition 4. *The right-shifted Grünwald-Letnikov approximation of a function for is [25]

Oldham and Spanier obtain a second-order Grünwald-Letnikov approximation simply by noting that the Riemann sum present in the derivative would converge faster if a simple modification was made.

*Definition 5. *The Grünwald-Letnikov derivative of order of a function with a central difference modification is [21] where .

The primary difference between Definition 5 and the present work is the derivation and the choice of , the upper limit of the Riemann sum. Although these discrepancies are slight we also address the issue of the asymptotic behavior of the fractional derivative as approaches zero, a difficult behavior for a difference scheme to capture.

#### 3. Method

In much the same way the Grünwald-Letnikov derivative is generalized from the first-order backward difference scheme, as described by Podlubny in [22], we derive a new approximate fractional derivative in the form of a difference scheme which is derived from the second-order central difference scheme. Yuste and Acedo [24] state that the standard Grünwald-Letnikov derivative is accurate to , while we derive a scheme based on the central-difference scheme; we hope to achieve an approximation that is . Podlubny [26] shows, by example, that the Grünwald-Letnikov derivative of a function or is accurate to and hence conclude that the Grünwald-Letnikov derivative of any function that can be written as a power series will be accurate to . The central difference approximation to the first-order derivative is given by or Similarly the second-order derivative can be approximated by Following this trend we may write Generalizing the order of the derivative to be a positive real value and substituting into (10) we obtain a general formula: where the binomial coefficient is calculated by the use of the Gamma function that generalizes the factorial operator: It is well known [26, 27] that the fractional differentiation operator is a nonlocal operator. This effect is seen in the discretization being affected by every point in the function and not just a neighborhood of the point . This insight guides our choice of as the upper limit of the summation. We choose so that the last term in the series corresponds with the terminating point of the function , and in terms of fractional differential equations and fractional partial differential equations this would correspond with a boundary or initial condition: where represents the ceiling function. Due to the singularity of a fractional derivative when it is convenient to select such that to alleviate the asymptotic behavior. However we do not want to restrict ourselves only to functions that obey this property. Hence we define so that and then since the fractional derivative is a linear operator. We may now exploit Property 2.2 in [23] which yields which then leads to Finally (15) becomes While correcting for the asymptotic behavior at produces sufficiently accurate results, we may go one step further and substitute the exact solution of the linear component of . This is done to ensure that the accuracy of the scheme is at least . Letthen

##### 3.1. Accuracy

As described by Oldham and Spanier in [21] the error in the above scheme for functions , , and is given by , , and , respectively. Since we make the exact substitutions for the linear components we are able to increase the accuracy of the scheme to for any analytic function.

Theorem 6. *The error of , denoted by is at least accurate for the function for all given that the exact derivative of order of is *

*Proof. *Let . Then (see [21], Chapter 8) (while Oldham and Spanier cite the relative error it is easy to show that the more restrictive, absolute error is also . More information can be found in [22], Chapter 7)Assume that the order of error is at least for ; that is, Now let ; then by the product rule for fractional derivatives (see [21], Chapter 5) we have Then using the recursive property of the Gamma function we have

*Using Theorem 6 we can deduce that the accuracy of for any analytic function is provided that the linear and constant components of that function are dealt with accordingly.*

*4. Results*

*This section presents some numerical results for different functions. Each example shows the absolute error over the domain of the function, the consistency of each derivative as tends to zero, the absolute error of the truncated series, and the absolute error over a range of . We compare the aforementioned criteria for four different methods: the standard Grünwald-Letnikov approximation , the second-order Grünwald-Letnikov approximation , the second-order Grünwald-Letnikov approximation with constant correction described in Section 3 , and finally the right-shifted Grünwald-Letnikov approximation .*

*4.1. Example 1*

*Let and let the order of derivative be , with time step . Since this choice of satisfies the approximation does not suffer from poor accuracy near and hence the second-order Grünwald-Letnikov and the second-order Grünwald-Letnikov with constant correction approximations are similar.*

*4.2. Example 2*

*Letthe order of derivative , and . As mentioned by Oldham and Spanier [21] this method is robust for all ; this example shows highly accurate results for .*

*4.3. Example 3*

*Letthe order of derivative , and again . This function suffers from the asymptotic behavior of the fractional derivative near and hence the constant correction method enjoys a vastly superior accuracy over the other methods.*

*4.4. Example 4*

*Letthe order of derivative , and .*

*4.5. Example 5*

*Letthe order of derivative , and .*

*The above results indicate that the present method enjoys a high accuracy. The nature of the fractional derivative near the lower limit of differentiation results in a high gradient. Difference schemes with fixed interval steps are unable to accurately capture very steep gradients, similar to those at the asymptote of the fractional derivative. The use of the constant correction allows one to avoid the asymptotic behavior near and the provable accuracy.*

*Figures 1, 5, 9, 13, and 17 show the increasing accuracy of the approximations as deviates from 0, as well as the robustness of our method to the asymptotic behavior. The consistency of the methods is evident in Figures 2, 6, 10, 14, and 18. These figures indicate that as tends to 0 the absolute error of the approximations also tends to 0. As mentioned in Section 1 the “short-memory” principle allows for the truncation of the approximating series in order to alleviate computational cost; our present method is highly accurate even for small values of as is evident in Figures 3, 7, 11, 15, and 19. The final set of figures, Figures 4, 8, 12, 16, and 20, indicate the robustness of our method to varying values of . It is important to note here the log scale of the graph.*