\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2014 (2014), No. 97, pp. 1--17.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2014 Texas State University - San Marcos.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2014/97\hfil Well-posedness]
{Well-posedness of fractional parabolic differential and difference
equations \\ with  Dirichlet-Neumann conditions}

\author[A. Ashyralyev, N. Emirov, Z. Cakir \hfil EJDE-2014/97\hfilneg]
{Allaberen Ashyralyev, Nazar Emirov, Zafer Cakir}  % in alphabetical order

\address{Allaberen Ashyralyev \newline
Department of Mathematics, Fatih University,
Buyukcekmece, Istanbul, Turkey}
\email{aashyr@fatih.edu.tr}

\address{Nazar Emirov \newline
Department of Mathematics, Fatih University,
Buyukcekmece, Istanbul, Turkey}
\email{nazaremirov@gmail.com}

\address{Zafer Cakir \newline
Department of Mathematical Engineering,
Gumushane University, Gumushane, Turkey}
\email{zafer@gumushane.edu.tr}

\thanks{Submitted December 26, 2013. Published April 10, 2014.}
\subjclass[2000]{35R11, 35B35, 47B39, 47B48}
\keywords{Fractional parabolic equations; Dirichlet-Neumann conditions;
\hfill\break\indent positive operator; difference schemes; stability}

\begin{abstract}
 We study  initial-boundary value problems for  fractional
 parabolic equations with the Dirichlet-Neumann conditions.
 We obtain a stable difference schemes for this problem, and
 obtain theorems on coercive stability estimates for the
 solution of the first order of accuracy difference scheme.
 A procedure of modified Gauss elimination method is
 applied for the solution of the first and second order of accuracy
 difference schemes of one-dimensional fractional parabolic
 differential equations.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\allowdisplaybreaks

\section{Introduction}

Theory, applications and methods of solutions of problems for
fractional differential equations have been studied extensively 
by many researchers (see, e.g., \cite{4}--\cite{7},
\cite{34}--\cite{5}, \cite{21}, \cite{27}, \cite{44},
\cite{41}--\cite{8}, \cite{10}--\cite{1}, \cite{43}--\cite{35} and
the references given therein). In this article, we study the
initial-boundary value problem
\begin{equation}
\begin{gathered}
D_t^{\alpha }u(t,x)-a(x)u_{xx}(t,x)+\sigma u(t,x)=f(t,x),\quad
0<x<l,\; 0<t<T, \\
u(t,0)=0,\quad u_{x}(t,l)=0, \quad 0\leq t\leq T, \\
u(0,x)=0,\quad 0\leq x\leq l
\end{gathered}  \label{1}
\end{equation}
for the fractional parabolic equation with the Dirichlet-Neumann
conditions. Here $D_t^{\alpha }=D_{0+}^{\alpha }$ is
the standard Riemann-Louville's derivative of order $\alpha \in [ 0,1)$.
 Here $a(x)(x\in (0,l))$ and $f(t,x)(t\in (0,T),x\in (0,l))$ are given
smooth functions, $a(x) \geq a>0$, $\sigma >0$. Theorem on coercive
stability estimates for the solution of the initial-boundary value problem 
\eqref{1} is established. Stable difference schemes for the
approximate solution of problem \eqref{1} are considered. Theorem on
coercive stability estimates for the solution of the first order of
accuracy in $t$ difference scheme is proved. A procedure of modified
Gauss elimination method is applied for the solution of the first
and second order of accuracy difference schemes for the fractional
parabolic equations.

 The organization of the present
paper as follows. The first section is introduction where we provide
the history and formulation of the problem. In Section 2, theorem on
coercivity stability of problem \eqref{1} is established. In Section
3, stable difference schemes for the approximate solution of problem
\eqref{1} are considered. Theorem on coercivity stability for the
first order of accuracy in $t$ difference scheme is proved. In
Section 4, the numerical application is given. Finally, Section 5 is
conclusion.

\section{Theorems on coercive stability}

We will give some statements which will be useful in the sequel.

Let $E$ be a Banach space, and $A:D(A)\subset E\to E$ be a
linear unbounded operator densely defined in $E$. We call $A$
$strongly$ $positive$ in the Banach space $E$, if its spectrum
$\sigma _{A}$ lies in the interior of the sector of angle $\phi $,
$0<2\phi <\pi $, symmetric with respect to the real axis, and if on
the edges of this sector, $S_1(\phi )=\{\rho e^{i\phi }:0\leq \rho
\leq \infty \}$ and $S_2(\phi )=\{\rho e^{-i\phi
}:0\leq \rho \leq \infty \}$, and outside of the sector the resolvent 
$(\lambda -A)^{-1}$ is subject to the bound
\begin{equation}
\| (\lambda -A)^{-1}\| _{E\to E}\leq \frac{M}{1+|\lambda |}.  \label{2}
\end{equation}
The infimum of such angles is called spectral angle $\varphi (A,E)$ of $A$.

 Throughout this article, positive constants have different
values in time and they will be indicated with $M$ On the other hand $
M\left( \alpha ,\beta ,\cdots \right) $ is used to focus on the fact
that the constant depends only on $\alpha ,\beta ,\cdots $. 

For a positive operator $A$ in the Banach space $E$, let us
introduce the fractional spaces 
$E_{\beta }=E_{\beta }(E,A)(0<\beta <1)$ consisting of those $\nu \in E$ 
for which the norm
\[
\| \nu \| _{E_{\beta }}=\sup_{\lambda >0}\lambda ^{\beta}
\| A(\lambda +A)^{-1}\nu \| _{E}+\| \nu\| _{E}
\]
is finite.

\begin{theorem}[\cite{18,13}]\label{thm1}
  Let $A$ and $B$ be two commutative positive operators with 
$\varphi(A,E)+\varphi(B,E)<\pi$. Then it follows that there exists the 
bounded operator $(A+B)^{-1}$ defined on whole space $E$. 
Moreover, for every $\beta\in(0,1)$ and $f$, there exists a unique solution 
$u=u(f)$ of the problem
 \[
   Au+Bu=f
 \]
and the following estimates hold
\begin{gather*}
   \| Au \|_{E_{\beta}(E,B)}+\| Bu \|_{E_{\beta}(E,B)}+\| Bu \|_{E_{\beta}(E,A)}
\leq M(\beta)\| f \|_{E_{\beta}(E,B)}, \\
 \| Au \|_{E_{\beta}(E,A)}+\| Bu \|_{E_{\beta}(E,A)}+\| Au \|_{E_{\beta}(E,B)}\leq M(\beta)\| f \|_{E_{\beta}(E,A)}.
\end{gather*}
 \end{theorem}

 \begin{theorem}[\cite{13}]\label{thm2}
  Let $A$ be the positive operator with $\varphi(A,E)<\pi$. 
Then for $\beta\leq\frac{1}{2},A^{\beta}$ is a positive operator with 
$\varphi(A^{\beta},E)<\frac{\pi}{2}$.
 \end{theorem}

\begin{theorem}[\cite{12}]\label{thm3}
  Let $A$ be the operator acting in $E=C[0,T]$ defined by the formula 
$Av(t)=v'(t)$, with the domain $D(A)=\{v(t):v'(t)\in C[0,T], v(0)=0\}$.
 Then $A$ is a positive operator in the Banach space $E=C[0,T]$ and
\[
    A^{\beta}f(t)=D_t^{\beta}f(t)
\]
for all $f(t)\in D(A)$.
 \end{theorem}

 From the above theorems it follows the following theorem.

\begin{theorem}\label{thm4}
 Let $A$ and $B$ be the positive operators with $\varphi(A,E)<\pi$ and 
$\varphi(B,E)\leq\frac{\pi}{2}$. Then for $\beta\leq\frac{1}{2}$ it follows 
that there exists bounded $(D^{\beta}+B)^{-1}$ defined on whole space $E$. 
Moreover, for every $f$, there exists a unique solution $u=u(f)$ of the problem
\[
    D^{\beta}u+Bu=f
\]
and the following estimate holds
\[
   \| D^{\beta}u \|_{E_{\beta}(E,B)}+\| Bu \|_{E_{\beta}(E,B)}
\leq M(\beta)\| f \|_{E_{\beta}(E,B)}.
\]
 \end{theorem}
 Now, we consider the second order differential operator
\begin{equation}
B^{x}u(x)=-a(x)u_{xx}(x)+\sigma u(x)  \label{3}
\end{equation}
with the domain $D(B^{x})=\{u;u,u',u''\in C[0,l],u(0)=0,u'(l)=0\}$.

Let us introduce the Banach space $C^{\gamma }[0,l]$,
$\gamma \in (0,1] $ of all continuous function $\varphi (x)$ defined on
$[0,l]$ and satisfying a H\"{o}lder condition for which the
following norm is finite
\[
\| \varphi \| _{C^{\gamma }[0,l]}=\| \varphi \|
_{C[0,l]}+\sup_{x_1\neq x_2}\frac{|\varphi (x_1)-\varphi (x_2)|}{
|x_1-x_2|^{\gamma }},
\]
where $C[0,l]$ is the Banach space of all continuous functions
$\varphi (x)$ defined on $[0,l]$ with the norm
\[
\| \varphi \| _{C[0,l]}=\max_{x\in [ 0,l]}|\varphi (x)|.
\]
The positivity of the operator $B^{x}$ in the Banach space $C[0,l]$
was established (see, \cite{45,46}). Moreover, we have that for any 
$\beta \in (0,1/2)$ the norms in the spaces $E_{\beta }(E,B)$ and
$C^{2\beta }[0,l]$ are equivalent.

\begin{theorem}\label{thm5}
 For $\beta \in (0,1/2)$, the norms of the space 
$E_{\beta}(C[0,l],B^{x})$ and the H\"{o}lder space $C^{2\beta}[0,l]$ are equivalent.
 \end{theorem}
 The proof of Theorem \ref{thm5} is based on the following estimates
\[
|G^{x}(x,x_0;\lambda )|\leq \frac{M(\sigma ,a)}{\sqrt{\sigma +\lambda }}
\begin{cases}
e^{-\frac{1}{2}\sqrt{\frac{\sigma +\lambda }{a}}(x-s)}, & 0\leq x_0 \leq x, \\
e^{-\frac{1}{2}\sqrt{\frac{\sigma +\lambda }{a}}(x_0-x)}, & x\leq x_0 \leq l ,
\end{cases}
\]
\[
|G_{x}^{x}(x,x_0;\lambda )|\leq M(\sigma ,a)
\begin{cases}
e^{-\frac{1}{2}\sqrt{\frac{\sigma +\lambda }{a}}(x-x_0)}, & 
0\leq x_0 \leq x, \\
e^{-\frac{1}{2}\sqrt{\frac{\sigma +\lambda }{a}}(x_0-x)}, & 
x\leq x_0 \leq l 
\end{cases}
\]
for the Green's function of the differential operator $B^{x}$
defined by the formula \eqref{3} and it follows the scheme of the
proof of the Theorem of paper \cite{20}.

\begin{theorem}\label{thm6}
 For the solution of problem \eqref{1} the  coercive stability 
estimate
\[
\max_{0\leq t \leq T} \|u_{xx}(t,.)\|_{C^{\beta}[0,l]}\leq M(\beta) \|
f(t,.)\|_{C^{\beta}[0,l]}
\]
holds, where $ M(\beta) $ does not depend on $ f(t,x)$ 
$(0\leq t\leq T$, $x\in [0,l])$ and $0<\beta<1$.
\end{theorem}

The proof of Theorem \ref{thm6} is based on the positivity of
differential operator $B^{x}$ defined by formula \eqref{3}, on the
Theorem \ref{thm3} on connection of fractional derivatives with
fractional powers of positive operators, on the Theorem \ref{thm2}
on spectral angle of fractional powers of positive operators, and on
the Theorem \ref{thm4} on fractional powers of coercively positive
sums two operators.

\section{Difference schemes and stability estimates}

The discretization of problem \eqref{1} is carried out in two steps.
In the first step, let us define the grid space
\[
[ 0,l]_h=(x_{n}=nh,\; 0\leq n\leq M,\; Mh=l)
\]
To the differential space operator $B^{x}$ generated by formula
\eqref{3}, we assign the difference operator $B_h^{x}$ by the
formula
\begin{equation}
B_h^{x}u^{h}=-a(x)u_{{x}_{n}\bar{x}_{n}}^{h}+\sigma u(x)^{h}
\label{3a}
\end{equation}
acting in the space of grid functions $u^{h}(x)$, satisfying the conditions
 $u^{h}(x)=0$ for all $x=0$ and $D^{h}u^{h}(x)=0$ for $x=l$. Here 
$D^{h}u^{h}(x)$ is the approximation of $u_{x}$. With the help of
$B_h^{x}$ we arrive at the initial boundary value problem
\begin{equation}
\begin{gathered}
D_t^{\alpha }v^{h}(t,x)+B_h^{x}v^{h}(t,x)=f^{h}(t,x),\quad
0<t<T,\quad x\in [ 0,l]_h, \\
v^{h}(0,x)=0,\quad x\in [ 0,l]_h
\end{gathered}  \label{1a}
\end{equation}
for a finite system of ordinary fractional differential
equations.

In the second step, applying the first order of approximation
formula (see \cite{12})
\[
D_{\tau }^{\alpha }u_{k}=\frac{1}{\Gamma (1-\alpha )}\sum_{r=1}^{k}
\frac{\Gamma (k-r-\alpha +1)}{(k-r)!}\frac{u_{r}-u_{r-1}}{\tau ^{\alpha }}
, 1\leq k\leq N
\]
for
\[
D_{\tau }^{\alpha }u(t_{k})=\frac{1}{\Gamma (1-\alpha )}\int
_0^{t_{k}}(t_{k}-s)^{-\alpha }u'(s)ds
\]
and using the first order of accuracy stable difference scheme for
parabolic equations, one can present the first order of accuracy
difference scheme with respect to $t$,
\begin{equation}
\begin{gathered}
\frac{1}{\Gamma (1-\alpha )}\sum_{r=1}^{k}
\frac{\Gamma(k-r-\alpha +1)}{(k-r)!}\frac{u_{r}^{h}(x)-u_{r-1}^{h}(x)}
{\tau ^{\alpha }} +B_h^{x}u_{k}^{h}(x)=f_{k}^{h}(x), 
\\
f_{k}^{h}(x)=f^{h}(t_{k},x),\text{ }t_{k}=k\tau , \quad 
1\leq k\leq N,\; N\tau =T,\;  x\in [ 0,l]_h, 
\\
u_0^{h}(x)=0,\quad x\in [ 0,l]_h
\end{gathered}  \label{2a}
\end{equation}
for the approximate solution of problem \eqref{1}. Moreover,
applying the second order of approximation formula:
for $k=1$,
\[ 
D_{\tau }^{\alpha }u_{k}
=-d\frac{2^{\alpha -1}}{(2-\alpha )(1-\alpha )}u_0+d\frac{2^{\alpha -1}}{
(2-\alpha )(1-\alpha )}u_1,
\]
for $k=2$,
\begin{align*}
D_{\tau }^{\alpha }u_{k}
&=d\Bigl[\frac{3^{5-\alpha }}{2^{4-\alpha }}\frac{1}{(1-\alpha)(2-\alpha )(3-\alpha )} 
-7\frac{3^{2-\alpha }}{2^{3-\alpha }}\frac{1}{(1-\alpha)(2-\alpha )}\Bigr]u_0 \\
&\quad+d\Bigl[-\frac{3^{4-\alpha }}{2^{2-\alpha }}\frac{1}{(1-\alpha
)(2-\alpha)(3-\alpha )}  
+\frac{3^{2-\alpha }}{2^{-\alpha }}\frac{1}{(1-\alpha )(2-\alpha )}\Bigr]u_1 \\
&\quad +d\Bigl[\frac{3^{4-\alpha }}{2^{4-\alpha }}\frac{1}{(1-\alpha
)(2-\alpha )(3-\alpha )} 
-\frac{3^{2-\alpha }}{2^{3-\alpha }}\frac{1}{(1-\alpha
)(2-\alpha )}\Bigr]u_2, 
\end{align*}
for $3\leq k\leq N$,
\begin{equation}
\begin{aligned}
D_{\tau }^{\alpha }u_{k}
&= d\sum_{m=2}^{k-1}\Bigr\{\Bigl[\frac{(k-m)}{1-\alpha }\xi (k-m)-
\frac{\eta (k-m)}{2-\alpha }\Bigr]u_{m-2}  \\
&\quad +\Bigl[\frac{(2m-2k-1)}{1-\alpha }\xi (k-m)+\frac{2 \eta (k-m)}{2-\alpha }
\bigr]u_{m-1}  \\
&\quad +\bigl[\frac{(k-m+1)}{1-\alpha }\xi (k-m)-\frac{\eta (k-m)}{2-\alpha }
\bigr]u_{m}\Bigr\}   \\
&\quad +d\Bigl[-\frac{2^{\alpha -2}}{2-\alpha }u_{k-2}
 -\Bigl(\frac{2^{\alpha -1} }{1-\alpha }-\frac{2^{\alpha -1}}{2-\alpha }
 \Bigr)u_{k-1} 
+\Bigl(\frac{2^{\alpha -1}}{1-\alpha }-\frac{2^{\alpha -2}}{2-\alpha }
\Bigr)u_{k}\Bigr]
\end{aligned} \label{fr}
\end{equation}
for
\[
D_t^{\alpha }u(t_{k}-\tau /2)=\frac{1}{\Gamma (1-\alpha )}
\int_0^{t_{k}-\tau /2}(t_{k}-\tau /2-s)^{-\alpha }u'(s)ds,
\]
and using a Crank-Nicholson difference scheme for parabolic equations,
one can present the second order of accuracy difference scheme with
respect to $t$ and $x$,
\begin{equation}
\begin{gathered}
D_{\tau }^{\alpha }u_{k}^{h}(x)+\frac{1}{2}B_h^{x}\Bigl(
u_{k}^{h}(x)+u_{k-1}^{h}(x)\Bigr)=f_{k}^{h}(x),\quad x\in [ 0,l]_h,
\\
f_{k}^{h}(x)=f(t_{k}-\tau /2,x),\quad t_{k}=k\tau ,\; 1\leq k\leq
N,\; N\tau =T, \\
u_0^{h}(x)=0,\quad x\in [ 0,l]_h
\end{gathered}  \label{1aaaa}
\end{equation}
for the approximate solution of problem \eqref{1}. Here,
\begin{gather*}
d=\frac{\tau ^{-\alpha }}{\Gamma (1-\alpha )}\text{,~}\xi (r)=\bigl(r+1/2
\bigr)^{1-\alpha }-\bigl(r-1/2\bigr)^{1-\alpha }, \\
\eta (r)=\bigl(r+1/2
\bigr)^{2-\alpha }-\bigl(r-1/2\bigr)^{2-\alpha }\text{.}
\end{gather*}
Now, we consider the equation
\begin{equation}
B_h^{x}u^{h}+\lambda u^{h}=f^{h}  \label{1aa}
\end{equation}
in the case $a(x)=1$.

\begin{lemma}\label{lem7}
Let $\lambda>0$. Then  \eqref{1aa} is uniquely solvable,
and the  formula
\begin{equation}\label{100}
    u^{h}=(B_h^{x}+\lambda)^{-1}f^{h}=\Bigl\{\sum_{i=1}^{M-1}
G(k,i;\lambda+\sigma)f_ih\Bigr\}_0^{M}
\end{equation}
is valid, where
\[
G(k,i;\lambda+\sigma)=\frac{h(R^{M-i}-R^{M+i})(R^{M-k}-R^{M+k})}
{(1-R^2)(1+R^{2M-1})}+\frac{h(R^{|k-i|+1}-R^{k+i+1})}{(1-R^2)}
\]
for $1\leq i\leq M-1$, and $1\leq k\leq M$,
\[
    R=(1+\delta h)^{-1}, \delta=\frac{1}{2}
\bigl(h(\lambda+\sigma)+\sqrt{(\lambda+\sigma)(4+h^2(\lambda+\sigma))}\bigr).
\]
\end{lemma}

The grid function $G(k,i;\lambda +\sigma )$ is called the Green
function of equation \eqref{1aa} and by the formulas for $R$ and
$\delta $, we get
\begin{equation}
\sum_{i=1}^{M-1}G(k,i;\lambda +\sigma )h=\frac{1}{\lambda +\sigma }-
\frac{1}{\lambda +\sigma }\frac{R^{k}+R^{2M-k-1}}{1+R^{2M-1}},\quad 
1\leq k\leq M. \label{ssss}
\end{equation}

To prove the positivity on $B_h^{x}$ in the Banach space $C_h$,
we need the following auxiliary lemmas  \cite{19}.

\begin{lemma} \label{8}
The following estimate holds
\begin{equation}\label{120}
    |\delta|\geq \max\Bigl\{\frac{|\lambda+\sigma|h}{2},
\sqrt{|\lambda+\sigma|}\Bigr\}.
\end{equation}
\end{lemma}

\begin{lemma}\label{9} The following estimate
\begin{equation}\label{130}
    |R|\leq \frac{1}{1+\sqrt{|\lambda+\sigma|}h\cos\theta}<1
\end{equation}
is valid, where $|\theta|<\pi/2$.
\end{lemma}

\begin{theorem} \label{thm10}
For all $\lambda$ in the sector
$\Sigma_{\theta}=\{\lambda:|\arg\lambda|\leq\theta,
0\leq\theta<\pi/2\}$ the resolvent $(\lambda I+B_h^{x})^{-1}$
defined by \eqref{100} satisfies the following estimate
\begin{equation}\label{140}
    \|(\lambda I+B_h^{x})^{-1}\|_{C_h\to C_h}
\leq \frac{M(\mu,\theta,\sigma)}{1+|\lambda|}.
\end{equation}
\end{theorem}

\begin{proof} 
First, we  consider the operator $B_h^{x}$ defined by formula \eqref{3a} 
in the case $a(x)=1$. Let us set $k=M$. Since
\begin{align*}
    u_{M} &=\frac{h^2R(1-R^{M-1})(1+R^{M-1})}{(1-R)(1+R^{2M-1})}f_{M-1} \\
&\quad +\frac{1}{(1-R)(1+R^{2M-1})}\sum_{i=1}^{M-2}
\Bigl(R^{M-i}-R^{M+i}\Bigr)h^2f_i,
\end{align*}
we have that
\begin{align*}
    \bigl|u_{M} \bigr| 
&\leq  2\Bigl|\frac{R}{1-R}\Bigr|h^2|f_{M-1}|+\frac{1}{\bigl(1-|R|\bigr)}\sum_{i=1}^{M-2}
    \Bigl(|R|^{M-i}+|R|^{M+i}\Bigr)h^2\bigl|f_i\bigr| \\
&\leq 2h^2\|f^{h}\|_{C_h}\Bigl\{|\frac{R}{1-R}|
+\frac{|R^2|}{\bigl(1-|R|\bigr)^2}\Bigr\}.
\end{align*}
Now, let us $1\leq k\leq M-1$. Then by formula
\eqref{100} and the triangle inequality, we obtain
\begin{align*}
|u_k|  &\leq \frac{\bigl(|R|^{M-k}
+|R|^{M+k}\bigr)}{|1-R^2|
\bigl|1+R^{2M-1}\bigr|}\sum_{i=1}^{M-1}\Bigl(|R|^{M-i}
+|R|^{M+i}\Bigr)h^2\bigl|f_i\bigr|  \\
&\quad +\frac{1}{|1-R^2|}\sum_{i=1}^{M-1}\Bigl
(|R|^{|k-i|+1}+|R|^{k+i+1}\Bigr)h^2\bigl|f_i\bigr|  \\
&\leq \frac{2}{|1-R^2|}
\sum_{i=1}^{M-1}\Bigl(|R|^{M-i+1}
+|R|^{M+i+1}\Bigr)h^2\bigl|f_i\bigr|  \\
&\quad +\frac{1}{|1-R^2|}\sum_{i=1}^{M-1}
\Bigl(|R|^{|k-i|+1}+|R|^{k+i+1}
\Bigr) h^2\bigl|f_i\bigr|  \\
&\leq \frac{4h^2}{|1-R^2|}\|f^{h}\|_{C_h}
\sum_{i=1}^{M-1}|R|^{M-i+1}  \\
&\quad +\frac{2h^2}{|1-R^2|}\|f^{h}\|_{C_h}\Bigl
\{\sum_{i=1}^{k-1}|R|^{k-i+1}+|R|
+\sum_{i=k+1}^{M-1}|R|^{i-k+1}\Bigr\}  \\
&\leq \frac{2h^2}{|1-R^2|}\|f^{h}\|_{C_h}\Bigl
\{\frac{2|R|^2}{1-|R|}
+\frac{2|R|^2}{1-|R|}+ |R|\Bigr\}  \\
&\leq M\Bigl\{\frac{|R|^2}{\bigl(1-|R|\bigr)^2}
\frac{h^2}{\bigl|1+R\bigr|}+\Bigl|\frac{R}{1-R}\Bigr|\Bigl|\frac{1}
{1+R}\Bigr|h^2\Bigr\}.
\end{align*}

From estimate \eqref{130} it follows that
\begin{equation}\label{150}
\frac{|R|^2}{\bigl(1-|R|\bigr)^2}\leq\Bigl(\frac{\frac{1}
{1+\sqrt{|\lambda+\sigma|}hcos\theta}}
{1-\frac{1}{1+\sqrt{|\lambda+\sigma|}hcos\theta}}\Bigr)^2
=\Bigl(\frac{1}{\sqrt{|\lambda+\sigma|}hcos\theta}\Bigr)^2.
\end{equation}
Clearly, we have that
\begin{align*}
|\lambda+\sigma|
&=|\rho\cos\theta+i\rho\sin\theta+\sigma|
&=\sqrt{\rho^2+2\rho\sigma\cos\theta+\sigma^2}
\\
&\geq\sqrt{\rho^2\cos^2\theta+2\rho\sigma\cos\theta+\sigma^2}
=|\lambda|\cos\theta+\sigma.
\end{align*}
Thus
\begin{equation} \label{160}
\begin{aligned}
 \frac{1}{|\lambda+\sigma|}
&\leq\frac{1}{|\lambda|\cos\theta+\sigma}\leq\frac{1}
    {|\lambda|\cos\theta+\sigma\cos\theta}\\
&=\frac{\frac{1}{\cos\theta}}{|\lambda|+\sigma}
    =\frac{\frac{1}{\sigma\cos\theta}}{1+\frac{1}{\sigma}|\lambda|}\\
&\leq\frac{M(\sigma,\theta)}{1+|\lambda|}.
\end{aligned}
\end{equation}
Combining estimates  \eqref{150} and \eqref{160}, we obtain that
\begin{equation}\label{170}
    \frac{h^2|R|^2}{\bigl(1-|R|\bigr)^2}
\leq\frac{\frac{1}{\cos^2\theta}}{|\lambda+\sigma|}
\leq\frac{M(\sigma,\theta)}{1+|\lambda|}.
\end{equation}
From the definition of $R$ and estimate \eqref{120}, it follows that
\begin{equation}\label{180}
    \Bigl|\frac{R}{1-R}\Bigr|h^2=\frac{h}{|\delta|}\leq\frac{2}{1+|\lambda|}.
\end{equation}
Combining estimates \eqref{170} and \eqref{180}, we obtain
\[
    \|u^{h}\|_{C_h}\leq\frac{M(\mu,\sigma,\theta)}{1+|\lambda|}\|f^{h}\|_{C_h}.
\]
This concludes the proof of Theorem \ref{thm10} in the case
$a(x)=1$.
Second, noted that the proof of this statement is based on
estimates for the Green's function. Under one more assumption that
$\sigma>0$ is sufficiently large number, applying a fixed point
Theorem, same estimates for the Green's function can be obtained.
Therefore, this statement of theorem is true also for difference
operator $B_h^{x}$ defined by formula \eqref{3a}. Theorem
\ref{thm10} is proved.
\end{proof}

\begin{theorem}\label{thm11}
Let $0<\beta<\frac{1}{2}$. Then, the norms of spaces
$E_{\beta}(C_h,B_h^{x}$) and $C_h^{2\beta}$ are equivalent
uniformly in $h$, $0<h<h_0$.
\end{theorem}

\begin{proof}
From \eqref{100} and \eqref{ssss} it follows that
\begin{align*}
    \bigl(\lambda^{\beta}B_h^{x}(B_h^{x}+\lambda)^{-1}f^{h}\bigr)_{k}
&= \frac{\sigma\lambda^{\beta}}{\lambda+\sigma}f_{k}
    +\frac{\lambda^{\beta+1}}{\lambda+\sigma}
 \frac{R^{k}+R^{2M-k-1}}{1+R^{2M-1}}f_{k}  \\
&\quad +\lambda^{\beta+1}\sum_{i=1}^{M-1}G(k,i;\lambda+\sigma)h(f_{k}-f_{j}).
\end{align*}
Applying the triangle inequality, we obtain
\begin{align*}
&\bigl|\bigl(\lambda^{\beta}B_h^{x}(B_h^{x}+\lambda)^{-1}f^{h}\bigr)_{k} \bigr|\\
&\leq \frac{\sigma\lambda^{\beta}}{\lambda+\sigma}\bigl|f_{k}\bigr|
  +\frac{\lambda^{\beta+1}}{\lambda+\sigma}\bigl|f_{k}\bigr| 
  +\lambda^{\beta+1}\sum_{i=1}^{M-1}|G(k,i;\lambda+\sigma)|h|f_{k}-f_{j}|  \\
&\leq \Bigl[\frac{\sigma\lambda^{\beta}}{\lambda+\sigma}+\frac{\lambda^{\beta+1}}
    {\lambda+\sigma} 
   +M(\sigma)\frac{\lambda^{\beta+1}}{\sqrt{\lambda+\sigma}}
    \sum_{i=1}^{M-1}R^{|k-i|}|(k-i)h|^{2\beta}h\Bigr]\|f^{h}\|_{C_h^{2\beta}}  \\
&\leq  M_1(\sigma)\|f^{h}\|_{C_h^{2\beta}}
\end{align*}
for any $\lambda>0$ and $x\in[0,l]$. Therefore, 
$f^{h}\in E_{\beta}(C_h,B_h^{x})$ and 
\[
    \|f^{h}\|_{E_{\beta}(C_h,B_h^{x})}\leq M_1(\sigma)\|f^{h}\|_{C_h^{2\beta}}.
\]
Now, we  prove the reverse inequality. For any positive operator
$B_h^{x}$, we can write
\[
    v=\int_0^{\infty}\sum_{i=1}^{M-1}G(k,i;\lambda+\sigma)B_h^{x}(B_h^{x}
+\lambda)^{-1}f_ih_1dt.
\]
Consequently,
\[
    f_{k}-f_{k+r}=\int_0^{\infty}\sum_{i=1}^{M-1}\lambda^{-\beta}[G(k+r,i;\lambda
    +\sigma)-G(k,i;\lambda+\sigma)]
\lambda^{\beta}A_h^{x}(A_h^{x}+\lambda)^{-1}f_ih_1dt,
\]
hence
\[
    |f_{k}-f_{k+r}|\leq\int_0^{\infty}\lambda^{-\beta}
\sum_{i=1}^{M-1}|G(k+r,i;\lambda
    +\sigma)-G(k,i;\lambda+\sigma)|h_1dt\|f^{h}\|_{E_{\beta}(C_h,B_h^{x})}.
\]
Let
\[
    T_h=|rh_1|^{-2\beta}\int_0^{\infty}\lambda^{-\beta}
    \sum_{i=1}^{M-1}|G(k+r,i;\lambda+\sigma)-G(k,i;\lambda+\sigma)|h_1dt.
\]
The proof of estimate
\[
    \frac{|f_{k}-f_{k+r}|}{|rh_1|^{2\beta}}
\leq T_h\|f^{h}\|_{E_{\beta}(C_h,B_h^{x})}
\]
 is based on the Lemmas \ref{8} and \ref{9}.  
Thus, for any $1\leq k <k+r\leq N-1$ we have established the inequality
\[
    \frac{|f_{k}-f_{k+r}|}{|rh_1|^{2\beta}}
\leq \frac{M}{\beta(1-2\beta)}\|f^{h}\|_{E_{\beta}(C_h,B_h^{x})}.
\]
This means that 
\[
    \|f^{h}\|_{C_h^{2\beta}}\leq \frac{M}{\beta(1-2\beta)}
\|f^{h}\|_{E_{\beta}(C_h,B_h^{x})}.
\]
Theorem \ref{thm11} in the case $a(x)=1$ is proved. Now, let $a(x)$
be continuous functions and let $x,x_0\in [0,1]$ be arbitrary
fixed points. Clearly, we have that
\[
    \|(B_h^{x}-B_h^{x_0})(B_h^{x_0})\|_{C_h\to C_h}\leq M.
\]
From the formula
\begin{align*}
    B_h^{x}(B_h^{x}+\lambda)^{-1}f^{h}
&=  B_h^{x_0}(B_h^{x_0}+\lambda)^{-1}f^{h}    \\
&\quad +\lambda(\lambda+B_h^{x})^{-1}[B_h^{x}
 -B_h^{x_0}](B_h^{x_0})^{-1}B_h^{x_0}(B_h^{x_0}+\lambda)^{-1}f^{h}
\end{align*}
it follows that
\begin{align*}
&\bigl|\lambda^{\beta}B_h^{x}(B_h^{x}+\lambda)^{-1}f^{h}\bigr|\\
&\leq \|f^{h}\|_{E_{\beta}(C_h,B_h^{x_0})}
 +M\lambda\|(\lambda+B_h^{x})^{-1}\|_{C_h\to C_h}\|f^{h}\|_{E_{\beta}(C_h,B_h^{x_0})}  \\
&\leq M_1\|f^{h}\|_{E_{\beta}(C_h,B_h^{x_0})}.
\end{align*}
Then
\[
    \|f^{h}\|_{E_{\beta}(C_h,B_h^{x})}\leq M_1\|f^{h}\|_{E_{\beta}(C_h,B_h^{x_0})}.
\]
Theorem \ref{thm11} is proved.
\end{proof}

\begin{theorem}[\cite{12}]\label{thm12}
 Let $A_{\tau}$ be the operator acting in $E_{\tau}=C[0,T]_{\tau}$ defined 
by the formula $A_{\tau}v^{\tau}=\{\frac{v_{k}-v_{k-1}}{\tau}\}_1^{N}$, with
$v_0=0$. Then $A_{\tau}$ is a positive operator in the Banach
space $E_{\tau}=C[0,T]_{\tau}$ and
\[
    A_{\tau}^{\beta}f^{\tau}=\Bigl\{\frac{1}{\Gamma(1-\beta)}\sum_{r=1}^{k}
    \frac{\Gamma(k-m-\beta+1)}{(k-m)!}\frac{f_{m}-f_{m-1}}{\tau^\beta}\Bigr\}_1^{N}.
\]
By the definition of fractional difference derivative
\[
    D_{\tau}^{\beta}f^{\tau}:=\Bigl\{\frac{1}{\Gamma(1-\beta)}\sum_{r=1}^{k}
    \frac{\Gamma(k-m-\beta+1)}{(k-m)!}\frac{f_{m}-f_{m-1}}{\tau^\beta}\Bigr\}_1^{N}.
\]
\end{theorem}

\begin{theorem}\label{thm13}
Let $A_{\tau}$ be the operator acting in $E_{\tau}=C[0,T]_{\tau}$
defined by the formula 
$A_{\tau}v^{\tau}(t)=\{\frac{v_{k}-v_{k-1}}{\tau}\}_1^{N}$ with the
domain
\[
D(A_{\tau})=\{v^{\tau}:\frac{v_{k}-v_{k-1}}{\tau}\in C[0,T]_{\tau},
v_0=0\}.
\]
Then $A$ is a positive operator in the Banach space
$E_{\tau}=C[0,T]_{\tau}$, and
\[
     A_{\tau}^{\beta}f^{\tau}(t)= D_{\tau}^{\beta}f^{\tau}(t)
\]
for all $f^{\tau}(t)\in D(A_{\tau})$.
\end{theorem}

Thus, we have the following result on coercive stability of
difference scheme \eqref{1aaaa}.

\begin{theorem}\label{thm14}
Let $ \tau $ and $h$ be sufficiently small positive numbers and
$0<\beta<1$. Then the solution of difference scheme \eqref{1aaaa}
satisfies the following coercive stability estimate:
\[
\max_{1\leq k \leq N} Big\| \Bigl \{
\frac{u_{n+1}^{k}-2u_{n}^{k}+u_{n-1}^{k}}{h^2}\Bigr
\}_{n=1}^{M-1}\|_{C^{\beta}[0,l]_h}\leq
M(\beta)\max_{1\leq k \leq N} \Big\|
f_{k}^{h}\|_{C^{\beta}[0,l]_h}.
\]
Here, $ M(\beta) $ does not depend on $ \tau, h $ and $
f_{k}^{h},1\leq k\leq N. $
\end{theorem}

The proof of Theorem \ref{thm14} is based on the Theorem \ref{thm10}
on positivity of difference space operator $B_h^{x}$ defined by formula (\eqref
{3a}, on the Theorem \ref{thm11} on the structure of fractional
space $E_{\beta }(C_h,B_h^{x_0})$, on the Theorem \ref{thm3}
on connection of fractional derivatives with fractional powers of
positive operators, on the Theorem \ref{thm2} on spectral angle of
fractional powers of positive operators, and on the Theorem
\ref{thm1} on fractional powers of coercively positive sums two
operators.

\section{A numerical application}

For numerical results, we consider the  example
\begin{equation}
\begin{gathered}
D_t^{\alpha }u(t,x)-u_{xx}(t,x)+u(t,x)=f(t,x), \\
f(t,x)=\frac{6\sin ^2(\pi x)t^{3-\alpha }}{\Gamma (4-\alpha
)}-2\pi^2t^{3}\cos (2\pi x) +t^{3}\sin ^2(\pi x),\\
0<t<1,\; 0<x<1, \\
u(0,x)=0, 0\leq x\leq 1, \\
u(t,0)=u_{x}(t,1)=0,\quad  0\leq t\leq 1
\end{gathered}  \label{ea17}
\end{equation}
for the one-dimensional fractional parabolic partial differential
equation with 
$0<\alpha <1$. The exact solution of problem \eqref{ea17} is 
$u(t,x)=t^{3}\sin^2\pi x$. Note that this function is independent of
$\alpha $, but $f(t,x)$ depends on $\alpha$.

 Applying the difference scheme \eqref{2a} for the numerical solution of 
\eqref{ea17}, we obtain
\begin{equation}
\begin{gathered}
\frac{1}{\Gamma (1-\alpha )}\sum_{m=1}^{k}\frac{\Gamma
(k-m-\alpha +1)}{(k-m)!}\frac{u_{n}^{m}-u_{n}^{m-1}}{\tau ^{\alpha}} 
-\frac{u_{n+1}^{k}-2u_{n}^{k}+u_{n-1}^{k}}{h^2}+u_{n}^{k}=\phi_{n}^{k},\\
\phi _{k}^{n}=f(t_{k},x_{n}),\quad t_{k}=k\tau ,\;
1\leq k\leq N,\quad N\tau =T, \\
x_{n}=nh,\quad 1\leq n\leq M-1, \\
u_{n}^{0}=0,\quad 0\leq n\leq M, \\
u_0^{k}=0,\quad u_{M-1}^{k}=u_{M}^{k},\quad 0\leq k\leq N.
\end{gathered}  \label{2b}
\end{equation}
We get the system of equations in the matrix form
\begin{equation}
\begin{gathered}
AU_{n+1}+BU_{n}+CU_{n-1}=D\phi _{n},\quad 1\leq n\leq M-1, \\
U_0=\widetilde{0},\quad U_{M-1}=U_{M},
\end{gathered} \label{2aa}
\end{equation}
where
\begin{gather*}
A= \bordermatrix{&&&   &&\cr &0&0&0&\dots&0&0\cr &0&a_n&0&\dots&0&0\cr
&0&0&a_n&\dots&0&0\cr &\dots&\dots&\dots&\dots&\dots&\dots\cr &0&0&0&\dots&a_n&0\cr
&0&0&0&\dots&0&a_n\cr}_{(N+1)x(N+1)},
\\
B= \bordermatrix{&&&   &&\cr &b_{11}&0&0&\dots&0&0\cr
&b_{21}&b_{22}&0&\dots&0&0\cr &b_{31}&b_{32}&b_{33}&\dots&0&0\cr
&\dots&\dots&\dots&\dots&\dots&\dots\cr &b_{N1}&b_{N2}&b_{N3}&\dots&b_{NN}&0\cr
&b_{N+1,1}&b_{N+1,2}&b_{N+1,3}&\dots&b_{N+1,N}&b_{N+1,N+1}\cr}_{(N+1)x(N+1)},
\\
C= \bordermatrix{&&&   &&\cr &0&0&0&\dots&0&0\cr &0&c_n&0&\dots&0&0\cr
&0&0&c_n&\dots&0&0\cr &\dots&\dots&\dots&\dots&\dots&\dots\cr &0&0&0&\dots&c_n&0\cr
&0&0&0&\dots&0&c_n\cr}_{(N+1)x(N+1)},
\\
D= \bordermatrix{&&&   &&\cr &0&0&0&\dots&0&0\cr &0&1&0&\dots&0&0\cr
&0&0&1&\dots&0&0\cr &\dots&\dots&\dots&\dots&\dots&\dots\cr &0&0&0&\dots&1&0\cr
&0&0&0&\dots&0&1\cr}_{(N+1)x(N+1)},
\\
\phi_{n}= \bordermatrix{&\cr &\phi_{n}^{0}\cr &\phi_{n}^{1}\cr
&\phi_{n}^2\cr &\vdots \cr&\phi_{n}^{N-1}\cr
&\phi_{n}^{N}\cr}_{(N+1)x(1)},\quad
 U_{n}= \bordermatrix{&\cr
&U_{q}^{0}\cr &U_{q}^{1}\cr &U_{q}^2\cr &\vdots \cr
&U_{q}^{N-1}\cr &U_{q}^{N}\cr}_{(N+1)x(1)}, q=\{n\pm 1,n\}, 
\\
a_{n}=-\frac{1}{h^2},~c_{n}=-\frac{1}{h^2},\quad 
b_{11}=1,\quad b_{21}=-\frac{1}{\tau ^{\alpha }},\quad
b_{22}=\frac{1}{\tau ^{\alpha }}+1+\frac{2}{h^2},
\\
b_{31}=-\frac{\Gamma (2-\alpha )}{\Gamma (1-\alpha )\tau ^{\alpha }},\quad
b_{32}=\frac{\Gamma (2-\alpha )}{\Gamma (1-\alpha )\tau ^{\alpha }}-
\frac{1}{\tau ^{\alpha }},\quad 
b_{33}=\frac{1}{\tau ^{\alpha }}+1+\frac{2}{h^2},
\end{gather*}
and
\begin{equation}
b_{ij}=\begin{cases}
-\frac{\Gamma (i-1-\alpha )}{\Gamma (1-\alpha )(i-2)!\tau ^{\alpha }}, &
j=1,
 \\
\frac{1}{\Gamma (1-\alpha )\tau ^{\alpha }}\Bigl[\frac{\Gamma
(i-j+1-\alpha )}{(i-j)!}-\frac{\Gamma (i-j-\alpha
)}{(i-j-1)!}\Bigr], &  2\leq j\leq i-2, \\
\frac{\Gamma (2-\alpha )-\Gamma (1-\alpha )}{\Gamma (1-\alpha )\tau
^{\alpha }}, & \ j=i-1,
  \\
\frac{1}{\tau ^{\alpha }}+1+\frac{2}{h^2}, &  j=i, \\
0, &  i<j\leq N+1
\end{cases}
  \label{22}
\end{equation}
for $i=4,5,\dots,N+1$ and
\[
\phi _{n}^{k}=\frac{6\sin ^2(\pi nh)(k\tau )^{3-\alpha }}{\Gamma
(4-\alpha )}-2\pi ^2(k\tau )^{3}\cos (2\pi nh)+(k\tau )^{3}\sin
^2(\pi nh).
\]

To solve the difference problem \eqref{2aa}, a procedure of modified
Gauss elimination method is applied. Hence, we seek a solution of
the matrix equation in the following form:
\[
U_{j}=\alpha _{j+1}U_{j+1}+\beta _{j+1},\,U_{M}=(I-\alpha
_{M})^{-1}\beta _{M},\,j=M-1,\dots,2,1
\]
where $\alpha _{j}$ $(j=1,2,\dots,M)$ are 
$(N+1)\times (N+1) $ square matrices, and \newline $\beta _{j}$
$(j=1,2,\dots,M)$ are $(N+1)\times 1$ column matrices defined by
\begin{gather*}
\alpha _{j+1}=-(B+C\alpha _{j})^{-1}A, \\
\beta_{j+1}=(B+C\alpha_{j})^{-1}(D\phi -C\beta_{j}),~\,j=1,2,\dots,M-1
\end{gather*}
where $j=1,2,\dots,M-1$, $\alpha _1$ is the $(N+1)\times (N+1)$ zero
matrix, and $\beta _1$ is the $(N+1)\times 1~$zero matrix.

Second, applying the difference scheme \eqref{1aaaa}, we obtain the
second order of accuracy difference scheme in $t$ and in $x$ and the
Crank-Nicholson difference scheme for parabolic equations, one can
represent the second order of accuracy difference scheme with
respect in $t$ and in $x$
\begin{equation}
\begin{gathered}
D_{\tau }^{\alpha }u_{n}^{k}-\frac{1}{2}\Bigl[\frac{
u_{n+1}^{k}-2u_{n}^{k}+u_{n-1}^{k}}{h^2}+\frac{
u_{n+1}^{k-1}-2u_{n}^{k-1}+u_{n-1}^{k-1}}{h^2}\Bigr]
+\frac{1}{2}\Bigr[u_{n}^{k}+u_{n}^{k-1}\Bigr]=\phi _{n}^{k}, 
\\
\phi _{n}^{k}=f(t_{k}-\frac{\tau }{2},x_{n}),\quad t_{k}=k\tau ,x_{n}=nh,
\\
1\leq k\leq N,\quad 1\leq n\leq M-1,
\\
u_{n}^{0}=0,\quad 0\leq n\leq M,
\\
u_0^{k}=0,\quad 3u_{M}^{k}-4u_{M-1}^{k}+u_{M-2}^{k}=0,\quad
 0\leq k\leq N.
\end{gathered}  \label{2aaa}
\end{equation}
Here $D_{\tau }^{\alpha }u_{n}^{k}$ is defined by \eqref{fr} for $u_{n}^{k}$. 
We get the system of equations in the matrix form
\begin{equation}
\begin{gathered}
AU_{n+1}+BU_{n}+CU_{n-1}=D\phi _{n},\quad 1\leq n\leq M-1, \\
U_0=\widetilde{0},\quad 3U_{M}-4U_{M-1}+U_{M-2}=0,
\end{gathered} \label{7}
\end{equation}
where
\begin{gather*}
A=\bordermatrix{&&&   &&\cr &0&0&0&\dots&0&0\cr &a_n&a_n&0&\dots&0&0\cr
&0&a_n&a_n&\dots&0&0\cr &\dots&\dots&\dots&\dots&\dots&\dots\cr
&0&0&0&\dots&a_n&0\cr &0&0&0&\dots&a_n&a_n\cr}_{(N+1)\text{x}(N+1)},
\\
B= \bordermatrix{&&&   &&\cr &b_{11}&0&0&\dots&0&0\cr
&b_{21}&b_{22}&0&\dots&0&0\cr &b_{31}&b_{32}&b_{33}&\dots&0&0\cr
&\dots&\dots&\dots&\dots&\dots&\dots\cr &b_{N1}&b_{N2}&b_{N3}&\dots&b_{NN}&0\cr
&b_{N+1,1}&b_{N+1,2}&b_{N+1,3}&\dots&b_{N+1,N}&b_{N+1,N+1}\cr}_{(N+1)\text{x}
(N+1)},
\\
C= \bordermatrix{&&&   &&\cr &0&0&0&\dots&0&0\cr &c_n&c_n&0&\dots&0&0\cr
&0&c_n&c_n&\dots&0&0\cr &\dots&\dots&\dots&\dots&\dots&\dots\cr
&0&0&0&\dots&c_n&0\cr &0&0&0&\dots&c_n&c_n\cr}_{(N+1)\text{x}(N+1)},
\\
D=\bordermatrix{&&&   &&\cr &0&0&0&\dots&0&0\cr &0&1&0&\dots&0&0\cr
&0&0&1&\dots&0&0\cr &\dots&\dots&\dots&\dots&\dots&\dots\cr &0&0&0&\dots&1&0\cr
&0&0&0&\dots&0&1\cr}_{(N+1)\text{x}(N+1)},
\\
\phi _{n}= \bordermatrix{&\cr &\phi_{n}^{0}\cr &\phi_{n}^{1}\cr
&\phi_{n}^2\cr &\vdots\cr &\phi_{n}^{N-1}\cr
&\phi_{n}^{N}\cr}_{(N+1)x(1)},\quad
U_{q}= \bordermatrix{&\cr &U_{q}^{0}\cr &U_{q}^{1}\cr &U_{q}^2\cr &\vdots\cr
 &U_{q}^{N-1}\cr &U_{q}^{N}\cr}_{(N+1)x(1)},\quad
q=\{n\pm 1,n\},\\
a_{n}=-\frac{1}{2h^2},\quad c_{n}=-\frac{1}{2h^2}, \\
b_{11}=1,\quad b_{21}=-d\frac{2^{\alpha -1}}{(2-\alpha )(1-\alpha )}+\frac{1}{
h^2}+\frac{1}{2}, \quad b_{22}=d\frac{2^{\alpha -1}}{(2-\alpha )(1-\alpha )}+
\frac{1}{h^2}+\frac{1}{2},
\\
b_{31}=d\Bigl[(3/2) ^{5-\alpha }\Big( \frac{1}{
1-\alpha }-\frac{2}{2-\alpha }+\frac{1}{3-\alpha }\Big) -7\frac{
3^{2-\alpha }}{2^{3-\alpha }}\frac{1}{(1-\alpha )(2-\alpha)}\Bigr],
\\
b_{32}=d\Bigl[-\frac{3^{4-\alpha }}{2^{3-\alpha }}\Big( \frac{1}{
1-\alpha }-\frac{2}{2-\alpha }+\frac{1}{3-\alpha }\Big) +\frac{
3^{2-\alpha }}{2^{-\alpha }}\frac{1}{(1-\alpha )(2-\alpha )}\Bigr]+\frac{1
}{h^2}+\frac{1}{2},
\\
b_{33}=d\Bigl[\frac{3^{4-\alpha }}{2^{5-\alpha }}\Big( \frac{1}{
1-\alpha }-\frac{2}{2-\alpha }+\frac{1}{3-\alpha }\Big) -\frac{
3^{2-\alpha }}{2^{3-\alpha }}\frac{1}{(1-\alpha )(2-\alpha )}\Bigr]+\frac{
1}{h^2}+\frac{1}{2},
\\
b_{41}=d\Bigl[\frac{1}{1-\alpha }\xi (1)-\frac{1}{2-\alpha }\eta (1)
\Bigr],
\\
b_{42}=d\Bigl[-\frac{5}{1-\alpha }\xi (1)+\frac{2}{2-\alpha }\eta (1)-
\frac{2^{\alpha -2}}{2-\alpha }\Bigr],
\\
b_{43}=d\Bigl[\frac{2}{1-\alpha }\xi (1)-\frac{1}{2-\alpha }\eta (1)-
\frac{2^{\alpha -1}}{1-\alpha }+\frac{2^{\alpha -1}}{2-\alpha }\Bigr]+
\frac{1}{h^2}+\frac{1}{2},
\\
b_{44}=d\Bigl[\frac{2^{\alpha -1}}{1-\alpha }-\frac{2^{\alpha -2}}{
2-\alpha }\Bigr]+\frac{1}{h^2}+\frac{1}{2},~b_{51}=d\Bigl[\frac{2}{
1-\alpha }\xi (2)-\frac{1}{2-\alpha }\eta (2)\Bigr],
\\
b_{52}=d\Bigl[-\frac{5}{1-\alpha }\xi (2)+\frac{2}{2-\alpha }\eta (2)+
\frac{1}{1-\alpha }\xi (1)-\frac{1}{2-\alpha }\eta (1)\Bigr],
\\
b_{53}=d\Bigl[-\frac{3}{1-\alpha }\xi (1)+\frac{2}{2-\alpha }\eta (1)+
\frac{3}{1-\alpha }\xi (2)-\frac{1}{2-\alpha }\eta (2)-\frac{2^{\alpha -2}
}{2-\alpha }\Bigr],
\\
b_{54}=d\Bigl[\frac{2}{1-\alpha }\xi (1)-\frac{1}{2-\alpha }\eta (1)-
\frac{2^{\alpha -1}}{1-\alpha }+\frac{2^{\alpha -1}}{2-\alpha }\Bigr]+
\frac{1}{h^2}+\frac{1}{2},
\\
b_{55}=d\Bigl[\frac{2^{\alpha -1}}{1-\alpha }-\frac{2^{\alpha -2}}{
2-\alpha }\Bigr]+\frac{1}{h^2}+\frac{1}{2},
\end{gather*}
and
\[
b_{ij}=\begin{cases}
d\Bigl[\frac{1}{1-\alpha }(i-3)\xi (i-3)-\frac{1}{2-\alpha }\eta (i-3)
\Bigr], & j=1, 
\\[4pt]
d\Bigl[\frac{1}{1-\alpha }(5-2i)\xi (i-3)+\frac{2}{2-\alpha }\eta(i-3) 
\\
+\frac{1}{1-\alpha }(i-4)\xi (i-4)-\frac{1}{2-\alpha }\eta (i-4)\Bigr], &
j=2, 
\\[4pt]
d\Bigl[\frac{1}{1-\alpha }(i-j+1)\xi (i-j)-\frac{1}{2-\alpha}\eta (i-j)  \\
+\frac{1}{1-\alpha }(2j-2i+1)\xi (i-j-1)+\frac{2}{2-\alpha }\eta (i-j-1)  \\
+\frac{1}{1-\alpha }(i-j-2)\xi (i-j-2)-\frac{1}{2-\alpha }\eta (i-j-2) \Bigr], 
& 3\leq j\leq i-3, 
\\[4pt]
d\Bigl[\frac{3}{1-\alpha }\xi (2)-\frac{1}{2-\alpha }\eta (2)
-\frac{3}{1-\alpha }\xi (1) \\
+\frac{2}{2-\alpha }\eta (1)-\frac{2^{\alpha -2}}{
2-\alpha }\Bigr], 
& j=i-2, 
\\[4pt]
d\Bigl[\frac{2\xi (1)}{1-\alpha }-\frac{\eta (1)}{2-\alpha }
-\frac{2^{\alpha -1}}{1-\alpha }+\frac{2^{\alpha -1}}{2-\alpha }\Bigr]+\frac{1}{
h^2}+\frac{1}{2}, 
& j=i-1,
\\[4pt]
d\Bigl[\frac{2^{\alpha -1}}{1-\alpha }-\frac{2^{\alpha -2}}{2-\alpha }
\Bigr]+\frac{1}{h^2}+\frac{1}{2}, &  j=i,  
  \\[4pt]
0, & i<j\leq N+1
\end{cases}
\]
for $i=6,7,\dots,N+1$ and
\[
\phi _{n}^{k}=\frac{6\sin ^2(\pi nh)(k\tau )^{3-\alpha }}{\Gamma
(4-\alpha )}-2\pi ^2(k\tau )^{3}\cos (2\pi nh)+(k\tau )^{3}\sin
^2(\pi nh).
\]

 For  solving of the matrix equation \eqref{7}, we use the
same algorithm as in the \eqref{2aa} with
\[
u_{M}=[3I-4\alpha _{M}+\alpha _{M-1}\alpha _{M}]^{-1}[(4I-\alpha
_{M-1})\beta _{M}-\beta _{M-1}].
\]
 Applying the difference schemes \eqref{2a} and \eqref{1aaaa} for
the numerical solution of \eqref{ea17}, we constructed first and
second order of accuracy difference schemes. The results of computer
calculations show that the Crank-Nicholson difference scheme is more
accurate than first order of accuracy difference scheme. Tables \ref{table1}1 and
\ref{table2} are constructed for $N=M=10, 20, 40, 80$, respectively.

\begin{table}[ht]
\caption{Error analysis of first and second order of accuracy
difference schemes for $\alpha =1/2$} \label{table1}
\begin{tabular}{|c|c|c|c|c|}
\hline Method & N=M=10 & N=M=20 & N=M=40 & N=M=80 \\ \hline $1^{st}$
order of accuracy & 1.1110 & 0.7049 & 0.3850 & 0.1998 \\ \hline
$2^{nd}$ order of accuracy & 0.0953 & 0.0111 & 0.0017 & $3.332\times10^{-4}$ \\
\hline
\end{tabular}
\end{table}

\begin{table}[ht]
\caption{Error analysis of first and second order of accuracy
difference schemes for $\alpha =1/3$}
\label{table2}
\begin{tabular}{|c|c|c|c|c|}
\hline Method & N=M=10 & N=M=20 & N=M=40 & N=M=80 \\ \hline $1^{st}$
order of accuracy & 1.1493 & 0.7333 & 0.4015 & 0.2086 \\ \hline
$2^{nd}$ order of accuracy & 0.1015 & 0.0121 & 0.0019 & $7.5456\times10^{-4}$ \\
\hline
\end{tabular}

\end{table}

\subsection*{Conclusion}

In \cite{33} the multidimensional fractional parabolic
equation with the Dirichlet-Neumann conditions was studied.
Stability estimates for the solution of the initial-boundary value
problem for this fractional parabolic equation were given without
proof. The stable difference schemes for this problem were
presented. Stability estimates for the solution of the first order
of accuracy difference scheme were given without proof. The
numerical result was given for the solution of first and second
order of accuracy difference schemes of one-dimensional fractional
parabolic differential equations without any discuss on the
realization.

In the present study, coercive stability estimates for the solution
of this initial-value problem for the fractional parabolic equation
with the Dirichlet-Neumann conditions are established. Stable the
first and second order of approximation in $t$ and first order of
approximation in $x$ difference schemes for this problem are
considered. Coercive stability estimates for the solution of the
first order of accuracy difference scheme are obtained. A procedure
of modified Gauss elimination method is applied for the solution of
the first and second order of accuracy difference schemes of
one-dimensional fractional parabolic differential equations.
Moreover, applying this approach we can constructed the first and
second of approximation in $t$ and a high order of approximation in
$x$ difference schemes. Of course, coercive stability estimates for
the solution of the first order of accuracy difference scheme can be
obtained.

\begin{thebibliography}{99}

\bibitem{4} R. P. Agarwal, B. de Andrade, C. Cuevas; 
\emph{Weighted pseudo-almost periodic solutions of a class of semilinear fractional
differential equations}. Nonlinear Analysis-Real World Applications
\textbf{11} (5) (2010), 3532--3554.

\bibitem{6} D. Amanov, A. Ashyralyev; 
\emph{Initial-boundary value
problem for fractional partial differential equations of higher
order}. Abstract and Applied Analysis (2012), Article ID 973102,
doi: 10.1155/2012/973102, 1--14.

\bibitem{12} A. Ashyralyev;
 \emph{A note on fractional derivatives and
fractional powers of operators}. Journal of Mathematical Analysis
and Applications \textbf{357} (1) (2009), doi:
10.1016/j.jmaa.2009.04.012, 232--236.

\bibitem{28} A. Ashyralyev; 
\emph{Well-posedness of fractional parabolic
equations}. Boundary Value Problems (2013), Article ID 31, doi:
10.1186/1687-2770-2013-31, 1--18.

\bibitem{29} A. Ashyralyev; 
\emph{Well-Posedness of Parabolic
Differential and Difference Equations with the Fractional
Differential Operator}. Malaysian Journal of Mathematical Sciences
\textbf{6(S)} (2012), 73--89.

\bibitem{11} A. Ashyralyev; 
\emph{Well-posedness of the Basset problem in
spaces of smooth functions}. Applied Mathematics Letters \textbf{24}
(7) (2011), doi: 10.1016/j.aml.2011.02.002, 1176--1180.

\bibitem{22} A. Ashyralyev, B. Hicdurmaz; 
\emph{On the numerical
solution of fractional Schrodinger differential equations with the
Dirichlet condition}. International Journal of Computer Mathematics
\textbf{89} (13-14) (2012), doi: 10.1080/00207160.2012.698841,
1927--1936.

\bibitem{7} A. Ashyralyev, B. Hicdurmaz; 
\emph{A note on the fractional
Schrodinger differential equations}. Kybernetes \textbf{40} (5-6)
(2011), doi: 10.1108/03684921111142287, 736--750.

\bibitem{20} A. Ashyralyev, D. Agirseven; 
\emph{Well-posedness of delay
parabolic difference equations}. Advances in Difference Equations
\textbf{2014} (18) (2014), doi: 10.1186/1687-1847-2014-18, 1--20.

\bibitem{34} A. Ashyralyev, F. Dal; 
\emph{Finite difference and iteration
methods for fractional hyperbolic partial differential equations
with the neumann condition}. Discrete Dynamics in Nature and Society
(2012), Article ID 434976, doi: 10.1155/2012/434976.

\bibitem{9} A. Ashyralyev, F. Dal, Z. Pinar; 
\emph{A note on the
fractional hyperbolic differential and difference equations}.
Applied Mathematics and Computation \textbf{217} (9) (2011), doi:
10.1016 /j.amc.2010.11.017, 4654--4664.

\bibitem{33} A. Ashyralyev, N. Emirov, Z. Cakir; 
\emph{Fractional
parabolic differential and difference equations with the
dirichlet-neumann conditions}. First International Conference on
Analysis and Applied Mathematics \textbf{1470} (2012), doi:
10.1063/1.4747641, 69--72.

\bibitem{19} A. Ashyralyev, N. Nalbant, Y. Sozen; 
\emph{Structure of
fractional spaces generated by second order difference operators}.
Journal of the Franklin Institute-Engineering and Applied
Mathematics \textbf{351} (2) (2014), doi:
10.1016/j.jfranklin.2013.07.009, 713--731.

\bibitem{23} A. Ashyralyev, Y. A. Sharifov; 
\emph{Existence and uniqueness
of solutions for the system of nonlinear fractional differential
equations with nonlocal and integral boundary conditions}. Abstract
and Applied Analysis (2012), doi: 10.1155/2012/594802, Article ID
594802.

\bibitem{25} A. Ashyralyev, Z. Cakir; 
\emph{On the numerical solution
of fractional parabolic partial differential equations with the
dirichlet condition}. Discrete Dynamics in Nature and Society
(2012), doi: 10.1155/2012/696179, Article ID 696179.

\bibitem{5} A. S. Berdyshev, A. Cabada, E. T. Karimov; 
\emph{On a non-local boundary problem for a parabolic-hyperbolic equation
involving a Riemann-Liouville fractional differential operator}.
Nonlinear Analysis--Theory Methods and Applications \textbf{75} (6)
(2012), doi: 10.1016/j.na.2011.12.033, 3268--3273.

\bibitem{18} S. Q. Bu, P. Clement, S. Guerre-Delabriere; 
\emph{Regularity of pairs of positive operators}. 
Illinois Journal of Mathematics \textbf{42} (3) (1998), 357--370.

\bibitem{21} Z. Cakir; 
\emph{Stability of difference schemes for
fractional parabolic pde with the dirichlet-neumann conditions}.
Abstract and Applied Analysis (2012), doi: 10.1155/2012/463746,
Article ID 463746, 1--17.

\bibitem{27} J. Y. Cao, C. J. Xu; 
\emph{A high order schema for the
numerical solution of the fractional ordinary differential
equations}. Journal of Computational Physics \textbf{238} (2013),
doi: 10.1016/j.jcp.2012.12.013, 154--168.

\bibitem{16} Ph. Clement; 
\emph{On ($L_{p}$-$L_{q}$) Coerciveness for a
Class of Integrodifferential Equation on the Line}. Two lectures at
the 23rd annual Voronezh Winter School of Mathematics (USSR), Tech.
Report 5-4-90, Voronezh, (1990).

\bibitem{44} E. Cuesta, C. Lubich, C. Palencia; 
\emph{Convolution
quadrature time discretization of fractional diffusion-wave
equations}. Mathematics of Computation \textbf{75} (254) (2006),
doi: 10.1090/S0025-5718-06-01788-1, 673--696.

\bibitem{14} G. Da Prato, P. Grisvard; 
\emph{Sums of linear-operators and
operational differential equations}. Journal De Mathematiques Pures
Et Appliquees \textbf{54} (1975), 305--387.

\bibitem{41} F. Dal; 
\emph{Application of variational iteration method to
fractional hyperbolic partial differential equations}. Mathematical
Problems in Engineering (2009), doi: 10.1155/2009/824385, Article ID
824385.

\bibitem{37} M. De la Sen; 
\emph{About robust stability of caputo linear
fractional dynamic systems with time delays through fixed point
theory}. Fixed Point Theory and Applications (2011), doi:
10.1155/2011/867932, Article ID 867932.

\bibitem{36} M. De la Sen; 
\emph{Positivity and stability of the
solutions of caputo fractional linear time-invariant systems of any
order with internal point delays}. Abstract and Applied Analysis
(2011), doi: 10.1155/2011/161246, Article ID 161246.

\bibitem{3} M. De la Sen, R. P. Agarwal, A. Ibeas, et al; 
\emph{On the
existence of equilibrium points, boundedness, oscillating behavior
and positivity of a sveirs epidemic model under constant and
impulsive vaccination}. Advances in Difference Equations (2011),
doi: 10.1155/2011/748608, Article ID 748608.

\bibitem{24} A. S. Erdogan, H. Uygun; 
\emph{A note on the inverse
problem for a fractional parabolic equation}. Abstract and Applied
Analysis (2012), doi: 10.1155/2012/276080, Article ID 276080.

\bibitem{31} J. Fan, J. H. He; 
\emph{Fractal derivative model for air
permeability in hierarchic porous media}. Abstract and Applied
Analysis (2012), doi: 10.1155/2012/354701, Article ID 354701.

\bibitem{2} A. A. Kilbas, H. M. Srivastava, J. J. Trujillo; 
\emph{Theory and Applications of Fractional Differential Equations}.
Elsevier, Amsterdam (2006).

\bibitem{8} M. Kirane, Y. Laskri;
 \emph{Nonexistence of global solutions
to a hyperbolic equation with a space-time fractional damping}.
Applied Mathematics and Computation \textbf{167} (2) (2005), doi:
10.1016/j.amc.2004.08.038, 1304--1310.

\bibitem{13} S. G. Krein; 
\emph{Linear Differential Equations in a Banach Space}. Nauka, Moscow, (1966).

\bibitem{10} V. Lakshmikantham, A. Vatsala, 
\emph{Basic theory of fractional differential equations}. Nonlinear Analysis-Theory
Methods and Applications \textbf{69} (8) (2008), doi:
10.1016/j.na.2007.08.042, 2677--2682.

\bibitem{42} C. G. Li, M. Kostic, M. Li, S. Piskarev; 
\emph{On a class of time-fractional differential equations}. 
Fractional Calculus and Applied Analysis \textbf{15} (4) (2012), doi:
10.2478/s13540-012-0044-x, 639--668.

\bibitem{1} I. Podlubny; 
\emph{Fractional Differential Equations}.
Mathematics in Science and Engineering \textbf{198}, Academic Press,
San Diego, California, (1999).

\bibitem{15} P. E. Sobolevskii; 
\emph{Fractional powers of coercively positive sums of operators}. 
Doklady Akademii Nauk (SSSR) \textbf{16} (1975), 1638--1641.

\bibitem{17} P. E. Sobolevskii; 
\emph{Fractional powers of coercive-positive sums of operators}. 
Siberian Mathematical Journal \textbf{18} (3) (1977), 637--657.

\bibitem{45} M. Z. Solomyak; 
\emph{Estimation of norm of resolvent of elliptic operator in spaces $L_{p}$}. 
Usp. Mat. Nauk \textbf{15} (6) (1960), 141--148.

\bibitem{46} H. B. Stewart; 
\emph{Generation of analytic semigroups by strongly elliptic operators}. 
Tran. Amer. Math. Soc. \textbf{199} (1974), doi: 10.2307/1996879, 141--162.

\bibitem{43} L. J. Su, P. Cheng; 
\emph{A high-accuracy moc/fd method for solving fractional advection-diffusion 
equations}. Journal of Applied Mathematics (2013), doi: 10.1155/2013/648595, 
Article ID 648595.

\bibitem{26} A. Yakar, M. E. Koksal; 
\emph{Existence results for solutions of nonlinear fractional differential equations}. Abstract and
Applied Analysis (2012), doi: 10.1155/2012/267108, Article ID
267108.

\bibitem{32} H. Yang; 
\emph{Existence of mild solutions for a class of
fractional evolution equations with compact analytic semigroup}.
Abstract and Applied Analysis (2012), doi: 10.1155/2012/903518,
Article ID 903518.

\bibitem{39} C. J. Yuan; 
\emph{Multiple positive solutions for $(n-1, 1)$-type semipositone conjugate 
boundary value problems for coupled
systems of nonlinear fractional differential equations}. Electronic
Journal of Qualitative Theory of Differential Equations \textbf{13}
(2011), 1--12.

\bibitem{40} C. J. Yuan; 
\emph{Multiple positive solutions for $(n-1, 1)$-type semipositone conjugate 
boundary value problems of nonlinear fractional differential equations}.
 Electronic Journal of Qualitative Theory of Differential Equations 
\textbf{36} (2010), 1--12.

\bibitem{38} C. J. Yuan; 
\emph{Multiple positive solutions for semipositone $(n, p)$-type boundary value 
problems of nonlinear fractional differential equations}. 
Analysis and Applications \textbf{9} (1) (2011),
 doi: 10.1142/S0219530511001753, 97--112.

\bibitem{30} C. J. Yuan; 
\emph{Two positive solutions for $(n-1,1)$-type
semipositone integral boundary value problems for coupled systems of
nonlinear fractional differential equations}. Communications in
Nonlinear Science and Numerical Simulation \textbf{17} (2) (2012),
doi: 10.1016/j.cnsns.2011.06.008, 930--942.

\bibitem{35} C. J. Yuan, D. Jiang, D. O'Regan, R. P. Agarwal; 
\emph{Multiple positive solutions to systems of nonlinear 
semipositone fractional differential equations with coupled boundary conditions}. 
Electronic Journal of Qualitative Theory of Differential Equations \textbf{13}
(2012), 1--17.

\end{thebibliography}

\end{document}