---
abstract: |
    In this paper we use the Fourier transform method and the expansion formula for the polylogarithm function to derive approximations for the Caputo derivative of order $2-{\alpha}$ and $2$. The approximations are applied for computing the numerical solutions of the fractional relaxation and subdiffusion equations. While the properties of the weights of the approximation of order $2-{\alpha}$ are similar to those of the $L1$ approximation, the corresponding numerical solution is more accurate for the examples discussed in the paper and many of the functions used in practice.

    [**2010 Math Subject Classification:**]{} 26A33, 34A08, 42A38, 65D30.\
    [**Key Words and Phrases:**]{} Caputo derivative, Fourier transform, asymptotic expansion, fractional differential equation.
author:
- |
    Yuri Dimitrov\
    Department of Applied Mathematics and Statistics\
    University of Rousse, Rousse 7017, Bulgaria\
    `ymdimitrov@uni-ruse.bg`
nocite: '[@*]'
title: '[**APPROXIMATIONS FOR THE CAPUTO DERIVATIVE (II)** ]{}'
---

Introduction
============

Approximations of fractional integrals and derivatives have recently been an active research topic [@Dimitrov2015; @DingLi2016; @GaoSunZhang2014; @Odibat2006; @PoosehAlmeidaTorres2012; @TianZhouDeng2012; @YanPalFord2014]. The fractional integral of order ${\alpha}>0$ and the Caputo derivative of order ${\alpha}$, where $0<{\alpha}<1$ are defined as $$I^{\alpha}y(x)={\dfrac{1}{{\Gamma}({\alpha})}}\int_0^x (x-t)^{{\alpha}-1} y(t)dt,$$ $$y^{({\alpha})}(x)=D^{{\alpha}} y(x)={\dfrac{1}{\Gamma (1-{\alpha})}}\int_0^x \dfrac{y'(t)}{(x-t)^{{\alpha}}}d t.$$ Let $h=x/n$, where $n$ is a positive integer, and $x_m=m h$, $y_m=y(x_m)$. The $L1$ approximation for the Caputo derivative is a commonly used approximation for numerical solution of fractional differential equations. $$\label{2_1}
y^{(\alpha)}_n  =\dfrac{1}{{\Gamma}(2-{\alpha})h^\alpha}\sum_{k=0}^{n} {\sigma}_k^{(\alpha)} y_{n-k}+O{\left (}h^{2-{\alpha}}{\right )},$$ where ${\sigma}_0^{(\alpha)}=1$, ${\sigma}_n^{(\alpha)}=(n-1)^{1-a}-n^{1-a}$ and $${\sigma}_k^{(\alpha)}=(k+1)^{1-\alpha}-2k^{1-\alpha}+(k-1)^{1-\alpha}, \quad (k=2,\cdots,n-1).$$ The weights ${\sigma}_k^{(\alpha)}$ of the $L1$ approximation have the following properties $$\begin{aligned}
\label{2_2}
&{\sigma}_0^{(\alpha)}>0,\; {\sigma}_1^{(\alpha)}<{\sigma}_2^{(\alpha)}<\cdots<{\sigma}_k^{(\alpha)}<\cdots<{\sigma}_{n-1}^{(\alpha)}<0,\;{\sigma}_{n}^{(\alpha)}<0,\nonumber\\
&\sum_{k= 0}^n {\sigma}_k^{(\alpha)} = 0,\quad \sum_{k= 1}^n k {\sigma}_k^{(\alpha)} = -n^{1-{\alpha}},\\
&{\sigma}_k^{(\alpha)}= {\dfrac{C_1}{k^{1+{\alpha}}}}+O{\left (}{\dfrac{1}{k^{2+{\alpha}}}} {\right )},\quad {\sigma}_n^{(\alpha)}= {\dfrac{C_2}{n^{\alpha}}}+O{\left (}{\dfrac{1}{n^{1+{\alpha}}}} {\right )},\nonumber\end{aligned}$$ where $C_1={\alpha}({\alpha}-1)$ and $C_2={\alpha}-1$. When the function $y$ has a continuous second derivative, the $L1$ approximation has accuracy $O(h^{2-\alpha})$ ([@LinXu2007]). The numerical solution of the fractional relaxation equation which uses the $L1$ approximation for the Caputo derivative is computed with [@Dimitrov2015] $$\begin{aligned}
\label{2_3}
 u_n={\dfrac{1}{{\sigma}_0^{({\alpha})}+{\Gamma}(2-{\alpha})h^{\alpha}}} {\left (}{\Gamma}(2-{\alpha})h^{\alpha}F_n  -\sum_{k=1}^{n}{\sigma}_k^{({\alpha})} u_{n-k}{\right )},\; u_0=y_0.\end{aligned}$$ In Table 1 we compute the error and the order of numerical solution for Equation I and ${\alpha}=0.25$, Equation II and ${\alpha}=0.5$ and Equation III, ${\alpha}=0.75$. In [@Dimitrov2015] we derived the second-order expansion of the $L1$ approximation $$\begin{aligned}
{\dfrac{1}{{\Gamma}(2-{\alpha})h^{\alpha}}}\sum_{k=0}^n {\sigma}_k^{({\alpha})} y(x-kh)=y^{({\alpha})}(x)+&{\dfrac{{\zeta}({\alpha}-1)}{{\Gamma}(2-{\alpha})}}y''(x)h^{2-{\alpha}}+O{\left (}h^2{\right )}.\end{aligned}$$ By approximating $y''(x)$ using second-order backward difference we obtain the second-order approximation for the Caputo derivative $$\label{2_4}
y^{({\alpha})}_n={\dfrac{1}{{\Gamma}(2-{\alpha})h^{\alpha}}}\sum_{k=0}^n {\delta}_k^{({\alpha})} y_{n-k}+O{\left (}h^{2}{\right )},$$ where ${\delta}_k^{({\alpha})}={\sigma}_k^{({\alpha})}$ for $2\leq k\leq n$ and $${\delta}_0^{({\alpha})}={\sigma}_0^{({\alpha})}-{\zeta}({\alpha}-1),\; {\delta}_1^{({\alpha})}={\sigma}_1^{({\alpha})}+2{\zeta}({\alpha}-1),\; {\delta}_2^{({\alpha})}={\sigma}_2^{({\alpha})}-{\zeta}({\alpha}-1).$$ The weights ${\delta}_k^{(\alpha)}$ of approximation satisfy $${\delta}_0^{(\alpha)}>0,\;{\delta}_1^{(\alpha)}<0,\;{\delta}_2^{(\alpha)}>0,\; {\delta}_3^{(\alpha)}<{\delta}_4^{(\alpha)}<\cdots<{\delta}_k^{(\alpha)}<\cdots<{\delta}_{n-1}^{(\alpha)}<0.$$ The asymptotic expansions of the trapezoidal approximation for the definite integral and the integral approximations for the fractional derivative involve the values of the [*Riemann zeta function*]{} defined as $${\zeta}({\alpha})=
    \displaystyle{\sum_{n=1}^\infty {\dfrac{1}{n^{\alpha}}}}, \;({\alpha}>1),
    \qquad  {\zeta}({\alpha})=\displaystyle{\dfrac{1}{1-2^{1-{\alpha}}}\sum_{n=1}^\infty {\dfrac{(-1)^{n-1}}{n^{\alpha}}}}, \;  ({\alpha}>0).$$ The Riemann zeta function is a special case $(x=1)$ of the [*polylogarithm function*]{} is defined as $$Li_{\alpha}(x)=\sum_{n=1}^{\infty}{\dfrac{x^n}{n^{\alpha}}}=x+{\dfrac{x^2}{2^{\alpha}}}+\cdots+{\dfrac{x^n}{n^{\alpha}}}+\cdots$$ The polylogarithm function has properties $$\begin{aligned}
&Li_{\alpha}(x)+Li_{\alpha}(-x)=2^{1-{\alpha}}Li_{\alpha}(x^2),\label{3_1}\\
&Li_{\alpha}(x)={\Gamma}(1-{\alpha}){\left (}\ln {\dfrac{1}{x}}{\right )}^{{\alpha}-1}+\sum_{n=0}^{\infty}{\dfrac{{\zeta}({\alpha}-n)}{n!}}{\left (}\ln x{\right )}^n,\label{3_2}\end{aligned}$$ where ${\alpha}\neq 1,2,3,\cdots$ and $|\ln x|<2\pi$. From with $x=e^{i {\omega}h}$ we obtain $$\begin{aligned}
\label{4_1}
Li_{\alpha}{\left (}e^{i w h}{\right )}={\Gamma}(1-&{\alpha})(-i w)^{{\alpha}-1}h^{{\alpha}-1}+{\zeta}({\alpha})-(-i w){\zeta}({\alpha}-1)h\\
&+(-i w)^2{\dfrac{{\zeta}({\alpha}-2)}{2}}h^2-(-i w)^3{\dfrac{{\zeta}({\alpha}-3)}{6}}h^3+O{\left (}h^4{\right )}.\nonumber\end{aligned}$$ In [@Dimitrov2016_1] we use the Fourier transform method to derive the asymptotic expansion formula of the trapezoidal approximation for the fractional integral $$\begin{aligned}
  h^{{\alpha}} \sum_{k=1}^{N-1} {\dfrac{y(x-k h)}{k^{1-{\alpha}}}}&=\int_0^x {\dfrac{y(t)}{(x-t)^{1-{\alpha}}}}dt +\sum_{k=0}^\infty (-1)^k{\dfrac{{\zeta}(1-{\alpha}-k)}{k!}}y^{(k)}(x) h^{k+{\alpha}}\\
-{\Gamma}({\alpha})&\sum_{k=0}^\infty {\dfrac{B_{k+1}}{(k+1)!}}{\left (} \sum_{m=0}^k (-1)^m\binom{k}{m}{\dfrac{x^{{\alpha}-m-1}}{{\Gamma}({\alpha}-m)}}y^{(m-k)}(0) {\right )}h^{k+1}.\end{aligned}$$ In the first part of this paper [@Dimitrov2016_2] we obtain approximations and for the Caputo derivative of order $2-{\alpha}$. $$\begin{aligned}
\label{4_2}
{\dfrac{1}{2{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}{\sigma}_k^{({\alpha})} y_{n-k}=y_n^{({\alpha})}+O{\left (} h^{2-{\alpha}} {\right )},\end{aligned}$$ where ${\sigma}_0^{({\alpha})}=1-2{\zeta}({\alpha}),\; {\sigma}_1^{({\alpha})}={\dfrac{1}{2^{\alpha}}}+2{\zeta}({\alpha})$ and $${\sigma}_k^{({\alpha})}={\dfrac{1}{(k+1)^{\alpha}}}-{\dfrac{1}{(k-1)^{\alpha}}},\quad (k=2,\cdots,n-2),$$ $${\sigma}_{n-1}^{({\alpha})}=-{\dfrac{1}{(n-2)^{\alpha}}}-2 \left(\sum _{k=1}^{n-1} \frac{1}{k^{\alpha}}-\frac{n^{1-{\alpha}}}{1-{\alpha}}-\zeta ({\alpha})\right),$$ $${\sigma}_{n}^{({\alpha})}=-{\dfrac{1}{(n-1)^{\alpha}}}+2 \left(\sum _{k=1}^{n-1} \frac{1}{k^{\alpha}}-\frac{n^{1-{\alpha}}}{1-{\alpha}}-\zeta ({\alpha})\right).$$ $$\label{4_3}
{\dfrac{1}{{\Gamma}(-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}{\sigma}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (}h^{2-{\alpha}}{\right )},$$ where ${\sigma}_0^{({\alpha})}={\zeta}({\alpha})-{\zeta}(1+{\alpha}),\; {\sigma}_1^{({\alpha})}=1-{\zeta}({\alpha})$ and $${\sigma}_k^{({\alpha})}={\dfrac{1}{k^{1+{\alpha}}}},\quad (k=2,\cdots,n-2),$$ $${\sigma}_{n-1}^{({\alpha})}={\dfrac{1}{(n-1)^{1+{\alpha}}}}-{\dfrac{n^{1-{\alpha}}}{{\alpha}(1-{\alpha})}}+n{\left (}{\zeta}(1+{\alpha})-\sum_{k=1}^{n-1}{\dfrac{1}{k^{1+{\alpha}}}} {\right )}-
{\left (}{\zeta}({\alpha})-\sum_{k=1}^{n-1}{\dfrac{1}{k^{{\alpha}}}} {\right )},$$ $${\sigma}_{n}^{({\alpha})}=(1-n){\left (}{\zeta}(1+{\alpha})-\sum_{k=1}^{n-1}{\dfrac{1}{k^{1+{\alpha}}}} {\right )}+{\left (}{\zeta}({\alpha})-\sum_{k=1}^{n-1}{\dfrac{1}{k^{{\alpha}}}} {\right )}+{\dfrac{n^{1-{\alpha}}}{{\alpha}(1-{\alpha})}}.$$ In section 2 we use the Fourier transform method to derive the asymptotic expansion formula of approximation . In section 3 we use we obtain approximation for the Caputo derivative of order $2-{\alpha}$ and the second-order approximation . $$\begin{aligned}
\label{5_1}
{\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n} {\sigma}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (} h^{2-{\alpha}} {\right )},\end{aligned}$$ where $$\begin{aligned}
{\sigma}_0^{({\alpha})}=2^{\alpha}-{\left (}2^{\alpha}-1{\right )}{\zeta}({\alpha}),\quad 
{\sigma}_1^{({\alpha})}={\left (}{\dfrac{2}{3}}{\right )}^{\alpha}-2^{\alpha}+{\left (}2^{\alpha}-1{\right )}{\zeta}({\alpha}),\end{aligned}$$ $${\sigma}_{k}^{({\alpha})}=2^{\alpha}{\left (}{\dfrac{1}{(2k+1)^{\alpha}}}-{\dfrac{1}{(2k-1)^{\alpha}}}{\right )},\quad (k=2,\cdots,n-2),$$ $${\sigma}_{n-1}^{({\alpha})}={\dfrac{2^{\alpha}}{(2n-1)^{\alpha}}}-{\dfrac{2^{\alpha}}{(2n-3)^{\alpha}}}+\frac{n^{1-{\alpha}}}{1-{\alpha}}-{\zeta}({\alpha})-2^{\alpha}{\left (}\sum _{k=1}^{n} {\dfrac{1}{(2k-1)^{\alpha}}}-\zeta ({\alpha}){\right )},$$ $${\sigma}_{n}^{({\alpha})}=-{\dfrac{2^{\alpha}}{(2n-1)^{\alpha}}}+{\zeta}({\alpha})-\frac{n^{1-{\alpha}}}{1-{\alpha}}+2^{\alpha}{\left (}\sum _{k=1}^{n} {\dfrac{1}{(2k-1)^{\alpha}}}-\zeta ({\alpha}){\right )}.$$ $$\begin{aligned}
\label{5_2}
{\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n} {\delta}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (} h^{2} {\right )},\end{aligned}$$ where ${\delta}_{k}^{({\alpha})}={\sigma}_{k}^{({\alpha})},\; (k=3,\cdots,n)$ and $${\delta}_0^{({\alpha})}=2^{\alpha}+{\dfrac{3}{2}}{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})+{\left (} 2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1),$$ $${\delta}_1^{({\alpha})}=2^{\alpha}{\left (}{\dfrac{1}{3^{\alpha}}}-1 {\right )}-2{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})-2{\left (} 2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1),$$ $${\delta}_2^{({\alpha})}=2^{\alpha}{\left (}{\dfrac{1}{5^{\alpha}}}-{\dfrac{1}{3^{\alpha}}}{\right )}+{\dfrac{1}{2}}{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})+{\left (} 2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1).$$ The weights of the $L1$ approximation satisfy . Approximations , and have similar properties of the weights. These properties are used in the proofs for the convergence of the numerical solutions. An important factor for the accuracy of the numerical solutions of order $2-{\alpha}$ which use approximations ,, and for the Caputo derivative is the coefficient of the term $y''(x)h^{2-{\alpha}}$ in the asymptotic expansion formulas of the approximations. Approximations ,, and have coefficients $$C_1({\alpha})=\frac{\zeta ({\alpha}-1)}{\Gamma (2-{\alpha})},C_8({\alpha})=\frac{\zeta ({\alpha})-2\zeta ({\alpha}-1)}{2\Gamma (1-{\alpha})},C_9({\alpha})=\frac{\zeta ({\alpha})-\zeta ({\alpha}-1)}{2\Gamma (-{\alpha})},$$ $$C_{10}({\alpha})=\frac{{\left (}2-2^{{\alpha}}{\right )}\zeta ({\alpha}-1)-{\left (}2^{{\alpha}}-1{\right )}\zeta ({\alpha})}{2\Gamma (1-{\alpha})}.$$

![Graph of the absolute value of the coefficients $C_1({\alpha})$(black), $C_8({\alpha})$(dashed), $C_9({\alpha})$(blue) and $C_{10}({\alpha})$(red) for $0<{\alpha}<1$.](graph1.eps){width="50.00000%"}

In Figure 1 we compare the absolute values of the coefficients $C_1({\alpha})$, $C_8({\alpha})$, $C_9({\alpha})$ and $C_{10}({\alpha})$. In the first part of this paper [@Dimitrov2016_2] we derive approximations and for the Caputo derivative and we apply the approximations for computing the numerical solution of the fractional relaxation equation. The accuracy of the numerical solution of Equation I, Equation II and Equation III which uses approximation is higher than the accuracy of the numerical solutions using and . In Table 1 and Table 3 we compute the errors of the numerical solutions of Equation I, Equation II and Equation III which use approximations and . The numerical solution which uses approximation for the Caputo derivative has a higher accuracy than the numerical solution using the $L1$ approximation . The improvement is $30\%$ for Equation I and ${\alpha}=0.25$, $15\%$ for Equation II and ${\alpha}=0.5$ and $4.7\%$ for Equation III with ${\alpha}=0.75$. In section 4 we construct a finite difference scheme for the fractional subdiffusion equation with accuracy $O{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}$, using approximation for the Caputo derivative and we analyse the convergence of the scheme.

Approximation for the Caputo Derivative of Order $\boldsymbol{4-{\alpha}}$
==========================================================================

The construction of approximation is based on the trapezoidal approximation for the fractional integral and the second-order backward difference approximation for the first derivative on each subinterval. The coefficient $C_8({\alpha})$ of the term $y''_n h^{2-{\alpha}}$ in the expansion of approximation is greater than the coefficient $C_{1}({\alpha})$ of the $L1$ approximation. The midpoint approximation for the definite integral is around twice more accurate than the trapezoidal approximation. In this section we use the Fourier transform method to derive approximations and for the Caputo derivative of order $2-{\alpha}$ and $2$ based on the midpoint approximation for the fractional integral. The coefficient $C_{10}({\alpha})$ of approximation is smaller than $C_1({\alpha})$. The accuracy of the numerical solutions of Equation I, Equation II and Equation III which use approximation for the Caputo derivative is higher than the accuracy of the numerical solutions which use the $L1$ approximation . The exponential Fourier Transform of the function $y$ is defined as $$\mathcal{F}[y(x)](w)=\hat{y}(w)=\int_{-\infty}^{\infty}e^{i w t}y(t)d t.$$ The Fourier transform has properties $\mathcal{F}[y(x-b)](w)=e^{i w b} \hat{y}(w)$ and $$\mathcal{F}[D^{\alpha}y(x)](w)=(-iw)^{\alpha}\hat{y}(w),\quad \mathcal{F}[I^{\alpha}y(x)](w)=(-iw)^{-{\alpha}} \hat{y}(w).$$ From the midpoint approximation for the fractional integral in the definition of the Caputo derivative $$y^{({\alpha})}(x)\approx {\dfrac{h}{{\Gamma}(1-{\alpha})}}\sum_{k=1}^{n}
{\dfrac{y'{\left (}x_{k-1/2} {\right )}}{(x-x_{k-1/2})^{\alpha}}}={\dfrac{h^{1-{\alpha}}}{{\Gamma}(1-{\alpha})}}\sum_{k=1}^{n}{\dfrac{y'_{k-1/2}}{(n-k+1/2)^{\alpha}}}.$$ By approximating $y'_{k-1/2}\approx (y_k-y_{k-1})/{h}$ we obtain $$y^{({\alpha})}(x)\approx {\dfrac{2^{\alpha}h^{1-{\alpha}}}{{\Gamma}(1-{\alpha})}}\sum_{k=1}^{n}{\dfrac{y'_{k-1/2}}{(2n-2k+1)^{\alpha}}}\approx 
{\dfrac{2^{\alpha}}{{\Gamma}(1-{\alpha}) h^{\alpha}}}\sum_{k=1}^{n}{\dfrac{y_k-y_{k-1}}{(2n-2k+1)^{\alpha}}}.$$ Substitute $K=n-k$ $${\Gamma}(1-{\alpha}){\left (}{\dfrac{h}{2}}{\right )}^{\alpha}y^{({\alpha})}(x)\approx  \sum_{K=0}^{n-1}{\dfrac{y_{n-K}-y_{n-K+1}}{(2K+1)^{\alpha}}}
=   \sum_{k=0}^{n-1}{\dfrac{y_{n-K}}{(2K+1)^{\alpha}}}-\sum_{k=0}^{n-1}{\dfrac{y_{n-K+1}}{(2K+1)^{\alpha}}},$$ $$y^{({\alpha})}(x)\approx {\dfrac{2^{\alpha}}{{\Gamma}(1-{\alpha}) h^{\alpha}}} {\left (}\sum_{k=0}^{n-1}{\dfrac{y_{n-k}}{(2k+1)^{\alpha}}}-\sum_{k=1}^{n}{\dfrac{y_{n-k}}{(2k-1)^{\alpha}}}{\right )}.$$ Denote $$\begin{aligned}
\label{8_1}
S_n[y]=\sum_{k=0}^{n-1}{\dfrac{y_{n-k}}{(2k+1)^{\alpha}}}-\sum_{k=1}^{n}{\dfrac{y_{n-k}}{(2k-1)^{\alpha}}}.\end{aligned}$$ $$\begin{aligned}
\label{8_2}
S_n[y]=y_n+\sum_{k=1}^{n-1}{\left (}{\dfrac{1}{(2k+1)^{\alpha}}}-{\dfrac{1}{(2k-1)^{\alpha}}}{\right )}y_{n-k}-{\dfrac{y_0}{(2n-1)^{\alpha}}}.\end{aligned}$$ The generating function of an approximation is related to the Fourier transform of the approximation. The Fourier transform/generating function method is used by Ding and Li [@DingLi2016], Lubich [@Lubich1986], Tian et. al. [@TianZhouDeng2012] for constructing approximations for the fractional derivative. By applying Fourier transform to and letting $n\rightarrow \infty$ $$\mathcal{F}[S_\infty[y]](w)={\left (}\sum_{k=0}^{\infty}{\dfrac{e^{iwkh}}{(2k+1)^{\alpha}}}-\sum_{k=1}^{\infty}{\dfrac{e^{iwkh}}{(2k-1)^{\alpha}}}  {\right )}\hat{y}(w).$$ Denote $W=e^{i w h/2}$. $$\mathcal{F}[S_\infty[y]](w)={\left (}\sum_{k=0}^{\infty}{\dfrac{W^{2k}}{(2k+1)^{\alpha}}}-\sum_{k=1}^{\infty}{\dfrac{W^{2k}}{(2k-1)^{\alpha}}}  {\right )}\hat{y}(w),$$ $$\mathcal{F}[S_\infty[y]](w)=\hat{y}(w){\left (}{\dfrac{1}{W}}-W{\right )}\sum_{k=1}^{\infty}{\dfrac{W^{2k-1}}{(2k-1)^{\alpha}}}.$$ From $$\begin{aligned}
2\sum_{k=1}^{\infty}{\dfrac{W^{2k-1}}{(2k-1)^{\alpha}}}&=Li_{\alpha}(W)-Li_{\alpha}(-W) =2Li_{\alpha}(W)-Li_{\alpha}(W)-Li_{\alpha}(-W)\\
&=2Li_{\alpha}(W)-2^{1-{\alpha}}Li_{\alpha}{\left (}W^2{\right )}.\end{aligned}$$ Then $$\mathcal{F}[S_\infty[y]](w)=\hat{y}(w){\left (}e^{-\frac{i w h}{2}}-e^{\frac{i w h}{2}}{\right )}{\left (}Li_{\alpha}(e^{\frac{i w h}{2}})-{\dfrac{1}{2^{\alpha}}}Li_{\alpha}{\left (}e^{i w h}{\right )}{\right )}.$$ From the Taylor expansion formula of the exponential function $$\begin{aligned}
\label{9_1}
e^{-\frac{i w h}{2}}-e^{\frac{i w h}{2}}=(-i w h)+{\dfrac{(-i w h)^3}{24}} +O{\left (}h^5{\right )}.\end{aligned}$$ From the function $Li_{\alpha}{\left (}e^{\frac{i w h}{2}}{\right )}$ has fourth-order expansion $$\begin{aligned}
\label{9_2}
Li_{\alpha}{\left (}e^{\frac{i w h}{2}}{\right )}={\dfrac{{\Gamma}(1-{\alpha})}{2^{{\alpha}-1}}}&(-i w h)^{{\alpha}-1}+{\zeta}({\alpha})-(-i w){\dfrac{{\zeta}({\alpha}-1)}{2}}h\\
+(-i w)^2&{\dfrac{{\zeta}({\alpha}-2)}{8}}h^2-(-i w)^3{\dfrac{{\zeta}({\alpha}-3)}{48}}h^3+O{\left (}h^4{\right )}.\nonumber\end{aligned}$$ By combining and we obtain $$\begin{aligned}
\label{9_3}
Li_{\alpha}{\left (}e^{\frac{i w h}{2}}{\right )}&-{\dfrac{1}{2^{\alpha}}}Li_{\alpha}{\left (}e^{i w h}{\right )}={\dfrac{{\Gamma}(1-{\alpha})}{2^{{\alpha}}}}(-i w)^{{\alpha}-1}h^{{\alpha}-1}
+{\left (}1-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha})\nonumber\\
-(-i& w){\left (}{\dfrac{1}{2}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha}-1)h
+(-i w)^2{\dfrac{1}{2}}{\left (}{\dfrac{1}{4}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha}-2)h^2\nonumber\\
-(-i& w)^3{\dfrac{1}{6}}{\left (}{\dfrac{1}{8}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha}-3)h^3+O{\left (}h^4{\right )}.\end{aligned}$$ From and $$\begin{aligned}
\mathcal{F}[S_\infty[y]](&w)/\hat{y}(w)={\dfrac{{\Gamma}(1-{\alpha})}{2^{{\alpha}}}}(-i w h)^{{\alpha}}
+(-i w){\left (}1-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha})h\\
&-(-i w)^2{\left (}{\dfrac{1}{2}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha}-1)h^2+{\dfrac{{\Gamma}(1-{\alpha})}{24.2^{\alpha}}}(-i w h)^{2+{\alpha}}\\
&+(-i w)^3{\left (}{\dfrac{1}{24}}{\left (}1-{\dfrac{1}{2^{\alpha}}}{\right )}+{\dfrac{1}{2}}{\left (}{\dfrac{1}{4}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\right )}{\zeta}({\alpha}-2)h^3+O{\left (}h^4{\right )}.\end{aligned}$$ By applying inverse Fourier transform we obtain the expansion formula $$\begin{aligned}
\label{9_4}
{\left (}{\dfrac{2}{h}} {\right )}^{\alpha}S_n[y]&= {\Gamma}(1-{\alpha})y_n^{({\alpha})}+2^{\alpha}{\left (}1-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha})y'_n h^{1-{\alpha}}+O{\left (}h^{4-{\alpha}}{\right )}\nonumber\\
&- 2^{\alpha}{\left (}{\dfrac{1}{2}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha}-1)y''_n h^{2-{\alpha}}
+{\dfrac{{\Gamma}(1-{\alpha})}{24}}{\dfrac{d^2}{dx^2}}y_n^{({\alpha})} h^2\\
&+ 2^{\alpha}{\left (}{\dfrac{1}{24}}{\left (}1-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha})+{\dfrac{1}{2}}{\left (}{\dfrac{1}{4}}-{\dfrac{1}{2^{\alpha}}}{\right )}{\zeta}({\alpha}-2){\right )}y'''_n h^{3-{\alpha}}.\nonumber\end{aligned}$$ Denote $${\omega}_0^{({\alpha})}=2^{\alpha},\;  {\omega}_k^{({\alpha})}=2^{\alpha}{\left (}{\dfrac{1}{(2k+1)^{\alpha}}}-{\dfrac{1}{(2k-1)^{\alpha}}}{\right )},\; {\omega}_n^{({\alpha})}=-{\dfrac{2^{\alpha}}{(2n-1)^{\alpha}}},$$ where $k=1,\cdots,n-1$. From $${\left (}{\dfrac{2}{h}}{\right )}^{\alpha}{\dfrac{S_n[y]}{{\Gamma}(1-{\alpha})}}={\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^n {\omega}_k^{({\alpha})} y_{n-k}\approx y^{({\alpha})}_n.$$ From and the properties of the inverse Fourier transform we obtain the expansion formula of approximation for the Caputo derivative.

Let $y(0)=y'(0)=y''(0)=0$. Then $$\begin{aligned}
\label{10_1}
{\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}&\sum_{k=0}^{n}{\omega}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+{\dfrac{1}{{\Gamma}(1-{\alpha})}}{\left (}2^{\alpha}-1{\right )}{\zeta}({\alpha})y'_n h^{1-{\alpha}}+O{\left (}h^{4-{\alpha}}{\right )}\nonumber\\
&- {\dfrac{1}{{\Gamma}(1-{\alpha})}} {\left (}2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1)y''_n h^{2-{\alpha}}
+{\dfrac{1}{24}}{\dfrac{d^2}{dx^2}}y_n^{({\alpha})} h^2\\
&+ {\dfrac{1}{{\Gamma}(1-{\alpha})}}{\left (}{\dfrac{1}{24}}{\left (}2^{\alpha}-1{\right )}{\zeta}({\alpha})+{\dfrac{1}{2}}{\left (}2^{{\alpha}-2}-1{\right )}{\zeta}({\alpha}-2){\right )}y'''_n h^{3-{\alpha}}.\nonumber\end{aligned}$$

Approximation for the Caputo derivative of order $1-{\alpha}$. $$\begin{aligned}
\label{10_2}
\mathcal{A}^{\omega}_n[y]={\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}{\omega}_k^{({\alpha})} y_{n-k}=y_n^{({\alpha})}+O{\left (} h^{1-{\alpha}} {\right )}.\end{aligned}$$

The function $y(x)=1$ has Caputo derivative $y^{({\alpha})}(x)=0$. Approximation satisfies $\mathcal{A}^{\omega}_n[1]=0$, because $\sum_{k=0}^n {\omega}_k^{({\alpha})}=0$. Denote by Equation I, Equation II and Equation III the fractional relaxation equations with exact solutions $y(x)=1+x+x^2+x^3+x^4,y(x)=e^x$ and $y(x)=\cos (2\pi x)$. $$\begin{aligned}
\tag{{\bf Equation I}}
y^{({\alpha})}(x)+y(x)=1&+x+x^2+x^3+x^4+{\dfrac{x^{1-{\alpha}}}{{\Gamma}(2-{\alpha})}}\nonumber\\
&+{\dfrac{2 x^{2-{\alpha}}}{{\Gamma}(3-{\alpha})}}
+{\dfrac{6 x^{3-{\alpha}}}{{\Gamma}(4-{\alpha})}}+{\dfrac{24 x^{4-{\alpha}}}{{\Gamma}(5-{\alpha})}},\; y(0)=1.\nonumber\end{aligned}$$ $$\begin{aligned}
 \tag{{\bf Equation II}}
y^{({\alpha})}(x)+y(x)=e^x+x^{1-{\alpha}}E_{1,2-{\alpha}}(x),\quad y(0)=1.\end{aligned}$$ $$\begin{aligned}
&y^{({\alpha})}(x)+y(x)=\cos (2\pi x)+i \pi x^{1-{\alpha}}
{\left (}E_{1,2-{\alpha}}(2\pi i x)-E_{1,2-{\alpha}}(-2\pi i x){\right )},\nonumber\\
& y(0)=1.\tag{{\bf Equation III}}\nonumber\end{aligned}$$

Approximations for the Caputo Derivative of Order $\boldsymbol{2-{\alpha}},\boldsymbol{2}$
==========================================================================================

In the present section we use approximation to obtain approximations for the Caputo derivative of order $2-{\alpha}$ and $2$ and we apply the approximations for computing the numerical solution of the fractional relaxation equation $$\begin{aligned}
\label{11_1}
y^{({\alpha})}(x)+y(x)=F(x),\quad y(0)=y_0.\end{aligned}$$ In [@Dimitrov2015] we showed that $$\bar{y}_1={\dfrac{y(0)+\Gamma(2-{\alpha})h^{\alpha}F(h)}{1+\Gamma(2-{\alpha})h^{\alpha}}}$$ is a second order approximation for the value of the solution $y(h)$. Let $${\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}{\lambda}_k^{({\alpha})} y_{n-k}\approx y_n^{({\alpha})}$$ be an approximation for the Caputo derivative. By approximating the Caputo derivative in equation we obtain $${\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}{\lambda}_k^{({\alpha})} y_{n-k}+y_n\approx F_n,$$ $${\left (}{\lambda}_0^{({\alpha})}+{\Gamma}(1-{\alpha})h^{\alpha}{\right )} y_{n}+\sum_{k=1}^{n}{\lambda}_k^{({\alpha})} y_{n-k}\approx {\Gamma}(1-{\alpha})h^{\alpha}F_n.$$ Let $u_n$ be an approximation for the value of the exact solution $y_n=y(x_n)$ of equation . The numbers $u_n$ are computed with $u_0=y(0),u_1=\bar{y}_1$ and $$\begin{aligned}
\tag{$\boldsymbol{NS[{\lambda}]}$}
 u_n={\dfrac{1}{{\lambda}_0^{({\alpha})}+{\Gamma}(1-{\alpha})h^{\alpha}}} {\left (}{\Gamma}(1-{\alpha})h^{\alpha}F_n  -\sum_{k=1}^{n}{\lambda}_k^{({\alpha})} u_{n-k}{\right )}.\end{aligned}$$

 Approximation for the Caputo derivative of order $\boldsymbol{2-{\alpha}}$
---------------------------------------------------------------------------

By approximating $y'_n$ in with first order backward difference $$y'_n={\dfrac{y_n-y_{n-1}}{h}}-{\dfrac{h}{2}}y''_n+O(h^2)$$ we obtain the approximation for the Caputo derivative $$\begin{aligned}
\label{12_1}
\mathcal{A}_n^{\bar{{\sigma}}}[y]={\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}\bar{{\sigma}}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (}h^{2-{\alpha}}{\right )},\end{aligned}$$ where $\bar{{\sigma}}_{k}^{({\alpha})}=\bar{{\omega}}_{k}^{({\alpha})}$ for $k=2,\cdots,n$ and $$\bar{{\sigma}}_0^{({\alpha})}={\omega}_0^{({\alpha})}-{\left (}2^{\alpha}-1{\right )}{\zeta}({\alpha}),\quad \bar{{\sigma}}_1^{({\alpha})}={\omega}_1^{({\alpha})}+{\left (}2^{\alpha}-1{\right )}{\zeta}({\alpha}).$$ Approximation has accuracy $O{\left (} h^{2-{\alpha}} {\right )}$ when the function $y$ satisfies $y(0)=y'(0)=0$, and $\mathcal{A}_n^{\bar{{\sigma}}}[1]=0$, because $\sum_{k=0}^{n}\bar{{\sigma}}_k^{({\alpha})}=0$. Denote $$W_n =\left({\zeta}({\alpha})-\frac{n^{1-{\alpha}}}{1-{\alpha}}+2^{\alpha}{\left (}\sum _{k=1}^{n} {\dfrac{1}{(2k-1)^{\alpha}}}-\zeta ({\alpha}){\right )}\right).$$

Let $y(x)=x$. Then $$\mathcal{A}_n^{\bar{{\sigma}}}[y]-y^{({\alpha})}(x)={\dfrac{h^{1-{\alpha}}W_n}{{\Gamma}(1-{\alpha})}}.$$

$${\Gamma}(1-{\alpha})h^{\alpha}\mathcal{A}_n^{\bar{{\sigma}}}[x]=\sum_{k=0}^{n}\bar{{\sigma}}_k^{({\alpha})}(x-k h) =x\sum_{k=0}^{n}\bar{{\sigma}}_k^{({\alpha})}-h\sum_{k=0}^{n}k\bar{{\sigma}}_k^{({\alpha})}=-h\sum_{k=1}^{n}k\bar{{\sigma}}_k^{({\alpha})},$$ $${\dfrac{{\Gamma}(1-{\alpha})}{h^{1-{\alpha}}}}\mathcal{A}_n^{\bar{{\sigma}}}[x]={\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})+2^{\alpha}\sum_{k=1}^{n-1}{\left (} {\dfrac{k}{(2k-1)^{\alpha}}}-  {\dfrac{k}{(2k+1)^{\alpha}}} {\right )} +{\dfrac{2^{\alpha}n}{(2 n-1)^{\alpha}}},$$ $$\begin{aligned}
\label{12_2}
\mathcal{A}_n^{\bar{{\sigma}}}[x]={\dfrac{h^{1-{\alpha}}}{{\Gamma}(1-{\alpha})}}{\left (}{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})+2^{\alpha}\sum_{k=1}^{n}{\dfrac{1}{(2 k-1)^{\alpha}}}{\right )}.\end{aligned}$$ The function $y(x)=x$ has fractional derivative $$\begin{aligned}
\label{13_1}
y^{({\alpha})}(x)={\dfrac{x^{1-{\alpha}}}{{\Gamma}(2-{\alpha})}}={\dfrac{n^{1-{\alpha}}h^{1-{\alpha}}}{(1-{\alpha}){\Gamma}(1-{\alpha})}}.\end{aligned}$$ The statement of Claim 3 follows from and .

The function $y(x)=x^{2}$ has fractional derivative $y^{({\alpha})}(x)=x^{2-{\alpha}}/{\Gamma}(3-{\alpha})$. Approximation for the Caputo derivative has accuracy $O{\left (}h^{2-{\alpha}}{\right )}$ for small values of $n$ when $\mathcal{A}_n^{{\sigma}}[1]=0$ and $\mathcal{A}_n^{{\sigma}}[x]=x^{1-{\alpha}}/{\Gamma}(2-{\alpha})$. Represent the function $y(x)$ as $$\begin{aligned}
\label{13_05}
y(x)=y(x)-y(0)-y'(0)x+y(0)+y'(0)x=z(x)+y(0)+y'(0)x.\end{aligned}$$ The function $z(x)=y(x)-y(0)-y'(0)x$ satisfies $z(0)=z'(0)=0$. Then $$\mathcal{A}_n^{\bar{{\sigma}}}[z(x)]=z^{({\alpha})}(x)+O{\left (} h^{2-{\alpha}} {\right )}=y^{({\alpha})}(x)-y'(0){\dfrac{x^{1-{\alpha}}}{{\Gamma}(2-{\alpha})}}+O{\left (} h^{2-{\alpha}} {\right )}.$$ From $$\mathcal{A}_n^{\bar{{\sigma}}}[y(x)]=\mathcal{A}_n^{\bar{{\sigma}}}[z(x)]+\mathcal{A}_n^{\bar{{\sigma}}}[y(0)+y'(0)x],$$ $$\mathcal{A}_n^{\bar{{\sigma}}}[y(x)]=y^{({\alpha})}(x)-y'(0){\dfrac{x^{1-{\alpha}}}{{\Gamma}(2-{\alpha})}}+y'(0)\mathcal{A}_n^{\bar{{\sigma}}}[x]+O{\left (} h^{2-{\alpha}} {\right )}.$$ From Claim 3 $$\begin{aligned}
{\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n}\bar{{\sigma}}_k y_{n-k}=y_n^{({\alpha})}+{\dfrac{y'(0) W_n h^{1-{\alpha}}}{{\Gamma}(1-{\alpha})}}+O{\left (} h^{2-{\alpha}} {\right )}.\end{aligned}$$ By approximating $y'_0$ with $y'_0=(y_1-y_{0})/h+O(h)$ we obtain $$\begin{aligned}
\label{13_2}
\mathcal{A}_n^{{\sigma}}[y]={\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n} {\sigma}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (} h^{2-{\alpha}} {\right )},\end{aligned}$$ where ${\sigma}^{({\alpha})}_k=\bar{{\sigma}}^{({\alpha})}_k$ for $k=0,\cdots,n-2$ and $${\sigma}^{({\alpha})}_{n-1}=\bar{{\sigma}}^{({\alpha})}_{n-1}-W_n,\quad
 {\sigma}^{({\alpha})}_{n}=\bar{{\sigma}}^{({\alpha})}_{n}+W_n.$$ The weights ${\sigma}_{k}^{({\alpha})}$ of approximation satisfy . When $n=2$ $${\sigma}_{0}^{({\alpha})}=1-2 \zeta ({\alpha}),{\sigma}_{1}^{({\alpha})}=4 \zeta ({\alpha})+\frac{2^{2-{\alpha}}}{1-{\alpha}}-2,{\sigma}_{2}^{({\alpha})}=-2 \zeta ({\alpha})-\frac{2^{2-{\alpha}}}{1-{\alpha}}+1.$$

Second-Order Approximation for the Caputo derivative 
-----------------------------------------------------

By approximating $y'_n$ and $y''_n$ in with $$y'_n={\dfrac{1}{h}}{\left (}{\dfrac{3}{2}}y_n-2y_{n-1}+{\dfrac{1}{2}}y_{n-2}   {\right )}
    +O{\left (}h^2{\right )},$$ $$y''_n={\dfrac{1}{h^2}}{\left (}y_n-2y_{n-1}+y_{n-2}   {\right )}
    +O{\left (}h{\right )},$$ we obtain the approximation for the Caputo derivative $$\begin{aligned}
\label{14_1}
{\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n} \bar{{\delta}}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (} h^{2} {\right )},\end{aligned}$$ where $\bar{{\delta}}_k^{({\alpha})}={\omega}_k^{({\alpha})}$ for $k=3,\cdots,n$ and $$\bar{{\delta}}_0^{({\alpha})}={\omega}_0^{({\alpha})}+{\dfrac{3}{2}}{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})+{\left (} 2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1),$$ $$\bar{{\delta}}_1^{({\alpha})}={\omega}_1^{({\alpha})}-2{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})-2{\left (} 2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1),$$ $$\bar{{\delta}}_2^{({\alpha})}={\omega}_2^{({\alpha})}+{\dfrac{1}{2}}{\left (}1-2^{\alpha}{\right )}{\zeta}({\alpha})+{\left (} 2^{{\alpha}-1}-1{\right )}{\zeta}({\alpha}-1).$$ Approximation has accuracy $O{\left (} h^{2} {\right )}$ when $y(0)=y'(0)=0$. From approximation we obtain the second order approximation approximation for the Caputo derivative $$\begin{aligned}
{\dfrac{1}{{\Gamma}(1-{\alpha})h^{\alpha}}}\sum_{k=0}^{n} {\delta}_k^{({\alpha})} y_{n-k}= y_n^{({\alpha})}+O{\left (} h^{2} {\right )},\end{aligned}$$ where ${\delta}_{k}^{({\alpha})}=\bar{{\delta}}^{({\alpha})}_{k}$ for $k=0,1,\cdots,n-2$ and $${\delta}^{({\alpha})}_{n-1}=\bar{{\delta}}^{({\alpha})}_{n-1}-W_n, \quad {\delta}^{({\alpha})}_{n}=\bar{{\delta}}^{({\alpha})}_{n}+W_n.$$ The formulas for the weights ${\delta}^{({\alpha})}_{n-1}$ and ${\delta}^{({\alpha})}_{n}$ are derived using the method for computing the weights ${\sigma}^{({\alpha})}_{n-1}$ and ${\sigma}^{({\alpha})}_{n}$ from the previous section. In Figure 2 we compare numerical solutions $NS[{\omega}],NS[{\sigma}]$ and $NS[{\delta}]$ of Equation III and ${\alpha}=0.6$. In Table 4 we compute the error and the order of numerical solution $NS[{\delta}]$ of Equation I and ${\alpha}=0.25$, Equation II and ${\alpha}=0.5$ and Equation III with ${\alpha}=0.75$.

![[Graph of the exact solution of Equation III and numerical solutions $NS[{\omega}]$(green), $NS[{\sigma}](red)$ and $NS[{\delta}]$(blue) for ${\alpha}=0.6,h=0.1$.]{}](graph2.eps){width="60.00000%"}

Numerical Solution of the Fractional Subdiffusion Equation
==========================================================

The analytical and numerical solutions of the fractional subdiffusion equation have been studied extensively [@Dimitrov2015; @JiSun2015; @JinLazarovZhou2016; @LinXu2007; @Ma2014; @SunWu2006; @ZengLiLiuTurner2015]. In this section we use approximation for the Caputo derivative to construct a finite-difference scheme for the fractional subdiffusion equation $$\label{16_1}
    \left\{
    \begin{array}{l l}
\dfrac{\partial^\alpha u(x,t)}{\partial t^\alpha}=D\dfrac{\partial^2 u(x,t)}{\partial x^2}+F(x,t),\quad (x,t)\in [0,1]\times[0,T],& \\
    u(x,0)=u_0(x),\;u(0,t)=u_L(t),\;u(1,t)=u_R(t).&  \\
    \end{array} 
        \right .$$ Let $h=1/N,{\tau}=T/M$, where $M$ and $N$ are positive integers, and $\mathcal{G}$ be a grid on the square $[0,1]\times[0,T]$ $$\mathcal{G}={\left \{}(n h,m {\tau}) \| 1\leq n\leq N, 1\leq m\leq M{\right \}}.$$ Denote $u_n^m=u(n h,m {\tau})$ and $F_n^m=F(n h,m{\tau}).$ By approximating the Caputo derivative in the time direction using and the second derivative in the space direction by second-order central difference approximation we obtain $${\dfrac{1}{{\Gamma}(1-{\alpha}){\tau}^{\alpha}}}\sum_{k=0}^n {\sigma}_k^{({\alpha})}u_n^{m-k} = D {\dfrac{u_{n-1}^m-2 u_n^m+u_{n+1}^m}{h^2}}+F_n^m+O{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}.$$ Let $\eta={\Gamma}(1-{\alpha})D{\tau}^{\alpha}/h^2.$ The numerical solution $\{U_n^m\}_{n=1}^{N-1}$ of equation satisfies the system of equations $$-\eta U_{n-1}^m+{\left (}{\sigma}_0^{({\alpha})}+2\eta{\right )} U_n^m-\eta U_{n+1}^m=-\sum_{k=1}^n {\sigma}_k^{({\alpha})}U_n^{m-k}+{\tau}^{\alpha}{\Gamma}(1-{\alpha}) F_n^m.$$ Let $\mathcal{K}$ be a tridiagonal matrix of dimension $N-1$ with values ${\sigma}_0^{({\alpha})}+2\eta$ on the main diagonal, and $-\eta$ on the diagonals above and below the main diagonal. The vector $\mathcal{U}^m={\left (}U_1^m,U_2^m,\cdots,U_{N-1}^m{\right )}^T$ of the numerical solution on the $m$-th layer of the grid $\mathcal{G}$ is a solution of the linear system $$\begin{aligned}
\label{16_2}
    \mathcal{K} \mathcal{U}^m=\mathcal{R}_1+\eta \mathcal{R}_2,
    \end{aligned}$$ where $\mathcal{R}_1$ and $\mathcal{R}_2$ are the column vectors of dimension $N-1$ $$\mathcal{R}_1^T={\left [}-\sum_{k=1}^n {\sigma}_k^{({\alpha})}U_n^{m-k}+{\Gamma}(1-{\alpha}){\tau}^{\alpha}F(n h,m{\tau}) {\right ]}_{n=1}^{N-1},$$ $$\mathcal{R}_2^T={\left [}u_L(m{\tau}),0,\cdots,0,u_R(m{\tau}){\right ]}^T.$$ The numerical solution of the fractional subdiffusion equation on the first layer of the grid $\mathcal{G}$ is computed with the approximation [@Dimitrov2015] $$\label{17_05}
y^{({\alpha})}({\tau})={\dfrac{y({\tau})-y(0)}{{\tau}^{\alpha}{\Gamma}(2-{\alpha})}}+O{\left (}{\tau}^{2-{\alpha}}{\right )}.$$ Let $\tilde{\eta}={\Gamma}(2-{\alpha})D{\tau}^{\alpha}/h^2$. The numerical solution on the first layer of the grid $\mathcal{G}$ satisfies the system of equations $$\left\{
    \begin{array}{l l}
U_0^1=u_L({\tau}),U_N^1=u_R({\tau}),\quad(n=1\cdots,N-1),& \\
-\tilde{\eta} U_{n-1}^1+(1+2\tilde{\eta}) U_n^1-\tilde{\eta} U_{n+1}^1=U_n^0+{\Gamma}(2-{\alpha}){\tau}^{\alpha}F_n^1.&  \\
    \end{array} 
        \right .$$ [**Example 1:**]{} The fractional subdiffusion equation $$\label{17_1}
    \left\{
    \begin{array}{l l}
\dfrac{\partial^\alpha v(x,t)}{\partial t^\alpha}=\dfrac{\partial^2 v(x,t)}{\partial x^2}+e^x{\left (}t^{1-{\alpha}}E_{1,2-{\alpha}}(t)-e^t  {\right )},& \\
    u(x,0)=e^x,\;u(0,t)=e^t, u(1,t)=e^{t+1},\quad (x,t)\in[0,1]\times[0,1],&  \\
    \end{array} 
        \right .$$ has the solution $u(x,t)=e^{x+t}$.

[**Example 2:**]{} The fractional subdiffusion equation $$\label{17_2}
    \left\{
    \begin{array}{l l}
\dfrac{\partial^\alpha u(x,t)}{\partial t^\alpha}={\dfrac{1}{\pi^2}}\dfrac{\partial^2 u(x,t)}{\partial x^2},\quad (x,t)\in[0,1]\times[0,1],& \\
    u(x,0)=\sin (\pi x),\;u(0,t)=0, u(1,t)=0,&  \\
    \end{array} 
        \right .$$ has the solution $u(x,t)=\sin (\pi x) E_{{\alpha}}{\left (}-t^{\alpha}{\right )}$. The first partial derivative of the solution is unbounded at $t=0$.

In Table 6 we compute the error and the order of numerical solution for equation . The numerical results from Table 6 suggest that the accuracy of numerical solution is $O(h)$. A numerical analysis of the scheme using the $L1$ approximation for the Caputo derivative is given by Jin, Lazarov and Zhou [@JinLazarovZhou2016]. We use the fractional Taylor polynomials of the solution in order to improve the differentiability class of the solution and the accuracy of the numerical solution. Denote by $D_t^{[{\beta}]} u(x,t)$ the Miller-Ross derivative of order ${\beta}$ of the solution $u(x,t)$ in the time direction. From equation $$D_t^{[{\alpha}]} u(x,t)={\dfrac{1}{\pi^2}}\dfrac{\partial^2 u(x,t)}{\partial x^2}.$$ By applying fractional differentiation of order ${\alpha}$ we obtain $$D_t^{[2{\alpha}]}u(x,t)=D_t^{{\alpha}}D_t^{{\alpha}}u(x,t)={\dfrac{1}{\pi^2}}
    D_t^{\alpha}\dfrac{\partial^2 u(x,t)}{\partial x^2} ={\dfrac{1}{\pi^4}}\dfrac{\partial^4 u(x,t)}{\partial x^4},$$ $$D_t^{[3{\alpha}]}u(x,t)=D_t^{{\alpha}}D_t^{[2{\alpha}]}u(x,t)={\dfrac{1}{\pi^4}}D_t^{\alpha}\dfrac{\partial^4 u(x,t)}{\partial x^4}={\dfrac{1}{\pi^6}}\dfrac{\partial^6 u(x,t)}{\partial x^6}.$$ By induction we obtain $$D_t^{[n{\alpha}]}u(x,t)={\dfrac{1}{\pi^{2n}}}\dfrac{\partial^{2n} u(x,t)}{\partial x^{2n}}.$$ Set $t=0$ $$D_t^{[n{\alpha}]}u(x,0)={\dfrac{1}{\pi^{2n}}}\dfrac{\partial^{2n} u(x,0)}{\partial x^{2n}}={\dfrac{1}{\pi^{2n}}}
        \dfrac{\partial^{2n} \sin x}{\partial x^{2n}}=(-1)^{n}\sin (\pi x).$$ The solution of the fractional subdiffusion equation has fractional Taylor polynomials in the time direction $$T^{({\alpha})}_m(x,t)=  \sum_{n=0}^m {\dfrac{t^{n{\alpha}}D_t^{[n{\alpha}]}u(x,0)}{{\Gamma}(n{\alpha}+1)}}=\sin (\pi x) \sum_{n=0}^m (-1)^n {\dfrac{h^{n{\alpha}}}{{\Gamma}(n{\alpha}+1)}}.$$ Substitute $$v(x,t)=u(x,t)-\sin (\pi x) \sum_{n=0}^m (-1)^n {\dfrac{t^{n{\alpha}}}{{\Gamma}(n{\alpha}+1)}}.$$ By differentiating $v(x,t)$ we obtain $$D_t^{\alpha}v(x,t)=D_t^{\alpha}u(x,t)+\sin (\pi x) \sum_{n=0}^{m-1} (-1)^{n+1} {\dfrac{t^{n{\alpha}}}{{\Gamma}(n{\alpha}+1)}},$$ $$\dfrac{\partial^{2}}{\partial x^{2}}v(x,t)=\dfrac{\partial^{2}}{\partial x^{2}}u(x,t)+\pi^2 \sin (\pi x) \sum_{n=0}^{m} (-1)^{n+1} {\dfrac{t^{n{\alpha}}}{{\Gamma}(n{\alpha}+1)}}.$$ The function $v(x,t)$ is a solution of the fractional subdiffusion equation $$\label{19_1}
    \left\{
    \begin{array}{l l}
    \dfrac{\partial^{{\alpha}} v(x,t)}{\partial t^{{\alpha}}}={\dfrac{1}{\pi^2}}\dfrac{\partial^2 v(x,t)}{\partial x^2}+(-1)^{m+1}\sin (\pi x){\dfrac{t^{m{\alpha}}}{{\Gamma}(m{\alpha}+1)}},&  \\
    v(x,0)=v(0,t)=v(\pi,t)=0.&  \\
    \end{array}
        \right .$$ The solution $v(x,t)$ of equation has a continuous second partial derivative in the time direction at the point $t=0$, when $m{\alpha}\geq 2$. In Table 7 we compute the error and the order of numerical solution of equation and ${\alpha}=0.25,=0.5,{\alpha}=0.75$. The accuracy of numerical solution of equation is $O{\left (}h^{2-{\alpha}}{\right )}$ and the error is significantly smaller than the error of the numerical solution of equation in Table 6.

In Theorem 5 we establish the convergence of numerical solution of the fractional subdiffusion equation. The proof uses estimate for ${\sigma}_m^{({\alpha})}$.

$$W_m={\dfrac{{\alpha}}{24 m^{1+{\alpha}}}}+O{\left (} {\dfrac{1}{m^{2+{\alpha}}}}{\right )}.$$

Denote $S_m=\sum_{k=1}^m 1/k^{\alpha}$. From the formula for the sum of powers $$S_{m-1}={\zeta}({\alpha})+{\dfrac{m^{1-{\alpha}}}{1-{\alpha}}}\sum_{m=0}^\infty \binom{1-{\alpha}}{m}{\dfrac{B_m}{n^m}},$$ $$\begin{aligned}
\label{20_1}
    S_m={\zeta}({\alpha})+\frac{m^{1-{\alpha}}}{1-{\alpha}}+{\dfrac{1}{2 m^{\alpha}}}-{\dfrac{{\alpha}}{12m^{1+{\alpha}}}}+O{\left (} {\dfrac{1}{m^{2+{\alpha}}}}{\right )}.
    \end{aligned}$$ We have that $$S_{2m}=\sum_{k=1}^{2m} 1/k^{\alpha}=\sum _{k=1}^{m} {\dfrac{1}{(2k-1)^{\alpha}}}+\sum _{k=1}^{m} {\dfrac{1}{(2k)^{\alpha}}}=\sum _{k=1}^{m} {\dfrac{1}{(2k-1)^{\alpha}}}+{\dfrac{1}{2^{\alpha}}}S_m.$$ From $$S_{2m}-{\dfrac{1}{2^{\alpha}}}S_m={\zeta}({\alpha}){\left (} 1-{\dfrac{1}{2^{\alpha}}}{\right )}+{\dfrac{m^{1-{\alpha}}}{(1-{\alpha})2^{\alpha}}}+{\dfrac{{\alpha}}{12 (2m)^{1+{\alpha}}}}+O{\left (} {\dfrac{1}{m^{2+{\alpha}}}}{\right )}.$$ Hence $$W_m={\zeta}({\alpha})-\frac{m^{1-{\alpha}}}{1-{\alpha}}+2^{\alpha}{\left (}S_{2m}-{\dfrac{1}{2^{\alpha}}}S_m-\zeta ({\alpha}){\right )}={\dfrac{{\alpha}}{24m^{1+{\alpha}}}}+O{\left (} {\dfrac{1}{m^{2+{\alpha}}}}{\right )}.$$

We can further show that $W_m<{\alpha}/{\left (}24 m^{1+{\alpha}}{\right )}$. From Claim 4 we obtain $$\begin{aligned}
\label{20_2}
    {\left |}{\sigma}_m^{({\alpha})}{\right |}>{\dfrac{2^{\alpha}}{(2m-1)^{\alpha}}}-W_m>{\dfrac{1}{m^{\alpha}}}-\frac{{\alpha}}{24 m^{1+{\alpha}}}>{\dfrac{23}{24 m^{\alpha}}}>{\dfrac{1}{2m^{\alpha}}}.
    \end{aligned}$$ The maximum (infinity) norm of the vector $\mathcal{V}={\left (}v_i{\right )}$ and the square matrix $\mathcal{L}={\left (}l_{ij}{\right )}$ of dimension $N-1$ are defined as $${\left \|} \mathcal{V}{\right \|}=\max_{1\leq i\leq N-1} |v_i|,\quad {\left \|} \mathcal{L} {\right \|}=\max_{1\leq j\leq N-1} \sum_{m=1}^{N-1}|l_{ij}|.$$ The matrix $\mathcal{K}$ is a diagonally dominant tridiagonal matrix with positive elements on the main diagonal and negative elements on the diagonals below and above the main diagonal. The matrix $\mathcal{K}^{-1}$ is a positive matrix. From the Ahlberg-Nilson-Varah bound [@AhlbergNilson1963; @Varah1975] $${\left \|} \mathcal{K}^{-1} {\right \|}\leq {\dfrac{1}{{\sigma}_0^{({\alpha})}}}.$$ The numbers ${\sigma}_k^{({\alpha})}$ satisfy $\sum_{k=1}^m {\left |} {\sigma}_k^{({\alpha})} {\right |}={\sigma}_0^{({\alpha})}$. From $$\begin{aligned}
\label{21_1}
\sum_{k=1}^{m-1} {\left |} {\sigma}_k^{({\alpha})} {\right |}={\sigma}_0^{({\alpha})}-{\left |}{\sigma}_m^{({\alpha})}{\right |}<{\sigma}_0^{({\alpha})}-{\dfrac{1}{2m^{\alpha}}}.\end{aligned}$$ Let $e_n^m=u_n^m-U_n^m$ be the error of numerical solution . The error vector $\mathcal{E}^m=(e_n^m)$ on the $m$-th layer of the grid $\mathcal{G}$ satisfies the system of equations $$\mathcal{K} \mathcal{E}^m=\mathcal{R}^m,$$ where $\mathcal{R}^m=(r^m_n)$ is an $N-1$ dimensional column vector with elements $$r_n^m=-\sum_{k=1}^{m-1} {\sigma}_k^{({\alpha})}e_n^{m-k}+{\tau}^{\alpha}{\left (}A_n^m {\tau}^{2-{\alpha}} +B_n^m h^2{\right )},$$ and $A_n^m {\tau}^{2-{\alpha}} +B_n^m h^2$ is the approximation error at the point $(n h,m{\tau})$. Let $A$ be a positive constant such that ${\left |}A_n^m{\right |}<A$ and ${\left |}B_n^m{\right |}<A$, for all $m$ and $n$. The numerical solution on the first layer of $\mathcal{G}$ is computed using approximation , which is obtained from the the $L1$ approximation when $n=1$. The convergence of the numerical solution of the time-fractional subdiffusion equation which uses the $L1$ approximation for the Caputo derivative is studied in [@LinXu2007; @Ma2014]. We can assume that the number $A$ is large enough, such that holds for $m=1$.

The error on the $m$-th layer of the grid $\mathcal{G}$ satisfies $$\begin{aligned}
\label{21_2}
{\left \|} \mathcal{E}^m{\right \|}\leq 2A m^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}.\end{aligned}$$

Induction on $m$. Suppose that holds for $k=1,\cdots,m-1$. $${\left |}r_n^m{\right |}\leq \sum_{k=1}^{m-1} {\left |}{\sigma}_k^{({\alpha})}{\right |}{\left |}e_n^{m-k}{\right |}+A {\tau}^{\alpha}{\left (} {\tau}^{2-{\alpha}} +h^2{\right )}.$$ From the induction assumption and $${\left |}r_n^m{\right |}\leq 2A m^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}\sum_{k=1}^{m-1} {\left |}{\sigma}_k^{({\alpha})}{\right |}+A {\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}} +h^2{\right )},$$ $${\left |}r_n^m{\right |}\leq 2A m^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}{\left (}{\sigma}_0^{({\alpha})}-{\dfrac{1}{2 m^{\alpha}}}{\right )}+A {\tau}^{\alpha}{\left (} {\tau}^{2-{\alpha}} +h^2{\right )},$$ $${\left |}r_n^m{\right |}\leq 2A{\sigma}_0^{({\alpha})} m^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )},$$ for all $n=1,\cdots,N-1$. Then $${\left \|} \mathcal{R}^m{\right \|}\leq 2A{\sigma}_0^{({\alpha})} m^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}.$$ We have that $\mathcal{E}^m=\mathcal{K}^{-1}\mathcal{R}^m$. Hence $${\left \|} \mathcal{E}^m{\right \|}\leq 
{\left \|} \mathcal{K}^{-1}{\right \|}{\left \|} \mathcal{R}^m{\right \|}\leq {\dfrac{1}{{\sigma}_0^{({\alpha})}}}{\left \|} \mathcal{R}^m{\right \|}\leq 2A m^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}.$$

From we obtain the estimate for the error on the grid $\mathcal{G}$ $${\left \|} \mathcal{E}^m{\right \|}\leq  2A M^{\alpha}{\tau}^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )}\leq  2A T^{\alpha}{\left (}{\tau}^{2-{\alpha}}+h^2{\right )},$$ for all $m=1,2,\cdots,M$.

[99]{} M. Abramowitz, I. A. Stegun: [*Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables,*]{} Dover, New York, 1964. J. H. Ahlberg, E. N. Nilson: [*Convergence Properties of the Spline Fit,*]{} Journal of the Society for Industrial and Applied Mathematics, [**11(1)**]{}, (1963), 95–104. D. Baleanu, K. Diethelm, E. Scalas and J. J. Trujillo: [*Fractional Calculus: Models and Numerical Methods,*]{} World Scientific, (2012). K. Diethelm: [*The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type,*]{} Springer; 2010. Y. Dimitrov: [*A second order approximation for the Caputo fractional derivative,*]{} [**arXiv:1502.00719**]{} (2015). Y. Dimitrov: [*Higher-Order Numerical Solutions of the Fractional Relaxation-Oscillation Equation using Fractional Integration,*]{} [**arXiv:1603.08733**]{} (2016). Y. Dimitrov: [*Approximations for the Caputo derivative (I)*]{}, In preparation (2016). H. Ding, C. Li: [*High-order algorithms for Riesz derivative and their applications (III)*]{}, Fract. Calc. Appl. Anal., 19(1), (2016), 19–55. G. H. Gao, Z. Z. Sun, H. W. Zhang: [*A new fractional numerical differ-entiation formula to approximate the Caputo fractional derivative and its applications,*]{} J. Comput. Phys., [**259**]{} (2014), 33–50. R. Gorenflo, F. Mainardi: [*FRACTIONAL CALCULUS: Integral and Differential Equations of Fractional Order*]{}, CISM Lecture notes (2008). M. Gülsu, Y. Öztürk, A. Anapal[i]{}: [*Numerical approach for solving fractional relaxation-oscillation equation,*]{} Applied Mathematical Modelling, [**37(8)**]{} (2013), 5927–5937. C.-c. Ji, Z.-z. Sun: [*A High-Order Compact Finite Difference Scheme for the Fractional Sub-diffusion Equation,*]{} Journal of Scientific Computing, [**64(3)**]{}, (2015), 959–985. B. Jin, R. Lazarov, J. Pasciak and Z. Zhou: [*Error analysis of semidiscrete FEM for inhomogeneous time-fractional diffusion,*]{} IMA J. Numer. Anal., [**35(2)**]{} (2015), 561–582 B. Jin, R. Lazarov and Z. Zhou: [*An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data,*]{} IMA J. Numer. Anal., [**36 (1)**]{} (2016), 197–221. C. Lubich: [*Discretized fractional calculus*]{}, SIAM J. Math. Anal., [**17**]{} (1986), 704–719. C. Li, R. Wu, H. Ding: [*High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations,*]{} Communications in Applied and Industrial Mathematics, [**6(2)**]{} (2014), e-536. Y. Lin, C. Xu: [*Finite difference/spectral approximations for the time-fractional diffusion equation,*]{} Journal of Computational Physics, [**225(2)**]{} (2007), 1533–1552. Y. Ma: [*Two implicit finite difference methods for time fractional diffusion equation with source term,*]{} Journal of Applied Mathematics and Bioinformatics, [**4(2)**]{} (2014), 125–145. Z. Odibat : [*Approximations of fractional integrals and Caputo fractional derivatives,*]{} Applied Mathematics and Computation, [178(2)]{} (2006), 527–533. I. Podlubny: [*Fractional Differential Equations,*]{} Academic Press, San Diego (1999). S. Pooseh, R. Almeida and D. F.M. Torres: [*Approximation of fractional integrals by means of derivatives,*]{} Computers & Mathematics with Applications, [**64(10)**]{} (2012), 3090–3100. Z.-z. Sun, X. Wu: [*A fully discrete scheme for a diffusion wave system,*]{} Appl. Numer. Math., 56(2) (2006), 193–209. C. Tadjeran, M. M. Meerschaert and H. P. Scheffer: [*A second-order accurate numerical approximation for the fractional diffusion equation,*]{} Journal of Computational Physics, [**213**]{} (2006), 205–213. W.Y. Tian, H. Zhou and W. Deng: [*A class of second order difference approximation for solving space fractional diffusion equations,*]{} Mathematics of Computation, [**84(294)**]{} (2012), 1703–1727. J.M. Varah: [*A lower bound for the smallest singular value of a matrix,*]{} Linear Algebra and its Applications, [**11(1)**]{}, (1975), 3–5. R. Wu, H. Ding and C. Li: [*Determination of coefficients of high-order schemes for Riemann-Liouville derivative,*]{} The Scientific World Journal (2014), Article ID 402373. Y. Yan, K. Pal and N. J. Ford: [*Higher order numerical methods for solving fractional differential equations,*]{} BIT Numerical Mathematics, [**54(2)**]{} (2014), 555–584. F. Zeng, C. Li, F. Liu and I. Turner: [*Numerical Algorithms for Time-Fractional Subdiffusion Equation with Second-Order Accuracy,*]{} SIAM Journal on Scientific Computing, 37(1), (2015), A55–A78.
