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.
5. Conclusion
This work has shown a novel derivation of a second-order accurate scheme for approximating the fractional order derivative of a function. To the best of the author’s knowledge performing a linear transformation on the function to avoid asymptotic behavior near the origin is also novel. Section 4 provides several examples with varying parameters and behaviors to illustrate the efficacy of the present method.
Using a higher-order accurate approximation leads to the potential for higher-order accurate numerical schemes for fractional differential equations and fractional partial differential equations. Tadjeran and Meerschaert [13] implement the right-shifted Grünwald-Letnikov derivative for fractional derivatives of order which lead to an unconditionally stable numerical scheme. According to Oldham and Spanier [21] the present method is robust for all values of alpha. Potential extensions to this work include constructing numerical schemes centered on this approximation and the fractional order derivative for nonlinear fractional differential equations as well as time-marching schemes for time or space fractional partial differential equations.
The accuracy of the method has been proven to be at least of , which suggests that the accuracy should be of order for . However in these examples the error is observed to be in the region of . Nevertheless the mathematical proof of accuracy as well as the compelling results presented substantiates the efforts toward this field.
Disclaimer
The opinions expressed and conclusions arrived at herein are those of the author and are not necessarily to be attributed to the CoE-MaSS.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The author would like to acknowledge and thank Professor Charis Harley for her helpful suggestions and useful insight. The DST-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS) is acknowledged for funding under Grant no. 94005. Thanks are due to the reviewers for their comments which helped to elevate the level of this work.