Numerical Approaches to Fractional Integrals and Derivatives: A Review

Abstract

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.

1. Introduction

1.1. Historical Review

The primary attempt, which was recorded in history to discuss the idea of generalizing the integer-order differentiation $\frac{{\mathrm{d}}^{n}f(t)}{\mathrm{d}{t}^n}$ to $\frac{{\mathrm{d}}^{\alpha }f(t)}{\mathrm{d}{t}^{\alpha }}$ with non-integer $\alpha$, 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

$${\int }_a^x\frac{\phi (t)\mathrm{d}t}{{(x-t)}^{\mu }}=f(x),x>a,0<\mu <1,$$

(1)

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-\mu$. 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, electrodynamics, quantum mechanics, plasma dynamics, and cosmic rays [27]. Achievements in this regard were also presented in Refs. [28,29,30,31].

1.2. Current Situations

For the time being, the most frequently utilized fractional integrals and derivatives in applications are the left- and right-sided Riemann–Liouville integrals [32] (fractional integrals for short),

$${}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)=\frac{1}{\mathsf{\Gamma }(\alpha )}{\int }_a^x\frac{f(t)\mathrm{d}t}{{(x-t)}^{1-\alpha }},\alpha >0,$$

(2)

$${}_{RL}{\mathrm{D}}_{x,b}^{-\alpha }f(x)=\frac{1}{\mathsf{\Gamma }(\alpha )}{\int }_x^b\frac{f(t)\mathrm{d}t}{{(t-x)}^{1-\alpha }},\alpha >0,$$

(3)

the left- and right-sided Riemann-Liouville derivatives [32],

$${}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)=\frac{1}{\mathsf{\Gamma }(m-\alpha )}\frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}{\int }_a^x\frac{f(t)\mathrm{d}t}{{(x-t)}^{\alpha +1-m}},m-1\le \alpha <m\in {\mathbb{Z}}^{+},$$

(4)

$${}_{RL}{\mathrm{D}}_{x,b}^{\alpha }f(x)=\frac{{(-1)}^m}{\mathsf{\Gamma }(m-\alpha )}\frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}{\int }_x^b\frac{f(t)\mathrm{d}t}{{(t-x)}^{\alpha +1-m}},m-1\le \alpha <m\in {\mathbb{Z}}^{+},$$

(5)

the left- and right-sided Caputo derivatives [32],

$${}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)=\frac{1}{\mathsf{\Gamma }(m-\alpha )}{\int }_a^x\frac{f^{(m)}(t)\mathrm{d}t}{{(x-t)}^{\alpha +1-m}},m-1\le \alpha <m\in {\mathbb{Z}}^{+},$$

(6)

$${}_C{\mathrm{D}}_{x,b}^{\alpha }f(x)=\frac{{(-1)}^m}{\mathsf{\Gamma }(m-\alpha )}{\int }_x^b\frac{f^{(m)}(t)\mathrm{d}t}{{(t-x)}^{\alpha +1-m}},m-1\le \alpha <m\in {\mathbb{Z}}^{+},$$

(7)

and Riesz derivative [33]

$${}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)={\mathsf{\Psi }}_{\alpha }\left[{}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)+{}_{RL}{\mathrm{D}}_{x,b}^{\alpha }f(x)\right],0<\alpha \ne 1,3,5,\dots ,$$

(8)

where ${\mathsf{\Psi }}_{\alpha }=-\frac{1}{2cos(\frac{\alpha \pi }{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 $\alpha$-th order Riemann–Liouville derivative of $sinx$ and $cosx$ with $\alpha >0$ are not $sin(x+\frac{\alpha \pi }{2})$ and $cos(x+\frac{\alpha \pi }{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=\frac{b-a}{N},N\in {\mathbb{Z}}^{+}$, and $x_j=a+jh,j=0,1,2,\dots ,N$.

2. 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.

2.1. Numerical Methods Based on Polynomial Interpolation

Assume that $f(x)$ is suitably smooth on $[a,b]$. Then the $\alpha$-th order fractional integral of $f(x)$ at $x=x_j$ with $1\le j\le N$ can be expressed as

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}=\frac{1}{\mathsf{\Gamma }(\alpha )}{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{\alpha -1}f(t)\mathrm{d}t.$$

(9)

It is reasonable to utilize an interpolate function $\tilde{f}(x)$ to approximate $f(x)$ on each subinterval, such that the integral ${\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{\alpha -1}\tilde{f}(t)\mathrm{d}t$ can be calculated exactly. This idea yields a series of numerical formulae in the form

$${\left[R_L{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx {\displaystyle \sum _{k=0}^j}{\omega }_{k}f(x_k),1\le j\le N,$$

(10)

where ${\omega }_k(k=0,1,\dots ,j)$ are the corresponding coefficients. To better understand this method, we retrospect some specific formulae with their brief derivations.

If $f(x)\in C[a,b]$ on the right-hand side of Equation (9) is approximated by a piecewise constant function

$$\tilde{f}(x)=f(x_k),x\in [x_k,x_{k+1}),0\le k\le j-1,$$

(11)

then there holds

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx \frac{1}{\mathsf{\Gamma }(\alpha )}{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{\alpha -1}f(x_k)\mathrm{d}t.$$

(12)

This yields the left fractional rectangular formula [38]

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx {\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}f(x_k),$$

(13)

where the convolution coefficients $b_k(0\le k\le j-1)$ are given by

$$b_k=\frac{1}{\mathsf{\Gamma }(\alpha )}{\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{\alpha -1}\mathrm{d}t=\frac{h^{\alpha }}{\mathsf{\Gamma }(\alpha +1)}\left[{(k+1)}^{\alpha }-k^{\alpha }\right].$$

(14)

Similarly, if the function $f(x)$ on the right-hand side of Equation (9) is replaced by

$$\tilde{f}(x)=f(x_{k+1}),x\in (x_k,x_{k+1}],$$

(15)

then we have the right fractional rectangular formula [38]

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx {\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}f(x_{k+1}),$$

(16)

with $b_k(0\le k\le j-1)$ given by Equation (14). Based on the left and right rectangular formulae, the weighted fractional rectangular formula [38] yields

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx {\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}\left[\theta f(x_k)+(1-\theta )f(x_{k+1})\right],0\le \theta \le 1,$$

(17)

or the similar form

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx {\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}f\left(x_k+(1-\theta )h\right),0\le \theta \le 1.$$

(18)

Remark 1.

(I) The left fractional rectangular formula (13) and the right fractional rectangular formula (16) will be recovered if $\theta =1$ and $\theta =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 $\alpha =1$ and $\theta =\frac{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 $f^{\prime }({\xi }_k)(t-x_k),0\le k\le j-1$ and $f^{\prime }({\eta }_{k+1})(t-x_{k+1}),0\le k\le j-1$, respectively. Therefore, the accuracy of fractional rectangular formulae are of first order accuracy for all $0\le \theta \le 1$. And all the above fractional rectangular formulae are of the first order accuracy.

Assume that $f(x)\in C[a,b]$. Replacing $f(x)$ in Equation (9) with the piecewise linear polynomial

$$\tilde{f}(x)=\frac{x_{k+1}-x}{x_{k+1}-x_k}f(x_k)+\frac{x-x_k}{x_{k+1}-x_k}f(x_{k+1}),x\in [x_k,x_{k+1}],$$

(19)

we obtain the fractional trapezoidal formula [38]

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx \frac{h^{\alpha }}{\mathsf{\Gamma }(\alpha +2)}{\displaystyle \sum _{k=0}^j}{a}_{k,j}f(x_k).$$

(20)

Here the coefficients $a_{k,j}(0\le k\le j)$ are given by

$$\left\{\begin{array}{cc}\hfill a_{0,j}=& {(j-1)}^{\alpha +1}-(j-1-\alpha )j^{\alpha },\hfill \\ \hfill a_{k,j}=& {(j-k+1)}^{\alpha +1}-2{(j-k)}^{\alpha +1}+{(j-k-1)}^{\alpha +1},1\le k\le j-1,\hfill \\ \hfill a_{j,j}=& 1.\hfill \end{array}\right.$$

(21)

Suppose that $f(x)\in C[a,b]$. For $0\le k\le j-1$, let $\{l_{k,i}(x)\}$ be Lagrangian functions defined on the grid points $\{x_{k+s},s\in S\}$ with $S=\{0,\frac{1}{2},1\}$, which are given by

$$l_{k,i}(x)=\prod _{s\in S,s\ne i}\frac{x-x_{k+s}}{x_{k+i}-x_{k+s}},i\in S.$$

(22)

Denote $x_{k+\frac{1}{2}}=\frac{x_k+x_{k+1}}{2}$ and utilize the piecewise quadratic polynomial

$$\tilde{f}(x)={\displaystyle \sum _{i\in S}}{l}_{k,i}(x)f(x_{k+i}),x\in [x_k,x_{k+1}].$$

(23)

Then we obtain the following fractional Simpson’s formula [38]

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx \frac{h^{\alpha }}{\mathsf{\Gamma }(\alpha +3)}\left[{\displaystyle \sum _{k=0}^j}{c}_{k,j}f(x_k)+{\displaystyle \sum _{k=0}^{j-1}}{\widehat{c}}_{k,j}f(x_{k+\frac{1}{2}})\right],$$

(24)

in which

$$\begin{array}{cc}\hfill {\widehat{c}}_{k,j}=& 4(\alpha +2)\left[{(j-k)}^{1+\alpha }+{(j-k-1)}^{1+\alpha }\right]-8\left[{(j-k)}^{2+\alpha }-{(j-k-1)}^{2+\alpha }\right],0\le k\le j-1,\hfill \end{array}$$

(25)

and

$$\left\{\begin{array}{cc}\hfill c_{0,j}& =4\left[j^{2+\alpha }-{(j-1)}^{2+\alpha }\right]-(\alpha +2)\left[3j^{1+\alpha }+{(j-1)}^{1+\alpha }\right]+(\alpha +2)(\alpha +1)j^{\alpha },\hfill \\ \hfill c_{k,j}& -(\alpha +2)=[{(j-k+1)}^{1+\alpha }+{(j-k-1)}^{1+\alpha }+6{(j-k)}^{1+\alpha }\hfill \\ & +4\left[{(j-k+1)}^{2+\alpha }-{(j-k-1)}^{2+\alpha }\right],1\le k\le j-1,\hfill \\ \hfill c_{j,j}& =2-\alpha .\hfill \end{array}\right.$$

(26)

Assume that $f(x)\in C[a,b]$. Let $f(x)$ be approximated by the following r-th degree polynomial on the grid points $\left\{x_k=x_0^{(k)},x_1^{(k)},\dots ,x_{r-1}^{(k)},x_r^{(k)}=x_{k+1}\right\}$,

$$p_{k,r}(x)={\displaystyle \sum _{i=0}^r}{l}_{k,i}(x)f(x_i^{(k)}),x\in [x_k,x_{k+1}],$$

(27)

with

$$l_{k,i}(x)={\displaystyle \prod _{\begin{array}{c}n=0,\\ n\ne i\end{array}}^r}\frac{x-x_n^{(k)}}{x_i^{(k)}-x_n^{(k)}}.$$

(28)

Then we obtain the fractional Newton–Cotes formula [38]

$${\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\right]}_{x=x_j}\approx {\left[{}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }{p}_{k,r}(x)\right]}_{x=x_j}={\displaystyle \sum _{k=0}^{j-1}}{\displaystyle \sum _{i=0}^r}{A}_{i,j}^{(k)}f(x_i^{(k)}),$$

(29)

with the coefficients being calculated by

$$A_{i,j}^{(k)}=\frac{1}{\mathsf{\Gamma }(\alpha )}{\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{\alpha -1}{l}_{k,i}(t)\mathrm{d}t.$$

(30)

Remark 2.

It has been demonstrated in Ref. [38] that the error estimate of Equation (29) is$\mathcal{O}(h^{r+1})$provided that$f\in C^{r+1}([a,b])$. The error estimate does not equal that for the classical one, which is$\mathcal{O}(h^{r+2})$. This inconsistency may be due to the asymmetry of the weakly singular kernel${(x_j-x)}^{\alpha -1}$, which leads to the non-symmetry of the remainder term${(x_j-x)}^{\alpha -1}{\displaystyle \prod _{i=0}^r}(x-x_i^{(k)})$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\ge 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

$$p_N(x)={\displaystyle \sum _{j=0}^N}{\tilde{p}}_j^{u,v}{P}_j^{u,v}(x),$$

(31)

based on the Jacobi polynomials ${\{P_j^{u,v}(x)\}}_{j=0}^N(u,v>-1)$. Here

$${\tilde{p}}_j^{u,v}=\frac{1}{{\delta }_j^{u,v}}{\displaystyle \sum _{k=0}^N}f(x_k)P_j^{u,v}(x_k){\omega }_k,j=0,1,\dots ,N,$$

(32)

with ${\{x_k\}}_{k=0}^N$ and ${\{{\omega }_k\}}_{k=0}^N$ being the collocation points and the corresponding quadrature weights [43]. The constants ${\delta }_j^{u,v}$ are given by

$${\delta }_j^{u,v}=\left\{\begin{array}{cc}{\gamma }_j^{u,v},\hfill & j=0,1,\dots ,N-1,\hfill \\ (2+\frac{u+v+1}{N}){\gamma }_N^{u,v},\hfill & j=N,\hfill \end{array}\right.$$

(33)

with ${\gamma }_j^{u,v}$ being defined by

$${\gamma }_j^{u,v}=\frac{2^{u+v+1}\mathsf{\Gamma }(j+u+1)\mathsf{\Gamma }(j+v+1)}{(2j+u+v+1)j!\mathsf{\Gamma }(j+u+v+1)}.$$

(34)

In this case, we have the following spectral approximation based on Jacobi polynomials

$${}_{RL}{\mathrm{D}}_{-1,x}^{-\alpha }f(x)\approx {}_{RL}{\mathrm{D}}_{-1,x}^{-\alpha }{p}_N(x)={\displaystyle \sum _{j=0}^N}{\tilde{p}}_j^{u,v}{\widehat{P}}_j^{u,v,\alpha }(x),x\in [-1,1].$$

(35)

Here ${\widehat{P}}_j^{u,v,\alpha }(x)=\frac{1}{\mathsf{\Gamma }(\alpha )}{\int }_{-1}^x{(x-t)}^{\alpha -1}{P}_j^{u,v}(t)\mathrm{d}t$ can be explicitly calculated by the recurrence formula

$$\left\{\begin{array}{cc}\hfill {\widehat{P}}_0^{u,v,\alpha }(x)& =\frac{{(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)},\hfill \\ \hfill {\widehat{P}}_1^{u,v,\alpha }(x)& =\frac{u+v+2}{2}\left(\frac{x{(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)}-\frac{\alpha {(x+1)}^{\alpha +1}}{\mathsf{\Gamma }(\alpha +2)}\right)+\frac{u-v}{2}{\widehat{P}}_0^{u,v,\alpha }(x),\hfill \\ \hfill {\widehat{P}}_{j+1}^{u,v,\alpha }(x)& =\frac{A_j^{u,v}x-B_j^{u,v}-\alpha A_j^{u,v}{\widehat{B}}_j^{u,v}}{1+\alpha A_j^{u,v}{\widehat{C}}_j^{u,v}}{\widehat{P}}_j^{u,v,\alpha }(x)\hfill \\ & +\frac{\alpha \left({\widehat{A}}_j^{u,v}{P}_{j-1}^{u,v}(-1)+{\widehat{B}}_j^{u,v}{P}_j^{u,v}(-1)+{\widehat{C}}_j^{u,v}{P}_{j+1}^{u,v}(-1)\right)}{\mathsf{\Gamma }(\alpha +1)(1+\alpha A_j^{u,v}{\widehat{C}}_j^{u,v})}{A}_j^{u,v}{(x+1)}^{\alpha }\hfill \\ & -\frac{C_j^{u,v}+\alpha A_j^{u,v}{\widehat{A}}_j^{u,v}}{1+\alpha A_j^{u,v}{\widehat{C}}_j^{u,v}}{\widehat{P}}_{j-1}^{u,v,\alpha }(x),j\ge 1,\hfill \end{array}\right.$$

(36)

which follows from the three-term recurrence relation of the Jacobi polynomials. Here the coefficients are given by

$$\left\{\begin{array}{cc}\hfill A_j^{u,v}=& \frac{(2j+u+v+1)(2j+u+v+2)}{2(j+1)(j+u+v+1)},\hfill \\ \hfill B_j^{u,v}=& \frac{(v^2-u^2)(2j+u+v+1)}{2(j+1)(j+u+v+1)(2j+u+v)},\hfill \\ \hfill C_j^{u,v}=& \frac{(j+u)(j+v)(2j+u+v+2)}{(j+1)(j+u+v+1)(2j+u+v)},\hfill \end{array}\right.$$

(37)

and

$$\left\{\begin{array}{cc}\hfill {\widehat{A}}_j^{u,v}=& \frac{-2(j+u)(j+v)}{(j+u+v)(2j+u+v)(2j+u+v+1)},\hfill \\ \hfill {\widehat{B}}_j^{u,v}=& \frac{2(u-v)}{(2j+u+v)(2j+u+v+2)},\hfill \\ \hfill {\widehat{C}}_j^{u,v}=& \frac{2(j+u+v+1)}{(2j+u+v+1)(2j+u+v+2)},\hfill \end{array}\right.$$

(38)

and ${\widehat{A}}_j^{u,v}=0$ if $j=1$. For $f(x)$ defined on an arbitrary interval $[a,b]$, it follows from the affine transformation $\widehat{x}=\frac{2x-a-b}{b-a}\in [-1,1]$ that

$${}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\approx {\left(\frac{b-a}{2}\right)}^{\alpha }{}_{RL}{\mathrm{D}}_{a,\widehat{x}}^{-\alpha }{p}_N(\widehat{x})={\left(\frac{b-a}{2}\right)}^{\alpha }{\displaystyle \sum _{j=0}^N}{\tilde{p}}_j^{u,v}{\widehat{P}}_j^{u,v,\alpha }(\widehat{x}).$$

(39)

Let $u=v=0$. Then the Jacobi polynomials ${\{P_j^{u,v}\}}_{j=0}^N$ reduce to the Legendre polynomials ${\{L_j(x)\}}_{j=0}^N$. Consequently, numerical scheme (35) becomes the spectral approximation based on Legendre polynomials

$${}_{RL}{\mathrm{D}}_{-1,x}^{-\alpha }f(x)\approx {\displaystyle \sum _{j=0}^N}{\tilde{l}}_j{\widehat{L}}_j^{\alpha }(x),x\in [-1,1].$$

(40)

Here the coefficients are given by

$${\tilde{l}}_j=\frac{1}{{\overline{\gamma }}_j}{\displaystyle \sum _{k=0}^N}f(x_k)L_j(x_k){\omega }_k,$$

(41)

with ${\overline{\gamma }}_j=\frac{2}{2j+1}$ for $0\le j\le N-1$, ${\overline{\gamma }}_N=\frac{2}{N}$, and ${\{{\omega }_k\}}_{k=0}^N$ being the corresponding quadrature weights. Recurrence formula for ${\widehat{L}}_j^{\alpha }(x)$ is in the form

$$\left\{\begin{array}{cc}\hfill {\widehat{L}}_0^{\alpha }(x)=& \frac{{(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)},\hfill \\ \hfill {\widehat{L}}_1^{\alpha }(x)=& \frac{x{(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)}-\frac{\alpha {(x+1)}^{\alpha +1}}{\mathsf{\Gamma }(\alpha +2)},\hfill \\ \hfill {\widehat{L}}_{j+1}^{\alpha }(x)=& \frac{(2j+1)x{\widehat{L}}_j^{\alpha }(x)-(j-\alpha ){\widehat{L}}_{j-1}^{\alpha }(x)}{j+1+\alpha },j\ge 1.\hfill \end{array}\right.$$

(42)

Correspondingly, for the case with arbitrary interval $[a,b]$, we also have

$${}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\approx {\left(\frac{b-a}{2}\right)}^{\alpha }{\displaystyle \sum _{j=0}^n}{\tilde{l}}_j{\widehat{L}}_j^{\alpha }(\widehat{x}),\widehat{x}=\frac{2x-a-b}{b-a}.$$

(43)

Let $u=v=-\frac{1}{2}$ in Equation (39). We obtain the spectral approximation based on Chebyshev polynomials

$${}_{RL}{\mathrm{D}}_{a,x}^{-\alpha }f(x)\approx {\left(\frac{b-a}{2}\right)}^{\alpha }\sum _{j=0}^N{\tilde{t}}_j{\widehat{T}}_j^{\alpha }(\widehat{x}),\widehat{x}=\frac{2x-a-b}{b-a}\in [-1,1],$$

(44)

which follows from the relation $P_j^{-\frac{1}{2},-\frac{1}{2}}(x)=\frac{\mathsf{\Gamma }\left(j+\frac{1}{2}\right)}{j!\sqrt{\pi }}{T}_j(x)$ with ${\{T_j(x)\}}_{j=0}^N$ being the Chebyshev polynomials. Here ${\widehat{T}}_j^{\alpha }(x)=\frac{1}{\mathsf{\Gamma }(\alpha )}{\int }_{-1}^x{(x-s)}^{\alpha -1}{T}_j(s)\mathrm{d}s$ can be computed by the recurrence formula

$$\left\{\begin{array}{cc}\hfill {\widehat{T}}_0^{\alpha }(x)=& \frac{{(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)},\hfill \\ \hfill {\widehat{T}}_1^{\alpha }(x)=& \frac{x{(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)}-\frac{\alpha {(x+1)}^{\alpha +1}}{\mathsf{\Gamma }(\alpha +2)},\hfill \\ \hfill {\widehat{T}}_2^{\alpha }(x)=& \frac{4x}{2+\alpha }{\widehat{T}}_1^{\alpha }(x)-\frac{2}{2+\alpha }{\widehat{T}}_0^{\alpha }(x)+\frac{\alpha {(x+1)}^{\alpha }}{(2+\alpha )\mathsf{\Gamma }(\alpha +1)},\hfill \\ \hfill {\widehat{T}}_{j+1}^{\alpha }(x)=& \frac{2(j+1)x}{j+1+\alpha }{\widehat{T}}_j^{\alpha }(x)-\frac{(j+1)(j-1-\alpha )}{(j+1+\alpha )(j-1)}{\widehat{T}}_{j-1}^{\alpha }(x)\hfill \\ & +\frac{2{(-1)}^j\alpha {(x+1)}^{\alpha }}{\mathsf{\Gamma }(\alpha +1)(j+1+\alpha )(j-1)},j\ge 2.\hfill \end{array}\right.$$

(45)

The coefficients ${\tilde{t}}_j(0\le j\le N)$ are determined by

$${\tilde{t}}_j=\frac{1}{{\sigma }_j}\sum _{k=0}^{N}f(x_k)T_j(x_k){\omega }_k,$$

(46)

with

$${\sigma }_j=\left\{\begin{array}{cc}{\gamma }_j^{-\frac{1}{2},-\frac{1}{2}},\hfill & j=0,1,\dots ,N-1,\hfill \\ 2{\gamma }_N^{-\frac{1}{2},-\frac{1}{2}},\hfill & j=N,\hfill \end{array}\right.$$

(47)

and ${\{{\omega }_k\}}_{k=0}^N$ 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.

Example 1.

Let$f(x)=x^{\mu }$,$x\in [0,1]$. Apply scheme (39) to evaluating its fractional integral.

Table 1. The absolute errors for Example 1.

$u=v=0,\mu =3.5$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
10	4.57 × ${10}^{-8}$	3.57 × ${10}^{-8}$	1.78 × ${10}^{-8}$	5.18 × ${10}^{-9}$	1.67 × ${10}^{-9}$	6.04 × ${10}^{-10}$
20	2.89 × ${10}^{-10}$	1.52 × ${10}^{-10}$	5.37 × ${10}^{-11}$	9.88 × ${10}^{-12}$	2.54 × ${10}^{-12}$	1.31 × ${10}^{-12}$
40	1.82 × ${10}^{-12}$	6.36 × ${10}^{-13}$	1.52 × ${10}^{-13}$	1.74 × ${10}^{-14}$	3.68 × ${10}^{-15}$	2.77 × ${10}^{-15}$
80	1.12 × ${10}^{-14}$	2.59 × ${10}^{-15}$	4.11 × ${10}^{-16}$	1.67 × ${10}^{-16}$	1.67 × ${10}^{-16}$	1.18 × ${10}^{-16}$
$u=v=-\frac{1}{2},\mu =3.5$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
10	5.49 × ${10}^{-8}$	4.59 × ${10}^{-8}$	2.54 × ${10}^{-8}$	8.33 × ${10}^{-9}$	3.09 × ${10}^{-9}$	1.67 × ${10}^{-9}$
20	3.08 × ${10}^{-10}$	1.96 × ${10}^{-10}$	7.62 × ${10}^{-11}$	1.70 × ${10}^{-11}$	4.97 × ${10}^{-12}$	2.93 × ${10}^{-12}$
40	1.81 × ${10}^{-12}$	7.79 × ${10}^{-13}$	2.14 × ${10}^{-13}$	3.23 × ${10}^{-14}$	8.09 × ${10}^{-15}$	5.61 × ${10}^{-15}$
80	1.06 × ${10}^{-14}$	3.05 × ${10}^{-15}$	5.66 × ${10}^{-16}$	3.33 × ${10}^{-16}$	1.80 × ${10}^{-16}$	1.73 × ${10}^{-16}$

shows the absolute maximum errors at the Jacobi–Gauss–Lobatto points. The spectral accuracy is visible.

Example 2.

Let$f(x)=sinx$,$x\in [0,1]$. Utilize scheme (39) to evaluate its fractional integral.

Table 2. The absolute errors for Example 2.

$u=v=0$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
4	4.46 × ${10}^{-6}$	5.40 × ${10}^{-6}$	3.88 × ${10}^{-6}$	1.71 × ${10}^{-6}$	6.94 × ${10}^{-7}$	2.90 × ${10}^{-7}$
8	4.79 × ${10}^{-12}$	4.72 × ${10}^{-12}$	2.73 × ${10}^{-12}$	8.94 × ${10}^{-13}$	2.88 × ${10}^{-13}$	6.96 × ${10}^{-14}$
16	6.66 × ${10}^{-16}$	2.22 × ${10}^{-16}$	2.22 × ${10}^{-16}$	1.67 × ${10}^{-16}$	1.67 × ${10}^{-16}$	8.33 × ${10}^{-17}$
$u=v=-\frac{1}{2}$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
4	6.08 × ${10}^{-6}$	7.91 × ${10}^{-6}$	6.26 × ${10}^{-6}$	3.41 × ${10}^{-6}$	1.87 × ${10}^{-6}$	1.27 × ${10}^{-6}$
8	6.58 × ${10}^{-12}$	6.63 × ${10}^{-12}$	3.96 × ${10}^{-12}$	1.38 × ${10}^{-12}$	4.97 × ${10}^{-13}$	1.36 × ${10}^{-13}$
16	3.33 × ${10}^{-16}$	3.33 × ${10}^{-16}$	2.22 × ${10}^{-16}$	5.55 × ${10}^{-17}$	2.22 × ${10}^{-16}$	5.55 × ${10}^{-17}$

displays the absolute maximum errors of the spectral approximations to the fractional integral with$u=v=0$and$u=v=-\frac{1}{2}$. We can observe that satisfactory results are obtained as well.

2.2. Fractional Linear Multistep Method

In Ref. [45], the convolution quadrature

$$I_h^{\alpha }f(x)=h^{\alpha }{\displaystyle \sum _{j=0}^n}{\omega }_{\ell ,n-j}f(jh)+h^{\alpha }{\displaystyle \sum _{j=0}^s}{w}_{n,j}f(jh),x=nh,$$

(48)

is utilized to evaluate fractional integrals (with $\alpha >0$) and derivatives (with $\alpha <0$). Here the convolution quadrature weights ${\omega }_{\ell ,j}(j\ge 0)$ and the starting quadrature weights $w_{n,j}(n\ge 0,j=0,\dots ,s;s\mathrm{fixed})$ do not depend on h.

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 $\alpha$, 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}{\mathrm{D}}_{0.x}^{-\alpha }f(x)$ with accuracy order $\mathcal{O}(h^{\ell })$, i.e.,

$${}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }f(x)=h^{\alpha }{\displaystyle \sum _{j=0}^n}{\omega }_{\ell ,n-j}f(jh)+h^{\alpha }{\displaystyle \sum _{j=0}^s}{w}_{n,j}f(jh)+\mathcal{O}(h^{\ell }),x=a+nh,$$

(49)

where s is a fixed integer with $s\le n$. Here the convolution coefficients ${\omega }_{\ell ,j}$ are respectively those of the Taylor series expansions of the corresponding generating functions

$$W_{\ell }^{(\alpha )}(z)={\left({\displaystyle \sum _{k=1}^{\ell }}\frac{1}{k}{(1-z)}^k\right)}^{\alpha }={\displaystyle \sum _{j=0}^{\infty }}{\omega }_{\ell ,j}{z}^j,\left|z\right|<1,$$

(50)

with ℓ being the order of consistency. Technically all the coefficients ${\omega }_{\ell ,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

$${}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }f(x)=h^{\alpha }{\displaystyle \sum _{j=0}^n}{\omega }_{\ell ,n-j}f(jh)+\mathcal{O}(h^{\ell }),x=a+nh,$$

(51)

when $f(0)=0$. For $s\ne 0$, the coefficients $w_{n,j}$ can be constructed such that Equation (49) exactly holds for power functions. Therefore, we recover

$${\displaystyle \sum _{j=0}^s}{w}_{n,j}{j}^q=\frac{\mathsf{\Gamma }(q+1)}{\mathsf{\Gamma }(q+\alpha +1)}{n}^{q+\alpha }-{\displaystyle \sum _{j=0}^n}{\omega }_{\ell ,n-j}{j}^q,q=0,\dots ,s,$$

(52)

and it makes sense to choose $s=\ell -1$.

In the case with the lower terminal $a\ne 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

$${\mathcal{W}}_2^{\alpha }(z)={\left(\frac{1}{2}\frac{1+z}{1-z}\right)}^{\alpha }$$

(53)

as the generating function, the fractional trapezoidal rule with second order accuracy for$\alpha \ge 0$is obtained.

Let${\gamma }_i(i=0,1,2,\dots )$denote the coefficients of

$${\displaystyle \sum _{i=0}^{\infty }}{\gamma }_i{(1-z)}^i={\left(\frac{lnz}{z-1}\right)}^{-\alpha },$$

(54)

and set

$${\tilde{\mathcal{W}}}_{\ell }^{\alpha }={(1-z)}^{-\alpha }\left[{\gamma }_0+{\gamma }_1(1-z)+\cdots +{\gamma }_{\ell -1}{(1-z)}^{\ell -1}\right],\ell =1,2,\dots$$

(55)

Then we obtain the generating function for the coefficients of the generalized Newton-Gregory formulae, which is convergent of order ℓ. Direct calculation gives${\gamma }_0=1$and${\gamma }_1=-\frac{\alpha }{2}$. Then the corresponding generating function for the second order generalized Newton–Gregory formula is given by

$${\tilde{\mathcal{W}}}_2^{\alpha }={(1-z)}^{-\alpha }\left[1-\frac{\alpha }{2}(1-z)\right].$$

(56)

More details for generating functions can be found in Refs. [45,49,50,51].

2.3. 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

$$\mathsf{\Gamma }(\alpha )={\int }_0^{\infty }{e}^{-z}{z}^{\alpha -1}\mathrm{d}z,$$

(57)

and

$$\mathsf{\Gamma }(1-\alpha )\mathsf{\Gamma }(\alpha )=\frac{\pi }{sin(\pi \alpha )},$$

(58)

the fractional integral with $0<\alpha <1$ can be rewritten as [52]

$${}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }f(x)=\frac{sin(\pi \alpha )}{\pi }{\int }_0^x\left[{\int }_0^{\infty }{e}^{-z}{\left(\frac{z}{x-t}\right)}^{1-\alpha }\frac{\mathrm{d}z}{z}\right]f(t)\mathrm{d}t.$$

(59)

Define the variable transformation $z=(x-t){\omega }^2,\omega \ge 0$. The Fubini’s Theorem yields

$${}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }f(x)=\frac{2sin(\pi \alpha )}{\pi }{\int }_0^{\infty }{\omega }^{1-2\alpha }\left({\int }_0^x{e}^{-(x-t){\omega }^2}f(t)\mathrm{d}t\right)\mathrm{d}\omega .$$

(60)

Introducing the auxiliary function

$$\varphi (\omega ,x)=\frac{2sin(\pi \alpha )}{\pi }{\omega }^{1-2\alpha }{\int }_0^x{e}^{-(x-t){\omega }^2}f(t)\mathrm{d}t,$$

(61)

we have

$${}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }f(x)={\int }_0^{\infty }\varphi (\omega ,x)\mathrm{d}\omega ,0<\alpha <1.$$

(62)

It follows from the definition of $\varphi (\omega ,x)$ that the auxiliary function satisfies

$$\left\{\begin{array}{c}\frac{\partial }{\partial x}\varphi (\omega ,x)=\frac{2sin(\pi \alpha )}{\pi }{\omega }^{1-2\alpha }f(x)-{\omega }^2\varphi (\omega ,x),\hfill \\ \varphi (\omega ,0)=0.\hfill \end{array}\right.$$

(63)

In this case, evaluating Riemann–Liouville integral ${}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }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^{\alpha -1}=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\infty }{e}^{-zx}{z}^{-\alpha }\mathrm{d}z.$$

(64)

Consequently, the Fubini’s Theorem gives

$$\begin{array}{cc}\hfill {}_{RL}{\mathrm{D}}_{0,x}^{-\alpha }f(x)=& \frac{1}{\mathsf{\Gamma }(\alpha )}\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\infty }\left({\int }_0^x{e}^{-z(x-t)}f(t)\mathrm{d}t\right)\frac{\mathrm{d}z}{z^{\alpha }}\hfill \\ & =\frac{sin(\pi \alpha )}{\pi }{\int }_0^{\infty }g(z,x)z^{-\alpha }\mathrm{d}z,\hfill \end{array}$$

(65)

where $g(z,x)$ is defined as

$$g(z,x)={\int }_0^x{e}^{-z(x-t)}f(t)\mathrm{d}t.$$

(66)

In order to generate nonreflecting boundary conditions [55] and accelerate convolutions with the heat kernel [56], literatures such as Ref. [57] usually recognize

$$g(z,x)=e^{-z\Delta x}g(z,x-\Delta x)+\mathsf{\Psi }(z,x,\Delta x),$$

(67)

where

$$\mathsf{\Psi }(z,x,\Delta x)={\int }_{x-\Delta x}^x{e}^{-z(x-t)}f(t)\mathrm{d}t.$$

(68)

Alternatively, other literatures have regarded $g(z,x)$ as the solution of a first order ordinary differential (ODE) equation [52,53,58,59],

$$\frac{\mathrm{d}g(z,x)}{\mathrm{d}x}=-zg(z,x)+f(x),g(z,0)=0.$$

(69)

Any approximate method for ODEs can be used to obtain $g(z,x),x=\Delta x,2\Delta x,\dots$, 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.

3. 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.

3.1. 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.

3.1.1. L1, L2, and L2C Methods

The well-known L1 method was originally introduced in Ref. [62] to evaluate Riemann–Liouville derivative with $0<\alpha <1$, which equivalently reads as

$${}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)=\frac{{(x-a)}^{-\alpha }}{\mathsf{\Gamma }(1-\alpha )}f(a)+\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_a^x{(x-t)}^{-\alpha }{f}^{\prime }(t)\mathrm{d}t.$$

(70)

Note that the second term on the right-hand side happens to be Caputo derivative with $0<\alpha <1$. That is the reason why we introduce the L1 method when considering numerical approximations to the Caputo derivative.

Let $f(x)\in C^2[a,b]$. On the setting of uniform grids ${\{x_k\}}_{k=0}^N$, utilizing the constant $\frac{f(x_{k+1})-f(x_k)}{h}$ to approximate $f^{\prime }(x)$ on each interval $[x_k,x_{k+1}]$ yields the following L1 method on uniform grids for Caputo derivative [62]

$$\begin{array}{cc}\hfill {\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}& =\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{-\alpha }\frac{f(x_{k+1})-f(x_k)}{h}\mathrm{d}t+\mathcal{O}(h^{2-\alpha })\hfill \\ & ={\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}\left[f(x_{k+1})-f(x_k)\right]+\mathcal{O}(h^{2-\alpha }),0<\alpha <1,1\le j\le N.\hfill \end{array}$$

(71)

Here the coefficients are given by

$$b_k=\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}\left[{(k+1)}^{1-\alpha }-k^{1-\alpha }\right],0\le k\le j-1.$$

(72)

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.

Theorem 2.

Let$0<\alpha <1$and$f(x)\in C^2[a,b]$. Denote by

$${\left[{\delta }_x^{\alpha }f(x)\right]}_{x=x_j}={\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}\left[f(x_{k+1})-f(x_k)\right],1\le j\le N.$$

(73)

Then it holds that

$$\left|{\left[{\delta }_x^{\alpha }f(x)\right]}_{x=x_j}-{\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}\right|\le C{h}^{2-\alpha },$$

(74)

where C is a positive constant given by

$$C=\frac{1}{\mathsf{\Gamma }(2-\alpha )}\left[\frac{1-\alpha }{12}+\frac{2^{2-\alpha }}{2-\alpha }-(2^{-\alpha }+1)\right]\underset{x_0\le x\le x_j}{max}\left|f^{″}(x)\right|.$$

(75)

Proof. 

Denote

$$A={\displaystyle \sum _{k=0}^{j-1}}\frac{f(x_{k+1})-f(x_k)}{h}{\int }_{x_k}^{x_{k+1}}\frac{\mathrm{d}t}{{(x_j-t)}^{\alpha }}-{\int }_{x_0}^{x_j}\frac{f^{\prime }(t)}{{(x_j-t)}^{\alpha }}\mathrm{d}t.$$

(76)

Then it immediately follows that

$$\frac{\left|A\right|}{\mathsf{\Gamma }(1-\alpha )}=|{\displaystyle \sum _{k=0}^{j-1}}{b}_{j-k-1}\left[f(x_{k+1})-f(x_k)\right]-{\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}|.$$

(77)

Using the Taylor expansion with integral remainder, we have for $t\in [x_k,x_{k+1}]$,

$$f^{\prime }(t)-\frac{f(x_{k+1})-f(x_k)}{h}=\frac{1}{h}\left[{\int }_{x_k}^t{f}^{″}(s)(s-x_k)\mathrm{d}s-{\int }_t^{x_{k+1}}{f}^{″}(s)(x_{k+1}-s)\mathrm{d}s\right],$$

(78)

which yields

$$\begin{array}{cc}\hfill A& ={\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}\left[f^{\prime }(t)-\frac{f(x_{k+1})-f(x_k)}{h}\right]{(x_j-t)}^{-\alpha }\mathrm{d}t\hfill \\ & =\frac{1}{h}{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}\left[{\int }_{x_k}^t{f}^{″}(s)(s-x_k)\mathrm{d}s-{\int }_t^{x_{k+1}}{f}^{″}(s)(x_{k+1}-s)\mathrm{d}s\right]\frac{\mathrm{d}t}{{(x_j-t)}^{\alpha }}.\hfill \end{array}$$

(79)

Exchanging the order of integration gives

$$\begin{array}{ccc}A\hfill & =\hfill & \frac{1}{h}{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}\left[{\int }_{x_k}^t{f}^{″}(s)(s-x_k)\mathrm{d}s-{\int }_t^{x_{k+1}}{f}^{″}(s)(x_{k+1}-s)\mathrm{d}s\right]\frac{\mathrm{d}t}{{(x_j-t)}^{\alpha }}\hfill \\ & =\hfill & \frac{1}{h}{\displaystyle \sum _{k=0}^{j-1}}\left[{\int }_{x_k}^{x_{k+1}}{f}^{″}(s)(s-x_k){\int }_s^{x_{k+1}}{(x_j-t)}^{-\alpha }\mathrm{d}t\mathrm{d}s-{\int }_{x_k}^{x_{k+1}}{f}^{″}(s)(x_{k+1}-s){\int }_{x_k}^s{(x_j-t)}^{-\alpha }\mathrm{d}t\mathrm{d}s\right]\hfill \\ & =\hfill & \frac{1}{1-\alpha }{\displaystyle \sum _{k=0}^{j-1}}[{\int }_{x_k}^{x_{k+1}}{f}^{″}(s)\frac{s-x_k}{h}\left[{(x_j-s)}^{1-\alpha }-{(x_j-x_{k+1})}^{1-\alpha }\right]\mathrm{d}s\hfill \\ & & -{\int }_{x_k}^{x_{k+1}}{f}^{″}(s)\frac{x_{k+1}-s}{h}\left[{(x_j-x_k)}^{1-\alpha }\right.-{(x_j-s)}^{1-\alpha }]\mathrm{d}s]\hfill \\ & =\hfill & \frac{1}{1-\alpha }{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}{f}^{″}(s)\{{(x_j-s)}^{1-\alpha }-\left[\frac{s-x_k}{h{(x_j-x_{k+1})}^{\alpha -1}}+\frac{x_{k+1}-s}{h{(x_j-x_k)}^{\alpha -1}}\right]\}\mathrm{d}s.\hfill \end{array}$$

(80)

In the following, we show that when $0<\alpha <1$,

$${\int }_{x_k}^{x_{k+1}}\left\{{(x_j-s)}^{1-\alpha }-\left[\frac{s-x_k}{h{(x_j-x_{k+1})}^{\alpha -1}}+\frac{x_{k+1}-s}{h{(x_j-x_k)}^{\alpha -1}}\right]\right\}\mathrm{d}s\ge 0$$

(81)

for $0\le k\le j-1$, and

$${\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}\{{(x_j-s)}^{1-\alpha }-\left[\frac{s-x_k}{h{(x_j-x_{k+1})}^{\alpha -1}}+\frac{x_{k+1}-s}{h{(x_j-x_k)}^{\alpha -1}}\right]\}\mathrm{d}s<+\infty .$$

(82)

Denote $g(s)={(x_j-s)}^{1-\alpha }$. Then it holds for any $s\in (x_k,x_{k+1})$ that

$$g(s)-\left[\frac{s-x_k}{h}g(x_{k+1})+\frac{x_{k+1}-s}{h}g(x_k)\right]=\frac{(1-\alpha )(-\alpha )(s-x_k)(s-x_{k+1})}{2{(x_j-{\xi }_k)}^{\alpha +1}}\ge 0,$$

(83)

with certain ${\xi }_k\in (x_k,x_{k+1})$. As a result, inequality (81) holds. For the inequality (82), one has

$$\begin{array}{cc}& {\displaystyle \sum _{k=0}^{j-3}}{\int }_{x_k}^{x_{k+1}}\{g(s)-\left[\frac{s-x_k}{h}g(x_{k+1})+\frac{x_{k+1}-s}{h}g(x_k)\right]\}\mathrm{d}s\hfill \\ \hfill =& {\displaystyle \sum _{k=0}^{j-3}}{\int }_{x_k}^{x_{k+1}}\frac{\alpha (1-\alpha )(s-x_k)(x_{k+1}-s)}{2{(x_j-{\xi }_k)}^{\alpha +1}}\mathrm{d}s\hfill \\ \hfill \le & \frac{\alpha (1-\alpha )}{2}{\displaystyle \sum _{k=0}^{j-3}}{(x_j-x_{k+1})}^{-\alpha -1}{\int }_{x_k}^{x_{k+1}}(s-x_k)(x_{k+1}-s)\mathrm{d}s\hfill \\ \hfill \le & \frac{h^2}{12}\alpha (1-\alpha ){\displaystyle \sum _{k=0}^{j-3}}{\int }_{x_{k+1}}^{x_{k+2}}{(x_j-s)}^{-\alpha -1}\mathrm{d}s\hfill \\ \hfill \le & \frac{1-\alpha }{12}{h}^{2-\alpha },\hfill \end{array}$$

(84)

and

$$\begin{array}{cc}& {\displaystyle \sum _{k=j-2}^{j-1}}{\int }_{x_k}^{x_{k+1}}\{g(s)-\left[\frac{s-x_k}{h}g(x_{k+1})+\frac{x_{k+1}-s}{h}g(x_k)\right]\}\mathrm{d}s\hfill \\ \hfill =& {\int }_{x_{j-2}}^{x_j}g(s)\mathrm{d}s-\left[\frac{g(x_{j-2})}{2}+g(x_{j-1})+\frac{g(x_j)}{2}\right]h\hfill \\ \hfill =& {\int }_{x_{j-2}}^{x_j}g(s)\mathrm{d}s-\left[\frac{g(x_{j-2})}{2}+g(x_{j-1})\right]h\hfill \\ \hfill =& {\int }_{x_{j-2}}^{x_j}{(x_j-s)}^{1-\alpha }\mathrm{d}s-\left[\frac{{(x_j-x_{j-2})}^{1-\alpha }}{2}+{(x_j-x_{j-1})}^{1-\alpha }\right]h\hfill \\ \hfill =& \left[\frac{2^{2-\alpha }}{2-\alpha }-(2^{-\alpha }+1)\right]h^{2-\alpha }.\hfill \end{array}$$

(85)

The above two equalities yield that Equation (82) holds.

Combining the above analysis, one has

$$\begin{array}{cc}\hfill 0& \le {\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}\{{(x_j-s)}^{1-\alpha }-\left[\frac{s-x_k}{h{(x_j-x_{k+1})}^{\alpha -1}}+\frac{x_{k+1}-s}{h{(x_j-x_k)}^{\alpha -1}}\right]\}\mathrm{d}s\hfill \\ & \le \left[\frac{1-\alpha }{12}+\frac{2^{2-\alpha }}{2-\alpha }-(2^{-\alpha }+1)\right]h^{2-\alpha }.\hfill \end{array}$$

(86)

Inserting the above estimate into Equation (80) gives

$$\left|A\right|\le \frac{1}{1-\alpha }\left[\frac{1-\alpha }{12}+\frac{2^{2-\alpha }}{2-\alpha }-(2^{-\alpha }+1)\right]\underset{x_0\le x\le x_j}{max}\left|f^{″}(x)\right|h^{2-\alpha }.$$

(87)

All this ends the proof. □

Remark 4.

The idea of proving Theorem 2 is borrowed from Ref. [70] where the case with$\alpha \in (1,2)$was considered. Such an estimate was also considered in [71].

Let $\{{\tilde{x}}_i\}$ be any division of $[a,b]$ with $a={\tilde{x}}_0<{\tilde{x}}_1<\dots <{\tilde{x}}_{n-1}<{\tilde{x}}_N=b$. Then the classical L1 method is generalized into the L1 method on nonuniform grids for Caputo derivative [72]

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x={\tilde{x}}_j}={\displaystyle \sum _{k=0}^{j-1}}{b}_{k+1}^j\left[f({\tilde{x}}_{k+1})-f({\tilde{x}}_k)\right]+\mathcal{O}({\tilde{h}}_{max}^{2-\alpha }),$$

(88)

provided that $\frac{\underset{0\le k\le j-1}{max}{\tilde{h}}_k}{\underset{0\le k\le j-1}{min}{\tilde{h}}_k}\le C$ with C being a positive constant. Here ${\tilde{h}}_k={\tilde{x}}_{k+1}-{\tilde{x}}_k$, ${\tilde{h}}_{max}=\underset{0\le k\le j-1}{max}{\tilde{h}}_k$, and the coefficients are given by

$$b_{k+1}^j=\frac{1}{\mathsf{\Gamma }(2-\alpha ){\tilde{h}}_k}\left[{({\tilde{x}}_j-{\tilde{x}}_k)}^{1-\alpha }-{({\tilde{x}}_j-{\tilde{x}}_{k+1})}^{1-\alpha }\right].$$

(89)

In the special case of nonuniform grids with ${\tilde{x}}_0=x_0$, ${\tilde{x}}_j=x_{j-\frac{1}{2}}=\frac{x_j+x_{j-1}}{2}$, $j=1,2,\dots$, scheme (88) is reduced to

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x={\tilde{x}}_{j+1}}=b_{0}f({\tilde{x}}_{j+1})-{\displaystyle \sum _{k=1}^j}(b_{j-k}-b_{j-k+1})f({\tilde{x}}_k)-B_{j}f({\tilde{x}}_0)+\mathcal{O}(h^{2-\alpha }).$$

(90)

Here the coefficients are given by

$$\left\{\begin{array}{c}{b}_k=\frac{{(k+1)}^{1-\alpha }-k^{1-\alpha }}{\mathsf{\Gamma }(2-\alpha )h^{\alpha }},0\le k\le j,\hfill \\ B_j=\frac{2{\left(j+\frac{1}{2}\right)}^{1-\alpha }-2j^{1-\alpha }}{\mathsf{\Gamma }(2-\alpha )h^{\alpha }},0\le j\le N.\hfill \end{array}\right.$$

(91)

Replacing $f({\tilde{x}}_k)=f(\frac{x_{k-1}+x_k}{2})$ with $\frac{f(x_k)+f(x_{k-1})}{2}$ yields the following modified L1 method for Caputo derivative [38]

$$\begin{array}{ccc}\hfill {\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_{j+\frac{1}{2}}}& =& -\frac{1}{2}{\displaystyle \sum _{k=1}^j}(b_{j-k}-b_{j-k+1})\left[f(x_{k-1})+f(x_k)\right]\hfill \\ & & +\frac{b_0}{2}\left[f(x_{j+1})+f(x_j)\right]-B_{j}f(x_0)+\mathcal{O}(h^{2-\alpha }).\hfill \end{array}$$

(92)

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<\alpha <2$ and the lower terminal $a=0$, there holds

$${\left[{}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{\mathsf{\Gamma }(2-\alpha )}\sum _{k=0}^{j-1}{\int }_{x_k}^{x_{k+1}}{t}^{1-\alpha }{f}^{″}(x_j-t)\mathrm{d}t.$$

(93)

Suppose that $f(x)\in C^3[a,b]$. Utilizing the central difference scheme

$$\frac{f(x_j-x_{k+1})-2f(x_j-x_k)+f(x_j-x_{k-1})}{h^2},$$

(94)

to approximate $f^{″}(x_j-t)$ on each interval $[x_k,x_{k+1}]$, we have the following L2 method for Caputo derivative [62]

$${\left[{}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{\mathsf{\Gamma }(3-\alpha )h^{\alpha }}{\displaystyle \sum _{k=-1}^j}{W}_{k}f(x_{j-k})+\mathcal{O}(h^{3-\alpha }),$$

(95)

in which

$$\left\{\begin{array}{c}{W}_{-1}=1,W_0=2^{2-\alpha }-3,\hfill \\ W_k=\left[{(k+2)}^{2-\alpha }-3{(k+1)}^{2-\alpha }+3k^{2-\alpha }-{(k-1)}^{2-\alpha }\right],1\le k\le j-2,\hfill \\ W_{j-1}=-2j^{2-\alpha }+3{(j-1)}^{2-\alpha }-{(j-2)}^{2-\alpha },\hfill \\ W_j=j^{2-\alpha }-{(j-1)}^{2-\alpha }.\hfill \end{array}\right.$$

(96)

In Ref. [78], the integral ${\int }_{x_k}^{x_{k+1}}{t}^{1-\alpha }{f}^{″}(x_j-t)\mathrm{d}t$ was evaluated in a more symmetric form. For $t\in [x_{k-1},x_k]$, if we replace $f^{″}(x_j-t)$ with the difference

$$\frac{f(x_j-x_{k+2})-f(x_j-x_{k+1})+f(x_j-x_{k-1})-f(x_j-x_k)}{2h^2},$$

(97)

then the L2C method for Caputo derivative

$${\left[{}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{2\mathsf{\Gamma }(3-\alpha )h^{\alpha }}{\displaystyle \sum _{k=-1}^{j+1}}{\widehat{W}}_{k}f(x_{j-k})+\mathcal{O}(h^{3-\alpha })$$

(98)

is obtained. Here the coefficients are given by

$$\left\{\begin{array}{c}{\widehat{W}}_{-1}=1,{\widehat{W}}_0=2^{2-\alpha }-2,{\widehat{W}}_1=3^{2-\alpha }-2^{3-\alpha },\hfill \\ {\widehat{W}}_k=\left[{(k+2)}^{2-\alpha }-2{(k+1)}^{2-\alpha }+2{(k-1)}^{2-\alpha }-{(k-2)}^{2-\alpha }\right],2\le k\le j-2,\hfill \\ {\widehat{W}}_{j-1}=2{(j-2)}^{2-\alpha }-j^{2-\alpha }-{(j-3)}^{2-\alpha },\hfill \\ {\widehat{W}}_j=2{(j-1)}^{2-\alpha }-j^{2-\alpha }-{(j-2)}^{2-\alpha },\hfill \\ {\widehat{W}}_{j+1}=j^{2-\alpha }-{(j-1)}^{2-\alpha }.\hfill \end{array}\right.$$

(99)

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^{\prime }(0)=0$ is met. For the case with lower terminal $a\ne 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$\alpha =1$. If$\alpha =2$, the L2 method reduces to the central difference method for the second order derivative and the L2C method reduces to

$$\frac{d^{2}f(x_k)}{d{x}^2}\approx \frac{f(x_{k+2})-f(x_k)+f(x_{k-1})-f(x_{k+1})}{2h^2}$$

(100)

with the first order accuracy. As a matter of fact, the error bound for the L2 method is$\mathcal{O}(h^{3-\alpha })$. Numerical experiments indicate that the L2 method is more accurate than the L2C method for$1.5<\alpha <2$, while the opposite result appears when$1<\alpha <1.5$. And these two methods behave in almost the same way near$\alpha =1.5$[78].

3.1.2. 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.

(I) $(3-\alpha )$-th order approximations

Let $f(x)\in C^3[a,b]$. For $0\le k\le j$ and $0<x-x_k<h$, it follows from Taylor expansions of $f(x_{k+1})$, $f(x_k)$, and $f(x_{k-1})$ at x that

$$\begin{array}{ccc}\hfill f^{\prime }(x)& =& \frac{f(x_{k+1})-f(x_{k-1})}{2h}+\frac{f(x_{k+1})-2f(x_k)+f(x_{k-1})}{h^2}(x-x_k)\hfill \\ & & -\frac{f^{(3)}(x_k)}{3!}{h}^2+\frac{f^{(3)}(x_k)}{2!}{(x-x_k)}^2+\mathcal{O}\left({(x-x_k)}^3\right).\hfill \end{array}$$

(101)

In this case, we have the following $(3-\alpha )$-th order approximation [79],

$$\begin{array}{cc}& {\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_{x_0}^{x_j}{(x_j-t)}^{-\alpha }{f}^{\prime }(t)\mathrm{d}t\hfill \\ =& \frac{h^{-\alpha }}{\mathsf{\Gamma }(3-\alpha )}{\displaystyle \sum _{k=0}^{j-1}}\left\{{\omega }_{1,j-k}\left[f(x_{k+1})-f(x_{k-1})\right]+{\omega }_{2,j-k}\left[f(x_{k+1})-2f(x_k)+f(x_{k-1})\right]\right\}+R^j,\hfill \end{array}$$

(102)

where $0<\alpha <1$, $R^j$ denotes the truncated error, and the coefficients are given by

$$\left\{\begin{array}{c}{\omega }_{1,j-k}=\frac{2-\alpha }{2}\left[{(j-k)}^{1-\alpha }-{(j-k-1)}^{1-\alpha }\right],\hfill \\ {\omega }_{2,j-k}={(j-k)}^{2-\alpha }-{(j-k-1)}^{2-\alpha }-(2-\alpha ){(j-k-1)}^{1-\alpha },\hfill \end{array}\right.$$

(103)

with $0\le k\le j-1$ and $1\le j\le N$.

Since the above $(3-\alpha )$-th order method is also widely used, we estimate its truncated error in detail.

Theorem 3.

([80]) Let $0<\alpha <1$ and $f(x)\in C^3[a,b]$. For the truncated error $R^j$ of approximation (102), it holds that

$$\left|R^j\right|\le c{h}^{3-\alpha },1\le j\le N,$$

(104)

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

$$\begin{array}{ccc}\hfill R^j& =& -\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{k=0}^{j-1}}{\int }_{x_k}^{x_{k+1}}{(x_j-t)}^{-\alpha }\left\{\frac{1}{2!·2h}\left[{\int }_t^{x_{k+1}}{(x_{k+1}-s)}^2{f}^{(3)}(s)\mathrm{d}s\right.\right.\hfill \\ & & \left.-{\int }_t^{x_{k-1}}{(x_{k-1}-s)}^2{f}^{(3)}(s)\mathrm{d}s\right]\hfill \\ & & +\frac{(t-x_k)}{2h^2}\left[{\int }_t^{x_{k+1}}{(x_{k+1}-s)}^2{f}^{(3)}(s)\mathrm{d}s-2{\int }_t^{x_k}{(x_k-s)}^2{f}^{(3)}(s)\mathrm{d}s\right.\hfill \\ & & \left.\left.+{\int }_t^{x_{k-1}}{(x_{k-1}-s)}^2{f}^{(3)}(s)\mathrm{d}s\right]\right\}\mathrm{d}t.\hfill \end{array}$$

(105)

Interchanging the order of integrations yields

$$\begin{array}{ccc}\hfill R^j& =& \frac{1}{2h\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=0}^{j-1}}\{{\int }_{x_k}^{x_{k+1}}{(x_{k+1}-s)}^2{f}^{(3)}(s)\left[\frac{{(x_j-s)}^{1-\alpha }-{(x_j-x_k)}^{1-\alpha }}{2}\right.\hfill \\ & & \left.+\frac{{(x_j-s)}^{1-\alpha }(s-x_k)}{h}+\frac{{(x_j-s)}^{2-\alpha }-{(x_j-x_k)}^{2-\alpha }}{h(2-\alpha )}\right]\mathrm{d}s\hfill \\ & & +\frac{2}{h}{\int }_{x_k}^{x_{k+1}}{(x_k-s)}^2{f}^{(3)}(s)[{(x_j-x_{k+1})}^{1-\alpha }h-{(x_j-s)}^{1-\alpha }(s-x_k)\hfill \\ & & +\frac{{(x_j-x_{k+1})}^{2-\alpha }-{(x_j-s)}^{2-\alpha }}{2-\alpha }]\mathrm{d}s\hfill \\ & & +{\int }_{x_k}^{x_{k+1}}{(x_{k-1}-s)}^2{f}^{(3)}(s)\left[\frac{{(x_j-x_{k+1})}^{1-\alpha }-{(x_j-s)}^{1-\alpha }}{2}-{(x_j-x_{k+1})}^{1-\alpha }\right.\hfill \\ & & \left.+\frac{{(x_j-s)}^{1-\alpha }(s-x_k)}{h}-\frac{{(x_j-x_{k+1})}^{2-\alpha }-{(x_j-s)}^{2-\alpha }}{h(2-\alpha )}\right]\mathrm{d}s\hfill \\ & & +{\int }_{x_{k-1}}^{x_k}{(x_{k-1}-s)}^2{f}^{(3)}(s)[\frac{{(x_j-x_{k+1})}^{1-\alpha }-{(x_j-x_k)}^{1-\alpha }}{2}\hfill \\ & & -\frac{{(x_j-x_{k+1})}^{2-\alpha }-{(x_j-x_k)}^{2-\alpha }}{h(2-\alpha )}-{(x_j-x_{k+1})}^{1-\alpha }\left]\mathrm{d}s\right\}\hfill \\ & & =\frac{1}{2h\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=0}^{j-1}}{S}_k.\hfill \end{array}$$

(106)

For $k=0,1,\dots ,j-1$, denote

$$\begin{array}{ccc}\hfill B_k& =& {\int }_{x_{k-1}}^{x_k}{(x_{k-1}-s)}^2{f}^{(3)}(s)[\frac{{(x_j-x_{k+1})}^{1-\alpha }-{(x_j-x_k)}^{1-\alpha }}{2}\hfill \\ & & -\frac{{(x_j-x_{k+1})}^{2-\alpha }-{(x_j-x_k)}^{2-\alpha }}{h(2-\alpha )}-{(x_j-x_{k+1})}^{1-\alpha }]\mathrm{d}s,\hfill \end{array}$$

(107)

and

$$A_k=S_k-B_k,$$

(108)

where the expression of $A_k$ can be derived from Equations (106) and (107) so is left out due to lengthiness.

Let $l=j-k,k=0,1,\dots ,j-1$. The affine transformation $s=x_{k-1}+\xi h$ with $\xi \in [0,1]$ yields

$$\begin{array}{c}{B}_{j-l}=h^{4-\alpha }{\int }_0^1{\xi }^2{f}^{(3)}(x_{j-l-1}+\xi h)\left[\frac{{(l-1)}^{1-\alpha }-l^{1-\alpha }}{2}-{(l-1)}^{1-\alpha }-\frac{{(l-1)}^{2-\alpha }-l^{2-\alpha }}{2-\alpha }\right]\mathrm{d}\xi \hfill \\ =h^{4-\alpha }{b}_l{\int }_0^1{\xi }^2{f}^{(3)}(x_{j-l-1}+\xi h)\mathrm{d}\xi ,1\le l\le j.\hfill \end{array}$$

(109)

It is evident that $b_1=\frac{1}{2-\alpha }-\frac{1}{2}$, and for $l\ge 2$,

$$b_l=l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}\left(\frac{1}{2n!}-\frac{1}{(n+1)!}\right)(-\alpha +1)\alpha (\alpha +1)\cdots (\alpha +n-2)\ge 0.$$

(110)

Thus, it holds for $l\ge 2$ that

$$\begin{array}{cc}& \left|B_{j-l}\right|=h^{4-\alpha }\left|b_l{\int }_0^1{\xi }^2{f}^{(3)}(x_{j-l-1}+\xi h)\mathrm{d}\xi \right|\hfill \\ \le & \frac{h^{4-\alpha }}{3}{\displaystyle \underset{x\in [x_{j-l-1},x_{j-l}]}{max}}\left|f^{(3)}(x)\right|l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}\left(\frac{1}{2n!}-\frac{1}{(n+1)!}\right)(1-\alpha )\alpha (\alpha +1)\cdots (\alpha +n-2)\hfill \\ \hfill \le & \frac{h^{4-\alpha }}{3}{\displaystyle \underset{x\in [x_{j-l-1},x_{j-l}]}{max}}\left|f^{(3)}(x)\right|l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}\left(\frac{1}{2}-\frac{1}{n+1}\right)\frac{1-\alpha }{n}\hfill \\ \hfill \le & \frac{h^{4-\alpha }}{3}{\displaystyle \underset{x\in [x_{j-l-1},x_{j-l}]}{max}}\left|f^{(3)}(x)\right|l^{-1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^{n-2}}·\frac{1}{2}·\frac{1-\alpha }{2}\hfill \\ \hfill \le & \frac{h^{4-\alpha }(1-\alpha )}{12}{l}^{-1-\alpha }\frac{l}{l-1}{\displaystyle \underset{x\in [x_{j-l-1},x_{j-l}]}{max}}\left|f^{(3)}(x)\right|\hfill \\ \hfill \le & \frac{h^{4-\alpha }}{6}{\displaystyle \underset{x\in [x_{j-l-1},x_{j-l}]}{max}}\left|f^{(3)}(x)\right|\frac{1-\alpha }{l^{1+\alpha }}.\hfill \end{array}$$

(111)

As a result,

$$\begin{array}{ccc}\hfill \left|{\displaystyle \sum _{k=0}^{j-1}}{B}_k\right|=\left|{\displaystyle \sum _{l=1}^j}{B}_{j-l}\right|& \le & h^{4-\alpha }\left\{\left|b_1{\int }_0^1{\xi }^2{f}^{(3)}(x_{j-2}+\xi h)\mathrm{d}\xi \right|+{\displaystyle \sum _{l=2}^j}\left|B_l\right|\right\}\hfill \\ & \hfill \le & h^{4-\alpha }{\displaystyle \underset{x\in [x_{-1},x_{j-1}]}{max}}\left|f^{(3)}(x)\right|\left[\frac{1}{3}\left(\frac{1}{2-\alpha }-\frac{1}{2}\right)+\frac{1-\alpha }{12}{\displaystyle \sum _{l=2}^j}{l}^{-1-\alpha }\right]\hfill \\ & \hfill \le & C_2{\displaystyle \underset{x\in [x_{-1},x_{j-1}]}{max}}\left|f^{(3)}(x)\right|h^{4-\alpha }\hfill \end{array}$$

(112)

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+\xi h,\xi \in [0,1]$ and $l=j-k,k=0,1,\dots ,j-1$ yield

$$\begin{array}{ccc}\hfill A_{j-l}& =& h^{4-\alpha }{\int }_0^1{f}^{(3)}(x_{j-l}+\xi h)\left\{-\frac{2}{2-\alpha }\left[{(l-1)}^{2-\alpha }-{(l-\xi )}^{2-\alpha }\right]\right.\hfill \\ & & +\frac{{(1-\xi )}^2}{2-\alpha }\left[{(l-1)}^{2-\alpha }-l^{2-\alpha }\right]+2\left[{(l-\xi )}^{1-\alpha }\xi -{(l-1)}^{1-\alpha }\right]\hfill \\ & & +{(1-\xi )}^2{(l-1)}^{1-\alpha }+\frac{{(1-\xi )}^2}{2}\left[{(l-\xi )}^{1-\alpha }-l^{1-\alpha }\right]\hfill \\ & & +\frac{{(\xi +1)}^2}{2}\left[{(l-1)}^{1-\alpha }-{(l-\xi )}^{1-\alpha }\right]\}\mathrm{d}\xi \hfill \\ & \hfill =& h^{4-\alpha }{\int }_0^1{f}^{(3)}(x_{j-l}+\xi h)a_l(\xi )\mathrm{d}\xi .\hfill \end{array}$$

(113)

Rewrite $a_l(\xi )$ in the form

$$a_l(\xi )=l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}{\tilde{a}}_n(\xi )(-\alpha +1)\alpha (\alpha +1)\cdots (\alpha +n-2),$$

(114)

with

$${\tilde{a}}_n(\xi )=\frac{2{\xi }^{n+1}-2+{(1-\xi )}^2}{(n+1)!}+\frac{2-{(1-\xi )}^2-\frac{1}{2}{(1+\xi )}^2}{n!},n\ge 2.$$

(115)

For $n\ge 2$, we have ${\tilde{a}}_n(\xi )\ge 0$ for arbitrary $\xi \in [0,1]$. To see this, recall that

$$\left\{\begin{array}{c}{\tilde{a}}_n^{\prime }(\xi )=\frac{2{\xi }^n}{n!}-\frac{2(1-\xi )}{(n+1)!}+\frac{1-3\xi }{n!},\hfill \\ {\tilde{a}}_n^{″}(\xi )=\frac{2{\xi }^{n-1}}{(n-1)!}+\frac{2}{(n+1)!}-\frac{3}{n!}.\hfill \end{array}\right.$$

(116)

When ${\xi }_0={\left(\frac{1}{2n}+\frac{1}{n+1}\right)}^{\frac{1}{n-1}}\in (0,1)$, there hold

$${\tilde{a}}_n^{″}({\xi }_0)=0,$$

(117)

and

$${\tilde{a}}_n^{″}(\xi )\left\{\begin{array}{c}<0,\xi \in [0,{\xi }_0),\hfill \\ \ge 0,\xi \in [{\xi }_0,1],\hfill \end{array}\right.$$

(118)

Note that

$$\left\{\begin{array}{c}{\tilde{a}}_n^{\prime }(1)=0,\hfill \\ {\tilde{a}}_n^{\prime }(0)=\frac{1}{n!}-\frac{2}{(n+1)!}>0.\hfill \end{array}\right.$$

(119)

One has ${\tilde{a}}_n^{\prime }({\xi }_0)<{\tilde{a}}_n^{\prime }(1)=0$, and there exits ${\xi }_1\in (0,{\xi }_0)$ such that ${\tilde{a}}_n^{\prime }({\xi }_1)=0$ since ${\tilde{a}}_n^{\prime }(0)>0$. Therefore,

$$\begin{array}{c}\hfill {\tilde{a}}_n^{\prime }(\xi )\left\{\begin{array}{c}>0,\xi \in [0,{\xi }_1),\hfill \\ \le 0,\xi \in [{\xi }_1,1].\hfill \end{array}\right.\end{array}$$

(120)

Since

$$\left\{\begin{array}{c}{\tilde{a}}_n(1)=0,\hfill \\ {\tilde{a}}_n(0)=\frac{1}{2n!}-\frac{1}{(n+1)!}>0,\hfill \end{array}\right.$$

(121)

it holds that ${\tilde{a}}_n(\xi )\ge 0$ for arbitrary $\xi \in [0,1]$ when $n\ge 2$. As a result,

$$a_l(\xi )=l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}{\tilde{a}}_n(\xi )(1-\alpha )\alpha (\alpha +1)\cdots (\alpha +n-1)\ge 0,2\le l\le j.$$

(122)

Furthermore,

$$\begin{array}{cc}\hfill a_l(\xi )& =l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}\left[\frac{2{\xi }^{n+1}-2+{(1-\xi )}^2}{(n+1)!}+\frac{2-{(1-\xi )}^2-\frac{1}{2}{(1+\xi )}^2}{n!}\right](1-\alpha )\alpha (\alpha +1)\cdots (\alpha +n-1)\hfill \\ & \le l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}\left\{\frac{2{\xi }^3-2+{(1-\xi )}^2}{n+1}+\left[2-{(1-\xi )}^2-\frac{1}{2}{(1+\xi )}^2\right]\right\}\frac{1-\alpha }{n}\hfill \\ & \le l^{1-\alpha }{\displaystyle \sum _{n=2}^{\infty }}\frac{1}{l^n}\left(\frac{1}{n+1}\frac{35+13\sqrt{13}}{54}+\frac{2}{3}\right)\frac{1-\alpha }{n}\hfill \\ & \le l^{-1-\alpha }(1-\alpha )\frac{143+13\sqrt{13}}{324}\frac{1}{1-\frac{1}{l}}.\hfill \end{array}$$

(123)

Especially, for $l\ge 2$,

$$a_l(\xi )\le l^{-1-\alpha }(1-\alpha )\frac{143+13\sqrt{13}}{162}.$$

(124)

As a result, it holds that

$$\begin{array}{cc}& \left|{\displaystyle \sum _{l=1}^j}{A}_l\right|\le h^{4-\alpha }\left\{\left|{\int }_0^1{a}_1(\xi )f^{(3)}(x_{j-1}+\xi h)\mathrm{d}\xi \right|+{\displaystyle \sum _{l=2}^j}\left|A_l\right|\right\}\hfill \\ \hfill \le & h^{4-\alpha }{\displaystyle \underset{x\in [x_0,x_j]}{max}}\left|f^{(3)}(x)\right|\left\{{\int }_0^1{a}_1(\xi )\mathrm{d}\xi +{\displaystyle \sum _{l=2}^j}{\int }_0^1{a}_l(\xi )\mathrm{d}\xi \right\}\hfill \\ \hfill \le & h^{4-\alpha }{\displaystyle \underset{x\in [x_0,x_j]}{max}}\left|f^{(3)}(x)\right|\left[\frac{2}{(2-\alpha )(3-\alpha )}-\frac{1}{3(2-\alpha )}-\frac{1}{6}+{\displaystyle \sum _{l=2}^j}{l}^{-1-\alpha }(1-\alpha )\frac{31\sqrt{13}+125}{162}\right]\hfill \\ \hfill \le & C_1{\displaystyle \underset{x\in [x_0,x_j]}{max}}\left|f^{(3)}(x)\right|h^{4-\alpha },\hfill \end{array}$$

(125)

with $C_1>0$ being a constant. Consequently, the truncated error has the bound

$$\begin{array}{cc}\hfill \left|R^j\right|& =\frac{1}{2h\mathsf{\Gamma }(2-\alpha )}\left|{\displaystyle \sum _{k=0}^{j-1}}\left(A_k+B_k\right)\right|\hfill \\ & \le \frac{1}{2h\mathsf{\Gamma }(2-\alpha )}\left\{\left|{\displaystyle \sum _{l=1}^j}{A}_{j-l}\right|+\left|{\displaystyle \sum _{l=1}^j}{B}_{j-l}\right|\right\}\hfill \\ & \le \frac{h^{3-\alpha }}{2\mathsf{\Gamma }(2-\alpha )}\left\{C_1\underset{x\in [x_0,x_j]}{max}\left|f^{(3)}(x)\right|+C_2\underset{x\in [x_{-1},x_{j-2}]}{max}\left|f^{(3)}(x)\right|\right\}.\hfill \end{array}$$

(126)

Note that the derivative $f^{(3)}(x)$ with $x\in [x_{-1},x_0]$ is needed in the above inequality. In this case, $f^{(3)}(x_k)$ with $k\ge 0$ can be utilized to approximate $f^{(3)}(x)$ when $x\in [x_{-1},x_0]$ and then $f^{(3)}(x)$ is bounded on $[x_{-1},x_0]$. Consequently, the desired estimate is obtained. □

Remark 7.

In formula (102), $f(x_{-1})$ is defined outside of $[a,b]$. In numerical calculation, we can approximate $f(x_{-1})$ based on the relation $f(x_{-1})=f(a)-h{f}^{\prime }(a)+\frac{h^2}{2}{f}^{″}(a)+\mathcal{O}(h^3)$. When $f^{\prime }(a)=f^{″}(a)=0$, then $f(x_{-1})=f(a)+\mathcal{O}(h^3)$, and we have $R^j=\mathcal{O}(h^{3-\alpha })$. When $f^{\prime }(a)=0$ and $f^{″}(a)\ne 0$, then $f(x_{-1})=f(a)+\frac{h^2}{2}{f}^{″}(a)+\mathcal{O}(h^3)$, and $R^j=\mathcal{O}(h^2)$. If $f^{\prime }(a)\ne 0$, then $R^j=\mathcal{O}(h)$.

Example 3.

([79]) Consider the function $f(x)=x^4,x\in [0,1]$. Evaluate its Caputo derivative at $x=1$ by formula (102). Absolute error (AE) and convergence order (CO) are shown in

Table 3. Numerical results for Example 3.

$\mathit{\alpha }$	h	AE	CO	$\mathit{\alpha }$	h	AE	CO
0.2	$\frac{1}{10}$	0.0015	-	0.6	$\frac{1}{10}$	0.0139	-
	$\frac{1}{40}$	3.7575 × ${10}^{-5}$	2.6809		$\frac{1}{40}$	5.4517 × ${10}^{-4}$	2.3606
	$\frac{1}{160}$	8.6640 × ${10}^{-7}$	2.7289		$\frac{1}{160}$	2.0146 × ${10}^{-5}$	2.3846
0.4	$\frac{1}{10}$	0.0052	-	0.8	$\frac{1}{10}$	0.0331	-
	$\frac{1}{40}$	1.6158 × ${10}^{-4}$	2.5282		$\frac{1}{40}$	0.0017	2.1414
	$\frac{1}{160}$	4.6455 × ${10}^{-6}$	2.5676		$\frac{1}{160}$	8.1011 × ${10}^{-5}$	2.1910

. It is obvious that the convergence order is $(3-\alpha )$, which is in line with the theoretical analysis.

In Ref. [81], another $(3-\alpha )$-th order approximation was proposed. Denote

$$\left\{\begin{array}{c}{\delta }_x{f}_{j-\frac{1}{2}}=\frac{f(x_j)-f(x_{j-1})}{h},j\ge 1\hfill \\ {\delta }_x^2{f}_j=\frac{1}{h}\left({\delta }_x{f}_{j+\frac{1}{2}}-{\delta }_x{f}_{j-\frac{1}{2}}\right),j\ge 1.\hfill \end{array}\right.$$

(127)

Let $f(x)\in C^3[a,b]$ and $0<\alpha <1$. We utilize the linear interpolation

$$P_{1,k}(x)=f(x_{k-1})\frac{x_k-x}{h}+f(x_k)\frac{x-x_{k-1}}{h},$$

(128)

on the first interval $[x_0,x_1]$, and the quadratic interpolation

$$P_{2,k}(x)={\displaystyle \sum _{l=0}^2}f(x_{k-l}){\displaystyle \prod _{\begin{array}{c}i=0\\ i\ne l\end{array}}^2}\frac{x-x_{k-i}}{x_{k-l}-x_{k-i}}=P_{1,k}(x)+\frac{1}{2}\left({\delta }_x^2{f}_{k-1}\right)(x-x_{k-1})(x-x_k)$$

(129)

on the remaining intervals $[x_{k-1},x_k](k\ge 2)$ to approximate $f(x)$. Denote $x_{j+\frac{1}{2}}=\frac{x_{j+1}+x_j}{2},j\ge 0$. We obtain the following L1-2 formula [81]

$$\begin{array}{cc}& {\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{\mathsf{\Gamma }(1-\alpha )}\left[{\int }_{x_0}^{x_1}\frac{{\left(P_{1,1}(t)\right)}^{\prime }}{{(x_j-t)}^{\alpha }}\mathrm{d}t+{\displaystyle \sum _{k=2}^j}{\int }_{x_{k-1}}^{x_k}\frac{{\left(P_{2,k}(t)\right)}^{\prime }}{{(x_j-t)}^{\alpha }}\mathrm{d}t\right]+R^j\hfill \\ \hfill =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}\left\{{\displaystyle \sum _{k=2}^j}{\int }_{x_{k-1}}^{x_k}\frac{{\delta }_x{f}_{k-\frac{1}{2}}+\left({\delta }_x^2{f}_{k-1}\right)(t-x_{k-\frac{1}{2}})}{{(x_j-t)}^{\alpha }}\mathrm{d}t+{\delta }_x{f}_{\frac{1}{2}}{\int }_{x_0}^{x_1}\frac{\mathrm{d}t}{{(x_j-t)}^{\alpha }}\right\}+R^j\hfill \\ \hfill =& \frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}\left[c_0^{(\alpha )}f(x_j)-{\displaystyle {\displaystyle \sum _{k=1}^{j-1}}}\left(c_{j-k-1}^{(\alpha )}-c_{j-k}^{(\alpha )}\right)f(x_k)-c_{j-1}^{(\alpha )}f(x_0)\right]+R^j,\hfill \end{array}$$

(130)

with the truncated errors $R^1=\mathcal{O}(h^{2-\alpha })$ and $R^j=\mathcal{O}(h^{3-\alpha }),j\ge 2$. The coefficient $c_0^{(\alpha )}=1$ when $j=1$. For $j\ge 2$, the coefficients are give by

$$c_k^{(\alpha )}=\left\{\begin{array}{cc}{a}_0^{(\alpha )}+b_0^{(\alpha )},\hfill & k=0,\hfill \\ a_k^{(\alpha )}+b_k^{(\alpha )}-b_{k-1}^{(\alpha )},\hfill & 1\le k\le j-2,\hfill \\ a_k^{(\alpha )}-b_{k-1}^{(\alpha )},\hfill & k=j-1\hfill \end{array}\right.$$

(131)

with

$$a_k^{(\alpha )}={(k+1)}^{1-\alpha }-k^{1-\alpha },0\le k\le j-1,$$

(132)

and

$$b_k^{(\alpha )}=\frac{{(k+1)}^{2-\alpha }-k^{2-\alpha }}{2-\alpha }-\frac{{(k+1)}^{1-\alpha }+k^{1-\alpha }}{2},0\le k\le j-2.$$

(133)

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.

Modifying the above L1-2 formula, Alikhanov proposed an overall $(3-\alpha )$-th order approximation. Let $\sigma =1-\frac{\alpha }{2}$ with $0<\alpha <1$, then the Caputo derivative of $f(x)\in C^3[a,b]$ at $x_{j+\sigma }=a+(j+\sigma )h$ with $0\le j\le N-1$ can be expressed by

$${\left[{}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)\right]}_{x=x_{j+\sigma }}=\frac{1}{\mathsf{\Gamma }(1-\alpha )}\left[{\displaystyle \sum _{k=1}^j}{\int }_{x_{k-1}}^{x_k}\frac{f^{\prime }(t)dt}{{(x_{j+\sigma }-t)}^{\alpha }}+{\int }_{x_j}^{x_{j+\sigma }}\frac{f^{\prime }(t)dt}{{(x_{j+\sigma }-t)}^{\alpha }}\right].$$

(134)

Applying the quadratic interpolation

$$\begin{array}{ccc}\hfill {\mathsf{\Pi }}_{2,k}f(x)& =& f(x_{k-1})\frac{(x-x_k)(x-x_{k+1})}{2h^2}-f(x_k)\frac{(x-x_{k-1})(x-x_{k+1})}{h^2}\hfill \\ & & +f(x_{k+1})\frac{(x-x_{k-1})(x-x_k)}{2h^2},x\in [x_{k-1},x_k],1\le k\le j,\hfill \end{array}$$

(135)

which is different from the one defined in Equation (129) to approximating $f(x)$, and utilizing the expression $f^{\prime }(t)\approx \frac{f(x_{j+1})-f(x_j)}{h}$ on the interval $[x_j,x_{j+\sigma }]$, we obtain the L2-1${}_{\sigma }$ formula [82]

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_{j+\sigma }}=\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=0}^j}{c}_{j-k}^{(\alpha ,\sigma )}\left[f(x_{k+1})-f(x_k)\right]+\mathcal{O}(h^{3-\alpha })$$

(136)

with $0\le j\le N-1$. Here $c_0^{(\alpha ,\sigma )}=a_0^{(\alpha ,\sigma )}$ when $j=0$, and for $j\ge 1$,

$$c_k^{(\alpha ,\sigma )}=\left\{\begin{array}{cc}{a}_0^{(\alpha ,\sigma )}+b_1^{(\alpha ,\sigma )},\hfill & k=0,\hfill \\ a_k^{(\alpha ,\sigma )}+b_{k+1}^{(\alpha ,\sigma )}-b_k^{(\alpha ,\sigma )},\hfill & 1\le k\le j-1,\hfill \\ a_j^{(\alpha ,\sigma )}-b_j^{(\alpha ,\sigma )},\hfill & k=j,\hfill \end{array}\right.$$

(137)

with $a_k^{(\alpha ,\sigma )}$ and $b_k^{(\alpha ,\sigma )}$ given by

$$\left\{\begin{array}{c}{a}_0^{(\alpha ,\sigma )}={\sigma }^{1-\alpha },\hfill \\ a_k^{(\alpha ,\sigma )}={(k+\sigma )}^{1-\alpha }-{(k+\sigma -1)}^{1-\alpha },k\ge 1,\hfill \end{array}\right.$$

(138)

and

$$b_k^{(\alpha ,\sigma )}=\frac{{(k+\sigma )}^{2-\alpha }-{(k+\sigma -1)}^{2-\alpha }}{2-\alpha }-\frac{{(k+\sigma )}^{1-\alpha }+{(k+\sigma -1)}^{1-\alpha }}{2}.$$

(139)

The comparison between the L2-1${}_{\sigma }$ and L1-2 methods in Ref. [82] shows that the L2-1${}_{\sigma }$ formula refines the accuracy indeed.

Remark 8.

([80]) The L2-1${}_{\sigma }$ formula for the right-sided Caputo derivative can be derived in a similar manner. In this case, the parameter should be chosen as $\sigma =\frac{\alpha }{2}$, $\alpha \in (0,1)$. The corresponding approximation is given by

$${\left[{}_C{\mathrm{D}}_{x,b}^{\alpha }f(x)\right]}_{x=x_{j+\sigma }}=\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=j}^{N-1}}{\tilde{c}}_{k-j}^{(\alpha ,\sigma )}\left[f(x_k)-f(x_{k+1})\right]+\mathcal{O}(h^{3-\alpha })$$

(140)

with $0\le j\le N-1$. Here the coefficients are given by

$${\tilde{c}}_0^{(\alpha ,\sigma )}=-{\tilde{a}}_0^{(\alpha ,\sigma )},$$

(141)

if $j=N-1$, and for $0\le j<N-1$,

$${\tilde{c}}_k^{(\alpha ,\sigma )}=\left\{\begin{array}{cc}{\tilde{b}}_1^{(\alpha ,\sigma )},\hfill & k=0,\hfill \\ {\tilde{b}}_{k+1}^{(\alpha ,\sigma )}-{\tilde{b}}_k^{(\alpha ,\sigma )},\hfill & 1\le k\le N-j-2,\hfill \\ -{\tilde{a}}_{N-j-1}^{(\alpha ,\sigma )}-{\tilde{b}}_{N-j-1}^{(\alpha ,\sigma )},\hfill & k=N-j-1,\hfill \end{array}\right.$$

(142)

where

$${\tilde{a}}_k^{(\alpha ,\sigma )}={(k+1-\sigma )}^{1-\alpha },$$

(143)

and

$${\tilde{b}}_k^{(\alpha ,\sigma )}=\frac{{(k+1-\sigma )}^{1-\alpha }-{(k-\sigma )}^{1-\alpha }}{2}-\frac{{(k+1-\sigma )}^{2-\alpha }-{(k-\sigma )}^{2-\alpha }}{2-\alpha }.$$

(144)

For other $(3-\alpha )$-th order approximations to Caputo derivative based on interpolation, one may refer to Refs. [83,84].

(II) $(4-\alpha )$-th order approximation

Let $0<\alpha <1$ and $f(x)\in C^4[x_0,x_j]$. A linear interpolation of $f(x)$ on the first subinterval $[x_0,x_1]$ yields

$${\int }_{x_0}^{x_1}{(x_j-t)}^{-\alpha }{f}^{\prime }(t)\mathrm{d}t\approx \frac{f(x_1)-f(x_0)}{h}{\int }_{x_0}^{x_1}{(x_j-t)}^{-\alpha }\mathrm{d}t=\frac{a_{j-1}}{h^{\alpha }(1-\alpha )}\left[f(x_1)-f(x_0)\right]$$

(145)

with $a_{j-1}=j^{1-\alpha }-{(j-1)}^{1-\alpha }$. On the second subinterval $[x_1,x_2]$, we similarly obtain

$${\int }_{x_1}^{x_2}{(x_j-t)}^{-\alpha }{f}^{\prime }(t)\mathrm{d}t\approx \frac{h^{-\alpha }}{1-\alpha }\left[(a_{j-2}+b_{j-2})f(x_2)-(a_{j-2}+2b_{j-2})f(x_1)+b_{j-2}f(x_0)\right]$$

(146)

through the quadratic interpolation, where

$$\left\{\begin{array}{c}{a}_{j-2}={(j-1)}^{1-\alpha }-{(j-2)}^{1-\alpha },\hfill \\ b_{j-2}=\frac{{(j-1)}^{2-\alpha }-{(j-2)}^{2-\alpha }}{2-\alpha }-\frac{{(j-1)}^{1-\alpha }+{(j-2)}^{1-\alpha }}{2}.\hfill \end{array}\right.$$

(147)

For the remaining subintervals, we use the cubic interpolation function

$$p_3(x)={\displaystyle \sum _{l=0}^3}f(x_{k-l}){\displaystyle \prod _{i=0,i\ne l}^3}\frac{x-x_{k-i}}{x_{k-l}-x_{k-i}},x\in [x_{k-1},x_k],k\ge 3,$$

(148)

to approximate $f(x)$. Consequently, it holds that

$$\begin{array}{cc}& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{k=3}^j}{\int }_{x_{k-1}}^{x_k}\frac{f^{\prime }(t)}{{(x_j-t)}^{\alpha }}\mathrm{d}t\approx \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{k=3}^j}{\int }_{x_{k-1}}^{x_k}{(x_j-t)}^{-\alpha }{p}_3^{\prime }(t)\mathrm{d}t\hfill \\ \hfill =& \frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=3}^j}\left[{\omega }_{1,j-k}f(x_k)+{\omega }_{2,j-k}f(x_{k-1})+{\omega }_{3,j-k}f(x_{k-2})+{\omega }_{4,j-k}f(x_{k-3})\right],j\ge 3,\hfill \end{array}$$

(149)

where the coefficients are given by

$$\begin{array}{cc}\hfill {\omega }_{1,j-k}& =\frac{2{(j-k+1)}^{1-\alpha }-11{(j-k)}^{1-\alpha }}{6}-\frac{2{(j-k)}^{2-\alpha }-{(j-k+1)}^{2-\alpha }}{2-\alpha }-\frac{{(j-k)}^{3-\alpha }-{(j-k+1)}^{3-\alpha }}{(2-\alpha )(3-\alpha )},\hfill \\ \hfill {\omega }_{2,j-k}& =\frac{6{(j-k)}^{1-\alpha }+{(j-k+1)}^{1-\alpha }}{2}+\frac{5{(j-k)}^{2-\alpha }-2{(j-k+1)}^{2-\alpha }}{2-\alpha }+\frac{3{(j-k)}^{3-\alpha }-3{(j-k+1)}^{3-\alpha }}{(2-\alpha )(3-\alpha )},\hfill \\ \hfill {\omega }_{3,j-k}& =-\frac{3{(j-k)}^{1-\alpha }+2{(j-k+1)}^{1-\alpha }}{2}-\frac{4{(j-k)}^{2-\alpha }-{(j-k+1)}^{2-\alpha }}{2-\alpha }-\frac{{(j-k)}^{3-\alpha }-{(j-k+1)}^{3-\alpha }}{(2-\alpha )(3-\alpha )},\hfill \\ \hfill {\omega }_{4,j-k}& =\frac{2{(j-k)}^{1-\alpha }+{(j-k+1)}^{1-\alpha }}{6}+\frac{{(j-k)}^{2-\alpha }}{2-\alpha }+\frac{{(j-k)}^{3-\alpha }-{(j-k+1)}^{3-\alpha }}{(2-\alpha )(3-\alpha )},3\le j\le N.\hfill \end{array}$$

In view of the above analysis, we obtain the numerical approximation [85]

$$\begin{array}{cc}\hfill {\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{k=1}^j}{\int }_{x_{k-1}}^{x_k}\frac{f^{\prime }(t)}{{(x_j-t)}^{\alpha }}\mathrm{d}t=\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=0}^j}{g}_{k}f(x_{j-k})+R^j,0<\alpha <1.\hfill \end{array}$$

(150)

The coefficients $g_k$ have different values for different j. When $j=1$,

$$g_0=a_0,g_1=-a_0.$$

(151)

When $j=2$,

$$\left\{\begin{array}{c}{g}_0=a_0+b_0,\hfill \\ g_1=a_1-a_0-2b_0,\hfill \\ g_2=b_0-a_1.\hfill \end{array}\right.$$

(152)

When $j=3$,

$$\left\{\begin{array}{c}{g}_0={\omega }_{1,0},g_1={\omega }_{2,0}+a_1+b_1,\hfill \\ g_2={\omega }_{3,0}+a_2-a_1-2b_1,\hfill \\ g_3={\omega }_{4,0}-a_2+b_1.\hfill \end{array}\right.$$

(153)

When $j=4$,

$$\left\{\begin{array}{c}{g}_0={\omega }_{1,0},g_1={\omega }_{1,1}+{\omega }_{2,0},\hfill \\ g_2={\omega }_{2,1}+{\omega }_{3,0}+a_2+b_2,\hfill \\ g_3={\omega }_{3,1}+{\omega }_{4,0}+a_3-a_2-2b_2,\hfill \\ g_4={\omega }_{4,1}-a_3+b_2.\hfill \end{array}\right.$$

(154)

When $j=5$,

$$\left\{\begin{array}{c}{g}_0={\omega }_{1,0},g_1={\omega }_{1,1}+{\omega }_{2,0},\hfill \\ g_2={\omega }_{1,2}+{\omega }_{2,1}+{\omega }_{3,0},\hfill \\ g_3={\omega }_{2,2}+{\omega }_{3,1}+{\omega }_{4,0}+a_3+b_3,\hfill \\ g_4={\omega }_{3,2}+{\omega }_{4,1}+a_4-a_3-2b_3,\hfill \\ g_5={\omega }_{4,2}-a_4+b_3.\hfill \end{array}\right.$$

(155)

When $j\ge 6$,

$$\left\{\begin{array}{c}{g}_0={\omega }_{1,0},g_1={\omega }_{1,1}+{\omega }_{2,0},\hfill \\ g_2={\omega }_{1,2}+{\omega }_{2,1}+{\omega }_{3,0},\hfill \\ g_k={\omega }_{1,k}+{\omega }_{2,k-1}+{\omega }_{3,k-2}+{\omega }_{4,k-3},3\le k\le j-3,\hfill \\ g_{j-2}=a_{j-2}+b_{j-2}+{\omega }_{2,j-3}+{\omega }_{3,j-4}+{\omega }_{4,j-5},\hfill \\ g_{j-1}={\omega }_{3,j-3}+{\omega }_{4,j-4}+a_{j-1}-a_{j-2}-2b_{j-2},\hfill \\ g_j={\omega }_{4,j-3}-a_{j-1}+b_{j-2}.\hfill \end{array}\right.$$

(156)

If $f(x)\in C^4[x_0,x_j]$ and $\alpha \in (0,1)$, the truncated error $R^j$ in Equation (150) satisfies

$$\left\{\begin{array}{c}\left|R^1\right|\le c_1{\displaystyle \underset{x_0\le x\le x_1}{max}}\left|f^{″}(x)\right|h^{2-\alpha },c_1>0,\hfill \\ \left|R^2\right|\le c_2{\displaystyle \underset{x_0\le x\le x_2}{max}}\left|f^{‴}(x)\right|h^{3-\alpha },c_2>0,\hfill \\ \left|R^j\right|\le \frac{1}{\mathsf{\Gamma }(1-\alpha )}\left\{\frac{2\alpha }{3}{\displaystyle \underset{x_0\le x\le x_2}{max}}\left|f^{‴}(x)\right|\right.{(x_j-x_2)}^{-\alpha -1}{h}^4\hfill \\ \left.+\left[\frac{1}{12}+\frac{3{\alpha }^2}{2(1-\alpha )(2-\alpha )}\right]{\displaystyle \underset{x_0\le x\le x_j}{max}}\left|f^{(4)}(x)\right|h^{4-\alpha }\right\},j\ge 3.\hfill \end{array}\right.$$

(157)

Numerical examples in Ref. [85] verify the above theoretical results.

Example 4.

Suppose that$0<\alpha <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. Numerical results for Example 4.

$\mathit{\alpha }$	h	ME	CO	$\mathit{\alpha }$	h	ME	CO
0.2	$\frac{1}{10}$	1.2176 × ${10}^{-4}$	-	0.6	$\frac{1}{10}$	1.0943 × ${10}^{-3}$	-
	$\frac{1}{40}$	6.9376 × ${10}^{-7}$	3.7336		$\frac{1}{40}$	9.9598 × ${10}^{-6}$	3.3918
	$\frac{1}{160}$	3.8404 × ${10}^{-9}$	3.7528		$\frac{1}{160}$	8.9946 × ${10}^{-8}$	3.3963
0.4	$\frac{1}{10}$	4.1401 × ${10}^{-4}$	-	0.8	$\frac{1}{10}$	2.6315 × ${10}^{-3}$	-
	$\frac{1}{40}$	2.9349 × ${10}^{-6}$	3.5741		$\frac{1}{40}$	3.1265 × ${10}^{-5}$	3.1982
	$\frac{1}{160}$	2.0437 × ${10}^{-8}$	3.5855		$\frac{1}{160}$	3.7065 × ${10}^{-7}$	3.1994

.

Example 5.

Let$f(x)=e^{2x}$. We evaluate Caputo derivative of$f(x)$at$x=1$by utilizing Equation (150). The maximum errors (ME) and convergence order (CO) are shown in

Table 5. Numerical results for Example 5.

$\mathit{\alpha }$	h	ME	CO	$\mathit{\alpha }$	h	ME	CO
0.2	$\frac{1}{10}$	4.9025 × ${10}^{-4}$	-	0.6	$\frac{1}{10}$	4.2309 × ${10}^{-3}$	-
	$\frac{1}{40}$	3.8638 × ${10}^{-6}$	3.5125		$\frac{1}{40}$	4.6839 × ${10}^{-5}$	3.2902
	$\frac{1}{160}$	3.1440 × ${10}^{-8}$	3.4447		$\frac{1}{160}$	4.5469 × ${10}^{-7}$	3.3536
0.4	$\frac{1}{10}$	1.6156 × ${10}^{-3}$	-	0.8	$\frac{1}{10}$	1.0190 × ${10}^{-2}$	-
	$\frac{1}{40}$	1.4478 × ${10}^{-5}$	3.4371		$\frac{1}{40}$	1.4521 × ${10}^{-4}$	3.1086
	$\frac{1}{160}$	1.1851 × ${10}^{-7}$	3.4669		$\frac{1}{160}$	1.8089 × ${10}^{-6}$	3.1747

.

(III) $(r+1-\alpha )$-th order approximation

Generalizing the above $(4-\alpha )$-th order approximation, an $(r+1-\alpha )$-th order approximation was proposed in Ref. [86] by virtue of the Lagrange polynomials of degree r. Let $f(x)\in C^r[a,b](r\ge 4)$ and $0<\alpha <1$. On the subintervals $[x_{k-1},x_k],j\ge k\ge r,N\ge j\ge r$, we utilize the Lagrange polynomial

$$p_r(x)={\displaystyle \sum _{l=0}^r}f(x_{k-l}){\displaystyle \prod _{i=0,i\ne l}^r}\frac{x-x_{k-i}}{x_{k-l}-x_{k-i}},x\in [x_{k-1},x_k],$$

(158)

to approximate $f(x)$. Denote

$$I_k[p_r(x)]=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_{x_{k-1}}^{x_k}{(x_j-t)}^{-\alpha }{p}_r^{\prime }(t)\mathrm{d}t.$$

(159)

Then it holds that

$$\begin{array}{cc}& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_{x_{k-1}}^{x_k}\frac{f^{\prime }(t)}{{(x_j-t)}^{\alpha }}\mathrm{d}t\approx I_k[p_r(x)]\hfill \\ \hfill =& \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{l=0}^r}\frac{{(-1)}^{l}f(x_{k-l})}{l!(r-l)!h^r}{\int }_{x_{k-1}}^{x_k}{(x_j-t)}^{-\alpha }{\left[{\displaystyle \prod _{i=0,i\ne l}^r}(t-x_{k-i})\right]}^{\prime }\mathrm{d}t\hfill \\ \hfill =& \frac{h^{-\alpha }}{\mathsf{\Gamma }(1-\alpha )}{\displaystyle \sum _{l=0}^r}{\omega }_{l,j-k}^{r}f(x_{k-l}).\hfill \end{array}$$

(160)

To compute the coefficients ${\omega }_{i,j-k}^r$, we denote by

$$a_k^s=s^{k-\alpha },b_k^s={(s+1)}^{k-\alpha },p_k=\prod _{l=1}^k(l-\alpha ),$$

(161)

and

$${\alpha }_{j,i}^k=\left\{\begin{array}{cc}{\varphi }_{j,i}^k,\hfill & k\ne 0,\hfill \\ 1,\hfill & k=0,\hfill \end{array}\right.{\beta }_{j,i}^k=\left\{\begin{array}{cc}{\psi }_{j,i}^k,\hfill & k\ne 0,\hfill \\ 1,\hfill & k=0.\hfill \end{array}\right.$$

(162)

Here ${\varphi }_{j,i}^k$ and ${\psi }_{j,i}^k$ are the sums of products of all different combinations of k elements in the sets $A_{j,i}=\{\overline{a}|\overline{a}\in [0,j-1],\overline{a}\ne i,\overline{a}\in \mathbb{Z}\}$, and $B_{j,i}=\{\overline{b}|\overline{b}\in [-1,j-2],\overline{b}\ne i-1,\overline{b}\in \mathbb{Z}\}$, respectively. Then

$${\omega }_{i,j-k}^r=\frac{{(-1)}^{i+1}}{i!(r-i)!}{\displaystyle \sum _{l=1}^r}\left[\frac{l!}{p_l}\left({\alpha }_{r+1,i}^{r-l}{a}_l^{j-k}-{\beta }_{r+1,i}^{r-l}{b}_l^{j-k}\right)\right],0\le i\le r-1.$$

(163)

On the subinterval $[x_{k-1},x_k],1<k<r,1\le j\le 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 $\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_{x_{k-1}}^{x_k}{(x_j-t)}^{-\alpha }{f}^{\prime }(t)\mathrm{d}t$. In summary, we obtain the following approximation [86]

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\left\{\begin{array}{c}{\displaystyle \sum _{k=1}^j}{I}_k[p_k(x)]+R_r^j,j<r,\hfill \\ {\displaystyle \sum _{k=1}^{r-1}}{I}_k[p_k(x)]+{\displaystyle \sum _{k=r}^j}{I}_k[p_r(x)]+R_r^j,r\le j\le N,\hfill \end{array}\right.$$

(164)

with $R_r^j$ being the truncated error. It has been proved that when $f(x)\in C^r[a,b](r\ge 4)$, the truncation error satisfies

 (1)

$\left|R_r^1\right|\le c_1{\displaystyle \underset{x_0\le x\le x_1}{max}}\left|f^{(2)}(x)\right|h^{2-\alpha },c_1>0;$

 (2)

$\left|R_r^2\right|\le c_2{\displaystyle \underset{x_0\le x\le x_2}{max}}\left|f^{(3)}(x)\right|h^{3-\alpha },c_2>0;$

 (3)

If $f^{(1)}(a)=f^{(2)}(a)=0$, then $\left|R_r^3\right|\le c_3{\displaystyle \underset{x_0\le x\le x_3}{max}}\left|f^{(4)}(x)\right|h^{4-\alpha },c_3>0$;

 (4)

Provided that $f^{(k)}(a)=0$ for $0<k\le j$, then

$$\left|R_r^j\right|\le c_j{\displaystyle \underset{x_0\le x\le x_j}{max}}\left|f^{(j+1)}(x)\right|h^{j+1-\alpha },2<j<r,c_j>0;$$

 (5)

Provided that $f^{(k)}(a)=0$ for $0\le k\le r-1$, then

$$\begin{array}{ccc}\hfill \left|R_r^j\right|& =& \frac{\alpha }{\mathsf{\Gamma }(1-\alpha )}\left[\frac{1}{r+1}\left(\frac{1}{\alpha }+\frac{1}{(1-\alpha )(2-\alpha )}\right)\right.{\displaystyle \underset{x_0\le x\le x_j}{max}}\left|f^{(r+1)}(x)\right|h^{r+1-\alpha }\hfill \\ & & +{\displaystyle \sum _{k=1}^{r-1}}\frac{{(x_j-x_{r-1})}^{-\alpha -1}{(r-1)}^{r-1}}{(r-k-1)!(k+1)}\left.{\displaystyle \underset{x_0\le x\le x_j}{max}}\left|f^{(r)}(x)\right|h^{r+1}\right],j\ge r.\hfill \end{array}$$

Numerical examples in Ref. [86] verify the above theoretical results.

Example 6.

Suppose$0<\alpha <1$, and let$f(x)=x^6$,$x\in [0,1]$. Use scheme (164) to compute Caputo derivative of$f(x)$at$x=1$with different stepsizes.

Table 6. Numerical results for Example 6.

$\mathit{\alpha }$	h	Errors ($\mathit{r}=4$)	Order	Errors ($\mathit{r}=5$)	Order
0.2	$\frac{1}{40}$	3.4260 × ${10}^{-7}$	-	7.5843 × ${10}^{-9}$	-
	$\frac{1}{50}$	1.2040 × ${10}^{-7}$	4.6900	2.0999 × ${10}^{-9}$	5.7583
	$\frac{1}{60}$	5.1128 × ${10}^{-8}$	4.7001	7.3082 × ${10}^{-10}$	5.8067
0.4	$\frac{1}{40}$	1.5377 × ${10}^{-6}$	-	3.3066 × ${10}^{-8}$	-
	$\frac{1}{50}$	5.5973 × ${10}^{-7}$	4.5319	9.5097 × ${10}^{-9}$	5.5862
	$\frac{1}{60}$	2.4468 × ${10}^{-7}$	4.5408	3.4281 × ${10}^{-9}$	5.6014
0.8	$\frac{1}{40}$	1.7551 × ${10}^{-5}$	-	3.7712 × ${10}^{-7}$	-
	$\frac{1}{50}$	6.9402 × ${10}^{-6}$	4.1604	1.1820 × ${10}^{-7}$	5.1992
	$\frac{1}{60}$	3.2473 × ${10}^{-6}$	4.1673	4.5817 × ${10}^{-8}$	5.1977

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-\alpha )$.

Example 7.

Suppose$0<\alpha <1$, and consider the function$f(x)=e^{2x}-2x-2x^2-\frac{4}{3}{x}^3-\frac{2}{3}{x}^4$,$x\in [0,1]$.

Table 7. Numerical results for Example 7.

$\mathit{\alpha }$	h	Errors ($\mathit{r}=4$)	Order	Errors ($\mathit{r}=5$)	Order
0.2	$\frac{1}{26}$	7.7151 × ${10}^{-7}$	-	4.7095 × ${10}^{-8}$	-
	$\frac{1}{28}$	5.4884 × ${10}^{-7}$	4.5977	3.1083 × ${10}^{-8}$	5.6001
	$\frac{1}{30}$	3.9948 × ${10}^{-7}$	4.6064	2.1014 × ${10}^{-8}$	5.7337
0.4	$\frac{1}{26}$	3.2710 × ${10}^{-6}$	-	1.9686 × ${10}^{-7}$	-
	$\frac{1}{28}$	2.3542 × ${10}^{-6}$	4.4402	1.3170 × ${10}^{-7}$	5.4275
	$\frac{1}{30}$	1.7322 × ${10}^{-6}$	4.4490	9.0617 × ${10}^{-8}$	5.4056
0.8	$\frac{1}{26}$	3.2210 × ${10}^{-5}$	-	1.9644 × ${10}^{-6}$	-
	$\frac{1}{28}$	2.3819 × ${10}^{-5}$	4.0750	1.3521 × ${10}^{-6}$	5.0431
	$\frac{1}{30}$	1.7973 × ${10}^{-5}$	4.0832	9.5398 × ${10}^{-7}$	5.0602

lists the numerical results with different values for α, and$r=4,5$. It is evident that scheme (164) can reach$(r+1-\alpha )$-th order accuracy.

Remark 9.

The$(3-\alpha )$-th,$(4-\alpha )$-th, and$(r+1-\alpha )$-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<\alpha <m\in {\mathbb{Z}}^{+}$ and $f(x)\in H^r[a,b]$ with $r\ge 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

$$p_N(x)={\displaystyle \sum _{j=0}^N}{\tilde{p}}_j^{u,v}{P}_j^{u,v}(x),x\in [-1,1]$$

(165)

be an approximation of $f(x)$ based on the Jacobi polynomials. Recall that

$$\left\{\begin{array}{c}\frac{{\mathrm{d}}^m}{\mathrm{d}{x}^m}{P}_j^{u,v}(x)=d_{j,m}^{u,v}{P}_{j-m}^{u+m,v+m}(x),j\ge m\in {\mathbb{Z}}^{+},\hfill \\ d_{j,m}^{u,v}=\frac{\mathsf{\Gamma }(j+m+u+v+1)}{2^m\mathsf{\Gamma }(j+u+v+2)}.\hfill \end{array}\right.$$

(166)

It holds that

$$\begin{array}{cc}\hfill {}_C{\mathrm{D}}_{-1,x}^{\alpha }{p}_N(x)& =\frac{1}{\mathsf{\Gamma }(m-\alpha )}{\int }_{-1}^x{(x-t)}^{m-\alpha -1}{\displaystyle \sum _{j=m}^N}{\tilde{p}}_j^{u,v}{d}_{j,m}^{u,v}{P}_{j-m}^{u+m,v+m}(t)\mathrm{d}t\hfill \\ & ={\displaystyle \sum _{j=m}^N}{\tilde{p}}_j^{u,v}{d}_{j,m}^{u,v}{\widehat{P}}_{j-m}^{u+m,v+m,m-\alpha }(x),\hfill \end{array}$$

(167)

where ${\tilde{p}}_j^{u,v}$ and ${\widehat{P}}_j^{u+m,v+m,m-\alpha }(x)$ are defined by Equations (32) and (36). Denote

$${\mathrm{D}}_j^{u,v,\alpha ,m}(x)=d_{j,m}^{u,v}{\widehat{P}}_{j-m}^{u+m,v+m,m-\alpha }(x)$$

(168)

with ${\mathrm{D}}_j^{u,v,\alpha ,m}(x)=0$ for $0\le j\le m-1$. Then it holds that

$${}_C{\mathrm{D}}_{-1,x}^{\alpha }f(x)\approx {}_C{\mathrm{D}}_{-1,x}^{\alpha }{p}_N(x)={\displaystyle \sum _{j=m}^n}{\tilde{p}}_j^{u,v}{\mathrm{D}}_j^{u,v,\alpha ,m}(x),x\in [-1,1].$$

(169)

The affine transformation $\widehat{x}=\frac{2x-a-b}{b-a}$ with $x\in [a,b]$ yields

$${}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\approx {\left(\frac{b-a}{2}\right)}^{-\alpha }{\displaystyle \sum _{j=m}^N}{\tilde{p}}_j^{u,v}{\mathrm{D}}_j^{u,v,\alpha ,m}(\widehat{x}).$$

(170)

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^{\mu }$,$x\in [0,1]$. We use formula (170) to compute ${}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)$.

Table 8. The absolute errors for Example 8.

$u=v=0,\mu =3.5$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
20	2.49 × ${10}^{-9}$	2.90 × ${10}^{-8}$	1.99 × ${10}^{-7}$	6.63 × ${10}^{-6}$	1.92 × ${10}^{-5}$	2.67 × ${10}^{-5}$
40	2.70 × ${10}^{-11}$	4.73 × ${10}^{-10}$	4.88 × ${10}^{-9}$	2.81 × ${10}^{-7}$	1.22 × ${10}^{-6}$	2.55 × ${10}^{-6}$
80	2.88 × ${10}^{-13}$	7.62 × ${10}^{-12}$	1.19 × ${10}^{-10}$	1.18 × ${10}^{-8}$	7.77 × ${10}^{-8}$	2.44 × ${10}^{-7}$
$u=v=-\frac{1}{2},\mu =3.5$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
20	2.12 × ${10}^{-9}$	2.11 × ${10}^{-8}$	1.62 × ${10}^{-7}$	5.37 × ${10}^{-6}$	1.86 × ${10}^{-5}$	3.95 × ${10}^{-5}$
40	2.15 × ${10}^{-11}$	3.22 × ${10}^{-10}$	3.77 × ${10}^{-9}$	2.17 × ${10}^{-7}$	1.14 × ${10}^{-6}$	3.66 × ${10}^{-6}$
80	2.20 × ${10}^{-13}$	5.00 × ${10}^{-12}$	8.89 × ${10}^{-11}$	8.92 × ${10}^{-9}$	7.10 × ${10}^{-8}$	3.45 × ${10}^{-7}$

shows the absolute maximum errors at the Jacobi-Gauss-Lobatto points. The spectral accuracy is obtained.

Example 9.

([42]) Let $f(x)=sinx$, $x\in [0,1]$. We utilized Equation (170) to evaluate ${}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)$.

Table 9. The absolute errors for Example 9.

$u=v=0$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
4	1.05 × ${10}^{-5}$	6.17 × ${10}^{-5}$	2.31 × ${10}^{-4}$	1.02 × ${10}^{-3}$	2.18 × ${10}^{-3}$	5.66 × ${10}^{-3}$
8	1.48 × ${10}^{-11}$	1.08 × ${10}^{-10}$	5.96 × ${10}^{-10}$	4.14 × ${10}^{-9}$	1.47 × ${10}^{-8}$	5.33 × ${10}^{-8}$
16	3.22 × ${10}^{-15}$	1.91 × ${10}^{-14}$	9.29 × ${10}^{-14}$	6.39 × ${10}^{-13}$	2.45 × ${10}^{-12}$	9.03 × ${10}^{-12}$
$u=v=-\frac{1}{2}$
n	$\alpha =0.2$	$\alpha =0.5$	$\alpha =0.8$	$\alpha =1.2$	$\alpha =1.5$	$\alpha =1.8$
4	1.34 × ${10}^{-5}$	5.66 × ${10}^{-5}$	1.95 × ${10}^{-4}$	9.93 × ${10}^{-4}$	2.05 × ${10}^{-3}$	5.40 × ${10}^{-3}$
8	1.98 × ${10}^{-11}$	1.04 × ${10}^{-10}$	3.92 × ${10}^{-10}$	3.34 × ${10}^{-9}$	1.15 × ${10}^{-8}$	4.34 × ${10}^{-8}$
16	4.44 × ${10}^{-16}$	1.22 × ${10}^{-15}$	7.55 × ${10}^{-15}$	4.40 × ${10}^{-14}$	2.32 × ${10}^{-13}$	9.82 × ${10}^{-13}$

presents the absolute maximum errors for the cases of $u=v=0$ and $u=v=-\frac{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,\dots ,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

$$f(x)\approx S(x)={\displaystyle \sum _{j=1}^N}{\lambda }_j\varphi (|x-x_j|).$$

(171)

In order to take the values $f(x_i),i=1,2,\dots ,N$, the expansion coefficients ${\lambda }_j$ are required to satisfy the matrix form

$$\mathbf{A}\overrightarrow{\lambda }=\overrightarrow{f}$$

(172)

with $\overrightarrow{\lambda }={({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_N)}^{\top }$, $\overrightarrow{f}={(f(x_1),f(x_2),\dots ,f(x_N))}^{\top }$, and ${\mathbf{A}}_{ij}=\varphi (|x_i-x_j|)$. Here $\varphi (·)$ is the radial basis function. Some popular choices of radial basis function are cubic ($\varphi (r)=r^3$), multiquadrics ($\varphi (r)=\sqrt{r^2+c^2}$), and Gaussian ($\varphi (r)=e^{-c{r}^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

$${\displaystyle \sum _{j=1}^N}{\lambda }_j{\left.{}_C{\mathrm{D}}_{a,x}^{\alpha }\varphi (|x-x_j|)\right|}_{x=x_i}\approx g(x_i),1\le i\le N,$$

(173)

which can be written in the matrix form

$$\mathbf{B}\overrightarrow{\lambda }\approx \overrightarrow{g}$$

(174)

with ${\mathbf{B}}_{i,j}={\left.{}_C{\mathrm{D}}_{a,x}^{\alpha }\varphi (|x-x_j|)\right|}_{x=x_i}$ and $\overrightarrow{g}={(g(x_1),g(x_2),\dots ,g(x_N))}^{\top }$. Here $g(x_i)$ is the value of ${}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)$ at the point $x=x_i$. Note that the collocation matrix $\mathbf{A}$ is unconditionally nonsingular [94]. Combing equations (172) and (174) gives

$$\overrightarrow{g}\approx \mathbf{B}{\mathbf{A}}^{-1}\overrightarrow{f}.$$

(175)

Therefore, the differential matrix $\mathbf{D}=\mathbf{B}{\mathbf{A}}^{-1}$ yields an radial basis function discretiazation of the operator ${}_C{\mathrm{D}}_{a,x}^{\alpha }$.

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.

3.1.3. 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

$$\left\{\begin{array}{c}{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)={}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)-{\displaystyle \sum _{k=0}^{m-1}}\frac{f^{(k)}(a){(x-a)}^{k-\alpha }}{\mathsf{\Gamma }(1+k-\alpha )},\hfill \\ {}_C{\mathrm{D}}_{x,b}^{\alpha }f(x)={}_{RL}{\mathrm{D}}_{x,b}^{\alpha }f(x)-{\displaystyle \sum _{k=0}^{m-1}}\frac{{(-1)}^k{f}^{(k)}(b){(b-x)}^{k-\alpha }}{\mathsf{\Gamma }(1+k-\alpha )},\hfill \end{array}\right.$$

(176)

where $m-1<\alpha <m\in {\mathbb{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.

Theorem 4.

Suppose that$f(x)\in C^{[\alpha ]+3}[a,b]$, and that the derivatives of$f(x)$up to order$[\alpha ]+4$belong to$L^1[a,b]$. Then there hold

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{j+p}}{\zeta }_{2,p,k}^{(\alpha )}f(x_{j-k+p})+\mathcal{O}(h^2),0<\alpha <1,$$

(177)

and

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{j+p}}{\xi }_{2,p,k}^{(\alpha )}f(x_{j-k+p})+\mathcal{O}(h^2),1<\alpha <2,$$

(178)

where the weights${\zeta }_{2,p,k}^{(\alpha )}$and${\xi }_{2,p,k}^{(\alpha )}(k=0,1,\dots )$are given by

$${\zeta }_{2,p,k}^{(\alpha )}=\left\{\begin{array}{cc}{\varpi }_{2,p,k}^{(\alpha )},\hfill & k=0,1,\dots ,j+p-1,\hfill \\ -{\displaystyle \sum _{l=0}^{j+p-1}}{\varpi }_{2,p,l}^{(\alpha )},\hfill & k=j+p,\hfill \end{array}\right.$$

(179)

and,

$${\xi }_{2,p,k}^{(\alpha )}=\left\{\begin{array}{cc}{\varpi }_{2,p,k}^{(\alpha )},\hfill & k=0,1,\dots ,j+p-3,\hfill \\ {\varpi }_{2,p,j+p-2}^{(\alpha )}+\frac{1}{2}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{2,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-2,\hfill \\ {\varpi }_{2,p,j+p-1}^{(\alpha )}-2{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{2,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-1,\hfill \\ -{\displaystyle \sum _{l=0}^{j+p-1}}{\varpi }_{2,p,l}^{(\alpha )}+\frac{3}{2}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{2,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p.\hfill \end{array}\right.$$

(180)

Here the shift $p\le \frac{3\alpha }{2}$, and the coefficients ${\varpi }_{2,p,k}^{(\alpha )}(0\le k\le j+p)$ are given by

$${\omega }_{2,p}^{(\alpha )}(z)={\left(\frac{3\alpha -2p}{2\alpha }-\frac{2(\alpha -p)}{\alpha }z+\frac{\alpha -2p}{2\alpha }{z}^2\right)}^{\alpha }={\displaystyle \sum _{k=0}^{\infty }}{\varpi }_{2,p,k}^{(\alpha )}{z}^k,\left|z\right|<1.$$

(181)

Theorem 5.

Suppose that$f(x)\in C^{[\alpha ]+4}[a,b]$, and that the derivatives of$f(x)$up to order$[\alpha ]+5$belong to$L^1[a,b]$. Then the third order schemes are given by

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{j+p}}{\zeta }_{3,p,k}^{(\alpha )}f(x_{j-k+p})+\mathcal{O}(h^3),0<\alpha <1,$$

(182)

and

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{j+p}}{\xi }_{3,p,k}^{(\alpha )}f(x_{j-k+p})+\mathcal{O}(h^3),1<\alpha <2,$$

(183)

where the weights${\zeta }_{3,p,k}^{(\alpha )}$and${\xi }_{3,p,k}^{(\alpha )}(k=0,1,\dots )$are

$${\zeta }_{3,p,k}^{(\alpha )}=\left\{\begin{array}{cc}{\varpi }_{3,p,k}^{(\alpha )},\hfill & k=0,1,\dots ,j+p-1,\hfill \\ -{\displaystyle \sum _{l=0}^{j+p-1}}{\varpi }_{3,p,l}^{(\alpha )},\hfill & k=j+p,\hfill \end{array}\right.$$

(184)

and,

$${\xi }_{3,p,k}^{(\alpha )}=\left\{\begin{array}{cc}{\varpi }_{3,p,k}^{(\alpha )},\hfill & k=0,1,\dots ,j+p-4,\hfill \\ {\varpi }_{3,p,j+p-3}^{(\alpha )}-\frac{1}{3}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{3,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-3,\hfill \\ {\varpi }_{3,p,j+p-2}^{(\alpha )}+\frac{3}{2}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{3,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-2,\hfill \\ {\varpi }_{3,p,j+p-1}^{(\alpha )}-3{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{3,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-1,\hfill \\ -{\displaystyle \sum _{l=0}^{j+p-1}}{\varpi }_{3,p,l}^{(\alpha )}+\frac{11}{6}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{3,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p.\hfill \end{array}\right.$$

(185)

Here the shift p satisfies$3p^2-12\alpha p+11{\alpha }^2\ge 0$, and${\varpi }_{3,p,k}^{(\alpha )}(0\le k\le j+p)$are given by

$${\omega }_{3,p}^{(\alpha )}(z)={\left(a_0+a_{1}z+a_2{z}^2+a_3{z}^3\right)}^{\alpha }={\displaystyle \sum _{k=0}^{\infty }}{\varpi }_{3,p,k}^{(\alpha )}{z}^k,\left|z\right|<1,$$

(186)

with

$$\left\{\begin{array}{c}{a}_0=\frac{11{\alpha }^2-12\alpha p+3p^2}{6{\alpha }^2},a_1=\frac{-18{\alpha }^2+30\alpha p-9p^2}{6{\alpha }^2},\hfill \\ a_2=\frac{9{\alpha }^2-24\alpha p+9p^2}{6{\alpha }^2},a_3=\frac{-2{\alpha }^2+6\alpha p-3p^2}{6{\alpha }^2}.\hfill \end{array}\right.$$

(187)

Theorem 6.

Suppose that$f(x)\in C^{[\alpha ]+5}[a,b]$, and that the derivatives of$f(x)$up to order$[\alpha ]+6$belong to$L^1[a,b]$. Then there hold

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{j+p}}{\zeta }_{4,p,k}^{(\alpha )}f(x_{j-k+p})+\mathcal{O}(h^4),0<\alpha <1,$$

(188)

and

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{j+p}}{\xi }_{4,p,k}^{(\alpha )}f(x_{j-k+p})+\mathcal{O}(h^4),1<\alpha <2,$$

(189)

where the weights${\zeta }_{4,p,k}^{(\alpha )}$ and ${\xi }_{4,p,k}^{(\alpha )}$ $(k=0,1,\dots )$ are defined by

$${\zeta }_{4,p,k}^{(\alpha )}=\left\{\begin{array}{cc}{\varpi }_{4,p,k}^{(\alpha )},\hfill & k=0,1,\dots ,j+p-1,\hfill \\ -{\displaystyle \sum _{l=0}^{j+p-1}}{\varpi }_{4,p,l}^{(\alpha )},\hfill & k=j+p,\hfill \end{array}\right.$$

(190)

and,

$${\xi }_{4,p,k}^{(\alpha )}=\left\{\begin{array}{cc}{\varpi }_{4,p,k}^{(\alpha )},\hfill & k=0,1,\dots ,j+p-5,\hfill \\ {\varpi }_{4,p,j+p-4}^{(\alpha )}+\frac{1}{4}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{4,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-4,\hfill \\ {\varpi }_{4,p,j+p-3}^{(\alpha )}-\frac{4}{3}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{4,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-3,\hfill \\ {\varpi }_{4,p,j+p-2}^{(\alpha )}+3{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{4,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-2,\hfill \\ {\varpi }_{4,p,j+p-1}^{(\alpha )}-4{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{4,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p-1,\hfill \\ -{\displaystyle \sum _{l=0}^{j+p-1}}{\varpi }_{4,p,l}^{(\alpha )}+\frac{25}{12}{\displaystyle \sum _{l=0}^{j+p}}{\varpi }_{4,p,l}^{(\alpha )}(j-l+p),\hfill & k=j+p.\hfill \end{array}\right.$$

(191)

Here the shift p satisfies$25{\alpha }^3-35{\alpha }^{2}p+15\alpha p^2-2p^3\ge 0$, and${\varpi }_{4,p,k}^{(\alpha )}(0\le k\le j+p)$are given by

$${\omega }_{4,p}^{(\alpha )}(z)={\left(b_0+b_{1}z+b_2{z}^2+b_3{z}^3+b_4{z}^4\right)}^{\alpha }={\displaystyle \sum _{k=0}^{\infty }}{\varpi }_{2,p,k}^{(\alpha )}{z}^k,\left|z\right|<1,$$

(192)

with

$$\left\{\begin{array}{c}{b}_0=\frac{25{\alpha }^3-35{\alpha }^{2}p+15\alpha p^2-2p^3}{12{\alpha }^3},\hfill \\ b_1=\frac{-48{\alpha }^3+104{\alpha }^{2}p-54\alpha p^2+8p^3}{12{\alpha }^3},\hfill \\ b_2=\frac{36{\alpha }^3-114{\alpha }^{2}p+72\alpha p^2-12p^3}{12{\alpha }^3},\hfill \\ b_3=\frac{-16{\alpha }^3+56{\alpha }^{2}p-42\alpha p^2+8p^3}{12{\alpha }^3},\hfill \\ b_4=\frac{3{\alpha }^3-11{\alpha }^{2}p+9\alpha p^2-2p^3}{12{\alpha }^3}.\hfill \end{array}\right.$$

(193)

Different from numerical algorithms based on the polynomial interpolation, in which the corresponding accuracy depends on the derivative order $\alpha$ 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.

Example 10.

([95]) Consider the function $f(x)=x^{6+\alpha },x\in [0,1]$. We utilize schemes (178), (183), and (189) to evaluate the α-th order Caputo derivative at $x=1$. The absolute errors and convergence orders for $p=0$ and $p=1$ are shown in

Table 10. The absolute errors and convergence orders of Example 10.

$p=0$
$\alpha$	h	$n=2$		$n=3$		$n=4$
		$E(h)$	rate	$E(h)$	rate	$E(h)$	rate
1.70	$\frac{1}{20}$	1.0791	-	1.5714 × ${10}^{-1}$	-	1.9056 × ${10}^{-2}$	-
	$\frac{1}{80}$	7.5767 × ${10}^{-2}$	1.94	2.8201 × ${10}^{-3}$	2.93	8.4873 × ${10}^{-5}$	3.94
	$\frac{1}{320}$	4.8717 × ${10}^{-3}$	1.99	4.5587 × ${10}^{-5}$	2.99	3.4203 × ${10}^{-7}$	3.96
1.80	$\frac{1}{20}$	1.4104	-	2.0548 × ${10}^{-1}$	-	2.4921 × ${10}^{-2}$	-
	$\frac{1}{80}$	9.9074 × ${10}^{-2}$	1.94	3.6878 × ${10}^{-3}$	2.93	1.1098 × ${10}^{-4}$	3.94
	$\frac{1}{320}$	6.3705 × ${10}^{-3}$	1.99	5.9612 × ${10}^{-5}$	2.98	4.4862 × ${10}^{-7}$	3.98
1.85	$\frac{1}{20}$	1.6112	-	2.3480 × ${10}^{-1}$	-	2.8478 × ${10}^{-2}$	-
	$\frac{1}{80}$	1.1321 × ${10}^{-1}$	1.94	4.2140 × ${10}^{-3}$	2.93	1.2682 × ${10}^{-4}$	3.94
	$\frac{1}{320}$	7.2796 × ${10}^{-3}$	1.99	6.8119 × ${10}^{-5}$	2.98	5.1220 × ${10}^{-7}$	3.98
$p=1$
$\alpha$	h	$n=2$		$n=3$		$n=4$
		$E(h)$	rate	$E(h)$	rate	$E(h)$	rate
1.70	$\frac{1}{20}$	3.0173 × ${10}^{-1}$	-	4.3788 × ${10}^{-2}$	-	5.0963 × ${10}^{-3}$	-
	$\frac{1}{80}$	1.9203 × ${10}^{-2}$	1.99	7.3362 × ${10}^{-4}$	2.97	2.1650 × ${10}^{-5}$	3.96
	$\frac{1}{320}$	1.2061 × ${10}^{-3}$	2.00	1.1665 × ${10}^{-5}$	2.99	8.6359 × ${10}^{-8}$	4.00
1.80	$\frac{1}{20}$	3.2844 × ${10}^{-1}$	-	5.0679 × ${10}^{-2}$	-	6.0965 × ${10}^{-3}$	-
	$\frac{1}{80}$	2.0842 × ${10}^{-2}$	1.99	8.4958 × ${10}^{-4}$	2.97	2.5938 × ${10}^{-5}$	3.96
	$\frac{1}{320}$	1.3081 × ${10}^{-3}$	2.00	1.3511 × ${10}^{-5}$	2.99	1.0348 × ${10}^{-7}$	3.99
1.85	$\frac{1}{20}$	3.3890 × ${10}^{-1}$	-	5.4135 × ${10}^{-2}$	-	6.6280 × ${10}^{-3}$	-
	$\frac{1}{80}$	2.1454 × ${10}^{-2}$	1.99	9.0750 × ${10}^{-4}$	2.97	2.8214 × ${10}^{-5}$	3.96
	$\frac{1}{320}$	1.3456 × ${10}^{-3}$	2.00	1.4432 × ${10}^{-5}$	2.99	1.1252 × ${10}^{-7}$	3.99

. The experiment convergence orders are consistent with theoretical analysis. Furthermore, the shifted numerical methods are more efficient than the unshifted one when $1<\alpha <2$.

Example 11.

([95]) Consider $f(x)=x^{6+\alpha }+{(x+1)}^2,x\in [0,1]$. In this case, $f(0)\ne 0$. We utilize schemes (178), (183), and (189) to evaluate its α-th order Caputo derivative at $x=1$. Numerical results are presented in

Table 11. The absolute errors and convergence orders of Example 11 with $p=1$.

$\alpha$	h	$n=2$		$n=3$		$n=4$
		$E(h)$	rate	$E(h)$	rate	$E(h)$	rate
1.70	$\frac{1}{20}$	3.0157 × ${10}^{-1}$	-	4.3799 × ${10}^{-2}$	-	5.0897 × ${10}^{-3}$	-
	$\frac{1}{80}$	1.9193 × ${10}^{-2}$	1.99	7.3380 × ${10}^{-4}$	2.97	2.1622 × ${10}^{-5}$	3.96
	$\frac{1}{320}$	1.2055 × ${10}^{-3}$	2.00	1.1668 × ${10}^{-5}$	2.99	8.6891 × ${10}^{-8}$	3.98
1.80	$\frac{1}{20}$	3.2833 × ${10}^{-1}$	-	5.0688 × ${10}^{-2}$	-	6.0949 × ${10}^{-3}$	-
	$\frac{1}{80}$	2.0836 × ${10}^{-2}$	1.99	8.4970 × ${10}^{-4}$	2.97	2.5935 × ${10}^{-5}$	3.96
	$\frac{1}{320}$	1.3077 × ${10}^{-3}$	2.00	1.3513 × ${10}^{-5}$	2.99	1.0350 × ${10}^{-7}$	3.99
1.85	$\frac{1}{20}$	3.3882 × ${10}^{-1}$	-	5.4142 × ${10}^{-2}$	-	6.6269 × ${10}^{-3}$	-
	$\frac{1}{80}$	2.1449 × ${10}^{-2}$	1.99	9.0760 × ${10}^{-4}$	2.97	2.8211 × ${10}^{-5}$	3.96
	$\frac{1}{320}$	1.3453 × ${10}^{-3}$	2.00	1.4433 × ${10}^{-5}$	2.99	1.1225 × ${10}^{-7}$	3.99

. These results imply that the numerical approximations can be used to compute Caputo derivatives of suitably smooth functions with inhomogeneous conditions at the initial time.

3.1.4. Diffusive Approximation

Recall that Caputo derivative is defined as the Riemann–Liouville integral of an integer-order derivative, i.e.,

$${}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)={}_{RL}{\mathrm{D}}_{0,x}^{-(m-\alpha )}{f}^{(m)}(x),m-1<\alpha <m\in {\mathbb{Z}}^{+}.$$

(194)

In this case, Equation (62) implies

$${}_C{\mathrm{D}}_{0,x}^{\alpha }f(x)={\int }_0^{\infty }\varphi (\omega ,x)\mathrm{d}\omega ,$$

(195)

with $m-1<\alpha <m\in {\mathbb{Z}}^{+}$ and $\varphi :(0,\infty )\times [0,x]\to \mathbb{R}$ being the auxiliary bivariate function defined by

$$\varphi (\omega ,x)={(-1)}^{m-1}\frac{2sin(\pi \alpha )}{\pi }{\omega }^{2\alpha -2m+1}{\int }_0^x{f}^{(m)}(t)e^{-(x-t){\omega }^2}\mathrm{d}t.$$

(196)

For fixed $\omega >0$, the function $\varphi (\omega ,·)$ satisfies the differential equation

$$\frac{\partial }{\partial x}\varphi (\omega ,x)={(-1)}^{m-1}\frac{2sin(\pi \alpha )}{\pi }{\omega }^{2\alpha -2m+1}{f}^{(m)}(x)-{\omega }^2\varphi (\omega ,x)$$

(197)

subject to the initial condition $\varphi (\omega ,0)=0$. Consequently, any implementation solving ODEs and suitable quadratures approximating the infinite integral (196) yield numerical approximations to Caputo derivative.

For more details of discussions, modifications, and applications of the diffusive approximation, see Refs. [60,61].

3.2. 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.

3.2.1. Numerical Methods Based on Linear Spline Interpolation

Let $f(x)\in C^4[a,b]$ and $1<\alpha <2$. For $j=1,2,\dots ,N-1$, there holds

$$\begin{array}{cc}\hfill {\left[{}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}& =\frac{1}{\mathsf{\Gamma }(2-\alpha )}{\left[\frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}{\int }_a^x{(x-t)}^{1-\alpha }f(t)\mathrm{d}t\right]}_{x=x_j}\hfill \\ & \approx \frac{h^{-2}}{\mathsf{\Gamma }(2-\alpha )}\left[{\mathcal{I}}_{\alpha }^l(x_{j-1})-2{\mathcal{I}}_{\alpha }^l(x_j)+{\mathcal{I}}_{\alpha }^l(x_{j+1})\right],\hfill \end{array}$$

(198)

where ${\mathcal{I}}_{\alpha }^l(x_j)={\int }_a^{x_j}{(x_j-t)}^{1-\alpha }f(t)\mathrm{d}t$. Approximate $f(x)$ with the linear spline

$$s_j^l(x)={\displaystyle \sum _{k=0}^j}f(x_k)s_{j,k}^l(x),$$

(199)

where

$$s_{j,k}^l(x)=\left\{\begin{array}{c}\frac{x-x_{k-1}}{x_k-x_{k-1}},x_{k-1}\le x\le x_k,\hfill \\ \frac{x_{k+1}-x}{x_{k+1}-x_k},x_k\le x\le x_{k+1},\hfill \\ 0,\mathrm{otherwise},\hfill \end{array}\right.$$

(200)

for $1\le k\le j-1$,

$$s_{j,0}^l(x)=\left\{\begin{array}{cc}\frac{x_1-x}{x_1-x_0},\hfill & \hfill x_0\le x\le x_1,\\ 0,\hfill & \hfill \mathrm{otherwise},\end{array}\right.$$

(201)

and

$$s_{j,j}^l(x)=\left\{\begin{array}{cc}\frac{x-x_{j-1}}{x_j-x_{j-1}},\hfill & x_{j-1}\le x\le x_j,\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(202)

Then we obtain an approximation to ${\mathcal{I}}_{\alpha }^l(x_j)$ given by

$$I_{\alpha }^l(x_j)={\int }_a^{x_j}{(x_j-t)}^{1-\alpha }{s}_j^l(t)\mathrm{d}t=\frac{h^{2-\alpha }}{(2-\alpha )(3-\alpha )}\sum _{k=0}^j{a}_{j,k}^{l}f(x_k)$$

(203)

with

$$a_{j,k}^l=\left\{\begin{array}{c}(3-\alpha )j^{2-\alpha }+{(j-1)}^{3-\alpha }-j^{3-\alpha },k=0,\hfill \\ {(j-k+1)}^{3-\alpha }-2{(j-k)}^{3-\alpha }+{(j-k-1)}^{3-\alpha },1\le k\le j-1,\hfill \\ 1,k=j.\hfill \end{array}\right.$$

(204)

As a result, there holds [96]

$$\begin{array}{cc}& {\left[{}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_j}\hfill \\ \hfill \approx & \frac{h^{-\alpha }}{\mathsf{\Gamma }(4-\alpha )}\left[{\displaystyle \sum _{k=0}^{j-1}}{a}_{j-1,k}^{l}f(x_k)-2{\displaystyle \sum _{k=0}^j}{a}_{j,k}^{l}f(x_k)+{\displaystyle \sum _{k=0}^{j+1}}{a}_{j+1,k}^{l}f(x_k)\right]\hfill \\ \hfill =& \frac{h^{-\alpha }}{\mathsf{\Gamma }(4-\alpha )}[{\displaystyle \sum _{k=0}^{j-1}}\left(a_{j-l,k}^l-2a_{j,k}^l+a_{j+l,k}^l\right)f(x_k)+\left(a_{j+1,j}^l-2a_{j,j}^l\right)f(x_j)+a_{j+1,j+1}^{l}f(x_{j+1})].\hfill \end{array}$$

(205)

Similarly, the right-sided Riemann–Liouville derivative can be approximated by [96]

$$\begin{array}{ccc}\hfill {\left[{}_{RL}{\mathrm{D}}_{x,b}^{\alpha }f(x)\right]}_{x=x_j}& \approx & \frac{h^{-\alpha }}{\mathsf{\Gamma }(4-\alpha )}\left[{\displaystyle \sum _{k=j+1}^N}\left(a_{j-1,k}^r-2a_{j,k}^r+a_{j,k}^r\right)f(x_k)\right.\hfill \\ & & +a_{j-1,j-1}^{r}f(x_{j-1})+\left(a_{j-1,j}^r-2a_{j,j}^r\right)f(x_j)],\hfill \end{array}$$

(206)

where

$$\left\{\begin{array}{c}{a}_{j,N}^r=(3-\alpha ){(N-j)}^{2-\alpha }+{(N-j-1)}^{3-\alpha }-{(N-j)}^{3-\alpha },\hfill \\ a_{j,k}^r={(k-j+1)}^{3-\alpha }-2{(k-j)}^{3-\alpha }+{(k-j-1)}^{3-\alpha },j+1\le k\le N-1,\hfill \\ a_{j,j}^r=1.\hfill \end{array}\right.$$

(207)

The truncated error of this approach has been proved in Ref. [96] to be $\mathcal{O}(h^2)$ provided that $f^{(4)}(x)$ has compact support on $[a,b]$.

In the particular cases with $a=-\infty$ and $b=+\infty$, Equations (205) and (206) can be written as [96,97,98]

$${\left[{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)\right]}_{x=x_j}\approx \frac{h^{-\alpha }}{\mathsf{\Gamma }(4-\alpha )}{\displaystyle \sum _{m=-1}^{\infty }}{q}_{m}f(x_{j-m}),$$

(208)

and

$${\left[{}_{RL}{\mathrm{D}}_{x,\infty }^{\alpha }f(x)\right]}_{x=x_j}\approx \frac{h^{-\alpha }}{\mathsf{\Gamma }(4-\alpha )}{\displaystyle \sum _{m=-1}^{\infty }}{q}_{m}f(x_{j+m}),$$

(209)

with

$$q_m=\left\{\begin{array}{cc}{a}_{m-1}-2a_m+a_{m+1},\hfill & m\ge 1,\hfill \\ -2a_0+a_1,\hfill & m=0,\hfill \\ a_0,\hfill & m=-1.\hfill \end{array}\right.$$

(210)

Here

$$a_m=\left\{\begin{array}{cc}1,\hfill & m=0,\hfill \\ {(m+1)}^{3-\alpha }-2m^{3-\alpha }+{(m-1)}^{3-\alpha },\hfill & m\ge 1.\hfill \end{array}\right.$$

(211)

Both series on the right-hand side of Equations (208) and (209) converge absolutely for $1<\alpha <2$ if $f(x)$ is bounded [96]. When $\alpha =1$, Equations (208) and (209) reduce to the second order finite difference formula for the first order derivative. When $\alpha =2$, Equations (208) and (209) are consistent with the central difference formula for the second order derivative.

3.2.2. Grünwald-Letnikov Type Approximations

Let $m-1<\alpha <m\in {\mathbb{Z}}^{+}$ and $f(x)$ be m times continuously differentiable. It is known that when $a=-\infty$ or $f(a)=0$, the Grünwald–Letnikov derivative

$${}_{GL}{\mathrm{D}}_{a,x}^{\alpha }f(x)=\underset{h\to 0}{lim}\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\left[\frac{x-a}{h}\right]}}{(-1)}^l\left(\genfrac{}{}{0pt}{}{\alpha }{l}\right)f(x-lh),$$

(212)

can approximate the $\alpha$-th order Riemann–Liouville derivative with first order accuracy [48], i.e.,

$${}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^j}{\omega }_l^{(\alpha )}f(x-lh)+\mathcal{O}(h),jh=x-a,$$

(213)

where ${\omega }_l^{(\alpha )}={(-1)}^l\left(\genfrac{}{}{0pt}{}{\alpha }{l}\right)$. Equation (213) can be verified through the Fourier transform. The above numerical approximation, which is called the classical Grünwald–Letnikov formula, is warmly applied to solving fractional differential equations. However, this approximation is not suitable for the discretization of fractional differential equations when $\alpha \in (1,2)$ since it leads to unstable numerical schemes [99].

One way to construct stable schemes for fractional differential equations is to make the corresponding coefficient matrix diagonally dominated via replacing $f(x-lh)$ in Equation (213) by $f(x-(l-p)h)$ with $p\in \mathbb{Z}$ being the shift. In this case, we obtain the shifted Grünwald–Letnikov formula

$${}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{[j+p]}}{\omega }_l^{(\alpha )}f(x-(l-p)h)+\mathcal{O}(h),jh=x-a.$$

(214)

It turns out that the best performances of the shifted Grünwald–Letnikov method come from minimizing $|p-\frac{\alpha }{2}|$. And it coincides with the central difference of the classical second order differentiation. If the shift is chosen to be non-integers, numerical method (214) may have superconvergent behaviors [100,101,102].

Introducing integer shifts to the classical Grünwald-Letnikov approximation may eliminate the instability indeed, while the truncated error is still $\mathcal{O}(h)$. A modification called the weighted and shifted Grünwald-Letnikov formulae was proposed in Ref. [103], in the spirit of eliminating the first order terms in the truncated errors of the shifted Grünwald–Letnikov formulae. In the following, we introduce the basic idea in details.

Let $f(x)\in L^1(\mathbb{R})$, ${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha +1}f$, and its Fourier transform belong to $L^1(\mathbb{R})$. It can be verified through the Fourier analysis that the shifted Grünwald–Letnikov difference operator

$$A_{h,p}^{\alpha }f(x)=h^{-\alpha }{\displaystyle \sum _{k=0}^{\infty }}{\omega }_k^{(\alpha )}f\left(x-(k-p)h\right),$$

(215)

with $p\in \mathbb{Z}$ approximates the left-sided Riemann-Liouville derivative ${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)$ with first order accuracy. Assume that the following weighted and shifted Grünwald–Letnikov difference (WSGD) operator

$$B_{h,p,q}^{\alpha }f(x)=a_1{A}_{h,p}^{\alpha }f(x)+b_1{A}_{h,q}^{\alpha }f(x),a_1,b_1\in \mathbb{R},$$

(216)

approximates the Riemann–Liouville derivative with second order accuracy. Then the Fourier transform gives

$$\begin{array}{cc}\hfill \mathcal{F}\left\{B_{h,p,q}^{\alpha }f(x);\xi \right\}& =\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{\infty }}{\omega }_k^{(\alpha )}\left(a_1{e}^{-i(k-p)h\xi }+b_1{e}^{-i(k-q)h\xi }\right)\mathcal{F}\{f(x);\xi \}\hfill \\ & =\frac{1}{h^{\alpha }}\left(a_1{(1-e^{-ih\xi })}^{\alpha }{e}^{iph\xi }+b_1{(1-e^{-ih\xi })}^{\alpha }{e}^{iqh\xi }\right)\mathcal{F}\{f(x);\xi \}\hfill \\ & ={\left(\frac{1-e^{ih\xi }}{ih\xi }\right)}^{\alpha }\left[a_1{e}^{iph\xi }+b_1{e}^{iqh\xi }\right]{(i\xi )}^{\alpha }\mathcal{F}\{f(x);\xi \}\hfill \end{array}$$

(217)

with $i^2=-1$. Note that

$$\mathcal{F}\{{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x);\xi \}={(i\xi )}^{\alpha }\mathcal{F}\{f(x);\xi \},$$

(218)

and

$${\left(\frac{1-e^z}{z}\right)}^{\alpha }{e}^{zr}=1+(r-\frac{\alpha }{2})z+{\mathcal{O}(|z|}^2).$$

(219)

Therefore, $a_1$ and $b_1$ need to satisfy

$$\left\{\begin{array}{c}{a}_1+b_1=1,\hfill \\ a_1(p-\frac{\alpha }{2})+b_1(q-\frac{\alpha }{2})=0,\hfill \end{array}\right.$$

(220)

to assure that $B_{h,p,q}^{\alpha }f(x)$ is of second order accuracy. In other words,

$$B_{h,p,q}^{\alpha }f(x)={}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)+\mathcal{O}(h^2)$$

(221)

holds uniformly when $a_1=\frac{\alpha -2q}{2p-2q}$ and $b_1=\frac{2p-\alpha }{2p-2q}$ [103].

Remark 11.

The relevant academic literature has revealed that numerical approximation (221) results in unstable schemes for fractional partial differential equations when the shift paring$(p,q)=(0,-1)$[99]. The corresponding schemes are stable when$(p,q)=(1,0)$and$(p,q)=(1,-1)$. Furthermore, the WSGD operator reduces to the centred difference approximation for the classical second order differentiation when$(p,q)=(1,0)$and$(p,q)=(1,-1)$in the case with$\alpha =2$, while for the classical first order one when$(p,q)=(1,0)$if$\alpha =1$.

A third order WSGD operator was also proposed in Ref. [103]. Nevertheless, it fails to obtain stable numerical schemes when $\alpha \in (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.

3.2.3. 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 $\alpha >0$ replaced by $-\alpha$. In fact, Lubich formulae (51) are applicable for evaluating Riemann–Liouville derivative indeed.

Let $f^{(k)}(a+)=0(k=0,1,\dots ,\ell -1)$. We have the classical Lubich formulae [106]

$${}_{RL}{\mathrm{D}}_{a,x}^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\left[\frac{x-a}{h}\right]}}{\varpi }_{\ell ,l}^{(\alpha )}f(x-lh)+\mathcal{O}(h^{\ell }),\alpha >0,$$

(222)

in which h is the stepsize. The convolution coefficients ${\varpi }_{\ell ,l}^{(\alpha )}$ are generated by $W_{\ell }^{(\alpha )}(z)$ defined in Equation (50). This can be readily verified by the Fourier transform.

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,

 (1)

If$0<\alpha <1$, then${\varpi }_{3,l}^{(\alpha )}\le {\varpi }_{3,l+1}^{(\alpha )}$for$l\ge 4$,${\varpi }_{4,l}^{(\alpha )}\le {\varpi }_{4,l+1}^{(\alpha )}$for$l\ge 7$, and${\varpi }_{5,l}^{(\alpha )}\le {\varpi }_{5,l+1}^{(\alpha )}$for$l\ge 12$;

 (2)

If$1<\alpha <2$, then${\varpi }_{3,l}^{(\alpha )}\ge {\varpi }_{3,l+1}^{(\alpha )}$for$l\ge 7$,${\varpi }_{4,l}^{(\alpha )}\ge {\varpi }_{4,l+1}^{(\alpha )}$for$l\ge 12$, and${\varpi }_{5,l}^{(\alpha )}\ge {\varpi }_{5,l+1}^{(\alpha )}$for$l\ge 16$.

Recently, the above conjecture for${\varpi }_{3,l}^{(\alpha )}$with$0<\alpha <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)\in C^{[\alpha ]+3}(\mathbb{R})$, and all the derivatives of$f(x)$up to order$[\alpha ]+4$belong to$L^1(\mathbb{R})$. Then we have

$${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{k}_{2,l}^{(\alpha )}f\left(x-(l-1)h\right)+\mathcal{O}(h^2),$$

(223)

and

$${}_{RL}{\mathrm{D}}_{x,+\infty }^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{k}_{2,l}^{(\alpha )}f\left(x+(l-1)h\right)+\mathcal{O}(h^2),$$

(224)

as$h\to 0$. Here$k_{2,l}^{(\alpha )}(l=0,1,\dots )$are generated by

$${\tilde{W}}_2(z)={\left(\frac{3\alpha -2}{2\alpha }-\frac{2(\alpha -1)}{\alpha }z+\frac{\alpha -2}{2\alpha }{z}^2\right)}^{\alpha }={\displaystyle \sum _{l=0}^{\infty }}{k}_{2,l}^{(\alpha )}{z}^l,\left|z\right|<1.$$

(225)

Theorem 8.

Let$p\ge 3$,$f(x)\in C^{[\alpha ]+p+1}(\mathbb{R})$, and all the derivatives of$f(x)$up to order$[\alpha ]+p+2$belong to$L^1(\mathbb{R})$. Then there hold

$${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{k}_{p,l}^{(\alpha )}f\left(x-(l-1)h\right)+\mathcal{O}(h^p),p\ge 3,$$

(226)

and

$${}_{RL}{\mathrm{D}}_{x,+\infty }^{\alpha }f(x)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{k}_{p,l}^{(\alpha )}f\left(x+(l-1)h\right)+\mathcal{O}(h^p),p\ge 3.$$

(227)

Here the generating functions of coefficients$k_{p.l}^{(\alpha )}(l=0,1,\dots )$with$p\ge 3$are

$${\tilde{W}}_p(z)=((1-z)+\frac{\alpha -2}{2\alpha }{(1-z)}^2+{\displaystyle \sum _{k=3}^p}\frac{{\lambda }_{k-1,k-1}^{(\alpha )}}{\alpha }{(1-z)}^k{)}^{\alpha },$$

(228)

in which the parameters${\lambda }_{k-1,k-1}^{(\alpha )}(k=3,4,\dots )$can be determined by the relation

$${\tilde{W}}_k(e^{-z})\frac{e^z}{z^{\alpha }}=1-{\displaystyle \sum _{l=k}^{\infty }}{\lambda }_{k,l}^{(\alpha )}{z}^l,k=2,3,\dots$$

(229)

Introducing suitable weights and shifts to the classical Lubich operators, we can obtain weighted and shifted Lubich formulae [105], which are not only of high-order accuracy but also stable when $\alpha \in (1,2)$.

Define the operator

$${\mathcal{A}}_p^{\alpha }=\frac{1}{h^{\alpha }}{\displaystyle \sum _{k=0}^{\infty }}{\varpi }_{2,k}^{(\alpha )}f(x-(k-p)h),$$

(230)

where the shift p is an integer. The coefficients ${\varpi }_{2,k}^{(\alpha )}$ can be calculated by Equation (50) with $\ell =2$. The weighted and shifted Lubich formulae, whose convergence and accuracy can be verified through the Fourier transform, are given as follows.

Theorem 9.

Let$f(x)$,${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha +1}f(x)$(or${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha +2}f(x)$) with$1<\alpha <2$and their Fourier transforms belong to$L^1(\mathbb{R})$. Then we have

$$\left\{\begin{array}{c}{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)={\mathcal{A}}_p^{\alpha }f(x)+\mathcal{O}(h),p\ne 0,\hfill \\ {}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)={\mathcal{A}}_p^{\alpha }f(x)+\mathcal{O}(h^2),p=0,\hfill \end{array}\right.$$

(231)

where${\mathcal{A}}_p^{\alpha }$is given by Equation (230)

Theorem 10.

When$f(x)$,${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha +2}f(x)$with$1<\alpha <2$, and their Fourier transforms belong to$L^1(\mathbb{R})$, there holds

$${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)={\mathcal{A}}_{p,q}^{\alpha }f(x)+\mathcal{O}(h^2).$$

(232)

Here

$${\mathcal{A}}_{p,q}^{\alpha }f(x)=W_p{\mathcal{A}}_p^{\alpha }f(x)+W_q{\mathcal{A}}_q^{\alpha }f(x)$$

(233)

with${\mathcal{A}}_p^{\alpha }$,${\mathcal{A}}_q^{\alpha }$being defined in Equation (230), $W_p=\frac{q}{q-p}$, $W_q=\frac{p}{p-q}$, $p\ne q$, and $p,q$ being integers.

It was proved in Ref. [105] that the approximation (232) with $1<\alpha <2$ works well for space fractional differential equations when the pair $(p,q)=(1,q)$ with $\left|q\right|\ge 2$.

Theorem 11.

Assume that$f(x)$,${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha +3}f(x)$with$1<\alpha <2$, and their Fourier transforms belong to$L^1(\mathbb{R})$. Then there holds

$${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)={\mathcal{A}}_{p,q,r,s}^{\alpha }f(x)+\mathcal{O}(h^3),$$

(234)

$${\mathcal{A}}_{p,q,r,s}^{\alpha }f(x)=W_{p,q}{\mathcal{A}}_{p,q}^{\alpha }f(x)+W_{r,s}{\mathcal{A}}_{r,s}^{\alpha }f(x),$$

(235)

where${\mathcal{A}}_{p,q}^{\alpha }$and${\mathcal{A}}_{r,s}^{\alpha }$are defined in Equation (233), $W_{p,q}=\frac{3rs+2\alpha }{3(rs-pq)}$, $W_{r,s}=\frac{3pq+2\alpha }{3(pq-rs)}$, $rs\ne pq$, and $p,q,r,s$ are integers.

When $(p,q,r,s)=(1,q,1,s),\left|q\right|\ge 2,\left|s\right|\ge 2,\mathrm{and}qs<0$, the approximation (234) with $1<\alpha <2$ works well for space fractional differential equations.

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,

$${}^L{\mathcal{B}}_2^{\alpha }f(x+sh)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{k}_{2,l}^{(\alpha )}f\left(x-(l-s-1)h\right),$$

(236)

and

$${}^R{\mathcal{B}}_2^{\alpha }f(x+sh)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{k}_{2,l}^{(\alpha )}f\left(x+(l-s-1)h\right),$$

(237)

where the coefficients$k_{2,l}^{(\alpha )}$are given by the function

$${\tilde{W}}_2(z)={\left(\frac{3\alpha -1}{2\alpha }-\frac{2(\alpha -1)}{\alpha }z+\frac{\alpha -2}{2\alpha }{z}^2\right)}^{\alpha }={\displaystyle \sum _{l=0}^{\infty }}{k}_{2,l}^{(\alpha )}{z}^l,\left|z\right|<1.$$

(238)

If we introduce the fractional-compact difference operator

$$\mathcal{L}f(x)=\left(1-\frac{2{\alpha }^2-6\alpha +3}{6\alpha }{\delta }_x^2\right)f(x),$$

(239)

with${\delta }_x^2$being a second order central difference operator defined by${\delta }_x^{2}f(x)=f(x+h)-2f(x)+f(x+h)$, then equalities

$${}^L{\mathcal{B}}_2^{\alpha }f(x)=\mathcal{L}{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)+\mathcal{O}(h^3)$$

(240)

and

$${}^R{\mathcal{B}}_2^{\alpha }f(x)=\mathcal{L}{}_{RL}{\mathrm{D}}_{x,\infty }^{\alpha }f(x)+\mathcal{O}(h^3)$$

(241)

hold uniformly for$x\in \mathbb{R}$, provided that$f(x)\in C^{[\alpha ]+4}(\mathbb{R})$and all derivatives of$f(x)$up to order$[\alpha ]+5$belong to$L^1(\mathbb{R})$.

Theorem 13.

Choose the generating function as

$${\tilde{\tilde{W}}}_2(z)={\left(\frac{3\alpha +2}{2\alpha }-\frac{2(\alpha +1)}{\alpha }z+\frac{\alpha +2}{2\alpha }{z}^2\right)}^{\alpha }={\displaystyle \sum _{l=0}^{\infty }}{\tilde{k}}_{2,l}^{(\alpha )}{z}^l,\left|z\right|<1,$$

(242)

and define the following difference operators

$${}^L{\tilde{\mathcal{B}}}_2^{\alpha }f(x+sh)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{\tilde{k}}_{2,l}^{(\alpha )}f\left(x-(l-s+1)h\right),$$

(243)

$${}^R{\tilde{\mathcal{B}}}_2^{\alpha }f(x+sh)=\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{\infty }}{\tilde{k}}_{2,l}^{(\alpha )}f\left(x+(l-s+1)h\right).$$

(244)

Then the equalities

$${}^L{\tilde{\mathcal{B}}}_2^{\alpha }f(x)=\tilde{\mathcal{L}}{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)+\mathcal{O}(h^3),$$

(245)

and

$${}^R{\tilde{\mathcal{B}}}_2^{\alpha }f(x)=\tilde{\mathcal{L}}{}_{RL}{\mathrm{D}}_{x,\infty }^{\alpha }f(x)+\mathcal{O}(h^3)$$

(246)

hold uniformly for$x\in \mathbb{R}$, provided that$f(x)\in C^{[\alpha ]+4}(\mathbb{R})$and all derivatives of$f(x)$up to order$[\alpha ]+5$belong to$L^1(\mathbb{R})$. Here the fractional-compact difference operator$\tilde{\mathcal{L}}$is given by

$$\tilde{\mathcal{L}}f(x)=\left(1-\frac{2{\alpha }^2+6\alpha +3}{6\alpha }{\delta }_x^2\right)f(x).$$

(247)

Theorem 14.

Let$f(x)\in C^{[\alpha ]+5}(\mathbb{R})$and all derivatives of$f(x)$up to order$[\alpha ]+6$belong to$L^1(\mathbb{R})$. Define the fractional-compact operator as

$$\mathcal{H}f(x)=\left[\left({\tilde{\sigma }}_{3,0}^{(\alpha )}-{\sigma }_{3,0}^{(\alpha )}\right)+\left({\sigma }_{2,0}^{(\alpha )}{\tilde{\sigma }}_{3,0}^{(\alpha )}-{\tilde{\sigma }}_{2,0}^{(\alpha )}{\sigma }_{3,0}^{(\alpha )}\right){\delta }_x^2\right]f(x),$$

(248)

where

$$\left\{\begin{array}{c}{\sigma }_{2,0}^{(\alpha )}=-\frac{2{\alpha }^2-6\alpha +3}{6\alpha },{\tilde{\sigma }}_{2,0}^{(\alpha )}=-\frac{2{\alpha }^2+6\alpha +3}{6\alpha },\hfill \\ {\sigma }_{3,0}^{(\alpha )}=\frac{3{\alpha }^3-11{\alpha }^2+12\alpha -4}{12{\alpha }^2},{\tilde{\sigma }}_{3,0}^{(\alpha )}=\frac{3{\alpha }^3+11{\alpha }^2+12\alpha +4}{12{\alpha }^2}.\hfill \end{array}\right.$$

(249)

Then

$${\tilde{\sigma }}_{3,0}^{(\alpha )}{}^L{\mathcal{B}}_2^{\alpha }f(x)-{\sigma }_{3,0}^{(\alpha )}{}^L{\tilde{\mathcal{B}}}_2^{\alpha }f(x)=\mathcal{H}{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)+\mathcal{O}(h^4),$$

(250)

and

$${\tilde{\sigma }}_{3,0}^{(\alpha )}{}^R{\mathcal{B}}_2^{\alpha }f(x)-{\sigma }_{3,0}^{(\alpha )}{}^R{\tilde{\mathcal{B}}}_2^{\alpha }f(x)=\mathcal{H}{}_{RL}{\mathrm{D}}_{x,\infty }^{\alpha }f(x)+\mathcal{O}(h^4)$$

(251)

hold uniformly on$\mathbb{R}$.

The idea of the above approximations can be applied to evaluating tempered fractional derivatives, see Ref. [110] for more details.

3.2.4. Fractional Average Central Difference Method

In Ref. [111], a shifted operator of the form

$${}_C{\Delta }_{-h}^{\alpha }f(x)={\displaystyle \sum _{k=0}^{\infty }}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)f\left(x-\left(k-\frac{\alpha }{2}\right)h\right)$$

(252)

with $h>0$ was proposed to approximate Riemann–Liouville derivative. This difference operator reduces to the standard central difference operator when $\alpha$ is a positive integer. It can be verified through the Fourier transform that the equality

$$\frac{1}{h^{\alpha }}{}_C{\Delta }_{-h}^{\alpha }f(x)={}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)+\mathcal{O}(h^2),$$

(253)

holds uniformly for $x\in \mathbb{R}$ as $h\to 0$, provided that $f\in C^{[\alpha ]+3}(\mathbb{R})$ and all derivative of f up to the order $[\alpha ]+4$ exist and belong to $L^1(\mathbb{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]

$${}_{AC}{\Delta }_{-h}^{\alpha }f(x)=\frac{1}{2}{\displaystyle \sum _{k=0}^{\infty }}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)\left[f(x-kh)+f(x-(k-\alpha )h)\right],$$

(254)

and

$${}_{AC}{\Delta }_{+h}^{\alpha }f(x)=\frac{1}{2}{\displaystyle \sum _{k=0}^{\infty }}{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)\left[f(x+kh)+f(x+(k-\alpha )h)\right].$$

(255)

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}{\mathrm{D}}_{-\infty ,x}^{\alpha +2}f(x)$and${}_{RL}{\mathrm{D}}_{x,+\infty }^{\alpha +2}f(x)$are in$L^1(\mathbb{R})$. Then the equalities

$${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)=\frac{{}_{AC}{\Delta }_{-h}^{\alpha }f(x)}{h^{\alpha }}+\mathcal{O}(h^2),$$

(256)

$${}_{RL}{\mathrm{D}}_{x,+\infty }^{\alpha }f(x)=\frac{{}_{AC}{\Delta }_{+h}^{\alpha }f(x)}{h^{\alpha }}+\mathcal{O}(h^2)$$

(257)

hold uniformly on$\mathbb{R}$.

Theorem 16.

When$f(x)$and the Fourier transforms of${}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha +4}f(x)$and${}_{RL}{\mathrm{D}}_{x,+\infty }^{\alpha +4}f(x)$are in$L^1(\mathbb{R})$, then the relations

$$\left[1+\frac{\alpha (3\alpha +1)}{24}{\delta }_x^2\right]{}_{RL}{\mathrm{D}}_{-\infty ,x}^{\alpha }f(x)=\frac{1}{h^{\alpha }}{}_{AC}{\Delta }_{-h}^{\alpha }f(x)+\mathcal{O}(h^4),$$

(258)

and

$$\left[1+\frac{\alpha (3\alpha +1)}{24}{\delta }_x^2\right]{}_{RL}{\mathrm{D}}_{x,+\infty }^{\alpha }f(x)=\frac{1}{h^{\alpha }}{}_{AC}{\Delta }_{+h}^{\alpha }f(x)+\mathcal{O}(h^4)$$

(259)

hold uniformly on$\mathbb{R}$. Here${\delta }_x^2$denotes the second order central difference operator defined by${\delta }_x^{2}f(x_j)=f(x_{j+1})-2f(x_j)+f(x_{j-1})$.

For functions defined on $[a,b]$, the fractional average central difference formulae can be modified through suitable extensions.

3.3. 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${}_{\sigma }$ 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.

3.3.1. Approximation Based on L2-1${}_{\sigma }$ Formulae

In Ref. [107], the L2-1${}_{\sigma }$ formulae are reformulated in the following forms,

$${\left[{}_C{\mathrm{D}}_{a,x}^{\alpha }f(x)\right]}_{x=x_{j+\sigma }}=\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=0}^j}{d}_{j-k}^{(\alpha ,\sigma )}\left[f(x_{k+1})-f(x_k)\right]+\mathcal{O}(h^{3-\alpha }),$$

(260)

and

$${\left[{}_C{\mathrm{D}}_{x,b}^{\alpha }f(x)\right]}_{x=x_{j+{\sigma }^{\prime }}}=-\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}{\displaystyle \sum _{k=j}^{N-1}}{\tilde{d}}_{k-j}^{(\alpha ,{\sigma }^{\prime })}\left[f(x_{k+1})-f(x_k)\right]+\mathcal{O}(h^{3-\alpha }),$$

(261)

with $0\le j\le N-1$, $\sigma =1-\frac{\alpha }{2}$, and ${\sigma }^{\prime }=\frac{\alpha }{2}$. When $j=0$, $d_0^{(\alpha ,\sigma )}=c_0^{(\alpha ,\sigma )}$. When $j=N-1$, ${\tilde{d}}_0^{(\alpha ,\sigma )}=c_0^{(\alpha ,{\sigma }^{\prime })}$. For $j\ge 1$, the coefficients are given by

$$d_k^{(\alpha ,\sigma )}=\left\{\begin{array}{cc}{c}_0^{(\alpha ,\sigma )}+{\tilde{c}}_1^{(\alpha ,\sigma )},\hfill & k=0,\hfill \\ c_k^{(\alpha ,\sigma )}+{\tilde{c}}_{k+1}^{(\alpha ,\sigma )}-{\tilde{c}}_k^{(\alpha ,\sigma )},\hfill & 1\le k\le j-1,\hfill \\ c_i^{(\alpha ,\sigma )}-{\tilde{c}}_i^{(\alpha ,\sigma )},\hfill & k=j,\hfill \end{array}\right.$$

(262)

in which the second case of the right-hand side of Equation (262) should be removed if $j=1$, and

$${\tilde{d}}_k^{(\alpha ,{\sigma }^{\prime })}=\left\{\begin{array}{c}{c}_0^{(\alpha ,{\sigma }^{\prime })}+{\tilde{c}}_2^{(\alpha ,{\sigma }^{\prime })},k=0,\hfill \\ c_{k+1}^{(\alpha ,{\sigma }^{\prime })}-{\tilde{c}}_{k+1}^{(\alpha ,{\sigma }^{\prime })}+{\tilde{c}}_{k+2}^{(\alpha ,{\sigma }^{\prime })},1\le k\le N-2-j,\hfill \\ c_{k+1}^{(\alpha ,{\sigma }^{\prime })}-{\tilde{c}}_{k+1}^{(\alpha ,{\sigma }^{\prime })},k=N-1-j,\hfill \end{array}\right.$$

(263)

in which the second case of the right-hand side of Equation (263) should be removed if $j=N-2$. Here

$$\left\{\begin{array}{ccc}\hfill c_0^{(\alpha ,\sigma )}& =& {\sigma }^{1-\alpha },\hfill \\ \hfill c_k^{(\alpha ,\sigma )}& =& {(k+\sigma )}^{1-\alpha }-{(k-1+\sigma )}^{1-\alpha },k\ge 1,\hfill \\ \hfill {\tilde{c}}_k^{(\alpha ,\sigma )}& =& \frac{1}{2-\alpha }\left[{(k+\sigma )}^{2-\alpha }-{(k-1+\sigma )}^{2-\alpha }\right]\hfill \\ & & -\frac{1}{2}\left[{(k+\sigma )}^{1-\alpha }+{(k-1+\sigma )}^{1-\alpha }\right],k\ge 1,\hfill \end{array}\right.$$

(264)

and

$$\left\{\begin{array}{ccc}\hfill c_0^{(\alpha ,{\sigma }^{\prime })}& =& {(1-{\sigma }^{\prime })}^{1-\alpha },\hfill \\ \hfill c_k^{(\alpha ,{\sigma }^{\prime })}& =& {(k-{\sigma }^{\prime })}^{1-\alpha }-{(k-1-{\sigma }^{\prime })}^{1-\alpha },k\ge 1,\hfill \\ \hfill {\tilde{c}}_k^{(\alpha ,{\sigma }^{\prime })}& =& \frac{1}{2-\alpha }\left[{(k-{\sigma }^{\prime })}^{2-\alpha }-{(k-1-{\sigma }^{\prime })}^{2-\alpha }\right]\hfill \\ & & -\frac{1}{2}\left[{(k-{\sigma }^{\prime })}^{1-\alpha }+{(k-1-{\sigma }^{\prime })}^{1-\alpha }\right],k\ge 1.\hfill \end{array}\right.$$

(265)

In order to combine Equations (260) and (261), $\sigma +{\sigma }^{\prime }=1$ should be satisfied such that $2\sigma -2+\alpha +2{\sigma }^{\prime }-\alpha =0$. Furthermore, assume that $\sigma ={\sigma }^{\prime }$, i.e., $\sigma ={\sigma }^{\prime }=\frac{1}{2}$, then $x_{j+\sigma }=x_{j+\frac{1}{2}}=x_{j+{\sigma }^{\prime }}$. One can get the following $(3-\alpha )$-th order scheme for Riesz derivatives at $x=x_{j+\frac{1}{2}}$ with $0\le j\le N-1$ [107],

$$\begin{array}{ccc}\hfill {\left[{}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)\right]}_{x=x_{j+\frac{1}{2}}}& =& -\frac{1}{2cos(\frac{\pi \alpha }{2})}\left[\frac{{(x_{j+\frac{1}{2}}-a)}^{-\alpha }f(a)}{\mathsf{\Gamma }(1-\alpha )}+\frac{{(b-x_{j+\frac{1}{2}})}^{-\alpha }f(b)}{\mathsf{\Gamma }(1-\alpha )}\right.\hfill \\ & & +\frac{h^{-\alpha }}{\mathsf{\Gamma }(2-\alpha )}\left({\displaystyle \sum _{k=0}^j}{d}_{j-k}^{(\alpha ,\frac{1}{2})}\left[f(x_{k+1})-f(x_k)\right]\right.\hfill \\ & & \left.\left.-{\displaystyle \sum _{k=j}^{N-1}}{\tilde{d}}_{k-j}^{(\alpha ,\frac{1}{2})}\left[f(x_{k+1})-f(x_k)\right]\right)\right]+\mathcal{O}(h^{3-\alpha }),\hfill \end{array}$$

(266)

in which $d_{j-k}^{(\alpha ,\frac{1}{2})}$ and ${\tilde{d}}_{k-j}^{(\alpha ,\frac{1}{2})}$ are defined by Equations (262) and (263) with $\sigma ={\sigma }^{\prime }=\frac{1}{2}$, respectively.

3.3.2. 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

$${\Delta }_h^{\alpha }f(x)={\displaystyle \sum _{k=-\infty }^{+\infty }}{g}_k^{(\alpha )}f(x-kh),$$

(267)

with the coefficients being given by

$$g_k^{(\alpha )}=\frac{{(-1)}^k\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\frac{\alpha }{2}-k+1)\mathsf{\Gamma }(\frac{\alpha }{2}+k+1)},k=0,\pm 1,\pm 2,\dots$$

(268)

It can be verified through the Fourier analysis that the relation [114]

$${}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)=-\frac{1}{h^{\alpha }}{\Delta }_h^{\alpha }f(x)+\mathcal{O}(h^2),1<\alpha \le 2,$$

(269)

holds uniformly for $x\in \mathbb{R}$ as $h\to 0^{+}$, provided that $f\in C^5(\mathbb{R})$ and all of its derivatives up to order five belong to $L^1(\mathbb{R})$. It has been pointed out by Ref. [115] that Equation (269) also holds for $0<\alpha \le 1$.

3.3.3. Weighted and Shifted Centred Difference Operators

Introducing shifts to the symmetric fractional centred difference operator in Equation (267), the shifted centred difference operators [112]

$${\mathcal{L}}_{\theta }f(x)={\displaystyle \sum _{k=-\infty }^{\infty }}{g}_k^{(\alpha )}f(x-(k+\theta )h),\left|\theta \right|=0,1,2,\dots$$

(270)

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(\mathbb{R})$ with all the derivatives up to order 7 in $L^1(\mathbb{R})$, then the relation

$${}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)=\frac{1}{h^{\alpha }}\left[\frac{\alpha }{24}{\mathcal{L}}_{-1}f(x)-\left(1+\frac{\alpha }{12}\right){\mathcal{L}}_{0}f(x)+\frac{\alpha }{24}{\mathcal{L}}_{1}f(x)\right]+\mathcal{O}(h^4)$$

(271)

holds uniformly for $x\in \mathbb{R}$.

Theorem 18.

([112]) Assume that $f(x)\in C^9(\mathbb{R})$ with all the derivatives up to order 9 in $L^1(\mathbb{R})$. Then

$$\begin{array}{cc}\hfill {}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)& =\frac{1}{h^{\alpha }}\left[{\mathcal{A}}_1{\mathcal{L}}_{-2}f(x)+{\mathcal{A}}_2{\mathcal{L}}_{-1}f(x)+{\mathcal{A}}_3{\mathcal{L}}_{0}f(x)+{\mathcal{A}}_2{\mathcal{L}}_{1}f(x)+{\mathcal{A}}_1{\mathcal{L}}_{2}f(x)\right]+\mathcal{O}(h^6)\hfill \end{array}$$

(272)

holds uniformly for $x\in \mathbb{R}$, in which

$$\left\{\begin{array}{c}{\mathcal{A}}_1=-\left(\frac{\alpha }{1152}+\frac{11}{2880}\right)\alpha ,\hfill \\ {\mathcal{A}}_2=\left(\frac{\alpha }{288}+\frac{41}{720}\right)\alpha ,\hfill \\ {\mathcal{A}}_3=-\left(\frac{{\alpha }^2}{192}+\frac{17\alpha }{160}+1\right).\hfill \end{array}\right.$$

(273)

Theorem 19.

([112]) Assume that $f(x)$ lies in $C^{11}(\mathbb{R})$ with all the derivatives up to order 11 in $L^1(\mathbb{R})$. Then

$$\begin{array}{ccc}\hfill {}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)& =& \frac{1}{h^{\alpha }}\left[{\mathcal{B}}_1{\mathcal{L}}_{-3}f(x)+{\mathcal{B}}_2{\mathcal{L}}_{-2}f(x)+{\mathcal{B}}_3{\mathcal{L}}_{-1}f(x)+{\mathcal{B}}_4{\mathcal{L}}_{0}f(x)\right.\hfill \\ & & +\left.{\mathcal{B}}_3{\mathcal{L}}_{1}f(x)+{\mathcal{B}}_2{\mathcal{L}}_{2}f(x)+{\mathcal{B}}_1{\mathcal{L}}_{3}f(x)\right]+\mathcal{O}(h^8)\hfill \end{array}$$

(274)

holds uniformly for $x\in \mathbb{R}$. Here

$$\left\{\begin{array}{c}{\mathcal{B}}_1=\left(\frac{{\alpha }^2}{82944}+\frac{11\alpha }{69120}+\frac{191}{362880}\right)\alpha ,{\mathcal{B}}_2=-\left(\frac{{\alpha }^2}{13824}+\frac{7\alpha }{3840}+\frac{211}{30240}\right)\alpha ,\hfill \\ {\mathcal{B}}_3=\left(\frac{5{\alpha }^2}{27648}+\frac{3\alpha }{512}+\frac{7843}{120960}\right)\alpha ,{\mathcal{B}}_4=-\left(\frac{5{\alpha }^3}{20736}+\frac{29{\alpha }^2}{3456}+\frac{5297\alpha }{45360}+1\right).\hfill \end{array}\right.$$

(275)

For much higher-order difference operators in this respect, see Ref. [112].

3.3.4. 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.

Theorem 20.

([117]) Suppose that $f(x)\in C^{2n+3}(\mathbb{R})$, and all the derivatives of $f(x)$ up to order $2n+4$ exist and belong to $L^1(\mathbb{R})$. Then

$$\left({\delta }_x^0-b_{n-1}{\delta }_x^{2n-2}\right){}_{RZ}{\mathrm{D}}_x^{\alpha }f(x)=\left({\displaystyle \sum _{l=0}^{n-2}}{b}_l{\delta }_x^{2l}\right)\left(-\frac{{\Delta }_h^{\alpha }f(x)}{h^{\alpha }}\right)+\mathcal{O}(h^{2n}),n\in {\mathbb{Z}}^{+},$$

(276)

where

$${\delta }_x^{2l}f(x_j)={\displaystyle \sum _{s=0}^{2l}}{(-1)}^s\left(\genfrac{}{}{0pt}{}{2l}{s}\right)f(x_{l+j-s}),l\ge 0.$$

(277)

Specifically, ${\delta }_x^0$ is the identity operator, i.e., ${\delta }_x^{0}f(x_j)=f(x_j)$. The coefficients $b_l(l=0,1,\dots ,n-2)$ satisfy the following equation

$$\begin{array}{cc}& {\displaystyle \sum _{l=0}^{n-2}}{b}_l\left(2{\displaystyle \sum _{s=0}^{l-1}}{\displaystyle \sum _{q=0}^{n-1}}{\displaystyle \sum _{p=0}^{n-1-q}}\frac{{(-1)}^{s+q}{(l-s)}^{2q}\left(\genfrac{}{}{0pt}{}{2l}{s}\right)a_p}{(2q)!}{\left|\omega h\right|}^{2(p+q)}+{(-1)}^l\left(\genfrac{}{}{0pt}{}{2l}{l}\right){\displaystyle \sum _{p=0}^{n-1}}{a}_p{\left|\omega h\right|}^{2p}\right)\hfill \\ \hfill =& 1-b_{n-1}\left({\displaystyle \sum _{s=0}^{n-2}}{(-1)}^s\left(\genfrac{}{}{0pt}{}{2n-2}{s}\right)\frac{2{(n-1-s)}^{2n-2}}{(2n-2)!}\right){(-1)}^{n-1}{\left|\omega h\right|}^{2n-2},\hfill \end{array}$$

(278)

and $a_p(p=0,1,\dots )$ satisfy the equation

$$\begin{array}{c}{\displaystyle \sum _{p=0}^{\infty }}{a}_p{\left|\omega h\right|}^{2p}={\left|\frac{2sin\left(\frac{\omega h}{2}\right)}{\omega h}\right|}^{\alpha }=[1-\frac{\alpha }{24}{\left|\omega h\right|}^2+\left(\frac{1}{1920}+\frac{\alpha -1}{1152}\right)\alpha {\left|\omega h\right|}^4\hfill \\ -\left(\frac{1}{322560}+\frac{\alpha -1}{46080}+\frac{(\alpha -1)(\alpha -2)}{82944}\right)\alpha {\left|\omega h\right|}^6+\cdots ].\hfill \end{array}$$

(279)

Remark 14.

In view of proofs in Refs. [111,118,119], conditions in Theorem 20 can be weakened as$f(x)\in {\mathcal{L}}^{2n+\alpha }(\mathbb{R})$, where

$$\begin{array}{cc}\hfill {\mathcal{L}}^{2n+\alpha }(\mathbb{R})=& \left\{f\left|f\in L^1(\mathbb{R}),\right.\mathrm{and}{\int }_{\mathbb{R}}{(1+|\omega |)}^{2n+\alpha }\left|\widehat{f}(\omega )\right|\mathrm{d}\omega <+\infty \right\}.\hfill \end{array}$$

(280)

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.

For function $f(x)$ defined on the bounded interval $[a,b]$ with $f(a)=f(b)=0$, we can extend $f(x)$ by zero outside of the domain. When the conditions in Theorem 20 or Remark 14 are satisfied, the compact centred difference formula can be written in the following form,

$$({\delta }_x^0-b_{n-1}{\delta }_x^{2n-2})\frac{{\partial }^{\alpha }f(x)}{{\partial \left|x\right|}^{\alpha }}=\left(-\frac{1}{h^{\alpha }}{\displaystyle \sum _{l=0}^{n-2}}{b}_l{\delta }_x^{2l}{\displaystyle \sum _{k=[-\frac{b-x}{h}]+l}^{[\frac{x-a}{h}]-l}}{g}_k^{(\alpha )}f(x-kh)\right)+\mathcal{O}(h^{2n}).$$

(281)

The following numerical examples demonstrate the accuracy of the fractional-compact centred formula (281) and the assertion in Remark 15.

Example 12.

([117]) Consider the function $f_n(x)=x^{2n}{(1-x)}^{2n},x\in [0,1],n=2,3,4,5$. Utilize numerical scheme (281) with $n=2,3,4,5$ to compute the Riesz derivative of $f(x)$ at $x=0.5$. The absolute error (AE) and experimental convergence order (CO) displayed in

Table 12. The absolute error and the experimental convergence order of function $f_2(x)$ in Example 12 by numerical scheme (281) with $n=2,3,4,5$.

$n=2$
$\alpha$	h	AE	CO	$\alpha$	h	AE	CO
1.1	$\frac{1}{20}$	1.985528 × ${10}^{-6}$	-	1.7	$\frac{1}{20}$	9.316621 × ${10}^{-6}$	-
	$\frac{1}{80}$	7.806460 × ${10}^{-9}$	3.9981		$\frac{1}{80}$	3.661963 × ${10}^{-8}$	3.9982
	$\frac{1}{320}$	3.050456 × ${10}^{-11}$	4.0000		$\frac{1}{320}$	1.431370 × ${10}^{-10}$	3.9995
1.3	$\frac{1}{20}$	3.418165 × ${10}^{-6}$	-	1.9	$\frac{1}{20}$	1.486627 × ${10}^{-5}$	-
	$\frac{1}{80}$	1.343913 × ${10}^{-8}$	3.9981		$\frac{1}{80}$	5.841643 × ${10}^{-8}$	3.9983
	$\frac{1}{320}$	5.251979 × ${10}^{-11}$	3.9998		$\frac{1}{320}$	2.284402 × ${10}^{-10}$	3.9988
$n=3$
1.1	$\frac{1}{20}$	3.120201 × ${10}^{-8}$	-	1.7	$\frac{1}{20}$	2.053203 × ${10}^{-7}$	-
	$\frac{1}{28}$	4.223802 × ${10}^{-9}$	5.9537		$\frac{1}{28}$	2.780411 × ${10}^{-8}$	5.9527
	$\frac{1}{36}$	9.422903 × ${10}^{-10}$	5.9735		$\frac{1}{36}$	6.204232 × ${10}^{-9}$	5.9727
1.3	$\frac{1}{20}$	6.008620 × ${10}^{-8}$	-	1.9	$\frac{1}{20}$	3.675466 × ${10}^{-7}$	-
	$\frac{1}{28}$	8.135715 × ${10}^{-9}$	5.9531		$\frac{1}{28}$	4.976658 × ${10}^{-8}$	5.9530
	$\frac{1}{33}$	1.815230 × ${10}^{-9}$	5.9730		$\frac{1}{36}$	1.110461 × ${10}^{-8}$	5.9727
$n=4$
1.1	$\frac{1}{30}$	3.442344 × ${10}^{-11}$	-	1.7	$\frac{1}{30}$	2.889640 × ${10}^{-10}$	-
	$\frac{1}{38}$	5.303531 × ${10}^{-12}$	7.9206		$\frac{1}{38}$	4.449502 × ${10}^{-11}$	7.9249
	$\frac{1}{46}$	1.164521 × ${10}^{-12}$	7.9393		$\frac{1}{46}$	9.752953 × ${10}^{-12}$	7.9502
1.3	$\frac{1}{30}$	7.195825 × ${10}^{-11}$	-	1.9	$\frac{1}{30}$	5.601931 × ${10}^{-10}$	-
	$\frac{1}{38}$	1.108608 × ${10}^{-11}$	7.9213		$\frac{1}{38}$	8.624075 × ${10}^{-11}$	7.9260
	$\frac{1}{46}$	2.433728 × ${10}^{-12}$	7.9402		$\frac{1}{46}$	1.890060 × ${10}^{-11}$	7.9506
$n=5$
1.1	$\frac{1}{30}$	2.669378 × ${10}^{-12}$	-	1.7	$\frac{1}{30}$	2.756366 × ${10}^{-11}$	-
	$\frac{1}{38}$	2.057061 × ${10}^{-13}$	9.5521		$\frac{1}{38}$	2.007171 × ${10}^{-12}$	9.8190
	$\frac{1}{46}$	2.983513 × ${10}^{-14}$	8.8684		$\frac{1}{46}$	2.368039 × ${10}^{-13}$	10.2781
1.3	$\frac{1}{30}$	6.076893 × ${10}^{-12}$	-	1.9	$\frac{1}{30}$	5.640743 × ${10}^{-11}$	-
	$\frac{1}{38}$	4.563785 × ${10}^{-13}$	9.6687		$\frac{1}{38}$	4.097488 × ${10}^{-12}$	9.8250
	$\frac{1}{46}$	6.127250 × ${10}^{-14}$	9.3784		$\frac{1}{46}$	4.868649 × ${10}^{-13}$	10.2198

are in line with the theoretical analysis.

Example 13.

([117]) Consider the function $f(x)=x(1-x),x\in [0,1]$. This function fails to meet prerequisites for Equation (281) and Remark 14. We numerically compute its Riesz derivative at $x=0.5$ by using scheme (281) with $n=2$. The absolute error (AE) and experimental convergence order (CO) are displayed in

Table 13. Numerical results of Example 13 by using scheme (281) with $n=2$.

$\mathit{\alpha }$	h	AE	CO	$\mathit{\alpha }$	h	AE	CO
1.1	$\frac{1}{10}$	5.900848 × ${10}^{-4}$	-	1.7	$\frac{1}{10}$	6.101600 × ${10}^{-4}$	-
	$\frac{1}{40}$	3.674047 × ${10}^{-5}$	2.0011		$\frac{1}{40}$	3.779590 × ${10}^{-5}$	2.0026
	$\frac{1}{160}$	2.295736 × ${10}^{-6}$	2.0001		$\frac{1}{160}$	2.360928 × ${10}^{-6}$	2.0002
1.3	$\frac{1}{10}$	6.828118 × ${10}^{-4}$	-	1.9	$\frac{1}{10}$	2.863584 × ${10}^{-4}$	-
	$\frac{1}{40}$	4.245071 × ${10}^{-5}$	2.0015		$\frac{1}{40}$	1.770032 × ${10}^{-5}$	2.0032
	$\frac{1}{160}$	2.652296 × ${10}^{-6}$	2.0001		$\frac{1}{160}$	1.105502 × ${10}^{-6}$	2.0002

. One can see that the expected fourth order accuracy is not achieved, which verifies the assertion in Remark 15.

4. 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.