On the order reduction of approximations of fractional derivatives: an explanation and a cure

Finite-difference based approaches are common for approximating the Caputo fractional derivative. Often, these methods lead to a reduction of the convergence rate that depends on the fractional order. In this note, we approximate the expressions in the fractional derivative components using a separate quadrature rule for the integral and a separate discretization of the derivative in the integrand. By this approach, the error terms from the Newton–Cotes quadrature and the differentiation are isolated and it is possible to conclude that the order dependent error is inevitable when the end points of the interval are included in the quadrature. Furthermore, we show experimentally that the theoretical findings carries over to quadrature rules without the end points included. Finally we show how to increase accuracy for smooth functions, and compensate for the order dependent loss.


Introduction
In comparison to normal integer derivatives, fractional derivatives include a memory effect, which is advantageous when modeling real materials in mechanics and electronics, as well as in the description of rheological properties of solids [26]. Other applications where this memory effect is important include modeling of speech [6], bioengineering [20], fluid dynamics [16] and finance [29]. There are several formulations of fractional derivatives, and the most common ones are Riemann-Lioville's, Grünwald-Letnikov's and Caputo's formulation [25].
Fractional differential equations (FDEs) extend ordinary differential equations by introducing or replacing the integer order derivatives with fractional derivatives. The Caputo derivative is suitable when modeling real physical problems, since it requires integer (and not fractional) order of the derivative as initial condition in FDEs. Together with a closed quadrature rule (including both end-points), the finite-difference approximation is a common approach for evaluating the Caputo fractional derivative.
High-order methods for the Caputo fractional derivative have been developed in [8,19] where the order is 3−α and in [22] of order 4−α. Another approach for the Caputo fractional derivative based on quadratic interpolation and compact operators, which obtain convergence rates of 3 − α, is presented in [12]. In [10], a similar procedure for the diffusion-wave equation (1 < α < 2) obtained a convergence rate of 4 − α.
In [24], Odibat presents a modification of the Trapezoidal rule for approximating the fractional integral and Caputo derivative to second order. A convergent scheme for the diffusion-wave equation where the truncation error depends on the fraction of the derivative (usually denoted by α) was derived in [31]. Similar error estimates are found in [18,23]. A comprehensive review of finite difference methods and the relation between accuracy, stability and treatment of boundary conditions is provided in [32], where order reduction in numerical approximations is observed at the boundaries.
Fractional derivatives can be used to form FDEs but also stand alone, for example when used in a P I λ D μ controller [26]. FDEs are a special type of Volterra equations with weakly singular kernels [13]. To treat the singularity, different strategies can be employed such as graded meshes or transformations [7,14,27] Other numerical methods include non-polynomial spline methods [17], spectral methods [11,33], piecewise constant approximation or integrabilization methods [3,4] and the discrete time random walk approach [1,2,5]. Yet, another strategy, which we use in this article, is integration by parts. This bypasses the singularity already in the continuous setting, but requires additional smoothness of the integrand.
High-order methods are beneficial when accurate solutions to partial differential equations is the target [15]. They heavily rely on stable and accurate weak boundary and interface conditions, which summation-by-parts operators (SBP) together with the Simultaneous-Approximation-Term technique (SAT) provide for finite differences [9]. We will in this work illustrate that the approximation of the individual components of the fractional derivative (integral, convolution kernel and integrand) by SBP operators yields a concise numerical approximation of the fractional derivative. This enables results from the extensive SBP literature [9] to be employed to fractional order problems.
In this note, we examine a methodology for numerically computing the Caputo fractional derivative of smooth functions with non-singular derivatives by approximating the expression component-wise using a quadrature for the integral and an appropriate discretization of the integrand using SBP operators. Of particular interest is the limitations of the convergence rates for closed quadratures and the influence of the parameter α as observed in the references above. We do not intend to develop a new numerical method, but rather aim to edify the source of the reduced accuracy obtained from traditional quadrature methods. Moreover, we illustrate the dependence of the convergence rate on the position within the domain, most appreciably so when the quadrature domain approaches the singularity in the integration kernel. Finally we show how to increase accuracy for smooth functions, and compensate for the order dependent loss. The insights provided herein may lead to the development of robust and more accurate quadrature methods for integrals of this nature.

Fractional calculus preliminaries
All preliminaries presented here can be found in [26]. The fractional integral is defined as is the Gamma function. The Caputo derivative is given by where the first equality is introduced to ease the notation for the upcoming analysis.
In this work we only consider the case of m = 1. This is justified by the possibility of decomposing a higher order fractional derivative into a system of first and αth (0 < α < 1) order differential equations.

The reduction of accuracy in quadrature
This section illustrates the common reduction of order in quadrature methods when approximating the convolution integral of the fractional derivative. Considering the domain τ ∈ [0, t], with uniform discretisation τ = {0, τ, . . . , t − τ, t}, we show that O( τ 1−α ) terms cannot be avoided or cancelled by any Newton-Cotes quadrature. This claim is proved in Sect. 4.
We begin by considering the convolution integral We parameterise the integration limits and consider an arbitrary 'sub' integral subject to the general integration limits τ ∈ [a, b], such that b − a = τ . This enables the change in behaviour as b → t to be quantified.
. We now replace the integrand g(τ ) with the Taylor series of this function expanded around the midpoint of the integral domain,τ = a+b 2 , yielding This expansion is a valid use of the Taylor series since we never explicitly evaluate g(τ ) (or any derivatives of g(τ )) at τ = t where we have a singularity. We now compute the integral in noting that the integral of odd derivative terms is zero, due to the expansion around the midpoint of the domain. We may construct a quadrature (trapezoidal) to cancel the leading term in Eq. (3.4) using noting that the error, E, is then given by, Now one might be inclined to conclude that E = O( τ 3 ) in all subintervals. However since, will include increasing derivatives of the kernel K (t − τ ), due to repeated application of the product rule. Note that, Then the first term in Eq. (3.7) is then given by This process can be repeated for every error term in Eq. (3.7) giving the general form, The form in (3.11) illustrates that there are two sources of truncation error, and the dominance of one over the other is predicated by the position in the integral, or the value of b. Examining (3.11), when b = τ (corresponding to the first sub integral), or when the quadrature is sufficiently far from the singularity the order, (3.11), is dominated by the τ n+1 since provided τ t. This phenomenon is also observed numerically where, away from the singularity, our quadrature recovers O( τ 3 ) per subinterval, and when N of these errors are summed, we observe error like O( τ 2 ). When b = t, or when b is sufficiently close to the singularity, then we have relation (3.11) By observing the convergence behaviour at the lower and upper integration terminals we are able to establish bounds on the order of convergence of the total quadrature scheme. However, we are also interested in how the order of convergence changes across the domain. To assess this, we visualise the convergence in (3.7) as a function of b. We note that the expansions used in (3.5) and (3.6) are centered in the midpoint of each domain, so as to avoid evaluating the integrand at the point of singularity. Figure 1 illustrates the convergence order of each subinterval, τ ∈ [b − τ, b], as a function of b. As b → 0 the integrand exhibits no singular behaviour and recovers the expected convergence rate. However, as b → t the degeneration of convergence rate from O( τ 3 ) to O( τ 1−α ) occurs gradually as the quadrature approaches the singularity at τ = t. More importantly, the leading error for any Newton-Cotes quadrature at b = t is given by (3.13) and is bound by τ 1−α regardless of the value of n. The natural approach in improving the accuracy of a quadrature is to increase the number of quadrature nodes to cancel increasing error terms. However, as was shown above, this procedure does not reduce the order of leading error term, regardless of how many quadrature nodes are used.
While Fig. 1 illustrates the order of the truncation error of each subdomain, τ ∈ [b − τ, b], parameterised by b, we need the accumulated truncation error over all subdomains to quantify the truncation error of the total quadrature.
Considering the sum over all sub-intervals we have, where τ = t/N , and p indicates the error behaviour for the pth term in the Taylor series. The function E tot is sampled with decreasing values of τ and various values of α, shown in Fig. 2. From this figure we can see that the dynamics of the total quadrature follow O( τ 1−α ), an inevitable limitation using this approach. This poor accuracy has been previously observed and can be improved upon using several methodologies [8,10,22], however these approaches are not directly related to the Newton-Cotes type quadrature present in this work. Subsequent sections in this work propose a methodology for improving the accuracy.

Generalisations
This section generalises the results shown in the previous section to an arbitrary order Newton-Cotes quadrature as well as the expected convergence ceiling imposed by a generalised weakly singular kernel in the form Typically Newton-Cotes quadratures with n nodes in a given stencil will exhibit a convergence rate of O( τ n ) for non-singular integrands.

Proposition 1 A Newton-Cotes quadrature of order n for integrals containing a weakly singular kernel of the form
have a convergence rate of Proof Consider a discretisation of the domain [a, b] with equally spaced nodes, so that {τ 0 , τ 1 , . . . , τ n } = {a, a + τ, . . . , b}. Given g(τ ) = K (t − τ ) f (τ ) and provided f (τ ) is sufficiently smooth, then we have an nth order Newton-Cotes quadrature defined as where ω i is obtained from computing (3.5) and (3.6) on each subinterval and collecting the equally sized coefficients. The values for ω i 's are computed to match and cancel the first n terms of (3.4) Hence, this quadrature has error The leading error term in (4.2) contains the nth derivative of g(τ ), g (n) (τ ). Moreover, generalising (3.9) using the generalised Leibniz Rule (product rule of high order derivatives) we have, with The term in the above summation limits the convergence rate for j = n, is K (n) (t − τ ) f (τ ). We can now more precisely define (4.2) as The leading cause of error occurs in the final term of the sum, i.e. for j = n, where τ = t − τ 2 is sufficiently close to the singular boundary and the nth derivative of the convolution kernel is given by (4.6)

Remark 4.1
In summary, we highlight that the error result obtained in (4.5) is independent of n, indicating that the source of error inherent in this class of quadrature for weakly singular kernels is independent of the order of the specific quadrature used.

Raising the accuracy to compensate for the previously described loss
To obtain a numerical approximation of (2. is used, where ω i > 0 are weights, then the singularity of the integrand in (2.3) at τ = t will lead to a blow-up.
To cure this anomaly, (2.2) is rewritten by using integration by parts which leads to the modified formula and removes the singularity. Approximating (5.2) numerically yields In (5.3), D 1 f ≈ f and D 2 f ≈ f are finite-difference approximations of f and f evaluated on the grid. The matrix D 1 approximates the first derivative [30], while D 2 generates an approximation of the second derivative [21]. Furthermore, T = diag(t, t − h, . . . , 0) and the exponent in T 1−α should be interpreted element-wise.
To verify the accuracy of (5.3), we compare the numerical result to a known analytical value. For polynomials of degree n, an analytic expression can be obtained by repeatedly using integration-by-parts: D α t n = n!t n−α / (n + 1 − α). Given two step-sizes, h 1 , h 2 , the error in , and the convergence rate is r h 1 = log(e h 1 /e h 2 )/ log(h 1 /h 2 ). To start, we let P 2 contain the weights of the second order Trapezoidal rule and D 1,21 , D 2,21 be the matrices containing the second order accurate central finite-difference schemes in the interior and appropriate first order forward and backward stencils at the boundaries. The error and The Trapezoidal rule approximates the integral and the second order finite difference operators, D 1,21 , D 2,21 , approximate the derivatives A quadrature rule of order seven is used and differentiation matrices D 1,84 , D 2,84 of order eight in the interior and four at the boundaries convergence rate is presented in Table 1 for f (t) = t 3 . We find that the convergence rate approaches 2 − α, thus we have raised the order of the accuracy by one.
In an ambition to further increase the convergence rate, we use a quadrature matrix P and differentiation matrix D of higher order. Table 2 presents the result for SBP84, where P 7 now is a diagonal and positive definite matrix forming a quadrature of order seven and D 1,84 , D 2,84 are eighth order in the interior and fourth order at the boundaries. These operators are sufficiently accurate such that numerical integration and differentiation errors are negligible [21,30,32]. However, as for the lower order case, the numerical experiments again suggest a convergence rate of 2 − α, with no improvement in accuracy. Similar results are obtained for other smooth functions f in the integrand.
The results in Tables 1 and 2 and the previous analysis in Sect. 4 illustrate that there are two factors limiting the accuracy: (1) the natural error source that is related to the fractional kernel and present for all closed quadratures and (2) the order of the specific quadrature. The most significant error source determines the convergence rate.
To verify our claim, we start by integrating (2.2) by parts to get (5.4) The quadrature matrix P 7 and the differentiation matrices D 1,84 , D 2,84 have been used in the computations In (5.4), D 3 f = D 2 D 1 f ≈ f and the quadrature matrices are the same as in Sect. 5. For the fractional kernel (t −τ ) 1−α , the convergence order of the numerical approximation (5.3) was 2 − α. Raising the power of the kernel to 2 − α as in (5.4), we expect a high-order version of the numerical approximation in (5.5) to generate a convergence of order 3 − α. The error and convergence rates for the scheme in (5.5) together with P 7 and D 1,84 , D 2,84 and f (t) = t 3 are presented in Table 3. As expected, the results suggest a convergence rate of 3 − α, meaning that (1) is the limiting factor. If P or D are of lower convergence order, then (2) will dominate instead, as seen in Table 4. In a similar manner, repeatedly using integrating-by-parts results in which can be used to increase the order. Discretizing (5.6) leads to  where D i f ≈ f (i) is the i-th derivative evaluated on the grid. The obvious downsides of (5.6) and (5.7) are that the function must be smooth enough with non-singular derivatives and secondly, the necessity of an high-order approximation of a high derivative.

The extension to open quadratures
The analysis above carries over to open quadratures as well, which we will illustrate numerically by considering the composite 2nd order midpoint rule and the composite 4th order Milne's rule [28], denoted by P O 2 and P O 4 , respectively. For the approximation t 0 (t −τ ) 2−α dτ ≈ 1 PT 2−α 1, we expect that P O 2 leads to 2nd order convergence and P O 4 leads to 3 − α. Tables 5 and 6 presents the results, which agree with our theoretical speculation.

Conclusions
We have investigated the numerical approximations of Caputo's formulation for fractional derivatives for closed quadrature rules. We showed that the convergence rate for such discretization depend naturally on the fractional parameter α. Two different error sources determine the convergence rate: (1) the inherent error associated with the numerical integration of the kernel and (2) the order accuracy of the numerical method used. The most significant one determines the convergence rate. Moreover, we have shown in the case of a general Newton-Cotes quadrature that the leading term in the error expression is limited not by the order of the quadrature but by the degree of the singularity present in the integrand. We have also shown how to increase the order of accuracy by using integration by parts. The present work is general for smooth integrands. In future work we aim to extend the analysis to non-smooth integrands.