Approximation of Caputo Fractional Derivative and Numerical Solutions of Fractional Differential Equations

1. Introduction

Fractional differential equations are an essential means of modeling complex processes in different fields of science [3,10,35,50,51,53]. The Riemann-Liouville fractional integral of order $\alpha >0$ is defined as

$${}_a{I}_x^{\alpha }f\left(x\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_a^x{(x-t)}^{\alpha -1}f\left(t\right)dt.$$

When n is a positive integer and $0<\alpha <1$, the Caputo and Riemann-Liouville fractional derivatives of order $\alpha$ are defined as

$$\begin{array}{cc}\hfill {}_a^C{D}_x^{n+\alpha }f\left(x\right)& =I^{1-\alpha }{D}^{n+1}f\left(x\right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_a^x{\displaystyle \frac{f^{(n+1)}\left(t\right)}{{(x-t)}^{\alpha }}}dt,\hfill \\ \hfill {}_a^{RL}{D}_x^{n+\alpha }f\left(x\right)& =D^{n+1}{I}^{1-\alpha }f\left(x\right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\displaystyle \frac{d^{n+1}}{d{x}^{n+1}}}{\int }_a^x{\displaystyle \frac{f\left(t\right)}{{(x-t)}^{\alpha }}}dt.\hfill \end{array}$$

The Caputo and Riemann-Liouville fractional derivatives are related as

$${}_a^{RL}{D}_x^{n+\alpha }f\left(x\right){=}_a^C{D}_x^{n+\alpha }f\left(x\right)+\sum _{k=0}^n{\displaystyle \frac{f^{\left(k\right)}\left(0\right)x^{k-n-\alpha }}{\Gamma (k-\alpha -n+1)}}.$$

Without loss of generality the lower limit of the integral in the definition of fractional derivative is zero. When $0<\alpha <1$ the Caputo derivative is defined as

$$f^{\left(\alpha \right)}\left(x\right)=D_x^{\alpha }f\left(x\right)={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\int }_0^x{\displaystyle \frac{f^{\prime }\left(t\right)}{{(x-t)}^{\alpha }}}dt.$$

Fractional derivatives and integrals are nonlocal and have a singularity at the endpoint. The properties of the fractional derivatives make them an important tool for describing memory processes. The exponential, sine and cosine functions have Caputo derivatives,

$$D_x^{\alpha }{e}^x=x^{1-\alpha }{E}_{1,2-\alpha }\left(x\right),D_x^{\alpha }sinx=x^{1-\alpha }{E}_{2,2-\alpha }(-x^2),D_x^{\alpha }cosx=-x^{2-\alpha }{E}_{2,3-\alpha }(-x^2),$$

where $E_{{\beta }_1,{\beta }_2}\left(x\right)$ is the Mittag-Lefler function

$$E_{{\beta }_1,{\beta }_2}\left(x\right)=\sum _{k=0}^{\infty }{\displaystyle \frac{x^k}{\Gamma ({\beta }_{1}k+{\beta }_2)}}.$$

Finite difference schemes for numerical solution of fractional differential equations use discretizations of the fractional derivative. Two important discretizations of the Caputo fractional derivative are the Grünwald difference approximation and L1 approximation [8,26,30,49]. Grünwald difference approximation has a first order accuracy and a generating function ${(1-x)}^{\alpha }$. Central difference approximations of integer order derivatives have a second order accuracy. Grünwald difference approximation is a second order shifted approximation of the fractional derivative with a shift parameter $\alpha /2$, and it is a generalization of finite difference approximations [37,58]. L1 approximation has an order $2-\alpha$ and a generating function $\frac{{(x-1)}^2{\mathrm{Li}}_{\alpha -1}\left(x\right)}{x\Gamma (2-\alpha )}$. In Apostolov [7] and Dimitrov [16] we use the L1 approximation and approximations of the second derivative for construction of a second order approximations of the fractional derivative. Alikhanov [4] and Gao et al. [22] construct approximations of the fractional derivative of order $3-\alpha$, which are called $L1-2$ and $L1-2_{\sigma }$ formulas. Finite difference schemes using L2 formula approximations of the fractional derivative of order $3-\alpha$, and their convergence are studied by Alikhanov [5], Lv and Xu [33], Wong and Ren [55]. High order approximations of the fractional derivative and finite difference schemes for fractional differential equations are constructed in [8,9,12,14,29,32,44,48,57]. Navot [38] uses Taylor polynomials to deal with the singularity of the fractional integral ${\int }_0^1{t}^{\alpha }f\left(t\right)dt$. He derives the asymptotic formula for Riemann sums on a uniform net, called extended Euler-Maclaurin summation formula. The formula is extended in Navot [39] to functions with a logarithmic singularity. In [18] we derive the asymptotic formula of the Riemann sum of the fractional integral using the series expansion of the generating function ${\mathrm{Li}}_{1-\alpha }\left(e^{-x}\right)$. Other types of numerical methods for fractional differential equations include fractional multistep methods and methods using spline interpolation and wavelets [1,20,31,34,46,47].

Let n be a positive integer, $h=x/n$ and $y_k=y\left(kh\right)$ for $k=0,1,\cdots ,n$. In [17] we derive an approximation of Caputo fractional derivative and its asymptotic formula

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{f_{n-k}}{k^{1+\alpha }}}=\Gamma (-\alpha )f_n^{\left(\alpha \right)}+{\displaystyle \frac{\zeta (1+\alpha )}{h^{\alpha }}}{f}_n-\zeta \left(\alpha \right)f_n^{\prime }{h}^{1-\alpha }+O\left(h^{2-\alpha }\right).\end{array}$$

(1)

Approximation (1) has a generating function ${\mathrm{Li}}_{1+\alpha }\left(x\right)$, the polylogarithmic function of order $1+\alpha$. By substituting the derivative in the right-hand side of (1) using first order backward difference $(f_n-f_{n-1})/h$ we find an approximation of the fractional derivative

$$\begin{array}{c}\hfill {\mathcal{L}}_{n}f\left(x\right)={\displaystyle \frac{h^{-\alpha }}{\Gamma (-\alpha )}}\sum _{k=0}^m{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+O\left(h^{2-\alpha }\right),\end{array}$$

(2)

$${\sigma }_0^{\left(\alpha \right)}=\zeta \left(\alpha \right)-\zeta (1+\alpha ),{\sigma }_1^{\left(\alpha \right)}=1-\zeta \left(\alpha \right),{\sigma }_k^{\left(\alpha \right)}={\displaystyle \frac{1}{k^{1+\alpha }}},(2\le k\le n).$$

Approximations (1) and (2) have an order $2-\alpha$ when the function f satisfies the condition $f\left(0\right)=f^{\prime }\left(0\right)=0$. In paper [17] we extend approximation (2) to all functions in the class $C^2[0,x]$ by assigning values of the last two weights.

$$\begin{array}{cc}\hfill & {\mathcal{A}}_{n}f\left(x\right)={\displaystyle \frac{h^{-\alpha }}{\Gamma (-\alpha )}}\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+O\left(h^{2-\alpha }\right),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & {\sigma }_0^{\left(\alpha \right)}=\zeta \left(\alpha \right)-\zeta (1+\alpha ),{\sigma }_1^{\left(\alpha \right)}=1-\zeta \left(\alpha \right),{\sigma }_k^{\left(\alpha \right)}={\displaystyle \frac{1}{k^{1+\alpha }}},(2\le k\le n-2),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & {\sigma }_{n-1}^{\left(\alpha \right)}={\displaystyle \frac{1}{{(n-1)}^{1+\alpha }}}-n{S}_n[1+\alpha ]+S_n\left[\alpha \right]-{\displaystyle \frac{n^{1-\alpha }}{\alpha (1-\alpha )}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & {\sigma }_n^{\left(\alpha \right)}=(n-1)S_n[1+\alpha ]-S_n\left[\alpha \right]+{\displaystyle \frac{n^{1-\alpha }}{\alpha (1-\alpha )}},\hfill \end{array}$$

(6)

where $S_n\left[\alpha \right]={\sum }_{k=1}^{n-1}{k}^{-\alpha }-\zeta \left(\alpha \right)$. In this paper we investigate the properties of approximations (2) and (3) and applications of the approximations for numerical solution of ordinary and partial fractional differential equations. The outline of the paper is as follows. In section 2 we study the properties of the weights of approximation (3) and we derive the inequality

$$\begin{array}{c}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}<\zeta (1+\alpha )n^{\alpha }-2.\end{array}$$

(7)

In section 3 we derive estimates for the errors of approximations (2) and (3). In sections 4 and 5 we construct numerical solutions of the two-term ordinary fractional differential equation and the time-fractional Black-Scholes equation for option pricing and we prove the convergence and order of the methods.

2. Properties of the Weights

L1 approximation and approximation (2) have an order $2-\alpha$ and similar performance and properties of the weights [17]

$${\sigma }_0^{\left(\alpha \right)}>0,{\sigma }_1^{\left(\alpha \right)}<{\sigma }_2^{\left(\alpha \right)}<\cdots <{\sigma }_{n-1}^{\left(\alpha \right)}<0,\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}=0.$$

The weights of ${\sigma }_{n-1}^{\left(\alpha \right)}$ and ${\sigma }_n^{\left(\alpha \right)}$ of approximation (2) are chosen such that the values of ${\mathcal{A}}_{n}1$ and ${\mathcal{A}}_{n}x$ are equal to the fractional derivatives.

$${\mathcal{A}}_{n}1=D_x^{\alpha }1=0,{\mathcal{A}}_{n}x=D_x^{\alpha }x={\displaystyle \frac{x^{1-\alpha }}{\Gamma (2-\alpha )}}.$$

In this section we prove inequality (7) and we derive an estimate for the weight ${\sigma }_{n-1}^{\left(\alpha \right)}$ of approximation (2). The proof uses the inequalities in Claim 1 and Claim 2.

Claim 1.

Let $0<\alpha <1$ and $1\le k\le n$ . Then

$$\begin{array}{c}\hfill {(n-k)}^{\alpha }<n^{\alpha }-{\displaystyle \frac{\alpha k}{n^{1-\alpha }}}.\end{array}$$

(8)

Proof. 

From the binomial formula

$$\begin{array}{c}\hfill {(n-k)}^{\alpha }=n^{\alpha }{\left(1-{\displaystyle \frac{k}{n}}\right)}^{\alpha }=n^{\alpha }\sum _{i=0}^{\infty }{(-1)}^i\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right){\left({\displaystyle \frac{k}{n}}\right)}^i.\end{array}$$

The numbers ${(-1)}^i\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)$ are negative for $k\ge 1$. Then

$${(n-k)}^{\alpha }-n^{\alpha }+{\displaystyle \frac{\alpha k}{n^{1-\alpha }}}=n^{\alpha }\sum _{i=2}^{\infty }{(-1)}^i\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right){\left({\displaystyle \frac{k}{n}}\right)}^i<0.$$

□

Claim 2.

Let $0<\alpha <1$. Then

$$\begin{array}{c}\hfill {\displaystyle \frac{1-\alpha +{\alpha }^2}{\alpha (1-\alpha )}}+{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{2^{1-\alpha }}}>2.\end{array}$$

(9)

Proof. 

The Riemann zeta function has a series expansion [36]

$$\begin{array}{c}\hfill \zeta \left(\alpha \right)={\displaystyle \frac{1}{\alpha -1}}+\gamma +\sum _{k=1}^{\infty }{\displaystyle \frac{{\gamma }_k}{k!}}{(1-\alpha )}^k,\end{array}$$

(10)

where $\gamma =0.5772$ is Euler’s constant and ${\gamma }_k$ are Stieltjes constants

$${\gamma }_1=-0.0728158,{\gamma }_2=-0.0096904,{\gamma }_3=0.0020538.$$

Then

$$\begin{array}{cc}\hfill \zeta \left(\alpha \right)& >{\displaystyle \frac{1}{\alpha -1}}+\gamma +(1-\alpha ){\gamma }_1+{\displaystyle \frac{{(1-\alpha )}^2}{2}}{\gamma }_2>{\displaystyle \frac{1}{\alpha -1}}+\gamma +(1-\alpha )\left({\gamma }_1+{\displaystyle \frac{{\gamma }_2}{2}}\right)\hfill \\ \hfill & >{\displaystyle \frac{1}{\alpha -1}}+\gamma -(1-\alpha )0.08>{\displaystyle \frac{1}{\alpha -1}}+0.49+0.08\alpha .\hfill \end{array}$$

(11)

The value of $ln2=0.693147<0.7$. From Taylor’s Theorem

$$\begin{array}{cc}\hfill 2^{-(1-\alpha )}=e^{-(1-\alpha )ln2}& <1-(1-\alpha )ln2+{(1-\alpha )}^2{\displaystyle \frac{{ln}^{2}2}{2}}\hfill \\ \hfill & <1-\left(ln2-{\displaystyle \frac{{ln}^{2}2}{2}}\right)(1-\alpha )\hfill \\ \hfill & <1-0.45(1-\alpha )=0.55+0.45\alpha .\hfill \end{array}$$

(12)

From (11) and (12)

$$\begin{array}{cc}\hfill {\displaystyle \frac{1-\alpha +{\alpha }^2}{\alpha (1-\alpha )}}+{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{2^{1-\alpha }}}& >{\displaystyle \frac{1-\alpha +{\alpha }^2}{\alpha (1-\alpha )}}+\alpha \left({\displaystyle \frac{1}{\alpha -1}}+0.497+0.08\alpha \right)(0.55+0.45\alpha )\hfill \\ \hfill & >{\displaystyle \frac{1}{\alpha }}+0.036a^3+0.26765a^2+0.72335\alpha .\hfill \end{array}$$

Inequality (9) holds when

$${\displaystyle \frac{1}{\alpha }}+0.036{\alpha }^3+0.26765{\alpha }^2+0.72335\alpha >2,$$

$$0.036a^4+0.26765a^3+0.72335{\alpha }^2-2\alpha +1>0.$$

The function $f\left(\alpha \right)=0.036a^4+0.26765a^3+0.72335{\alpha }^2-2\alpha +1$ has a first derivative

$$f^{\prime }\left(\alpha \right)=0.144{\alpha }^3+0.80295{\alpha }^2+1.4467\alpha -2$$

and a minimum value $f\left(0.882\right)=0.00414$. Therefore $f\left(\alpha \right)$ is positive on $[0,1]$. □

The asymptotic formula of ${\sum }_{k=1}^{n-1}{(n-k)}^{\alpha }{k}^{-1-\alpha }$ is obtained from the asymptotic formula of the fractional integral [18,38]

$$\begin{array}{cc}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}=\Gamma (1+\alpha )\Gamma (-\alpha )& +\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right){\displaystyle \frac{\zeta (1+\alpha -k)}{n^{k-\alpha }}}\hfill \\ \hfill & +\sum _{k=0}^{\infty }{(-1)}^k\left(\genfrac{}{}{0pt}{}{-1-\alpha }{k}\right){\displaystyle \frac{\zeta (-\alpha -k)}{n^{\alpha +k+1}}}.\hfill \end{array}$$

(13)

Let $\beta =min\{2+\alpha ,3-\alpha \}$. From the properties of gamma function

$$\Gamma (1+\alpha )\Gamma (-\alpha )=-{\displaystyle \frac{\pi }{sin\left(\pi \alpha \right)}}.$$

From (13) we find the asymptotic formula of ${\sum }_{k=1}^{n-1}{(n-k)}^{\alpha }{k}^{-1-\alpha }$ of order $\beta$

$$\begin{array}{c}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}=\zeta (1+\alpha )n^{\alpha }-{\displaystyle \frac{\pi }{sin\left(\pi \alpha \right)}}-{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{n^{1-\alpha }}}+{\displaystyle \frac{\zeta (-\alpha )}{n^{1+\alpha }}}-\alpha (1-\alpha ){\displaystyle \frac{\zeta (\alpha -1)}{2n^{2-\alpha }}}+O\left(n^{-\beta }\right).\end{array}$$

Lemma 3.

Let $n\ge 2$. Then

$$\begin{array}{c}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}<\zeta (1+\alpha )n^{\alpha }-2.\end{array}$$

(14)

Proof. 

From inequality (8)

$$\begin{array}{c}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}<\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}\left(n^{\alpha }-{\displaystyle \frac{\alpha k}{n^{1-\alpha }}}\right)<n^{\alpha }\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}-{\displaystyle \frac{\alpha }{n^{1-\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}.\end{array}$$

Sum of the $\alpha$-th powers of the first $n-1$ integers has an asymptotic formula [21]

$$\begin{array}{c}\hfill \sum _{k=1}^{n-1}{k}^{\alpha }=\zeta (-\alpha )+{\displaystyle \frac{n^{1+\alpha }}{1+\alpha }}\sum _{p=1}^{\infty }\left(\genfrac{}{}{0pt}{}{1+\alpha }{p}\right){\displaystyle \frac{b_p}{n^p}},\end{array}$$

(15)

where $b_p$ are the Bernoulli numbers. Formula (15) is derived using Euler-Maclaurin formula and the actual value of ${\sum }_{k=1}^{n-1}{k}^{\alpha }$ lies between two consecutive partial sums [21]. Then

$$\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}<\zeta (1+\alpha )-{\displaystyle \frac{1}{\alpha n^{\alpha }}}-{\displaystyle \frac{1}{2n^{1+\alpha }}}-{\displaystyle \frac{1+\alpha }{12n^{2+\alpha }}}+{\displaystyle \frac{(1+\alpha )(2+\alpha )(3+\alpha )}{720n^{4+\alpha }}},$$

$$\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}>{\displaystyle \frac{n^{1-\alpha }}{1-\alpha }}+\zeta \left(\alpha \right)-{\displaystyle \frac{1}{2n^{\alpha }}}-{\displaystyle \frac{\alpha }{12n^{1+\alpha }}}$$

and

$$\begin{array}{cc}\hfill & n^{\alpha }\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}<\zeta (1+\alpha )n^{\alpha }-{\displaystyle \frac{1}{\alpha }}-{\displaystyle \frac{1}{2n}}-{\displaystyle \frac{1+\alpha }{12n^2}}+{\displaystyle \frac{1}{30n^4}},\hfill \\ \hfill & {\displaystyle \frac{\alpha }{n^{1-\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}>{\displaystyle \frac{\alpha }{1-\alpha }}+{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{n^{1-\alpha }}}-{\displaystyle \frac{\alpha }{2n}}-{\displaystyle \frac{{\alpha }^2}{12n^2}}.\hfill \end{array}$$

Therefore

$$\begin{array}{cc}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}& <n^{\alpha }\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}-{\displaystyle \frac{\alpha }{n^{1-\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}\hfill \\ \hfill <\zeta (1& +\alpha )n^{\alpha }-{\displaystyle \frac{1-\alpha +{\alpha }^2}{\alpha (1-\alpha )}}-{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{n^{1-\alpha }}}-{\displaystyle \frac{1-\alpha }{2n}}-{\displaystyle \frac{1+\alpha -{\alpha }^2}{12n^2}}+{\displaystyle \frac{1}{30n^4}}\hfill \\ \hfill <\zeta (1& +\alpha )n^{\alpha }-{\displaystyle \frac{1-\alpha +{\alpha }^2}{\alpha (1-\alpha )}}-{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{n^{1-\alpha }}}-{\displaystyle \frac{1}{12n^2}}+{\displaystyle \frac{1}{30n^4}}\hfill \\ \hfill <\zeta (1& +\alpha )n^{\alpha }-{\displaystyle \frac{1-\alpha +{\alpha }^2}{\alpha (1-\alpha )}}-{\displaystyle \frac{\alpha \zeta \left(\alpha \right)}{n^{1-\alpha }}}.\hfill \end{array}$$

From Claim 2 it follows that

$$\begin{array}{c}\hfill \sum _{k=1}^{n-1}{\displaystyle \frac{{(n-k)}^{\alpha }}{k^{1+\alpha }}}<\zeta (1+\alpha )n^{\alpha }-2.\end{array}$$

□

Now we derive an estimate for the weight ${\sigma }_{n-1}^{\left(\alpha \right)}$ of approximation (3). Denote

$${\overline{\sigma }}_{n-1}^{\left(\alpha \right)}={\sigma }_{n-1}^{\left(\alpha \right)}-{\displaystyle \frac{1}{{(n-1)}^{1+\alpha }}}=S_n\left[\alpha \right]-n{S}_n[1+\alpha ]-{\displaystyle \frac{n^{1-\alpha }}{\alpha (1-\alpha )}}.$$

Claim 4.

Let $n\ge 2$. Then

$$\left|{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}\right|<{\displaystyle \frac{1+2\alpha }{12n^{1+\alpha }}}.$$

Proof. 

From (15)

$$\begin{array}{c}\hfill S_n\left[\alpha \right]={\displaystyle \frac{n^{1-\alpha }}{1-\alpha }}-{\displaystyle \frac{1}{2n^{\alpha }}}+R_1,n{S}_n[\alpha +1]=-{\displaystyle \frac{n^{1-\alpha }}{\alpha }}-{\displaystyle \frac{1}{2n^{\alpha }}}+R_2,\end{array}$$

where

$$|R_1|<{\displaystyle \frac{1+\alpha }{12n^{1+\alpha }}},\left|R_2\right|<{\displaystyle \frac{\alpha }{12n^{1+\alpha }}}.$$

Then

$$\begin{array}{cc}\hfill \left|{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}\right|& =\left|S_n\left[\alpha \right]-n{S}_n[1+\alpha ]-{\displaystyle \frac{n^{1-\alpha }}{\alpha (1-\alpha )}}\right|\hfill \\ \hfill & <|R_1|+|R_2|<{\displaystyle \frac{1+2\alpha }{12n^{1+\alpha }}}\hfill \end{array}$$

□

3. Estimate for the Error

In this section we use the method from [38] to derive estimates for the errors of approximations (2) and (3). Let $0<t<x$ and

$$\begin{array}{c}\hfill g\left(t\right)=f\left(t\right)-f\left(x\right)+(x-t)f^{\prime }\left(x\right)-{\displaystyle \frac{{(x-t)}^2}{2}}{f}^{\prime \prime }\left(x\right).\end{array}$$

(16)

The function $g\left(t\right)$ satisfies $g\left(x\right)=g^{\prime }\left(x\right)=g^{\prime \prime }\left(x\right)=0$. From Taylor’s Theorem

$$g\left(t\right)={\displaystyle \frac{{(t-x)}^3}{6}}{f}^{\prime \prime \prime }\left({\xi }_{0,t}\right),g^{\prime }\left(t\right)={\displaystyle \frac{{(t-x)}^2}{2}}{f}^{\prime \prime \prime }\left({\xi }_{1,t}\right),g^{\prime \prime }\left(t\right)=(t-x)f^{\prime \prime \prime }\left({\xi }_{2,t}\right)$$

where $t<{\xi }_{i,t}<x$ for $i=0,1,2$. Let

$$M_i=\underset{0<t<x}{max}\left|f^{\left(i\right)}\left(t\right)\right|,(i=0,1,2,3).$$

Then

$$\left|g\left(t\right)\right|<{\displaystyle \frac{M_3}{6}}{(x-t)}^3,|g^{\prime }\left(t\right)|<{\displaystyle \frac{M_3}{2}}{(x-t)}^2,\left|g^{\prime \prime }\left(t\right)\right|<M_3(x-t).$$

Denote $G\left(t\right)=g\left(t\right)/{(x-t)}^{1+\alpha }$. The function $G\left(t\right)$ satisfies $G\left(x\right)=G^{\prime }\left(x\right)=0$ and its second derivative of is not defined at the point x.

Claim 5.

Let $0<t<x$. Then

$$|G^{\prime }\left(t\right)|<{\displaystyle \frac{5}{6}}{M}_3{(x-t)}^{1-\alpha },\left|G^{\prime \prime }\left(t\right)\right|<{\displaystyle \frac{4M_3}{{(x-t)}^{\alpha }}}.$$

Proof. 

$$G^{\prime }\left(t\right)={\displaystyle \frac{(x-t)g^{\prime }\left(x\right)+(1+\alpha )g\left(t\right)}{{(x-t)}^{2+\alpha }}},$$

$$G^{\prime \prime }\left(t\right)={\displaystyle \frac{{(x-t)}^2{g}^{\prime \prime }\left(t\right)+2(\alpha +1)(x-t)g^{\prime }\left(t\right)+(\alpha +1)(\alpha +2)g\left(t\right)}{{(x-t)}^{3+\alpha }}}.$$

Therefore

$$\begin{array}{cc}\hfill |G^{\prime }\left(t\right)|& <{\displaystyle \frac{1}{2}}{M}_3{(x-t)}^{1-\alpha }+{\displaystyle \frac{1+\alpha }{6}}{M}_3{(x-t)}^{1-\alpha }<{\displaystyle \frac{5}{6}}{M}_3{(x-t)}^{1-\alpha },\hfill \\ \hfill |G^{\prime \prime }\left(t\right)|& <\left((2+\alpha )(1+\alpha ){\displaystyle \frac{M_3}{6}}+(1+\alpha )M_3+M_3\right){(x-t)}^{-\alpha }<{\displaystyle \frac{4M_3}{{(x-t)}^{\alpha }}}.\hfill \end{array}$$

□

Lemma 6.

Let $f\left(0\right)=f^{\prime }\left(0\right)=0$. Then

$${\int }_0^{x}G\left(t\right)dt=\Gamma (-\alpha )f^{\left(\alpha \right)}\left(x\right)+{\displaystyle \frac{f\left(x\right)}{\alpha x^{\alpha }}}+{\displaystyle \frac{f^{\prime }\left(x\right)x^{1-\alpha }}{1-\alpha }}-{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^{2-\alpha }}{2(2-\alpha )}}.$$

Proof. 

$$\begin{array}{cc}\hfill {\int }_0^{x}G\left(t\right)dt& ={\int }_0^x{\displaystyle \frac{g\left(t\right)}{{(x-t)}^{1+\alpha }}}dt={\int }_0^{x}g\left(t\right)d{\displaystyle \frac{{(x-t)}^{-\alpha }}{\alpha }}\hfill \\ \hfill & =g\left(t\right){\left.{\displaystyle \frac{{(x-t)}^{-\alpha }}{\alpha }}\right|}_0^x-{\displaystyle \frac{1}{a}}{\int }_0^x{\displaystyle \frac{g^{\prime }\left(t\right)}{{(x-t)}^{\alpha }}}dt\hfill \\ \hfill & =-{\displaystyle \frac{g\left(0\right)}{\alpha x^{\alpha }}}+{\displaystyle \frac{\Gamma (-\alpha )}{\Gamma (1-\alpha )}}{\int }_0^x{\displaystyle \frac{g^{\prime }\left(t\right)}{{(x-t)}^{\alpha }}}dt.\hfill \end{array}$$

The function $g\left(t\right)$ satisfies

$$g\left(0\right)=-f\left(x\right)+f^{\prime }\left(x\right)x-{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{x}^2.$$

Hence

$$\begin{array}{c}\hfill {\int }_0^x{\displaystyle \frac{g\left(t\right)}{{(x-t)}^{1+\alpha }}}dt=\Gamma (-\alpha )g^{\left(\alpha \right)}\left(x\right)+{\displaystyle \frac{1}{\alpha }}\left({\displaystyle \frac{f\left(x\right)}{x^{\alpha }}}-x^{1-\alpha }{f}^{\prime }\left(x\right)+{\displaystyle \frac{x^{2-\alpha }}{2}}{f}^{\prime \prime }\left(x\right)\right).\end{array}$$

(17)

The fractional derivative of the function $g\left(x\right)$ satisfies

$$\begin{array}{cc}\hfill g^{\left(\alpha \right)}\left(x\right)& =f^{\left(\alpha \right)}\left(x\right)-{\displaystyle \frac{f^{\prime }\left(x\right)}{\Gamma (1-\alpha )}}{\int }_0^x{(x-t)}^{-\alpha }dt+{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2\Gamma (1-\alpha )}}{\int }_0^{x}2{(x-t)}^{1-\alpha }dt\hfill \\ \hfill & =f^{\left(\alpha \right)}\left(x\right)-{\displaystyle \frac{x^{1-\alpha }{f}^{\prime }\left(x\right)}{(1-\alpha )\Gamma (1-\alpha )}}+{\displaystyle \frac{x^{2-\alpha }{f}^{\prime \prime }\left(x\right)}{(2-\alpha )\Gamma (1-\alpha )}},\hfill \end{array}$$

$$\begin{array}{c}\hfill \Gamma (-\alpha )g^{\left(\alpha \right)}\left(x\right)=\Gamma (-\alpha )f^{\left(\alpha \right)}\left(x\right)+{\displaystyle \frac{x^{1-\alpha }{f}^{\prime }\left(x\right)}{\alpha (1-\alpha )}}-{\displaystyle \frac{x^{2-\alpha }{f}^{\prime \prime }\left(x\right)}{\alpha (2-\alpha )}}.\end{array}$$

(18)

From (17) and (18) it follows that

$${\int }_0^x{\displaystyle \frac{g\left(t\right)}{{(x-t)}^{1+\alpha }}}dt=\Gamma (-\alpha )f^{\left(\alpha \right)}\left(x\right)+{\displaystyle \frac{f\left(x\right)}{\alpha x^{\alpha }}}+{\displaystyle \frac{f^{\prime }\left(x\right)x^{1-\alpha }}{1-\alpha }}-{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^{2-\alpha }}{2(2-\alpha )}}.$$

□

The trapezoidal rule of a function $F\left(t\right)$ has a second order accuracy when $F\in C^2[0,x]$. The error of the trapezoidal rule $E_T^F$ satisfies [56]

$$|E_T^F|<{\displaystyle \frac{1}{12}}{h}^2{x}^3\underset{0<t<x}{max}\left|F^{\prime \prime }\left(t\right)\right|.$$

The second derivative of the function $G\left(t\right)$ is undefined at the point x, which leads to a lower order of accuracy of the trapezoidal rule in the interval $[0,x]$. Now we estimate the error of the trapezoidal rule of the function $G\left(t\right)$. Let $h=x/n$ and $T_n^1$ be the error of trapezoidal rule of $G\left(t\right)$ in the interval $[0,x-h]$.

$$\begin{array}{c}\hfill h\left({\displaystyle \frac{G\left(0\right)}{2}}+\sum _{k=1}^{n-2}G\left(kh\right)+{\displaystyle \frac{G(x-h)}{2}}\right)={\int }_0^{x-h}G\left(t\right)dt+T_n^1.\end{array}$$

(19)

Claim 7.

$$|T_n^1|<{\displaystyle \frac{1}{3}}{M}_3{x}^3{h}^{2-\alpha }.$$

Proof. 

The error $T_n^1$ of trapezoidal rule (19) satisfies

$$\begin{array}{cc}\hfill \left|T_n^1\right|& \le {\displaystyle \frac{1}{12}}{h}^2{(x-h)}^3\underset{0<t<x-h}{max}\left|G^{\prime \prime }\left(t\right)\right|\hfill \\ \hfill & \le {\displaystyle \frac{1}{12}}{h}^2{x}^3{\displaystyle \frac{4M_3}{h^{\alpha }}}\le {\displaystyle \frac{1}{3}}{M}_3{x}^3{h}^{2-\alpha }.\hfill \end{array}$$

□

Let $T_n^2$ be the error of the trapezoidal rule of $G\left(t\right)$ in the interval $[x-h,x]$.

$$\begin{array}{c}\hfill {\displaystyle \frac{hG(x-h)}{2}}={\int }_{x-h}^{x}G\left(t\right)dt+T_n^2.\end{array}$$

(20)

Claim 8.

$$|T_n^2|<{\displaystyle \frac{5}{24}}{M}_3{h}^{3-\alpha }.$$

Proof. 

$${\int }_{x-h}^{x}G\left(t\right)dt={\int }_{x-h}^{x}G\left(t\right)d\left(t+{\displaystyle \frac{h}{2}}-x\right)={\displaystyle \frac{h}{2}}G(x-h)-{\int }_{x-h}^x\left(t+{\displaystyle \frac{h}{2}}-x\right)G^{\prime }\left(t\right)dt.$$

The error $T_n^2$ satisfies

$$T_n^2={\int }_{x-h}^x\left(t+{\displaystyle \frac{h}{2}}-x\right)G^{\prime }\left(t\right)dt,$$

$$|T_n^2|<{\int }_{x-h}^x\left|t+{\displaystyle \frac{h}{2}}-x\right|\left|G^{\prime }\left(t\right)\right|dt<{\displaystyle \frac{5}{6}}{M}_3{h}^{1-\alpha }{\int }_{x-h}^x\left|t+{\displaystyle \frac{h}{2}}-x\right|dt<{\displaystyle \frac{5}{24}}{M}_3{h}^{3-\alpha }.$$

□

Let $T_n=T_n^1+T_n^2$ be the error of the trapezoidal rule of $G\left(t\right)$ in the interval $[0,x]$.

Lemma 9.

$$|T_n|\le {\displaystyle \frac{1}{3}}{M}_{3}x(x^2+1)h^{2-\alpha }.$$

Proof. 

From Claim 7 and Claim 8

$$\begin{array}{cc}\hfill |T_n|\le |T_n^1|+|T_n^2|& <{\displaystyle \frac{1}{3}}{M}_3{x}^3{h}^{2-\alpha }+{\displaystyle \frac{5}{24}}{M}_3{h}^{3-\alpha }\hfill \\ \hfill & <{\displaystyle \frac{1}{3}}{M}_3\left(x^3+h\right)h^{2-\alpha }<{\displaystyle \frac{1}{3}}{M}_{3}x(x^2+1)h^{2-\alpha }.\hfill \end{array}$$

□

Let $E_n$ be the error of approximation (1).

$${\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{f_{n-k}}{k^{1+\alpha }}}=\Gamma (-\alpha )f_n^{\left(\alpha \right)}+{\displaystyle \frac{\zeta (1+\alpha )}{h^{\alpha }}}{f}_n-\zeta \left(\alpha \right)f_n^{\prime }{h}^{1-\alpha }+E_n.$$

Theorem 10.

Let $f\left(0\right)=f^{\prime }\left(0\right)=0$. Then

$$|E_n|<{\displaystyle \frac{2M_{3}x(x^2+1)+M_2}{6}}{h}^{2-\alpha }.$$

Proof. 

The trapezoidal rule of the function $G\left(t\right)$ on the interval $[0,x]$ satisfies

$${\displaystyle \frac{1}{h^{\alpha }}}\left(\sum _{k=1}^{n-1}{\displaystyle \frac{g_{n-k}}{k^{1+\alpha }}}+{\displaystyle \frac{g_0}{2n^{1+\alpha }}}\right)={\int }_0^x{\displaystyle \frac{g\left(t\right)}{{(x-t)}^{1+\alpha }}}dt+T_n$$

and the function $g\left(t\right)$ has a value at zero

$$\begin{array}{cc}\hfill {\displaystyle \frac{g\left(0\right)}{2n^{1+\alpha }{h}^{\alpha }}}& =-{\displaystyle \frac{f\left(x\right)}{2n^{1+\alpha }{h}^{\alpha }}}+{\displaystyle \frac{f^{\prime }\left(x\right)x}{2n^{1+\alpha }{h}^{\alpha }}}-{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^2}{4n^{1+\alpha }{h}^{\alpha }}}={\displaystyle \frac{h}{2}}\left(-{\displaystyle \frac{f\left(x\right)}{x^{1+\alpha }}}+{\displaystyle \frac{f^{\prime }\left(x\right)}{x^{\alpha }}}-{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^{1-\alpha }}{2}}\right).\hfill \end{array}$$

Then

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{g_{n-k}}{k^{1+\alpha }}}={\int }_0^x{\displaystyle \frac{g\left(t\right)}{{(x-t)}^{1+\alpha }}}d& t+\left({\displaystyle \frac{f\left(x\right)}{x^{1+\alpha }}}-{\displaystyle \frac{f^{\prime }\left(x\right)}{x^{\alpha }}}+{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^{1-\alpha }}{2}}\right){\displaystyle \frac{h}{2}}+T_n\hfill \\ \hfill =\Gamma (-\alpha )f^{\left(\alpha \right)}\left(x\right)+& {\displaystyle \frac{f\left(x\right)}{\alpha x^{\alpha }}}+{\displaystyle \frac{f\left(x\right)h}{2x^{1+\alpha }}}+{\displaystyle \frac{x^{1-\alpha }{f}^{\prime }\left(x\right)}{1-\alpha }}-{\displaystyle \frac{f^{\prime }\left(x\right)h}{2x^{\alpha }}}\hfill \\ \hfill +& {\displaystyle \frac{x^{1-\alpha }{f}^{\prime \prime }\left(x\right)h}{4}}-{\displaystyle \frac{x^{2-\alpha }{f}^{\prime \prime }\left(x\right)}{2(2-\alpha )}}+T_n.\hfill \end{array}$$

From the definition (16) of the function $g\left(t\right)$

$${\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{g_{n-k}}{k^{1+\alpha }}}={\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{f_{n-k}}{k^{1+\alpha }}}-{\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}+f^{\prime }\left(x\right)h^{1-\alpha }\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}-{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }\sum _{k=1}^{n-1}{k}^{1-\alpha }.$$

From (15)

$$\begin{array}{cc}\hfill & \sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}=\zeta (1+\alpha )-{\displaystyle \frac{1}{\alpha n^{\alpha }}}-{\displaystyle \frac{1}{2n^{1+\alpha }}}+R_1,\hfill \\ \hfill & \sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}={\displaystyle \frac{n^{1-\alpha }}{1-\alpha }}+\zeta \left(\alpha \right)-{\displaystyle \frac{1}{2n^{\alpha }}}+R_2,\hfill \\ \hfill & \sum _{k=1}^{n-1}{k}^{1-\alpha }={\displaystyle \frac{n^{2-\alpha }}{2-\alpha }}-{\displaystyle \frac{n^{1-\alpha }}{2}}+\zeta (\alpha -1)+R_3,\hfill \end{array}$$

where

$$|R_1|<{\displaystyle \frac{1+\alpha }{12n^{2+\alpha }}},|R_2|<{\displaystyle \frac{\alpha }{12n^{1+\alpha }}},\left|R_3\right|<{\displaystyle \frac{1-\alpha }{12n^{\alpha }}}.$$

Then

$$\begin{array}{cc}\hfill {\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{1+\alpha }}}& =\zeta (1+\alpha ){\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}-{\displaystyle \frac{f\left(x\right)}{\alpha n^{\alpha }{h}^{\alpha }}}-{\displaystyle \frac{f\left(x\right)}{2n^{1+\alpha }{h}^{\alpha }}}+R_1{\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}\hfill \\ \hfill & =\zeta (1+\alpha ){\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}-{\displaystyle \frac{f\left(x\right)}{\alpha x^{\alpha }}}-{\displaystyle \frac{f\left(x\right)h}{2x^{1+\alpha }}}+R_1{\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}.\hfill \end{array}$$

Similarly

$$\begin{array}{c}\hfill f^{\prime }\left(x\right)h^{1-\alpha }\sum _{k=1}^{n-1}{\displaystyle \frac{1}{k^{\alpha }}}=f^{\prime }\left(x\right)h^{1-\alpha }\zeta \left(\alpha \right)+{\displaystyle \frac{f^{\prime }\left(x\right)x^{1-\alpha }}{1-\alpha }}-{\displaystyle \frac{f^{\prime }\left(x\right)h}{2x^{\alpha }}}+f^{\prime }\left(x\right)h^{1-\alpha }{R}_2,\end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }\sum _{k=1}^{n-1}{k}^{1-\alpha }=\zeta (\alpha -1){\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }+{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^{2-\alpha }}{2-\alpha }}-{\displaystyle \frac{f^{\prime \prime }\left(x\right)x^{1-\alpha }h}{4}}+R_3{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }.\end{array}$$

Hence

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{h^{\alpha }}}\sum _{k=1}^{n-1}{\displaystyle \frac{f_{n-k}}{k^{1+\alpha }}}=\Gamma (-\alpha )f^{\left(\alpha \right)}\left(x\right)+{\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}& \zeta (1+\alpha )-f^{\prime }\left(x\right)h^{1-\alpha }\zeta \left(\alpha \right)\hfill \\ \hfill & +{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }\zeta (\alpha -1)+E_n,\hfill \end{array}$$

where

$$E_n=T_n+R_1{\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}-f^{\prime }\left(x\right)h^{1-\alpha }{R}_2+R_3{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha },$$

$$\begin{array}{c}\hfill |E_n|<|T_n|+|R_1{\displaystyle \frac{f\left(x\right)}{h^{\alpha }}}|+|f^{\prime }\left(x\right)h^{1-\alpha }{R}_2|+|R_3{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }|.\end{array}$$

(21)

From Taylor’s Theorem

$$f\left(x\right)={\displaystyle \frac{1}{2}}{f}^{\prime \prime }\left({\xi }_1\right)x^2,f^{\prime }\left(x\right)=f^{\prime \prime }\left({\xi }_2\right)x,$$

where $0<{\xi }_2<x$ for $i=1,2$. Therefore

$$M_0\le {\displaystyle \frac{1}{2}}{M}_2{x}^2,M_1\le M_{2}x.$$

The terms on the right hand side of (21) satisfy the estimates

$$\begin{array}{cc}\hfill & |R_{1}f\left(x\right)h^{-\alpha }|<{\displaystyle \frac{(1+\alpha )M_0{h}^{-\alpha }}{12n^{2+\alpha }}}={\displaystyle \frac{(1+\alpha )M_0{h}^2}{12x^{2+\alpha }}}<{\displaystyle \frac{(1+\alpha )M_2{h}^{2-\alpha }}{24}},\hfill \\ \hfill & |R_2{f}^{\prime }\left(x\right)h^{1-\alpha }|<{\displaystyle \frac{\alpha M_1{h}^{1-\alpha }}{12m^{1+\alpha }}}={\displaystyle \frac{\alpha M_1{h}^2}{12x^{1+\alpha }}}<{\displaystyle \frac{\alpha M_2{h}^2}{12x^{\alpha }}}<{\displaystyle \frac{\alpha M_2{h}^{2-\alpha }}{12}},\hfill \\ \hfill & |R_3{\displaystyle \frac{f^{\prime \prime }\left(x\right)}{2}}{h}^{2-\alpha }|<{\displaystyle \frac{(1-\alpha )M_2{h}^{2-\alpha }}{24n^{\alpha }}}={\displaystyle \frac{(1-\alpha )M_2{h}^2}{24x^{\alpha }}}<{\displaystyle \frac{(1-\alpha )M_2{h}^{2-\alpha }}{24}}.\hfill \end{array}$$

Hence

$$\begin{array}{cc}\hfill |E_n|& <{\displaystyle \frac{1}{3}}{M}_{3}x(x^2+1)h^{2-\alpha }+{\displaystyle \frac{(1+\alpha )M_2{h}^{2-\alpha }}{12}}\hfill \\ \hfill & <{\displaystyle \frac{1}{3}}{M}_{3}x(x^2+1)h^{2-\alpha }+{\displaystyle \frac{M_2{h}^{2-\alpha }}{6}}={\displaystyle \frac{2M_{3}x(x^2+1)+M_2}{6}}{h}^{2-\alpha }.\hfill \end{array}$$

□

Let $L_n$ be the error of approximation (2).

$${\mathcal{L}}_{n}f\left(x\right)={\displaystyle \frac{h^{-\alpha }}{\Gamma (-\alpha )}}\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+L_n{h}^{2-\alpha },$$

$${\sigma }_0^{\left(\alpha \right)}=\zeta \left(\alpha \right)-\zeta (1+\alpha ),{\sigma }_1^{\left(\alpha \right)}=1-\zeta \left(\alpha \right),{\sigma }_k^{\left(\alpha \right)}={\displaystyle \frac{1}{k^{1+\alpha }}},(2\le k\le n).$$

Approximation (2) is constructed from approximation (1) by substituting $f_n^{\prime }$ with first order backward difference approximation $(f_n-f_{n-1})/h$.

Claim 11.

Let $f\in C^3[0,x]$ and $f\left(0\right)=f^{\prime }\left(0\right)=0$. Then

$$|L_n|<{\displaystyle \frac{2M_{3}x(x^2+1)+(1+3|\zeta \left(\alpha \right)\left|\right)M_2}{6|\Gamma (-\alpha \left)\right|}}{h}^{2-\alpha }.$$

Proof. 

From Taylor’s Theorem

$$f_{n-1}=f_n-h{f}_n^{\prime }+{\displaystyle \frac{h^2}{2}}{f}^{\prime \prime }\left(\xi \right),$$

$$\left|f_n^{\prime }-{\displaystyle \frac{f_n-f_{n-1}}{h}}\right|<{\displaystyle \frac{h}{2}}\left|f^{\prime \prime }\left(\xi \right)\right|,$$

where $\xi \in (x_{n-1},x_n)$. The error $L_n$ of approximation (2) satisfies

$$L_n={\displaystyle \frac{1}{\Gamma (-\alpha )}}\left(E_n-{\displaystyle \frac{1}{2}}\zeta \left(\alpha \right)f^{\prime \prime }\left(\xi \right)h^{2-\alpha }\right),$$

$$\begin{array}{cc}\hfill |L_n|& <{\displaystyle \frac{1}{|\Gamma (-\alpha \left)\right|}}\left({\displaystyle \frac{2M_{3}x(x^2+1)+M_2}{6}}{h}^{2-\alpha }+{\displaystyle \frac{\left|\zeta \left(\alpha \right)\right|M_2}{2}}{h}^{2-\alpha }\right)\hfill \\ \hfill & <{\displaystyle \frac{2M_{3}x(x^2+1)+(1+3|\zeta \left(\alpha \right)\left|\right)M_2}{6|\Gamma (-\alpha \left)\right|}}{h}^{2-\alpha }.\hfill \end{array}$$

□

Let $A_n$ be the error of approximation (3).

$$\begin{array}{c}\hfill {\mathcal{A}}_{n}f\left(x\right)={\displaystyle \frac{h^{-\alpha }}{\Gamma (-\alpha )}}\sum _{k=1}^n{\sigma }_k^{\left(\alpha \right)}{f}_{n-k}=f_n^{\left(\alpha \right)}+A_n.\end{array}$$

The weights of approximation (3) are defined with (4), (5) and (6). Denote

$$\tilde{f}\left(t\right)=f\left(t\right)-f\left(0\right)-f^{\prime }\left(0\right)t.$$

The function $\tilde{f}$ satisfies $\tilde{f}\left(0\right)={\tilde{f}}^{\prime }\left(0\right)=0$ and the second and third derivatives of f and $\tilde{f}$ are equal. Therefore

$$\underset{0\le t\le x}{max}|{\tilde{f}}^{\prime \prime }\left(t\right)|=M_2,\underset{0\le t\le x}{max}\left|{\tilde{f}}^{\prime \prime \prime }\left(t\right)\right|=M_3.$$

Lemma 12.

Let $f\in C^3[0,x]$. Then

$$|A_n|<{\displaystyle \frac{4M_{3}x(x^2+1)+3(1+2|\zeta \left(\alpha \right)\left|\right)M_2}{12|\Gamma (-\alpha \left)\right|}}{h}^{2-\alpha }.$$

Proof. 

$$\begin{array}{cc}\hfill {\mathcal{A}}_{n}f\left(x\right)& ={\mathcal{A}}_n\tilde{f}\left(x\right)+{\mathcal{A}}_n[f\left(x\right)-\tilde{f}\left(x\right)]\hfill \\ \hfill & ={\mathcal{L}}_n\tilde{f}\left(x\right)+{\displaystyle \frac{{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}{\tilde{f}}_1}{\Gamma (-\alpha )h^{\alpha }}}+D_x^{\alpha }[f\left(x\right)-\tilde{f}\left(x\right)].\hfill \end{array}$$

From Claim 4

$$\left|{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}\right|<{\displaystyle \frac{1+2\alpha }{12n^{1+\alpha }}}<{\displaystyle \frac{1}{4n}}<{\displaystyle \frac{1}{8}}.$$

From Taylor’s Theorem

$$|{\tilde{f}}_1|={\displaystyle \frac{h^2}{2}}\left|{\tilde{f}}^{\prime \prime }\left(\xi \right)\right|\le {\displaystyle \frac{h^2{M}_2}{2}}.$$

Then

$$\begin{array}{cc}\hfill |A_n|& <|L_n|+{\displaystyle \frac{1}{|\Gamma (-\alpha )|h^{\alpha }}}\left|{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}\right|\left|{\tilde{f}}_1\right|\hfill \\ \hfill & <\left({\displaystyle \frac{2M_{3}x(x^2+1)+(1+3|\zeta \left(\alpha \right)\left|\right)M_2}{6|\Gamma (-\alpha \left)\right|}}+{\displaystyle \frac{M_2}{16\Gamma (-\alpha )}}\right)h^{2-\alpha }\hfill \\ \hfill & <{\displaystyle \frac{4M_{3}x(x^2+1)+3(1+2|\zeta \left(\alpha \right)\left|\right)M_2}{12|\Gamma (-\alpha \left)\right|}}{h}^{2-\alpha }.\hfill \end{array}$$

□

4. Numerical Solution of Two-Term Equation

The two-term equation is an ordinary fractional differential equation in the form

$$\begin{array}{c}\hfill y^{\left(\alpha \right)}\left(t\right)+Dy\left(t\right)=F\left(t\right),y\left(0\right)=y_0.\end{array}$$

(22)

Numerical and analytical solutions of ordinary fractional differential equations are studied in [6,24,25,41,43,45]. In this section we construct the numerical solution of the two-term equation (22) which uses approximation (3) of the fractional derivative and prove its convergence and order. Let $h=x/N$ and $x_n=nh$ for $n=0,\cdots ,N$. By approximating the fractional derivative at the point $x_n$ we find

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\Gamma (-\alpha )h^{\alpha }}}\sum _{k=0}^n{\sigma }_k^{\left(\alpha \right)}{y}_{n-k}+D{y}_n=F_n+A_n,\end{array}$$

$$\begin{array}{c}\hfill \left({\sigma }_0^{\left(\alpha \right)}+\Gamma (-\alpha )D{h}^{\alpha }\right)y_n+\sum _{k=1}^n{\sigma }_k^{\left(\alpha \right)}{y}_{n-k}+D{y}_n=\Gamma (-\alpha )h^{\alpha }{F}_n+B_n{h}^2,\end{array}$$

where $B_n=\Gamma (-\alpha )A_n{h}^{\alpha -2}$.

$$\begin{array}{c}\hfill \left(\zeta \left(\alpha \right)-\zeta (1+\alpha )+\Gamma (-\alpha )D{h}^{\alpha }\right)y_n-\zeta \left(\alpha \right)y_{n-1}+\sum _{k=1}^n{\displaystyle \frac{y_{n-k}}{k^{1+\alpha }}}-{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}{u}_1=\Gamma (-\alpha )h^{\alpha }{F}_n+B_n{h}^2.\end{array}$$

The numerical solution ${\left\{u_n\right\}}_{n=2}^N$ of equation (22) is computed as

$$\begin{array}{c}\hfill \left(\zeta \left(\alpha \right)-\zeta (1+\alpha )+\Gamma (-\alpha )D{h}^{\alpha }\right)u_n-\zeta \left(\alpha \right)u_{n-1}+\sum _{k=1}^n{\displaystyle \frac{u_{n-k}}{k^{1+\alpha }}}+{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}{u}_1+D{u}_n=\Gamma (-\alpha )h^{\alpha }{F}_n,\end{array}$$

$$\begin{array}{c}\hfill u_n={\displaystyle \frac{1}{\zeta \left(\alpha \right)-\zeta (1+\alpha )+\Gamma (-\alpha )D{h}^{\alpha }}}\left(\Gamma (-\alpha )h^{\alpha }{F}_n+\zeta \left(\alpha \right)u_{n-1}-\sum _{k=1}^n{\displaystyle \frac{u_{n-k}}{k^{1+\alpha }}}-{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}{u}_1\right).\end{array}$$

(23)

The value of the numerical solution $u_1$ on the first step is computed with the following approximation [17]

$$\begin{array}{c}\hfill y_1^{\left(\alpha \right)}={\displaystyle \frac{y_1-y_0}{\Gamma (2-\alpha )h^{\alpha }}}+a_1{h}^{2-\alpha }.\end{array}$$

(24)

From the estimate for the error of L1 approximation [8]

$$\begin{array}{c}\hfill |a_1|<{\displaystyle \frac{1}{\Gamma (2-\alpha )}}\left({\displaystyle \frac{2+\alpha }{2^{\alpha }(2-\alpha )}}-{\displaystyle \frac{11+\alpha }{12}}\right)\underset{0\le t\le h}{max}\left|y^{\prime \prime }\left(t\right)\right|<{\displaystyle \frac{2M_2}{\Gamma (2-\alpha )}},\end{array}$$

(25)

where $M_2={max}_{0\le t\le x}\left|y^{\prime \prime }\left(t\right)\right|$. Numerical solution (23) has initial conditions [17]

$$u_0=y\left(0\right),u_1={\displaystyle \frac{y\left(0\right)+\Gamma (2-\alpha )F\left(h\right)h^{\alpha }}{1+D\Gamma (2-\alpha )h^{\alpha }}}.$$

When the solution of two-term equation (22) satisfies the condition $y\left(0\right)=y^{\prime }\left(0\right)=0$ both approximations (2) and (3) can be used for computation of numerical solution (23).

Example 1: Consider the following two-term equation

$$y^{\left(\alpha \right)}+Dy=-{\alpha }^2{t}^{2-\alpha }{E}_{2,3-\alpha }\left(-{\left(\alpha t\right)}^2\right)+D(cos\left(\alpha t\right)-1),y_0=0.$$

(26)

Two-term equation (26) has a solution $y\left(x\right)=cos\left(\alpha t\right)-1$ which satisfies $y\left(0\right)=y^{\prime }\left(0\right)=0$. Experimental results for the error and order of numerical solution (23) of equation (26) which uses approximation (2) of the fractional derivative are given in Table 1 and Table 2.

Example 2:

$$y^{\left(\alpha \right)}+Dy=\alpha t^{1-\alpha }{E}_{1,2-\alpha }\left(\alpha t\right)+D{e}^{\alpha t},y_0=1.$$

(27)

Equation (27) has a solution $y\left(t\right)=e^{\alpha t}$. Experimental results for the error and order of numerical solution (23) of equation (27) are given in Table 3 and Table 4. The experimental results in Table 1, Table 2, Table 3 and Table 4 suggest that numerical solution (23) of two-term equation (22) converges and has an order $2-\alpha$ for all positive values of the parameter D and negative values with small and large modulus.

Denote by $e_n=y_n-u_n$ the error of numerical solution (23) at the point $x_n=nh$, where $n=0,1,\cdots ,N$. The errors $e_n$ satisfy

$$\begin{array}{c}\hfill e_n={\displaystyle \frac{1}{\zeta \left(\alpha \right)-\zeta (1+\alpha )+\Gamma (-\alpha )D{h}^{\alpha }}}\left(\zeta \left(\alpha \right)e_{n-1}-\sum _{k=1}^{n-1}{\displaystyle \frac{e_{n-k}}{k^{1+\alpha }}}-{\overline{\sigma }}_{n-1}^{\left(\alpha \right)}{e}_1+B_n{h}^2\right),\end{array}$$

(28)

$$e_0=0,e_1={\displaystyle \frac{\Gamma (2-\alpha )a_1{h}^2}{1+D\Gamma (2-\alpha )h^{\alpha }}}.$$

In Theorem 13 and Theorem 14 we prove the convergence of numerical solution (23) of two-term equation (22) and derive estimates for the error, depending on the value of the parameter D. The proofs use induction on n, where $n=1,2,\cdots ,N$.

Theorem 13.

Let $D\in \left(-\infty ,-\frac{2}{h^{\alpha }}\frac{\zeta \left(\alpha \right)-\zeta (1+\alpha )}{\Gamma (-\alpha )}\right)\cup \left(0,\infty \right)$ and $y\in C^2[0,1]$. Then

$$\begin{array}{c}\hfill |e_n|<C{h}^2{n}^{\alpha },(1\le n\le N),\end{array}$$

(29)

where ${\displaystyle C=2M_2+{\displaystyle \frac{1}{2}}\underset{1\le n\le N}{max}\left|B_n\right|}$.

Proof. 

The value of the error $e_1$ on the first step satisfies

$$\begin{array}{c}\hfill |e_1|={\displaystyle \frac{\Gamma (2-\alpha )|a_1|h^2}{|1+D\Gamma (2-\alpha )h^{\alpha }|}}<{\displaystyle \frac{2M_2{h}^2}{|1+D\Gamma (2-\alpha )h^{\alpha }|}}.\end{array}$$

(30)

When $D>0$ the denominator $|1+D\Gamma (2-\alpha )h^{\alpha }|>1$ and

$$|e_1|<2M_2{h}^2<C{h}^2.$$

Now we estimate the error $e_1$ when the parameter $D<-\frac{2}{h^{\alpha }}\frac{\zeta \left(\alpha \right)-\zeta (1+\alpha )}{\Gamma (-\alpha )}$.

$$D\Gamma (-\alpha )h^{\alpha }>2\left(\right|\zeta \left(\alpha \right)|+\zeta (1+\alpha )),$$

$$\left|D\right|\Gamma (2-\alpha )h^{\alpha }>2\alpha (1-\alpha )\left(\right|\zeta \left(\alpha \right)|+\zeta (1+\alpha )).$$

From (10)

$$\zeta (1+\alpha )>{\displaystyle \frac{1}{\alpha }}+\gamma ,\left|\zeta \left(\alpha \right)\right|>{\displaystyle \frac{1}{1-\alpha }}-\gamma ,$$

$$\left|\zeta \left(\alpha \right)\right|+\zeta (1+\alpha )>{\displaystyle \frac{1}{1-\alpha }}+{\displaystyle \frac{1}{\alpha }}={\displaystyle \frac{1}{\alpha (1-\alpha )}}.$$

Hence $\left|D\right|\Gamma (2-\alpha )h^{\alpha }>2.$ In this case the denominator of (30) is again greater than one, $|1+D\Gamma (2-\alpha )h^{\alpha }|>1$ and the error $e_1$ satisfies

$$|e_1|<2M_2{h}^2<C{h}^2.$$

Suppose that inequality (29) holds for all $n=1,\cdots ,m-1$. The values of $\zeta \left(\alpha \right)$ and $\Gamma (-\alpha )$ are negative. In both cases $D>0$ and $D<-\frac{2}{h^{\alpha }}\frac{\zeta \left(\alpha \right)-\zeta (1+\alpha )}{\Gamma (-\alpha )}$ the following estimate holds

$$|\zeta \left(\alpha \right)-\zeta (1+\alpha )+\Gamma (-\alpha )D{h}^{\alpha }|>|\zeta \left(\alpha \right)|+\zeta (1+\alpha ).$$

From (28) and Claim 4

$$\begin{array}{cc}\hfill |e_m|& <{\displaystyle \frac{1}{\left|\zeta \right(\alpha \left)\right|+\zeta (1+\alpha )}}\left(\left|\zeta \left(\alpha \right)\right|e_{m-1}+\sum _{k=1}^{m-1}{\displaystyle \frac{\left|e_{m-k}\right|}{k^{1+\alpha }}}+\left|{\overline{\sigma }}_{m-1}^{\left(\alpha \right)}{e}_1\right|h^2+\left|B_m\right|h^2\right)\hfill \\ \hfill & <{\displaystyle \frac{1}{\left|\zeta \right(\alpha \left)\right|+\zeta (1+\alpha )}}\left(\left|\zeta \left(\alpha \right)\right|e_{m-1}+\sum _{k=1}^{m-1}{\displaystyle \frac{\left|e_{m-k}\right|}{k^{1+\alpha }}}+{\displaystyle \frac{M_2}{4}}{h}^2+\left|B_m\right|h^2\right)\hfill \\ \hfill & <{\displaystyle \frac{1}{\left|\zeta \right(\alpha \left)\right|+\zeta (1+\alpha )}}\left(\left|\zeta \left(\alpha \right)\right|e_{m-1}+\sum _{k=1}^{m-1}{\displaystyle \frac{\left|e_{m-k}\right|}{k^{1+\alpha }}}+2C{h}^2\right).\hfill \end{array}$$

By induction hypoyhesis

$$\begin{array}{cc}\hfill |e_m|& <{\displaystyle \frac{C{h}^2}{\left|\zeta \right(\alpha \left)\right|+\zeta (1+\alpha )}}\left({\left|\zeta \left(\alpha \right)\right|(m-1)}^{\alpha }+\sum _{k=1}^{m-1}{\displaystyle \frac{{(m-k)}^{\alpha }}{k^{1+\alpha }}}+2\right).\hfill \end{array}$$

From inequality (7)

$$\begin{array}{cc}\hfill |e_m|& <{\displaystyle \frac{C{h}^2}{\left|\zeta \right(\alpha \left)\right|+\zeta (1+\alpha )}}\left(\left|\zeta \left(\alpha \right)\right|m^{\alpha }+(\zeta (1+\alpha )m^{\alpha }-2)+2\right)=C{h}^2{m}^{\alpha }.\hfill \end{array}$$

□

Theorem 14.

Let $\frac{1}{\Gamma (-\alpha )}<D<0$ and $y\in C^2[0,1]$. Then

$$\begin{array}{c}\hfill |e_n|<K{h}^2{n}^{\alpha },(1\le n\le N),\end{array}$$

(31)

where ${\displaystyle K=3M_2+\underset{1\le n\le N}{max}\left|B_n\right|}$.

Proof. 

The error $e_1$ on the first step satisfies the bound

$$\begin{array}{c}\hfill |e_1|<{\displaystyle \frac{2M_2{h}^2}{1-\alpha (1-\alpha )D\Gamma (-\alpha )h^{\alpha }}}<{\displaystyle \frac{2M_2{h}^2}{1-\alpha (1-\alpha )}}\le {\displaystyle \frac{8}{3}}{M}_2{h}^2.\end{array}$$

Assume that (31) holds for all $n=1,\cdots ,m-1$. Then

$$\begin{array}{cc}\hfill |e_m|& <{\displaystyle \frac{1}{\left|\zeta \left(\alpha \right)\right|+\zeta (1+\alpha )-\Gamma (-\alpha )D{h}^{\alpha }}}\left(\left|\zeta \left(\alpha \right)\right|e_{m-1}+\sum _{k=1}^{m-1}{\displaystyle \frac{\left|e_{m-k}\right|}{k^{1+\alpha }}}+\left|{\overline{\sigma }}_{m-1}^{\left(\alpha \right)}{e}_1\right|h^2+\left|B_m\right|h^2\right)\hfill \\ \hfill & <{\displaystyle \frac{1}{\left|\zeta \left(\alpha \right)\right|+\zeta (1+\alpha )-h^{\alpha }}}\left(\left|\zeta \left(\alpha \right)\right|K{h}^2{(m-1)}^{\alpha }+K{h}^2\sum _{k=1}^{m-1}{\displaystyle \frac{{(m-k)}^{\alpha }}{k^{1+\alpha }}}+K{h}^2\right)\hfill \\ \hfill & <{\displaystyle \frac{K{h}^2}{\left|\zeta \left(\alpha \right)\right|+\zeta (1+\alpha )-h^{\alpha }}}\left(\left|\zeta \left(\alpha \right)\right|m^{\alpha }+\zeta (1+\alpha )m^{\alpha }-1\right)\hfill \\ \hfill & <{\displaystyle \frac{K{h}^2{m}^{\alpha }}{\left|\zeta \left(\alpha \right)\right|+\zeta (1+\alpha )-h^{\alpha }}}\left(\left|\zeta \left(\alpha \right)\right|+\zeta (1+\alpha )-{\displaystyle \frac{1}{m^{\alpha }}}\right)\le K{h}^2{m}^{\alpha }.\hfill \end{array}$$

□

From Theorem 13 and Theorem 14 the errors $e_n$ of numerical solution (23) of the two-term equation in the interval $[0,1]$, satisfy the estimate

$$|e_n|<K{h}^2{n}^{\alpha }\le K{h}^{2-\alpha }.$$

(32)

for all $n=1,\cdots ,N$. Now we generalize the results of Theorem 13 and Theorem 14 to an arbitrary interval $[a,b]$. Consider the two-term equation

$$y^{\left(\alpha \right)}\left(t\right)+Dy\left(t\right)=F\left(t\right),y\left(a\right)=y_0,t\in [a,b].$$

(33)

Substitute $t=a+(b-a)s$ and $z\left(s\right)=y\left(t\right)=y(a+(b-a\left)s\right).$ The function $z\left(s\right)$ satisfies a two-term equation

$$z^{\left(\alpha \right)}\left(s\right)+{(b-a)}^{\alpha }Dz\left(s\right)=G\left(s\right)={(b-a)}^{\alpha }F(a+(b-a)s),z\left(0\right)=y_0,s\in [0,1]$$

(34)

Let $\mathcal{E}={\left\{e_n\right\}}_{n=1}^N$ be an N-dimensional vector with the errors of numerical solution (23) of equation (33). The $L_{\infty }$ (maximum) norm of a vector is the maximum of the absolute values of its elements. From Theorem 13 and Theorem 14 we get conditions for the convergence of numerical solution (23) of two-term equations (33) and (34).

Corollary 15.

The maximum error of numerical solution (23) of two-term equation (33), satisfies

$$\parallel \mathcal{E}\parallel <K{(b-a)}^{\alpha }{h}^{2-\alpha }$$

when the solution $y\in C^2[a,b]$ and $h=(b-a)/N$,

$$D\in \left(-\infty ,-\frac{2\left(\zeta \right(\alpha )-\zeta (1+\alpha \left)\right)}{\Gamma (-\alpha )h^{\alpha }}\right)\cup \left(\frac{1}{\Gamma (-\alpha ){(b-a)}^{\alpha }},\infty \right)$$

□

5. Time Fractional Black–Scholes Equation

The time fractional Black–Scholes equation is a fractional partial differential equation which is used for modeling the prices of the options [1,11,13,28,42,54,59].

$$\begin{array}{c}\hfill {\displaystyle \frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}}+{\displaystyle \frac{1}{2}}{\sigma }^2{S}^2{\displaystyle \frac{{\partial }^{2}C(S,t)}{\partial S^2}}+rS{\displaystyle \frac{\partial C(S,t)}{\partial S}}-rC(S,t)=0,\end{array}$$

(35)

where C is the option price, T is the expiry time, r is the risk-free rate, $\sigma$ is the volatility and S is the underlying stock price. Fractional Black-Scholes equation (35) has terminal and boundary conditions

$$C(S,T)=C_0\left(S\right),C(0,t)=C_1\left(t\right),C(\infty ,t)=C_2\left(t\right)$$

and the fractional derivative is

$$\begin{array}{c}\hfill {\displaystyle \frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}}={\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\displaystyle \frac{d}{dt}}{\int }_t^T{\displaystyle \frac{C(S,\xi )-C(S,T)}{{(\xi -z)}^{\alpha }}}d\xi .\end{array}$$

(36)

The terminal condition is transformed to an initial condition with the substitution [59]

$$t=T-z,\eta =T-\xi ,v(S,\mu )=C(S,\xi )=C(S,t-\mu ).$$

By applying the substitution to (36) we find [59]

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\partial }^{\alpha }C(S,t)}{\partial t^{\alpha }}}& =-{\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\displaystyle \frac{d}{dz}}{\int }_z^T{\displaystyle \frac{C(S,\xi )-C(S,T)}{{(\xi -z)}^{\alpha }}}d\xi \hfill \\ \hfill & =-{\displaystyle \frac{1}{\Gamma (1-\alpha )}}{\displaystyle \frac{d}{dz}}{\int }_0^z{\displaystyle \frac{v(S,\mu )-v(S,0)}{{(z-\mu )}^{\alpha }}}d\mu =D_z^{\alpha }v(S,z).\hfill \end{array}$$

The function $v(S,z)$ satisfies the partial fractional differential equation

$$\begin{array}{c}\hfill D_z^{\alpha }v(S,z)={\displaystyle \frac{1}{2}}{\sigma }^2{S}^2{\displaystyle \frac{{\partial }^{2}v(S,z)}{\partial S^2}}+rS{\displaystyle \frac{\partial v(S,z)}{\partial S}}-rC(S,z).\end{array}$$

(37)

The substitution $u(x,t)=v(e^x,z)$ transforms (37) into a linear partial fractional differential equation [59]

$$D_t^{\alpha }u(x,t)={\displaystyle \frac{1}{2}}{\sigma }^2{\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial x^2}}+\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right){\displaystyle \frac{\partial u(x,t)}{\partial x}}-ru(x,t).$$

(38)

Finite difference schemes for the time fractional Black-Scholes equation are constructed in [19,23,52,59]. Now we construct an implicit finite difference scheme which uses approximation (3) of the fractional derivative and central difference approximations of the partial derivatives for the following time fractional Black-Scholes equation

$$D_t^{\alpha }u(x,t)={\displaystyle \frac{1}{2}}{\sigma }^2{\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial x^2}}+\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right){\displaystyle \frac{\partial u(x,t)}{\partial x}}-ru(x,t)+F(x,t),$$

(39)

where $(x,t)\in \mathcal{R}=[B_l,B_r]\times [0,T]$. Equation (39) has initial and boundary conditions

$$u(x,0)=u_0\left(x\right),u(B_l,t)=u_1\left(t\right),u(B_r,t)=u_2\left(t\right).$$

Let M and N be positive integers and $\mathcal{G}$ be a rectangular grid on $\mathcal{R}$ which has a step size $h=(B_r-B_l)/N$ in space and $\tau =T/M$ in time,

$$\mathcal{G}=\left\{\right(nh,m\tau \left)\right|0\le n\le N,0\le m\le M\}.$$

Denote $u_n^m=u(nh,m\tau )$. The central difference approximations of the partial derivatives of equation (39) have second order accuracy

$${\displaystyle \frac{\partial u(nh,m\tau )}{\partial x}}={\displaystyle \frac{u_{n+1}^m-u_{n-1}^m}{2h}}+e_{n,m}^1{h}^2,$$

$${\displaystyle \frac{{\partial }^{2}u(nh,m\tau )}{\partial x^2}}={\displaystyle \frac{u_{n+1}^m-2u_n^m+u_{n-1}^m}{h^2}}+e_{n,m}^2{h}^2,$$

where $\left|e_{n,m}^1\right|<{\overline{M}}_3/3$ and $\left|e_{n,m}^2\right|<{\overline{M}}_4/12$. The numbers ${\overline{M}}_3$ and ${\overline{M}}_4$ are bounds for the third and fourth order partial derivatives

$${\overline{M}}_3=\underset{(x,t)\in \mathcal{R}}{max}\left|{\displaystyle \frac{{\partial }^{3}u(x,t)}{\partial x^3}}\right|,{\overline{M}}_4=\underset{(x,t)\in \mathcal{R}}{max}\left|{\displaystyle \frac{{\partial }^{4}u(x,t)}{\partial x^4}}\right|.$$

The numerical solution ${\left\{U_n^m\right\}}_{n=1}^{N-1}$ of equation (39) on layer m of the grid $\mathcal{G}$ satisfies

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\tau }^{-\alpha }}{\Gamma (-\alpha )}}\sum _{k=0}^m{\sigma }_k^{\left(\alpha \right)}{U}_n^{m-k}=& {\displaystyle \frac{{\sigma }^2}{2h^2}}\left(U_{n+1}^m-2U_n^m+U_{n-1}^m\right)\hfill \\ \hfill & +{\displaystyle \frac{1}{2h}}\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)\left(U_{n+1}^m-U_{n-1}^m\right)-r{U}_n^m+F_n^m\hfill \end{array}$$

for $m=2,\cdots ,M$ and has has initial conditions

$$U_n^0=u_0\left(nh\right),U_0^m=u_1\left(m\tau \right),U_N^m=u_2\left(m\tau \right).$$

Denote $\eta =\Gamma (-\alpha ){\tau }^{\alpha }/h^2$ and

$$K_1=-\left({\displaystyle \frac{{\sigma }^2}{2}}-\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)h\right)\eta ,K_2={\sigma }_0^{\left(\alpha \right)}+{\sigma }^2\eta +r{h}^2\eta ,K_3=-\left({\displaystyle \frac{{\sigma }^2}{2}}+\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)h\right)\eta .$$

The numerical solution of time-fractional Black-Scholes equation (39) is a solution of the system of linear equations

$$\begin{array}{c}\hfill K_1{U}_{n-1}^m-K_2{U}_n^m+K_3{U}_{n+1}^m=-\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{U}_n^{m-k}+h^2\eta F_n^m.\end{array}$$

(40)

Denote by $E_n^m=U_n^m-u_n^m$ the error of finite difference scheme (40) at the point $(nh,m\tau )$ and by $A_n^m{\tau }^{2-\alpha }$ and $B_n^m{h}^2$ be the the truncation errors of the approximations of the left-hand and right-hand sides of (39). The errors $E_n^m$ on row $m\ge 2$ of the grid $\mathcal{G}$ satisfy

$$K_1{E}_{n-1}^m+K_2{E}_n^m+K_3{E}_{n+1}^m=-\sum _{k=1}^m{\sigma }_k^{\left(\alpha \right)}{E}_n^{m-k}+\Gamma (-\alpha ){\tau }^{\alpha }\left(A_n^m{\tau }^{2-\alpha }+B_n^m{h}^2\right).$$

The errors are equal to zero, $E_n^0=E_0^m=E_N^m=0$ at the boundary points of $\mathcal{G}$. Let $\mathcal{K}$ be the $(N-1)$-dimensional tridiagonal matrix with entries $K_2$ on the main diagonal and $K_1$ and $K_3$ bellow and above the main diagonal. The error vector ${\mathcal{E}}^m$ of row m is a solution of the matrix equation

$$\mathcal{K}{\mathcal{E}}^m={\mathcal{Q}}^m;{\mathcal{E}}^m={\mathcal{K}}^{-1}{\mathcal{Q}}^m,$$

where

$$\begin{array}{c}\hfill Q_n^m=\zeta \left(\alpha \right)E_n^{m-1}-\sum _{k=1}^m{\displaystyle \frac{E_n^{m-k}}{k^{1+\alpha }}}-{\overline{\sigma }}_{m-1}^{\left(\alpha \right)}{E}_n^1+\Gamma (-\alpha ){\tau }^{\alpha }\left(A_n^m{\tau }^{2-\alpha }+B_n^m{h}^2\right).\end{array}$$

(41)

The numerical solution on the first row of $\mathcal{G}$ is computed with approximation (24)

$$\begin{array}{cc}\hfill {\displaystyle \frac{U_n^1-U_n^0}{\Gamma (2-\alpha ){\tau }^{\alpha }}}=& {\displaystyle \frac{{\sigma }^2}{2h^2}}\left(U_{n+1}^1-2U_n^1+U_{n-1}^1\right)+{\displaystyle \frac{1}{2h}}\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)\left(U_{n+1}^1-U_{n-1}^1\right)-r{U}_n^1+F_n^1.\hfill \end{array}$$

Denote $\tilde{\eta }=\Gamma (2-\alpha ){\tau }^{\alpha }/h^2$. Then

$${\tilde{K}}_1{U}_{n-1}^1+{\tilde{K}}_2{U}_n^1+{\tilde{K}}_3{U}_{n+1}^1=U_n^0+h^2\eta F_n^1,$$

where $1\le n\le N$ and

$${\tilde{K}}_1=-\left({\displaystyle \frac{{\sigma }^2}{2}}-\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)h\right)\tilde{\eta },{\tilde{K}}_2=1+{\sigma }^2\tilde{\eta }+r{h}^2\tilde{\eta },{\tilde{K}}_3=-\left({\displaystyle \frac{{\sigma }^2}{2}}+\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)h\right)\tilde{\eta }.$$

The errors of the numerical solution of the fractional Black-Scholes equation (39) on the first row of $\mathcal{G}$ are the solutions of the system of linear equations

$${\tilde{K}}_1{E}_{n-1}^1+{\tilde{K}}_2{E}_n^1+{\tilde{K}}_3{E}_{n+1}^1=\Gamma (2-\alpha ){\tau }^{\alpha }\left(A_n^1{\tau }^{2-\alpha }+B_n^1{h}^2\right).$$

Let $\tilde{\mathcal{K}}$ be the $(N-1)$-dimensional tridiagonal matrix with elements ${\tilde{K}}_1,{\tilde{K}}_2$ and ${\tilde{K}}_3$. The error vector for the first row satisfies

$$\tilde{\mathcal{K}}{\mathcal{E}}^1={\mathcal{Q}}^1;{\mathcal{E}}^1={\tilde{\mathcal{K}}}^{-1}{\mathcal{Q}}^1,$$

where

$$Q_n^1=\Gamma (2-\alpha ){\tau }^{\alpha }\left(A_n^1{\tau }^{2-\alpha }+B_n^1{h}^2\right).$$

Let $C={max}_{n,m}\left\{\right|A_n^m|,|B_n^m\left|\right\}$ for all $n=1,\cdots ,N;m=1,\cdots ,M$ and $K=|\Gamma (-\alpha \left)C\right|$. The $L_{\infty }$ norm of a matrix is the maximum of the absolute row sums.

Theorem 16.

Let $u\in C^4[B_l,B_r]\times C^2[0,T]$ and $r<\frac{{\sigma }^2}{2h}$. Then

$$\begin{array}{c}\hfill \parallel {\mathcal{E}}^m\parallel <K{m}^{\alpha }{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)\end{array}$$

(42)

for all $m=1,\cdots M$.

Proof. 

When $0<\alpha <1$, the gamma function satisfies [15]

$$0.88<\Gamma (2-\alpha )<1,$$

$$|\Gamma (-\alpha )|={\displaystyle \frac{\Gamma (2-\alpha )}{\alpha (1-\alpha )}}>3.$$

Then

$$\parallel Q^1\parallel <C\Gamma (2-\alpha ){\tau }^{\alpha }\left({\tau }^{2-\alpha }+h^2\right)<C{\tau }^{\alpha }\left({\tau }^{2-\alpha }+h^2\right).$$

When $r<{\sigma }^2/\left(2h\right)$ we have that

$$\left|K_1\right|+|K_3|=\left|{\displaystyle \frac{{\sigma }^2}{2}}-\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)h\right|\left|\eta \right|+\left|{\displaystyle \frac{{\sigma }^2}{2}}+\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right)h\right|\left|\eta \right|\le {\sigma }^2\left|\eta \right|,$$

$$\left|K_2\right|-\left|K_1\right|-|K_3|>|{\sigma }_0^{\left(\alpha \right)}|.$$

Similarly

$$\left|{\tilde{K}}_2\right|-\left|{\tilde{K}}_1\right|-\left|{\tilde{K}}_3\right|>1.$$

The matrices $\mathcal{K}$ and $\tilde{\mathcal{K}}$ are M-matrices. From the Ahlberg-Nilson-Varah bound [27,40]

$$\parallel {\tilde{\mathcal{K}}}^{-1}\parallel <1,\parallel {\mathcal{K}}^{-1}\parallel <{\displaystyle \frac{1}{|{\sigma }_0^{\left(\alpha \right)}|}}={\displaystyle \frac{1}{\left|\zeta \right(\alpha \left)\right|+\zeta (1+\alpha )}}.$$

The $L_{\infty }$ norm of the error vector of the first row of $\mathcal{G}$ satisfies the bound

$$\parallel {\mathcal{E}}^1\parallel \le \parallel {\tilde{\mathcal{K}}}^{-1}\parallel \parallel {\mathcal{Q}}^1\parallel <C{\tau }^{\alpha }\left({\tau }^{2-\alpha }+h^2\right)<K{\tau }^{\alpha }\left({\tau }^{2-\alpha }+h^2\right).$$

Suppose that inequality (42) holds for all $m=1,\cdots ,p-1$. From (41)

$$\begin{array}{cc}\hfill \left|Q_n^p\right|& <\left|\zeta \left(\alpha \right)\right|\parallel {\mathcal{E}}^{p-1}\parallel +\sum _{k=1}^p{\displaystyle \frac{\parallel {\mathcal{E}}^{p-k}\parallel }{k^{1+\alpha }}}+\left|{\overline{\sigma }}_{p-1}^{\left(\alpha \right)}{E}_n^1\right|+\Gamma (-\alpha )C{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)\hfill \\ \hfill & <K{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)\left(\left|\zeta \left(\alpha \right)\right|p^{\alpha }+\sum _{k=1}^p{\displaystyle \frac{{(p-k)}^{\alpha }}{k^{1+\alpha }}}+{\displaystyle \frac{1}{8|\Gamma (-\alpha \left)\right|}}+1\right)\hfill \\ \hfill & <K{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)\left(\left|\zeta \left(\alpha \right)\right|p^{\alpha }+\zeta (1+\alpha )p^{\alpha }-1+{\displaystyle \frac{1}{8|\Gamma (-\alpha \left)\right|}}\right).\hfill \end{array}$$

Therefore

$$\parallel {\mathcal{Q}}^p\parallel <\left(\right|\zeta \left(\alpha \right)|+\zeta (1+\alpha ))K{p}^{\alpha }{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right),$$

$$\parallel {\mathcal{E}}^p\parallel <\parallel {\mathcal{K}}^{-1}\parallel \parallel {\mathcal{Q}}^p\parallel <K{p}^{\alpha }{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right).$$

□

Corollary 17.

Let $r<\frac{{\sigma }^2}{2h}$. Then

$$\parallel {\mathcal{E}}^m\parallel <K{T}^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)$$

for all $m=1,2,\cdots ,M$.

Proof. 

From Theorem 16

$$\parallel {\mathcal{E}}^m\parallel <K{m}^{\alpha }{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)\le K{M}^{\alpha }{\tau }^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right)=K{T}^{\alpha }\left(h^2+{\tau }^{2-\alpha }\right).$$

□

Example 3: Consider the following time fractional Black-Scholes equation

$$\left\{\begin{array}{c}{D}_t^{\alpha }u(x,z)={\displaystyle \frac{1}{2}}{\sigma }^2{\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial S^2}}+\left(r-{\displaystyle \frac{{\sigma }^2}{2}}\right){\displaystyle \frac{\partial u(x,t)}{\partial S}}-ru(x,t)+t^{1-\alpha }{e}^x{E}_{1,2-\alpha }\left(t\right),\hfill \\ u(x,0)=e^x,u(0,t)=e^t,u(1,t)=e^{1+t}.\hfill \end{array}\right.$$

(43)

Equation (43) has a solution $u(x,t)=e^{t+x}$. Experimental results of the numerical solution of equation (43) with parameters $\sigma =0.35,r=0.05$ are given in Table 5 and Table 6. The orders of convergence of the finite difference scheme are ${log}_{\frac{{\tau }_1}{{\tau }_2}}{\displaystyle \frac{{\mathrm{error}}_1}{{\mathrm{error}}_2}}$ in time and ${log}_{\frac{h_1}{h_2}}{\displaystyle \frac{{\mathrm{error}}_1}{{\mathrm{error}}_2}}$ in space.

Example 4: Consider the time fractional Black-Scholes equation (37) for pricing European call options, with a source term $F(x,t)=0$ and initial and boudary conditions

$${\displaystyle u_0\left(x\right)={(e^x-K)}^{+},u_1\left(t\right)=0,u_2\left(t\right)=e^{B_r}-K{e}^{-rt}.}$$

The European call option premium curves for different values of $\alpha$ are given in Figure 1, for values of the parameters $r=2\%$, $\sigma =30\%$, $B_l=0$, $B_r=200$, $T=1$ (year) and strike price $K=100$. Regarding near-the-money options, lower values of $\alpha$ prices them lower, and evaluates higher out-of-the-money and in-the-money options, compared to the classical Black-Scholes dynamics.

The graphs of call option and put option premiums for $\alpha =0.7$ and all t are given in Figure 2 and Figure 3.

6. Conclusions

In this paper we study the convergence and order of numerical solutions of ordinary and partial fractional differential equations which use approximation (3) of the fractional derivative. The results of the numerical experiments given in the paper illustrate the theoretical results. In future work we will use the method from [7] for construction of secon-order and high-order approximations of the fractional derivative and we will study the convergence of the finite difference schemes for numerical solution of fractional differential equations.