Numerical Approaches to Fractional Integrals and Derivatives: A Review

Fractional calculus, albeit a synonym of fractional integrals and derivatives which have two main characteristics—singularity and nonlocality—has attracted increasing interest due to its potential applications in the real world. This mathematical concept reveals underlying principles that govern the behavior of nature. The present paper focuses on numerical approximations to fractional integrals and derivatives. Almost all the results in this respect are included. Existing results, along with some remarks are summarized for the applied scientists and engineering community of fractional calculus.


Historical Review
The primary attempt, which was recorded in history to discuss the idea of generalizing the integer-order differentiation d n f (t) dt n to d α f (t) dt α with non-integer α, was contained in the correspondence of Leibniz [1]. Some remarks were made on the possibility of considering differentials and derivatives of order one-half. Then, the formulation for derivative of non-integer orders was considered by Euler [2] and Fourier [3]. At the end of the 19-th century, the theory of more-or-less complete form appeared for fractional calculus, primarily due to Liouville [4][5][6][7][8][9][10][11], Riemann [12], Grünwald [13], Letnikov [14][15][16][17], and Marchaud [18,19].
Theoretical analysis of fractional calculus has been booming since the 20-th century. Results in this respect are fruitful, for example, in mapping properties of fractional integration and integro-differentiation [20], Leibniz rule for the generalized differentiation [21], formulae for fractional integration by parts [22], and the Bernstein-type inequality for fractional integration and differentiation operators [23,24], et al.
It is believed that the proper history of fractional calculus began in the realm of physics, with the papers by Abel [25,26]. In those two papers the integral equation was solved in connection with the tautochrone problem. Although Abel did not intend to generalize differentiation, the left-hand side of the integral equation leads to the fractional integral operator of order 1 − µ. Fractional integro-differentiation in such a form was sharpened somewhat later. It was not until the recent few decades that scholars came to realize the importance of fractional calculus for applied sciences, such as rheology, continuum mechanics, porous media, thermodynamics, the left-and right-sided Riemann-Liouville derivatives [32], the left-and right-sided Caputo derivatives [32], and Riesz derivative [33] RZ D α x f (x) = Ψ α RL D α a,x f (x) + RL D α x,b f (x) , 0 < α = 1, 3, 5, . . . , (8) where Ψ α = − 1 2 cos( απ 2 ) . They are the subjects of this paper. Fractional integrals and derivatives of other kinds such as ones in [34] and the very newly defined ones in [35,36] and their approximations are omitted here.
Fractional calculus which has two main characteristics-singularity and nonlocality from its origin, is a generalization of the classical one to some extent. However, these two concepts are different. First of all, fractional integral and Riemann-Liouville derivatives coincide with the integer-order ones while Caputo derivative and Riesz derivative fail to be consistent with integer-order derivatives in general cases. Besides, semigroup property is valid for fractional integral while is invalid for the case with fractional derivatives. Fractional derivatives of periodic functions are not in the same form of those in the integer-order case, either. For example, the α-th order Riemann-Liouville derivative of sin x and cos x with α > 0 are not sin(x + απ 2 ) and cos(x + απ 2 ) unless one adopts the new axiom system proposed in Ref. [37], which differs from the commonly used one. Last but not least, definite conditions for fractional differential equations are in general different from the integer-order case. Especially, in the case with fractional derivatives, boundary and/or initial conditions usually contain fractional derivatives/integrals at the terminals or integer-order derivatives/integrals at points close to the terminals [38,39]. The behavior of the solutions to fractional differential equations may also differ from that of the solutions to the general class of difference equations presented in [40].
In light of potential applications of fractional integration and differentiation operators, there is a substantial demand for efficient algorithms for their numerical handling. Discretizing fractional integrals and derivatives gives a series of quadrature formulae. Different choices of nodes and coefficients give distinct accuracies. Numerical approximations to fractional integrals and derivatives are mainly derived from three distinct paths. Based on the polynomial interpolation, numerical schemes can be obtained with accuracy generally depending on the order of integration and differentiation, for example, the L1, L2, and L2C methods. Convolution quadratures, which preserve properties of fractional integrals and derivatives can be viewed as numerical evaluations of fractional integrals and derivatives with integer-order accuracy. For instance, the fractional multistep methods, among which the fractional backward difference formulae are mostly used, are of integer-order accuracy independent of the order of integration and differentiation. These methods can be verified through Fourier analysis. So do the Grünwald-Letnikov type approximations and fractional centered difference methods. Reformulating fractional integrals and derivatives as infinite integrals of solutions to integer-order ordinary differential equations, the diffusive approximation for fractional calculus can be obtained. Those numerical approximations are discussed in the coming sections. Without specific clarification, the introduced methods are in the setting of uniformed mesh with h = b−a N , N ∈ Z + , and x j = a + jh, j = 0, 1, 2, . . . , N.

Numerical Approximations to Fractional Integral
The weak singularities in Equations (2) and (3) often make it difficult to calculate fractional integrals directly. In the following, several kinds of numerical methods are introduced.

Numerical Methods Based on Polynomial Interpolation
Assume that f (x) is suitably smooth on [a, b]. Then the α-th order fractional integral of f (x) at x = x j with 1 ≤ j ≤ N can be expressed as It is reasonable to utilize an interpolate function f (x) to approximate f (x) on each subinterval, such that the integral x k+1 x k (x j − t) α−1 f (t)dt can be calculated exactly. This idea yields a series of numerical formulae in the form where ω k (k = 0, 1, . . . , j) are the corresponding coefficients. To better understand this method, we retrospect some specific formulae with their brief derivations. If f (x) ∈ C[a, b] on the right-hand side of Equation (9) is approximated by a piecewise constant function then there holds This yields the left fractional rectangular formula [38] where the convolution coefficients b k (0 ≤ k ≤ j − 1) are given by Similarly, if the function f (x) on the right-hand side of Equation (9) is replaced by then we have the right fractional rectangular formula [38] with b k (0 ≤ k ≤ j − 1) given by Equation (14). Based on the left and right rectangular formulae, the weighted fractional rectangular formula [38] yields or the similar form

Remark 1. (I)
The left fractional rectangular formula (13) and the right fractional rectangular formula (16) will be recovered if θ = 1 and θ = 0, respectively. In addition, the weighted rectangular formula (17) (or (18)) is reduced to the composite trapezoidal formula (or midpoint formula) [41] for the classical integral provided that α = 1 and θ = 1 2 . (II) Leading terms of remainders for left-and right-rectangular formulae generally can not be canceled out by introducing weights as the remainders depend on respectively. Therefore, the accuracy of fractional rectangular formulae are of first order accuracy for all 0 ≤ θ ≤ 1. And all the above fractional rectangular formulae are of the first order accuracy. (9) with the piecewise linear polynomial we obtain the fractional trapezoidal formula [38] Here the coefficients a k,j (0 ≤ k ≤ j) are given by Suppose that f (x) ∈ C[a, b]. For 0 ≤ k ≤ j − 1, let {l k,i (x)} be Lagrangian functions defined on the grid points {x k+s , s ∈ S} with S = {0, 1 2 , 1}, which are given by Denote x k+ 1 2 = x k +x k+1 2 and utilize the piecewise quadratic polynomial Then we obtain the following fractional Simpson's formula [38] in which and ]. Let f (x) be approximated by the following r-th degree polynomial on the grid points with Then we obtain the fractional Newton-Cotes formula [38] with the coefficients being calculated by i ) in the integrand. Note that formulae (13), (16), (20), and (24) are special cases of Equation (29). Therefore, they are of the first, second, and third-order accuracy, respectively.
Consider the function f (x) in the Sobolev space H r [a, b] with r ≥ 1 being an integer. Generalizing the above approaches, we can derive spectral approximations [42]. For f (x) defined on [−1, 1], we consider the interpolation function based on the Jacobi with {x k } N k=0 and {ω k } N k=0 being the collocation points and the corresponding quadrature weights [43]. The constants δ u,v j are given by with γ u,v j being defined by In this case, we have the following spectral approximation based on Jacobi polynomials Here P u,v,α dt can be explicitly calculated by the recurrence formula which follows from the three-term recurrence relation of the Jacobi polynomials. Here the coefficients are given by and , and Consequently, numerical scheme (35) becomes the spectral approximation based on Legendre polynomials Here the coefficients are given by withγ j = 2 2j+1 for 0 ≤ j ≤ N − 1,γ N = 2 N , and {ω k } N k=0 being the corresponding quadrature weights.
Correspondingly, for the case with arbitrary interval [a, b], we also have Let u = v = − 1 2 in Equation (39). We obtain the spectral approximation based on Chebyshev polynomials which follows from the relation P The coefficients t j (0 ≤ j ≤ N) are determined by with and {ω k } N k=0 being the corresponding quadrature weights. The above spectral approximations can be rewritten in matrix forms. For differential matrices for fractional integrals and derivatives, see Refs. [42,44] for more details. Here we present numerical examples given by Ref. [42] to verify the spectral accuracy of spectral approximations. 1]. Apply scheme (39) to evaluating its fractional integral. Table 1 shows the absolute maximum errors at the Jacobi-Gauss-Lobatto points. The spectral accuracy is visible.   (39) to evaluate its fractional integral. Table 2 displays the absolute maximum errors of the spectral approximations to the fractional integral with u = v = 0 and u = v = − 1 2 . We can observe that satisfactory results are obtained as well.
On the basis of Dahlquist's theorem on linear multistep methods [46], the proposed convolution quadrature was proved to be convergent of order if and only if it is stable and consistent of order . An easy way of obtaining such a convolution quadrature is by using an -th order linear multistep method to the power α, which gives fractional linear multistep methods. The widely used one is fractional backward difference formula (the fractional BDF), whose implementations are as follows.
Theorem 1 ([45,47,48]). The convolution quadrature (48) approximates the fractional integral RL where s is a fixed integer with s ≤ n. Here the convolution coefficients ω ,j are respectively those of the Taylor series expansions of the corresponding generating functions with being the order of consistency. Technically all the coefficients ω ,j can be computed by using any implementation of the fast Fourier transform. For the starting weights w n,j , we can consider the fixed s = 0.
In this case, we obtain the Lubich formulae for fractional integrals when f (0) = 0. For s = 0, the coefficients w n,j can be constructed such that Equation (49) exactly holds for power functions. Therefore, we recover s ∑ j=0 w n,j j q = Γ(q + 1) and it makes sense to choose s = − 1.
In the case with the lower terminal a = 0, we can readily adopt affine transform to modify the fractional linear multistep methods.

Remark 3.
Apart from the choice given by Equation (50), which corresponds to the fractional BDF, there are alternatives for the generating functions of the convolution coefficients. When we choose as the generating function, the fractional trapezoidal rule with second order accuracy for α ≥ 0 is obtained.
and set Then we obtain the generating function for the coefficients of the generalized Newton-Gregory formulae, which is convergent of order . Direct calculation gives γ 0 = 1 and γ 1 = − α 2 . Then the corresponding generating function for the second order generalized Newton-Gregory formula is given by More details for generating functions can be found in Refs. [45,49-51].

Diffusive Approximation
The above numerical methods may result in expensive computational costs. To eliminate this deficiency, the diffusive approximation reformulates the model containing the fractional integral as a system of differential equations.
Recalling the relations and the fractional integral with 0 < α < 1 can be rewritten as [52] Define the variable transformation z = (x − t)ω 2 , ω ≥ 0. The Fubini's Theorem yields Introducing the auxiliary function we have It follows from the definition of φ(ω, x) that the auxiliary function satisfies In this case, evaluating Riemann-Liouville integral RL D −α 0,x f (x) consists of two steps: solving the first order differential equation (63), and computing the infinite integral in Equation (62) via suitable quadratures.
Instead of utilizing the properties of the Gamma function, Chatterjee adopted a popular integral representation [53,54] x Consequently, the Fubini's Theorem gives where g(z, x) is defined as In order to generate nonreflecting boundary conditions [55] and accelerate convolutions with the heat kernel [56], literatures such as Ref. [57] usually recognize where Alternatively, other literatures have regarded g(z, x) as the solution of a first order ordinary differential (ODE) equation [52,53,58,59], Any approximate method for ODEs can be used to obtain g(z, x), x = ∆x, 2∆x, . . ., in an amount of work that is linear.
The principle difficulty of implementing both approaches lies in the discretization of the integrals on the right-hand side of Equations (62) and (65). The choices of quadrature nodes and corresponding weights have been investigated in several literatures, see Refs. [57,60,61] for more details.

Numerical Approximations to Fractional Derivatives
We introduce the existing numerical evaluations to Caputo, Riemann-Liouville, and Riesz derivatives in this section. The basic ideas of these methods are presented as well.

Numerical Caputo Differentiation
Caputo derivatives in Equations (6) and (7) can be viewed as Riemann-Liouville fractional integrals of integer-order derivatives. As a result, most of the numerical evaluations of Caputo derivatives follow from those of fractional integrals. We derive numerical evaluations of the Caputo derivative as follows.

L1, L2, and L2C Methods
The well-known L1 method was originally introduced in Ref. [62] to evaluate Riemann-Liouville derivative with 0 < α < 1, which equivalently reads as Note that the second term on the right-hand side happens to be Caputo derivative with 0 < α < 1. That is the reason why we introduce the L1 method when considering numerical approximations to the Caputo derivative.
Here the coefficients are given by Normally, the L1 method can lead to unconditionally stable algorithms [63][64][65][66][67][68][69]. Therefore, it is frequently used in the discretization of time fractional differential equations. Since the proof for this scheme available is not very direct or a little cryptic, it is necessary to present clear proof of its truncated error for reference as it is mostly used.
Then it holds that where C is a positive constant given by

Proof. Denote
Then it immediately follows that Using the Taylor expansion with integral remainder, we have for t ∈ [x k , x k+1 ], which yields Exchanging the order of integration gives In the following, we show that when 0 < α < 1, Denote g(s) = (x j − s) 1−α . Then it holds for any s ∈ (x k , x k+1 ) that with certain ξ k ∈ (x k , x k+1 ). As a result, inequality (81) holds. For the inequality (82), one has and The above two equalities yield that Equation (82) holds.
Combining the above analysis, one has Inserting the above estimate into Equation (80) gives All this ends the proof.

Remark 4.
The idea of proving Theorem 2 is borrowed from Ref. [70] where the case with α ∈ (1, 2) was considered. Such an estimate was also considered in [71].
Then the classical L1 method is generalized into the L1 method on nonuniform grids for Caputo derivative [72] provided that and the coefficients are given by In the special case of nonuniform grids with Here the coefficients are given by Replacing yields the following modified L1 method for Caputo derivative [38]

Remark 5. (I)
The modified L1 method (92) is useful to obtain the Crank-Nicolson scheme for the time-fractional subdiffusion equation [73,74], which can be regarded as a natural extension of the classical Crank-Nicolson scheme [75].
(II) The (weak) singularity makes it difficult to evaluate fractional derivatives. In this case, approximations such as Equations (88) and (92) on nonuniform meshes or graded meshes can be utilized. One can refer to [72,76,77] for more details in this respect.
For the case with 1 < α < 2 and the lower terminal a = 0, there holds In Ref. [78], the integral then the L2C method for Caputo derivative is obtained. Here the coefficients are given by Note that in the above two schemes the value of f (x −1 ) is needed. We can set f (x −1 ) = f (x 1 ) when the condition f (0) = 0 is met. For the case with lower terminal a = 0, we can utilize affine transformation before applying the L2 and L2C methods. Remark 6. The L2 and L2C methods reduce to the backward difference method and the central difference method for the first order derivative, respectively, when α = 1. If α = 2, the L2 method reduces to the central difference method for the second order derivative and the L2C method reduces to with the first order accuracy. As a matter of fact, the error bound for the L2 method is O(h 3−α ).
Numerical experiments indicate that the L2 method is more accurate than the L2C method for 1.5 < α < 2, while the opposite result appears when 1 < α < 1.5. And these two methods behave in almost the same way near α = 1.5 [78].

Numerical Methods Based on Polynomial Interpolation
It is evident that the higher-order accuracy can be achieved by utilizing the higher-order interpolation, provided that f (x) is suitably smooth. In the following, we introduce numerical approximations in this respect.
In this case, we have the following (3 − α)-th order approximation [79], where 0 < α < 1, R j denotes the truncated error, and the coefficients are given by with 0 ≤ k ≤ j − 1 and 1 ≤ j ≤ N.
Since the above (3 − α)-th order method is also widely used, we estimate its truncated error in detail.
with c being a positive constant and f (x −1 ) in Equation (102) being used.
Proof. It is clear that the truncated error is given by Interchanging the order of integrations yields For k = 0, 1, . . . , j − 1, denote and where the expression of A k can be derived from Equations (106) and (107) so is left out due to lengthiness.
It is evident that b 1 = 1 2−α − 1 2 , and for l ≥ 2, Thus, it holds for l ≥ 2 that As a result, with C 2 > 0 being a constant. Note that A k contains all the terms in Equation (106) with the form of integrals over [x k , x k+1 ]. Then the affine transformation s = x k + ξh, ξ ∈ [0, 1] and l = j − k, k = 0, 1, . . . , j − 1 yield Rewrite a l (ξ) in the form with a n (ξ) = For n ≥ 2, we have a n (ξ) ≥ 0 for arbitrary ξ ∈ [0, 1]. To see this, recall that When and a n (ξ) Note that One has a n (ξ 0 ) < a n (1) = 0, and there exits ξ 1 ∈ (0, ξ 0 ) such that a n (ξ 1 ) = 0 since a n (0) > 0. Therefore, Since it holds that a n (ξ) ≥ 0 for arbitrary ξ ∈ [0, 1] when n ≥ 2. As a result, Furthermore, Especially, for l ≥ 2, As a result, it holds that with C 1 > 0 being a constant. Consequently, the truncated error has the bound Note that the derivative x 0 ] is needed in the above inequality. In this case, Consequently, the desired estimate is obtained.  Table 3. It is obvious that the convergence order is (3 − α), which is in line with the theoretical analysis. In Ref. [81], another (3 − α)-th order approximation was proposed. Denote Let f (x) ∈ C 3 [a, b] and 0 < α < 1. We utilize the linear interpolation on the first interval [x 0 , x 1 ], and the quadratic interpolation with the truncated errors with a (α) Numerical results in Ref. [81] imply that the computational errors given by the L1-2 formula are obviously much smaller than those of the L1 formula.
Applying the quadratic interpolation which is different from the one defined in Equation (129) to approximating f (x), and utilizing the when j = 0, and for j ≥ 1, with a The comparison between the L2-1 σ and L1-2 methods in Ref. [82] shows that the L2-1 σ formula refines the accuracy indeed.

Remark 8 ([80]
). The L2-1 σ formula for the right-sided Caputo derivative can be derived in a similar manner. In this case, the parameter should be chosen as σ = α 2 , α ∈ (0, 1). The corresponding approximation is given by Here the coefficients are given by where a (α,σ) k For other (3 − α)-th order approximations to Caputo derivative based on interpolation, one may refer to Refs. [83,84].
(II) (4 − α)-th order approximation with On the second subinterval [x 1 , x 2 ], we similarly obtain through the quadratic interpolation, where For the remaining subintervals, we use the cubic interpolation function to approximate f (x). Consequently, it holds that where the coefficients are given by , , In view of the above analysis, we obtain the numerical approximation [85] The coefficients g k have different values for different j. When j = 1, When j = 5, When j ≥ 6, If f (x) ∈ C 4 [x 0 , x j ] and α ∈ (0, 1), the truncated error R j in Equation (150) satisfies Numerical examples in Ref. [85] verify the above theoretical results.

Example 4.
Suppose that 0 < α < 1 and f (x) = x 4 . Evaluate the α-th order Caputo derivative of f (x) at x = 1 by Equation (150). Maximum errors (ME) and convergence order (CO) are presented in Table 4.  The maximum errors (ME) and convergence order (CO) are shown in Table 5.  (III) (r + 1 − α)-th order approximation Generalizing the above (4 − α)-th order approximation, an (r + 1 − α)-th order approximation was proposed in Ref. [86] by virtue of the Lagrange polynomials of degree r. Let f (x) ∈ C r [a, b] (r ≥ 4) and 0 < α < 1. On the subintervals [x k−1 , x k ], j ≥ k ≥ r, N ≥ j ≥ r, we utilize the Lagrange polynomial to approximate f (x). Denote Then it holds that To compute the coefficients ω r i,j−k , we denote by and Here φ k j,i and ψ k j,i are the sums of products of all different combinations of k elements in the sets On the subinterval [x k−1 , x k ], 1 < k < r, 1 ≤ j ≤ N, there are no enough nodes to obtain an r-th degree Lagrange polynomial. In this case, we use I k [p k (x)] to approximate the integral In summary, we obtain the following approximation [86] with R j r being the truncated error. It has been proved that when f (x) ∈ C r [a, b] (r ≥ 4), the truncation error satisfies Numerical examples in Ref. [86] verify the above theoretical results.
Example 6. Suppose 0 < α < 1, and let f (x) = x 6 , x ∈ [0, 1]. Use scheme (164) to compute Caputo derivative of f (x) at x = 1 with different stepsizes. Table 6 lists the computational errors and convergence orders at x = 1 with different values for α, and r = 4, 5. It can be observed that the numerical convergence order of the utilized scheme is (r + 1 − α).  Table 7 lists the numerical results with different values for α, and r = 4, 5. It is evident that scheme (164) can reach (r + 1 − α)-th order accuracy. Remark 9. The (3 − α)-th, (4 − α)-th, and (r + 1 − α)-th order numerical schemes established in Refs. [79,85,86] are of unconditional stability in the practical sense when solving fractional differential equations. In other words, numerical schemes for fractional differential equations based on these approximations are stable only if α lies in their respective subsets of the interval (0, 1). On the other hand, there are some other interesting methods in this respect. See [87,88] for more details.
(IV) Spectral approximations Let m − 1 < α < m ∈ Z + and f (x) ∈ H r [a, b] with r ≥ 2m. Now we introduce spectral approximations to Caputo derivative [42,44]. Here we take the Jacobi approximation as a representative example since the others such as Chebyshev approximation are special cases of the Jacobi one. Let the polynomial be an approximation of f (x) based on the Jacobi polynomials. Recall that It holds that where p u,v j and P u+m,v+m,m−α j (x) are defined by Equations (32) and (36). Denote Then it holds that The affine transformation For the corresponding differential matrix, see Ref. [44] for more details.
The following numerical examples verify the efficiency of the spectral approximation.

Example 8 ([42]
). Let f (x) = x µ , x ∈ [0, 1]. We use formula (170) to compute C D α 0,x f (x). Table 8 shows the absolute maximum errors at the Jacobi-Gauss-Lobatto points. The spectral accuracy is obtained.   Table 9 presents the absolute maximum errors for the cases of u = v = 0 and u = v = − 1 2 . The expected results can be observed.   (V) Radial basis function discretization Being a natural generalization of univariate polynomial splines to a multivariate setting, radial basis functions work for arbitrary geometry with high dimensions and it does not require a mesh at all [89]. Numerically solving fractional differential equations based on radial basis functions has attracted sustained attention in engineering and science community. See [90][91][92][93] and references cited therein. In [94], radial basis functions are utilized to evaluate fractional differential operators. In the following, we introduce the basic idea of this method.
Take the one-dimensional case as an example. Let x j (j = 1, 2, . . . , N) be the collocation points in the interval [a, b]. An radial basis function interpolant of a given function f (x) is defined in the form In order to take the values f (x i ), i = 1, 2, . . . , N, the expansion coefficients λ j are required to satisfy the matrix form Here φ(·) is the radial basis function. Some popular choices of radial basis function are cubic (φ(r) = r 3 ), multiquadrics (φ(r) = √ r 2 + c 2 ), and Gaussian (φ(r) = e −cr 2 ), where the free parameter c is called the shape parameter for the radial basis function. The smooth radial functions (such as multiquadrics and Gaussian) give rise to spectrally accurate function representation while the piecewise smooth radial functions (such as cubic) only produce algebraically accurate representations [94]. Applying the Caputo differentiation operator to (171) yields which can be written in the matrix form a,x f (x) at the point x = x i . Note that the collocation matrix A is unconditionally nonsingular [94]. Combing equations (172) and (174) gives Therefore, the differential matrix D = BA −1 yields an radial basis function discretiazation of the operator C D α a,x .

Remark 10. (I)
The above procedure of deriving differential matrix based radial basis functions is applicable for other fractional differentiation operators as well.
(II) Finding a closed form analytical expression for the fractional derivative of a given function may be challenging. In practice, one has to represent the radial basis function in the form of Taylor series before applying fractional differentiation operator term by term. Then the infinite sum can be truncated once the terms are smaller in magnitude than machine precision.
(III) The standard radial basis function methods may result in ill-conditioning which often impairs the convergence. To offset this deficiency, the so-called RBF-QR method can be utilized instead of the standard one. See Ref. [94] for more details.

Fractional Backward Difference Formulae
It has been mentioned in Ref.
[45] that the fractional linear multistep method is applicable for numerical Riemann-Liouville differentiation, provided that the generating functions are properly chosen. We can therefore derive fractional backward multistep methods for Caputo derivative, on the basis of the relation where m − 1 < α < m ∈ Z + . In Ref. [95], shifted fractional backward difference formulae for Caputo derivative were derived through three steps. Shifted Lubich formulae for Riemann-Liouville derivative on bounded domain were first introduced for the case with homogeneous conditions. At that stage, generating functions of the coefficients were constructed to maintain high-order accuracy. By virtue of adopting suitable auxiliary functions, the shifted formulae were modified for the case with inhomogeneous conditions. Finally, the shifted formulae for Caputo derivative are obtained based on relation (176). Theoretical results which can be proved through Fourier analysis are as follows. and where the weights ζ and, Here the shift p ≤ 3α 2 , and the coefficients [a, b], and that the derivatives of f (x) up to order [α] + 5 belong to L 1 [a, b]. Then the third order schemes are given by and where the weights ζ and, 3,p,k , k = 0, 1, . . . , j + p − 4, Here the shift p satisfies 3p 2 − 12αp + 11α 2 ≥ 0, and (α) 3,p,k (0 ≤ k ≤ j + p) are given by with and where the weights ζ and, Here the shift p satisfies 25α 3 − 35α 2 p + 15αp 2 − 2p 3 ≥ 0, and (α) 4,p,k (0 ≤ k ≤ j + p) are given by Different from numerical algorithms based on the polynomial interpolation, in which the corresponding accuracy depends on the derivative order α and homogenous conditions are needed, the formulae presented in Theorems 4-6 for Caputo derivatives have no restriction on the initial conditions, and are of integer-order accuracy. Here we display two numerical examples to verify these arguments.

Numerical Riemann-Liouville Differentiation
Now we consider numerical approximations to Riemann-Liouville derivatives. The relation (176) indicates that numerical evaluations of Riemann-Liouville derivative can be readily obtained based on those of Caputo derivative. Numerical approximations derived from evaluations of Caputo derivative are therefore omitted in this section. Here we present alternative approaches.

Numerical Methods Based on Linear Spline Interpolation
Let f (x) ∈ C 4 [a, b] and 1 < α < 2. For j = 1, 2, . . . , N − 1, there holds where and Then we obtain an approximation to I l α (x j ) given by As a result, there holds [96] Similarly, the right-sided Riemann-Liouville derivative can be approximated by [96] where The truncated error of this approach has been proved in Ref. [96] to be O(h 2 ) provided that f (4) (x) has compact support on [a, b].
In the particular cases with a = −∞ and b = +∞, Equations (205) and (206) can be written as [96][97][98] and with Here Let f (x) ∈ L 1 (R), RL D α+1 −∞,x f , and its Fourier transform belong to L 1 (R). It can be verified through the Fourier analysis that the shifted Grünwald-Letnikov difference operator with p ∈ Z approximates the left-sided Riemann-Liouville derivative RL D α −∞,x f (x) with first order accuracy. Assume that the following weighted and shifted Grünwald-Letnikov difference approximates the Riemann-Liouville derivative with second order accuracy. Then the Fourier transform gives and Therefore, a 1 and b 1 need to satisfy to assure that B α h,p,q f (x) is of second order accuracy. In other words, holds uniformly when a 1 = α−2q 2p−2q and b 1 = 2p−α 2p−2q [103].
A third order WSGD operator was also proposed in Ref. [103]. Nevertheless, it fails to obtain stable numerical schemes when α ∈ (1, 2). To offset this situation, the compact-WSGD operator [104] was introduced through combining WSGD operators with Taylor expansions of the shifted Grünwald formula for sufficiently smooth function f (x) that satisfies homogeneous initial conditions.
The construction of the WSGD operators implies the possibility of deriving higher-order numerical approximations to Riemann-Liouville derivative by imposing various weights and shifts on higher-order Lubich formulae. In Ref. [105], numerical algorithms with second, third, and fourth order accuracy are proposed based on the second order Lubich formula.

Fractional Backward Difference Formulae and Their Modifications
It is evident that the classical Grünwald-Letnikov approximation coincides with the first order Lubich formula (51) with α > 0 replaced by −α. In fact, Lubich formulae (51) are applicable for evaluating Riemann-Liouville derivative indeed.
Remark 12. In Ref. [106], the coefficients of high-order approximations (till 10-th order) for Riemann-Liouville derivative were computed. Furthermore, a conjecture on coefficients of the third, fourth, and fifth order schemes was proposed by Li and Ding and was rephrased on Page 80 of Ref. [38], stated as follows, 4,l+1 for l ≥ 7, and (α) Recently, the above conjecture for (α) 3,l with 0 < α < 1 has been proved in Ref. [107].
Similarly to the case of the classical Grünwald-Letnikov approximation, the classical Lubich formulae may produce unstable numerical schemes for fractional differential equations due to the eigenvalue issue [105]. In this case, we often introduce shifts. To maintain the high-order accuracy, the corresponding generating functions need modifying. The shifted fractional backward difference formulae [108], which can be proved via the Fourier transform method, are presented as follows.
Theorem 7. Suppose that f (x) ∈ C [α]+3 (R), and all the derivatives of f (x) up to order [α] + 4 belong to L 1 (R). Then we have and as h → 0. Here k (α) 2,l (l = 0, 1, . . .) are generated by , and all the derivatives of f (x) up to order [α] + p + 2 belong to L 1 (R). Then there hold and For higher-order weighted and shifted Lubich formulae such as the fourth order one, see Ref. [105] for more details. Remark 13. All the above weighted and shifted Lubich formulae are applicable to the bounded domain (a, b) through performing zero extension, whenever the zero extended function satisfies the corresponding assumptions of the approximations.
An alternative approach modifying the Lubich formulae is to introduce compact operators, which gives the following fractional-compact formulae [109]. The corresponding accuracy can be proved by the Fourier transform.
Theorem 12. Define the following two difference operators, where the coefficients k (α) 2,l are given by the function If we introduce the fractional-compact difference operator with δ 2 x being a second order central difference operator defined by δ 2 hold uniformly for x ∈ R, provided that f (x) ∈ C [α]+4 (R) and all derivatives of f (x) up to order [α] + 5 belong to L 1 (R).
Theorem 13. Choose the generating function as and define the following difference operators Then the equalities hold uniformly for x ∈ R, provided that f (x) ∈ C [α]+4 (R) and all derivatives of f (x) up to order [α] + 5 belong to L 1 (R). Here the fractional-compact difference operator L is given by where Then and hold uniformly on R.
The idea of the above approximations can be applied to evaluating tempered fractional derivatives, see Ref. [110] for more details.

Fractional Average Central Difference Method
In Ref. [111], a shifted operator of the form with h > 0 was proposed to approximate Riemann-Liouville derivative. This difference operator reduces to the standard central difference operator when α is a positive integer. It can be verified through the Fourier transform that the equality holds uniformly for x ∈ R as h → 0, provided that f ∈ C [α]+3 (R) and all derivative of f up to the order [α] + 4 exist and belong to L 1 (R). Recently, the above shifted operator has been modified to obtain higher-order approximations, which are based on the fractional left and right average central difference operators [112] and The main results, which can be verified through Fourier transform, are presented as follows.
Theorem 15. Assume that f (x), the Fourier transform of RL D α+2 −∞,x f (x) and RL D α+2 x,+∞ f (x) are in L 1 (R). Then the equalities hold uniformly on R.
Theorem 16. When f (x) and the Fourier transforms of RL D α+4 −∞,x f (x) and RL D α+4 x,+∞ f (x) are in L 1 (R), then the relations and hold uniformly on R. Here δ 2 x denotes the second order central difference operator defined by δ 2 For functions defined on [a, b], the fractional average central difference formulae can be modified through suitable extensions.

Numerical Riesz Differentation
Since Riesz derivative can be viewed as a linear combination of the left-and right-sided Riemann-Liouville derivatives. Several numerical approximations to Riesz derivative can be readily obtained based on the aforementioned methods of Riemann-Liouville derivative. Here we only present the one based on L2-1 σ formulae when introducing indirect evaluations of Riesz derivative. For more details of these indirect approaches, one can refer to Refs. [108,109,112,113]. Then we focus on introducing some schemes evaluating Riesz derivative in direct ways.

Asymmetric Centred Difference Operators
Slightly different from the fractional average central difference operators in the previous section, the symmetric fractional centred difference operator is defined as with the coefficients being given by , k = 0, ±1, ±2, . . .
It can be verified through the Fourier analysis that the relation [114] RZ D α holds uniformly for x ∈ R as h → 0 + , provided that f ∈ C 5 (R) and all of its derivatives up to order five belong to L 1 (R). It has been pointed out by Ref. [115] that Equation (269) also holds for 0 < α ≤ 1.

Weighted and shifted centred difference operators
Introducing shifts to the symmetric fractional centred difference operator in Equation (267), the shifted centred difference operators [112] can be obtained. To achieve high-order accuracy, the following high-order approximations to Riesz derivative can be derived through combining these shifted operators with suitable weights.
Theorem 17 ([116]). If f (x) lies in C 7 (R) with all the derivatives up to order 7 in L 1 (R), then the relation holds uniformly for x ∈ R.
Theorem 18 ([112]). Assume that f (x) ∈ C 9 (R) with all the derivatives up to order 9 in L 1 (R). Then Theorem 19 ([112]). Assume that f (x) lies in C 11 (R) with all the derivatives up to order 11 in L 1 (R). Then For much higher-order difference operators in this respect, see Ref. [112].

Compact Centred Difference Operators
As another variant of the centred difference operator, compact centred difference operators are based on the idea of introducing compact operators to maintain even-order accuracy.
Remark 15. One should bear in mind that some suitable smooth conditions for f (x) are necessary and cannot be dropped. Once these conditions are violated, the expected accuracy cannot be achieved. Example 13 will verify this assertion latter.

Conclusions
In this paper we focus on numerical approximations to fractional integrals and derivatives, which are essential for solving fractional differential equations. This work is targeted at systematically clarifying basic ideas of the existing numerical evaluations, which provides the readers with comprehensive understanding of numerical methods for fractional calculus.
As the experimental advances further reveal nonlocality, memory, and hereditary properties of numerous materials and processes, the importance of the fractional calculus is becoming obvious. We hope that this work, which is designated to compressively review numerical approximations to fractional calculus, will become the first step in elucidating underlying principles and results of a wider variety of fractional dynamics.
Author Contributions: Investigation, M.C. and C.L.; writing-original draft preparation, M.C. and C.L.; writing-review and editing, M.C. and C.L. These two authors contributed equally to this paper. All authors have read and agreed to the published version of the manuscript.
Funding: The present job was partially supported by the National Natural Science Foundation under grant no. 11872234.