Higher-order approximations
for space-fractional diffusion equation

Abstract.

Second-order and third-order finite difference approximations for fractional derivatives are derived from a recently proposed unified explicit form. The Crank-Nicholson schemes based on these approximations are applied to discretize the space-fractional diffusion equation. We theoretically analyze the convergence and stability of the Crank-Nicholson schemes, proving that they are unconditionally stable. These schemes exhibit unconditional stability and convergence for fractional derivatives of order $\alpha$ in the range $\frac{4}{3}\le \alpha \le 2$. Numerical examples further confirm the convergence orders and unconditional stability of the approximations, demonstrating their effectiveness in practice.

1. Introduction

Fractional derivatives, for example, the Riemann-Liouville, Caputo, and Grünwald-Letnikov derivatives, have found numerous applications across various applied fields, including physics [10, 13, 27, 28], biology[7, 29, 30], finance[25, 26], and engineering [11, 12, 15, 16, 17, 18]. Their unique characteristics such as non-locality give more suitable descriptions for various phenomena, including anomalous diffusion, population dynamics, fractional Brownian motion, etc. compared to traditional derivatives. However, the non-local nature of fractional derivatives often leads to complex formulas, making it difficult to solve fractional-order differential equations, such as fractional-order diffusion equations, using both analytical and approximate methods [6, 8, 24].

The Grünwald difference (GD) approximation presents a finite difference technique for approximating fractional derivatives. It utilizes an infinite sum of terms derived from the power series expansion of the generator $W_1(z)={(1-z)}^{\alpha }$ about zero [6]. Despite its first-order accuracy, the GD approximation yields unstable results, even with stable methods like Euler and Crank-Nicholson, which are typically reliable for integer-order diffusion equations. To address this, a shifted version with shift $r=1$ was introduced [8], aiming to restore unconditional stability while preserving first-order accuracy for fractional derivatives.

The shifted Grünwald approximation has served as a cornerstone for constructing higher-order finite difference approximations for fractional order differential equations. Meerchaert et al. [20] employed extrapolation technique on the Crank-Nicholson scheme of the Shifted Grünwald approximation for the fractional diffusion equation to obtain the second-order accuracy for the space discretization. Nasir et al. [19], from the shifted Grünwal approximations with a non-integer shift $r=\frac{1}{2}$, obtained a second-order accuracy which displays super convergence. Convex combinations of the shifted form of the Grünwal approximation with various shifts were used to obtain some second-order approximations [21]. This technique is referred in the literature as the weighted and shifted Grünwald difference(WSGD) approximations. A third-order approximation through WSGD was not successful as it fails to give the desired stability for a fractional derivative order $\alpha$ in the range $1\le \alpha \le 2$. However, achieving stability for the third-order approximation is possible within a restricted range of the fractional order as shown in [23]. Hao [22] derived a fourth-order approximation using a quasi compact difference approximation technique on a WSGD approximation. Additionally, by utilizing super-convergent approximations for fractional derivatives, Zhao and Deng [34] proposed a series of higher-order difference schemes for the space fractional diffusion equation.

Moreover, Lubich [5] proposed generators in the form of power or rational polynomials to construct higher-order approximations for fractional derivatives. While these generators provide coefficients for higher-order accuracy without shifts, their shifted forms only yield first-order approximations regardless of the original accuracy orders.

Nasir and Nafa [4] introduced polynomial-type generators for higher-order approximations with shifts and derived a second-order finite-difference scheme for the one-dimensional fractional diffusion equation. In construction of this work, Nasir and Nafa [3] and Gunarathna et al. [14] developed quasi-compact schemes with third-and fourth-order accuracy, respectively, both derived from the second-order approximation and applied them to the one-dimensional fractional diffusion equation.

The generators for the Nasir and Nafa [4] approximations are usually obtained manually by hand calculations, solving a resulting system of linear equations or by symbolic computations and these processes are specific to the problem at hand. To alleviate those difficulties, Gunarathna et al. [9] have obtained an explicit form for generators that gives approximations for fractional derivatives with shifts retaining their higher orders. This form generalizes the Lubich form with shift and hence the Lubich form becomes a special case with no shift. Gunarathna et al. [31] then extended the explicit form developed in [9] to a more general unified explicit form that gives more new approximations for fractional derivatives and various finite difference formulas for any classical derivative.

In this paper, we apply the unified explicit form in [9] to the space fractional diffusion equation given by (1). We consider a second-order approximation derived from this unified form. Subsequently, a new quasi-compact third-order approximation is derived from this second-order approximation. Using these approximations, second-and third -order Crank-Nicholson (C-N) schemes are constructed for the space fractional diffusion equation. Theoretical analyses of stability and convergence are established for both the C-N schemes.

(1) $$\frac{\partial u(x,t)}{\partial t}=K_1{D}_{x-}^{\alpha }u(x,t)+K_2{D}_{x+}^{\alpha }u(x,t)+f(x,t),$$

with the initial and boundary conditions:

(2) $u(x,0)=s_0(x),x\in [a,b];u(a,t)=Z_1(t),u(b,t)=Z_2(t),t\in [0,T],$

where $u(x.t)$ is the unknown function to be determined; $K_1,K_2$ are non-negative constant diffusion coefficient with $K_1+K_2\ne 0$, i.e., not both are simultaneously zero; $f(x,t)$ is a known source term. The fractional derivatives $D_{x-}^{\alpha }$ and $D_{x+}^{\alpha }$, in Riemann-Liouville sense, are given in Definition 1.

The remaining sections of this paper are structured as follows: Section 2 presents essential preliminaries and terminologies. Section 3 applies the unified form to obtain new second-and third-order approximations for fractional derivatives and obtain their Crank-Nicholson (C-N) schemes with order $2$ and order $3$ to solve the space fractional diffusion equation. Section 4 analyses the C-N schemes derived in Section 3. Section 5 presents numerical examples. Section 6 concludes the paper.

2. Preliminaries and Terminologies

This section presents the requisite materials and definitions relevant to the subject of the paper. Let $f(x)$ be a sufficiently smooth function defined on a real domain $ℝ$.

Definition 1 ([6]).

The left(-) and right(+) Riemann-Liouville (R-L) fractional derivatives of a real order $\alpha >0$ are defined as

(3) $${}^{RL}{D}_{x-}^{\alpha }f(x)=\frac{1}{\mathrm{\Gamma }(n-\alpha )}\frac{d^n}{d{x}^n}{\int }_{-\mathrm{\infty }}^x\frac{f(\eta )}{{(x-\eta )}^{\alpha +1-n}}d\eta ,$$

and

(4) $${}^{RL}{D}_{x+}^{\alpha }f(x)=\frac{{(-1)}^n}{\mathrm{\Gamma }(n-\alpha )}\frac{d^n}{d{x}^n}{\int }_x^{\mathrm{\infty }}\frac{f(\eta )}{{(\eta -x)}^{\alpha +1-n}}d\eta$$

respectively, where $n=[\alpha ]+1$, an integer with and $\mathrm{\Gamma }(\cdot )$ denotes the gamma function.

Definition 2 ([4]).

Let $\{w_k^{(\alpha )}\}$ be a sequence of real numbers with generating function

$$W(z)=\sum _{k=0}^{\mathrm{\infty }}{w}_k^{(\alpha )}{z}^k.$$

Define a shifted difference formula with shift $r$ as

(5) $${\mathrm{\Delta }}_{\pm h,p,r}^{\alpha }f(x)=\frac{1}{h^{\alpha }}\sum _{k=0}^{\mathrm{\infty }}{w}_k^{(\alpha )}f(x\mp (k-r)h).$$

1. $W(z)$ is said to approximate the fractional derivatives $D_{x\mp }^{\alpha }$ if

   (6) $$\underset{h\to 0}{lim}{\mathrm{\Delta }}_{\pm h,p,r}^{\alpha }f(x)=D_{x\mp }^{\alpha }f(x).$$

2. $W(z)$ is said to approximate the fractional derivatives $D_{x\mp }^{\alpha }$ with order $p$ if

   (7) $${\mathrm{\Delta }}_{\pm h,p,r}^{\alpha }f(x)=D_{x\mp }^{\alpha }f(x)+𝒪(h^p).$$

Proposition 3 (Theorem 1 [4, 3]).

Let $\alpha >0,n=[\alpha ]+1,$ and a non-negative integer $m$ be given. Let a function $f(x)\in C^{m+n+1}(ℝ)$ and $D^{k}f(x)=\frac{d^k}{d{x}^k}f(x)\in L_1(ℝ)$ for $0\le k\le m+n+1$. Then, a generator $W(z)$ approximates the fractional derivatives $D_{x\pm }^{\alpha }f(x)$ with order $p$ and shift $r$, $1\le p\le m$, if and only if

(8) $$G(z)=\frac{1}{z^{\alpha }}W(e^{-z})e^{rz}=1+𝒪(z^p).$$

Moreover, if $G(z)=1+{\sum }_{l=p}^{\mathrm{\infty }}{a}_l{z}^l$, where $a_l\equiv a_l(\alpha ,r)$, then we have

${\mathrm{\Delta }}_{\pm h,p,r}^{\alpha }f(x)$ $=D_{x\pm }^{\alpha }f(x)+h^p{a}_p{D}_{x\pm }^{\alpha +p}f(x)+h^{p+1}{a}_{p+1}{D}_{x\pm }^{\alpha +p+1}f(x)+\mathrm{\dots }$

(9) $+h^m{a}_m{D}_{x\pm }^{\alpha +m}f(x)+𝒪(h^{m+1}).$

2.1. The unified explicit form

In this section, the unified form appearing in [31] is presented. This unified form extends the explicit form in [9] to a more general form that covers compact finite difference formulas for higher order classical derivatives as well as some new Lubich type generators for fractional derivatives. For this, we introduce a base differential order $d$, a positive integer, to express the fractional differential operator as

$$D_{x\pm }^{\alpha }={\left(D^d\right)}_{x\pm }^{\frac{\alpha }{d}},$$

and consider approximating the fractional derivative by a Lubich type generator of the form

(10) $$W(z)={\left({\beta }_0+{\beta }_{1}z+\mathrm{\dots }+{\beta }_{N-1}{z}^{N-1}\right)}^{\frac{\alpha }{d}}={\left(P(z)\right)}^{\frac{\alpha }{d}},$$

where $P(z)$ corresponds to the classical derivative operator $D^d$. The coefficients ${\beta }_j$ in (10) are to be determined based on the fractional order $\alpha$, the required approximation order $p$, and shift $r$. The degree $N-1$ of $P(z)$ is similarly determined based on $p$ and $d$. This setup leads to the following theorem.

Theorem 4.

With assumptions of Proposition 3, the generator of the form $W(z)=W_{p,r,d}(z)={\left({\beta }_0+{\beta }_{1}z+\mathrm{\dots }+{\beta }_{N-1}{z}^{N-1}\right)}^{\frac{\alpha }{d}}$, where $d$ is a positive integer, approximates the fractional derivatives $D_{x\mp }^{\alpha }f(x)$ at $x$ with order $p$ and shift $r$ if and only if the coefficients ${\beta }_j$ satisfy the linear system

(11) $$\sum _{j=0}^{N-1}{(\lambda -j)}^k{\beta }_j=d!{\delta }_{d,k},k=0,1,\mathrm{\dots },N-1,$$

where $\lambda =rd/\alpha ,N=p+d$ and ${\delta }_{d,k}$ is the Kronecker delta having value of one for $k=d$ and zero otherwise.

Proof.

In view of Proposition 3, we have $G(z)=\frac{1}{z^{\alpha }}W\left(e^{-z}\right)e^{rz}=1+𝒪(z^p)$. This gives

$G(z)$ $=\frac{1}{z^{\alpha }}{\left({\displaystyle \sum _{j=0}^{N-1}}{\beta }_j{e}^{-jz}\right)}^{\frac{\alpha }{d}}{e}^{rz}=\frac{1}{z^{\alpha }}{\left({\displaystyle \sum _{j=0}^{N-1}}{\beta }_j{e}^{(rd/\alpha -j)z}\right)}^{\frac{\alpha }{d}}$

$={\left(\frac{1}{z^d}{\displaystyle \sum _{j=0}^{N-1}}{\beta }_j{e}^{{\lambda }_{j}z}\right)}^{\frac{\alpha }{d}}={\left(\frac{1}{z^d}{\displaystyle \sum _{j=0}^{N-1}}{\beta }_j{\displaystyle \sum _{k=0}^{\mathrm{\infty }}}{\displaystyle \frac{1}{k!}}{\lambda }_j^k{z}^k\right)}^{\frac{\alpha }{d}}={\left(\frac{1}{z^d}{\displaystyle \sum _{k=0}^{\mathrm{\infty }}}{b}_k{z}^k\right)}^{\frac{\alpha }{d}}$

$={\left(\frac{b_0}{z^d}+\frac{b_1}{z^{d-1}}+\mathrm{\dots }+\frac{b_{d-1}}{z}+b_d+{\displaystyle \sum _{k=d+1}^{\mathrm{\infty }}}{b}_k{z}^{k-d}\right)}^{\frac{\alpha }{d}}$

$=1+𝒪(z^p),$

where ${\lambda }_j=\lambda -j,\lambda =\frac{rd}{\alpha }$ and

(12) $$b_k=\frac{1}{k!}\sum _{j=0}^{N-1}{\lambda }_j^k{\beta }_j,k=0,1,2,\mathrm{\dots }.$$

Since $G(z)$ does not have any pole singularities by virtue of (8), we have $b_k=0$ for $k=0,1,\mathrm{\dots },d-1$. Moreover, since $G(0)=1$, we have $b_d=1$. These are the consistency condition for the GTA with generator $W(z)$. Now, for order $p=1$, these conditions give the system (11) with $N=1+d$ and the proof ends. For $p>1$, $G(z)$ reduces to

$$G(z)={(1+\sum _{k=d+1}^{\mathrm{\infty }}{b}_k{z}^{k-d})}^{\frac{\alpha }{d}}=:{(1+X)}^{\gamma }=1+𝒪(z^p),$$

where $\gamma =\frac{\alpha }{d}$ and

(13) $$X=\sum _{k=d+1}^{\mathrm{\infty }}{b}_k{z}^{k-d}.$$

Expansion of ${(1+X)}^{\gamma }$ gives

(14) $$1+\gamma X+\frac{\gamma (\gamma -1)}{2!}{X}^2+\mathrm{\dots }=1+𝒪(z^p).$$

The term with $z$ appears in the term $\gamma X$ only on the left-hand side of (14). This gives $b_{d+1}=0$. The same is true for $b_k,k=d+1,d+2,\mathrm{\dots },p+d-1$, by successively comparing the coefficients of $z^{k-d}$ to gain $𝒪(z^p)$ in (14). Altogether, we have $b_k={\delta }_{d,k}$, $k=0,1,2,\mathrm{\dots },p+d-1$ which yield the linear system (11) with (12) and $N=p+d$. ∎

Theorem 5.

Let $\alpha >0$, a positive integer $d\ge 1$ and $f(x)$ be a sufficiently smooth function such that $D_{x\pm }^{\alpha }f(x)$ is defined. For an approximation (5) for $D_{x\pm }^{\alpha }f(x)$ of order $p$ and shift $r$ with the generator in the form (10), the coefficients ${\beta }_j$ are given by

(15) $${\beta }_j=\frac{N_j}{D_j},j=0,1,\mathrm{\dots },N-1,$$

where $N=p+d$ and

(16)

For the proof of Theorem 5, the interesting readers are referred to [31].

3. Applications of the unified explicit form

This section applies the unified explicit form to derive a second and third order approximations.

3.1. Second-order approximation

To derive the second order approximation, the following generating function form:

(17) $$W_{2,2}(z)={({\beta }_0+{\beta }_{1}z+{\beta }_2{z}^2+{\beta }_3{z}^3)}^{\alpha /2}$$

is considered. The coefficients ${\beta }_0,{\beta }_1,{\beta }_2$, and ${\beta }_3$ are computed using the unified form equation (15). The computed coefficients are: ${\beta }_0=-\lambda +2,{\beta }_1=3\lambda -5,{\beta }_2=-3\lambda +4,$ and ${\beta }_3=\lambda -1$, where $\lambda =\frac{2r}{\alpha }$.

3.2. Third-order approximation: Quasi-compact form

Now, we derive a new quasi-compact third-order approximation from the second-order approximation described in the Section 3.1. In view of Proposition 3, we have

$$H(z)=H_{r,2}(z)=\frac{1}{z^{\alpha }}{W}_{r,2}(e^{-z})e^{rz}=1+a_2(r)z^2+a_3(r)z^3+a_4(r)z^4+\mathrm{\dots },$$

where $a_2(r)=-\frac{1}{24\alpha }\left(11{\alpha }^2-36\alpha r+24r^2\right)$ and

$a_3(r)=\frac{1}{6{\alpha }^2}\left(3{\alpha }^3-13{\alpha }^{2}r+18\alpha r^2-8r^3\right)$. Also, Equation (3) gives

${\mathrm{\Delta }}_{\mp h,2,r}^{\alpha }f(x)$ $=D_{x\pm }^{\alpha }f(x)+a_2{h}^2{D}_{x\pm }^{2+\alpha }f(x)+𝒪(h^3)$

$=\left(1+a_2{h}^2{D}^2\right)D_{x\pm }^{\alpha }f(x)+𝒪(h^3)$

(18) $=P_x{D}_{x\pm }^{\alpha }f(x)+𝒪(h^3),$

where $P_x=\left(I+h^2{a}_2{D}_x^{\alpha }\right)$ and $I$ is the identity differential operator. The differential operator $D_x^2$ may be approximated by the standard central difference operator ${\delta }_h^{2}f(x)$ with $D_x^{2}f(x)={\delta }_h^{2}f(x)+𝒪(h^2)$. Using this, an approximation is obtained for $P_x$ such that

$P_{x}f(x)$ $=(I+a_2{h}^2{D}_x^2)f(x)$

$=(I+a_2{h}^2({\delta }_h^2+𝒪(h^2)))f(x)$

$=(I+a_2{h}^2{\delta }_h^2)f(x)+𝒪(h^4)$

(19) $=P_{h}f(x)+𝒪(h^4),$

where

(20) $$P_h=(I+a_2{h}^2{\delta }_h^2).$$

Then (3.2) and (3.2) reason to give

${\mathrm{\Delta }}_{\pm h,2,r}f(x)$ $=(P_h+𝒪(h^4))D_{x\pm }^{\alpha }f(x)+𝒪(h^3)$

(21) $=P_h{D}_{x\pm }^{\alpha }f(x)+𝒪(h^3).$

3.3. Discretization of the space fractional diffusion equation

The discretization of the space fractional diffusion equation (1) in the domain $[a,b]\times [0,T]$ is considered. The function $u(x,t)$ is zero-extended outside the space domain interval $[a,b]$ so that the left and right fractional derivatives are applicable. For a numerical scheme, the space domain $[a,b]$ is partitioned into a uniform mesh of size $N$ with sub-interval length $h=(b-a)/N$, and the time domain $[0,T]$ into a uniform partition of size $M$ with sub-interval length $\tau =T/M$. These form a uniform partition on the 2-D domain $[a,b]\times [0,T]$ with grid points $(x_i,t^m),0\le i\le N,0\le m\le M$, where $x_i=a+ih$ and $t^m=m\tau$. The following notations are also introduced for conciseness:

$$u_i^m=u(x_i,t^m),t^{m+1/2}=\frac{1}{2}\left(t^{m+1}+t^m\right),f_i^{m+1/2}=f(x_i,t^{m+1/2}),$$

$$U^m={(u_0^m,u_1^m,\mathrm{\dots },u_N^m)}^T,\text{and}{F}^{m+1/2}={(f_0^{m+1/2},f_1^{m+1/2},\mathrm{\dots },f_N^{m+1/2})}^T.$$

Furthermore, prior to the construction of the C-N schemes, it should be noted that the time derivative at $(x,t+\tau /2)$ may be approximated with order $2$ accuracy as follows:

(22) $\frac{\partial u(x,t+\tau /2)}{\partial t}$ $=\frac{1}{\tau }(u(x,t+\tau )-u(x,t))+𝒪({\tau }^2),$

and

(23) $u(x,t+\tau /2)$ $=\frac{1}{2}(u(x,t+\tau )+u(x,t))+𝒪({\tau }^2).$

3.4. Second-order Crank-Nicholson scheme

Using Equations (22) and (23) the FDE at $(x_i,t^{m+1/2})$ gives the C-N scheme:

(24) $$\frac{u_i^{m+1}-u_i^m}{\tau }=\frac{1}{2}{\mathrm{\Delta }}_{2,2,}(u_i^{m+1}+u_i^m)+f_i^{m+1/2}+𝒪({\tau }^2+h^2),$$

where ${\mathrm{\Delta }}_{2,2}=K_1{\mathrm{\Delta }}_{-h,2,1}^{\alpha }+K_2{\mathrm{\Delta }}_{+h,2,1}^{\alpha }$. Rearranging (24) gives

(25) $$u_i^{m+\mathrm{‘}1}-\frac{\tau }{2}{\mathrm{\Delta }}_{2,2}{u}_i^{m+1}=u_i^m+\frac{\tau }{2}{\mathrm{\Delta }}_{2,2}^{\alpha }{u}_i^m+\tau f_i^{m+1/2}+𝒪({\tau }^3+\tau h^2)$$

for all $i=0,1,2,\mathrm{\dots },N$ and $m=0,1,2,\mathrm{\dots },M-1$. The corresponding matrix form of (25) is given by

(26) $$(I-B_{\alpha })U^{m+1}=(I+B_{\alpha })U^m+\tau F^{m+1/2}+𝒪({\tau }^3+\tau h^2)$$

for all $m=0,1,2,\mathrm{\dots },M-1$, where $B_{\alpha }=\frac{\tau }{2}\left(K_1{A}_{2,1}+K_2{A}_{2,1}^{\text{T}}\right)$,

where $w_{i-j+1,2,1}$ are the power series coefficients of the generating function $W=W_{2,2}(z)$ as seen in (17). They can be computed using the J. C. P. Miller recurrence.

Now, let ${\widehat{U}}^m$ be the solution of (26) after neglecting the $𝒪({\tau }^3+\tau h^2)$ terms with ${\widehat{U}}^m={[{\widehat{u}}_1^m,{\widehat{u}}_2^m,\mathrm{\dots },{\widehat{u}}_{N-1}^m]}^T$, where its entries ${\widehat{u}}_i^m$ become the approximate values of the exact values $u_i^m$. Also, let ${\widehat{P}}_{\alpha }$ and ${\widehat{B}}_{\alpha }$ be the reduced matrix from $P_{\alpha }$ and $B_{\alpha }$, respectively by deleting the first and last rows and columns, and ${\widehat{F}}^{m+1/2}$ be the reduced vector obtained from $F^{m+1/2}$ by removing its first and last entries. After imposing the boundary conditions (2), Equation (26) becomes to be on the ready-to-solve form:

(27) $$(\widehat{I}-{\widehat{B}}_{\alpha }){\widehat{U}}^{m+1}=(\widehat{I}+{\widehat{B}}_{\alpha }){\widehat{U}}^m+\tau {\widehat{F}}^{m+1/2}+{\widehat{𝐛}}^m,m=1,2,\mathrm{\dots },M-1,$$

where ${\widehat{𝐛}}^m={\widehat{B}}_0(u_0^{m+1}+u_0^m)+{\widehat{B}}_N(u_N^{m+1}+u_N^m)$ and ${\widehat{B}}_0$ and ${\widehat{B}}_N$ are the first($0^{th}$) and last($N^{th}$) column vectors of the matrix $B_{\alpha }$ reduced again as before.

3.5. Third-order quasi-compact Crank-Nicholson scheme

The new order $3$ quasi compact approximation is applied to numerically solve the space fractional diffusion equation with third-order accuracy. Pre-operating (1) by $P_h$ which is given by (20) gives

(28) $$P_h\frac{\partial u(x,t)}{\partial t}=K_1{P}_h{D}_{x-}^{\alpha }u(x,t)+K_2{P}_h{D}_{x+}^{\alpha }u(x,t)+P_{h}f(x,t).$$

With the aid of the second-order approximations given by Equations (22) and (23), the FDE at $(x_i,t^{m+1/2})$ gives the C-N scheme

(29) $$P_h\frac{1}{\tau }(u_i^{m+1}-u_i^m)=\frac{1}{2}{C}_h(u_i^{m+1}+u_i^m)+P_h{f}_i^{m+1/2}+𝒪({\tau }^2+h^3),$$

where $C_h=K_1{\mathrm{\Delta }}_{-h,2,1}^{\alpha }+K_2{\mathrm{\Delta }}_{+h,2,1}^{\alpha }$ Rearranging (29) yields

(30) $$\left(P_h-\frac{\tau }{2}{C}_h\right)u_i^{m+1}=\left(P_h+\frac{\tau }{2}{C}_h\right)u_i^m+\tau P_h{f}_i^{m+\frac{1}{2}}+𝒪({\tau }^3+\tau h^3).$$

Consequently, in matrix language, the C-N scheme (30) can be read as

(31) $$\left(P_{\alpha }-C_{\alpha }\right)U^{m+1}=\left(P_{\alpha }+C_{\alpha }\right)U^m+\tau P_{\alpha }{F}^{m+1/2}+𝒪\left({\tau }^3+\tau h^3\right)$$

for $m=0,1,2,\mathrm{\dots },M-1$, where $P_{\alpha }=\text{Tri}[c_2,1-2c_2,c_2]$ is a tri-diagonal matrix with size $N+1$ and $C_{\alpha }=\frac{\tau }{2}\left(K_1{A}_{2,1}+K_2{A}_{2,1}^{\text{T}}\right)$. After imposing the boundary conditions (2), Equation (31) reduces to the following form:

(32) $$({\widehat{P}}_{\alpha }-{\widehat{C}}_{\alpha }){\widehat{U}}^{m+1}=({\widehat{P}}_{\alpha }+{\widehat{C}}_{\alpha }){\widehat{U}}^m+\tau {\widehat{P}}_{\alpha }{\widehat{F}}^{m+1/2}+{\widehat{𝐛}}^m,m=0,1,2,\mathrm{\dots },M-1,$$

where ${\widehat{𝐛}}^m={\widehat{C}}_0(u_0^{m+1}+u_0^m)+{\widehat{C}}_N(u_N^{m+1}+u_N^m)$ and ${\widehat{C}}_0$ and ${\widehat{C}}_N$ are the first($0^{th}$) and last($N^{th}$) column vectors of the matrix $C_{\alpha }$ reduced again as before.

4. Stability and convergence analysis

This section analyzes the stability and convergence of the C-N schemes presented in Section 3.4 and Section 3.5 for the fractional diffusion equation. The analysis also requires certain properties of definite matrices and equivalent norms, to which the reader is referred in the references [33, 34], in addition to the following useful results.

Lemma 6 ([33]).

Let $H=(A+A^{\ast })/2$ be the Hermiatian part of a complex matrix $A$. For any eigenvalue $\lambda (A)$ of $A$ with its real part $\mathrm{\Re }(\lambda )$ , we have

$${\lambda }_{min}(H)\le \mathrm{\Re }(\lambda (A))\le {\lambda }_{max}(H),$$

where ${\lambda }_{min}(H)$and ${\lambda }_{max}(H)$ are the minimum and maximum eigenvalues of $H$, respectively.

Definition 7 ([34]).

A function $G(x)={\sum }_{n=0}^{\mathrm{\infty }}{t}_n{x}^n$ is called the generator of a Toeplitz matrix $T=[t_{i-j}]$ if

$$t_n=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }G(x)e^{𝐢nx}𝑑x.$$

Lemma 8 (Grenander-Szego theorem, [32]).

Let the generator $G(x)$ of a Toeplitz matrix $T$ be a $2\pi$-periodic continuous real-valued function. Then

$$G_{min}\le {\lambda }_{min}(T)\le {\lambda }_{max}(T)\le G_{max},$$

where $G_{min},G_{max}$ denote the minimum and maximum values of $G(x)$ respectively in $[-\pi ,\pi ]$. Moreover, if , then all eigenvalues of $T$ satisfy and furthermore, if $G_{min}\ge 0$, then $T$ is positive definite.

Lemma 9.

If $G(x)$ is the generating function for a Toeplitz matrix $T=[t_{i-j}]$, then $G(x)e^{𝐢rx}$ is the generating function of the shifted Toeplitz matrix $T_r=[t_{i-j+r}]$.

Proof.

Let $T_r=s_{i-j}$ be the Toeplitz matrix for the generating function $G(x)e^{𝐢rx}$. Then

$$s_n=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }G(x)e^{𝐢rx}{e}^{𝐢nx}𝑑x=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }G(x)e^{𝐢(n+r)x}𝑑x=t_{n+r}.$$

The result follows with $n=i-j$. ∎

Lemma 10.

The generating functions of the matrices $A_{2,r}$ and $A_{2,r}^T$ corresponding to the approximation operators ${\mathrm{\Delta }}_{-h,r,2}$ and ${\mathrm{\Delta }}_{+h,r,2}$ of the second-order approximation with shift $r$ are given by $W_{2,r}(e^{-𝐢x})e^{𝐢rx}$ and its conjugate $W_{2,r}(e^{𝐢x})e^{-𝐢rx}$, respectively.

Proof.

The matrix $A_{2,r}$ is Toeplitz given by $A_{2,r}=[t_{i-j}]=[w_{i-j+r}]$. Now,

$$\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{W}_{2,r}(e^{-𝐢x})e^{𝐢rx}{e}^{𝐢nx}𝑑x=\frac{1}{2\pi }\sum _{k=0}^{\mathrm{\infty }}{\int }_{-\pi }^{\pi }{w}_k{e}^{-𝐢kx}{e}^{𝐢(n+r)x}𝑑x=w_{n+r}.$$

The result follows with $n=i-j.$ A similar argument follows for the $A_{2,r}^T$ as well. ∎

Note that the two generating functions are mutually conjugate. Furthermore, the following results are also required.

Let $V_h=\{v|v=(v_0,v_1,\mathrm{\dots },v_N),v_i\in ℝ,v_0=0=V_N)\}$ be the space of grid functions in the interval computational domain $[a,b]$ with $N$ uniform subintervals of length $h$. Associated with the analysis carried out in [22], for $𝐮,𝐯\in V_h$, the following discrete inner products and the corresponding norms are defined below:

$$\left(\text{u, v}\right)=h\sum _{i=1}^{N-1}{u}_i{v}_i,\Vert \text{u}\Vert =\sqrt{(𝐮,𝐮)},$$

$$⟨{\delta }_h\text{u},{\delta }_h\text{v}⟩=h\sum _{i=1}^{N-1}\left({\delta }_h{u}_{i-1/2}\right)\left({\delta }_h{v}_{i-1/2}\right),{\left|\text{u}\right|}_1=\sqrt{⟨{\delta }_h\text{u},{\delta }_h\text{v}⟩}.$$

Define difference operators on the component of $𝐯\in V_h$ as

$$v_{i+1/2}=\frac{1}{2}(v_{i+1}+v_i)\text{and}{\delta }_h{v}_{i-1/2}=\frac{1}{h}(v_i-v_{i-1}).$$

Theorem 11.

Let $P$ be a self adjoint operator defined on $V_h$ such that ${\mu }_1\Vert 𝐮\Vert \le {\Vert 𝐮\Vert }_P\le {\mu }_2\Vert 𝐮\Vert ,{\mu }_1,{\mu }_2>0$ and $A$ be a negative definite operator. Suppose that there exists a vector ${𝐯}^m=(v_0^m,v_1^m,\mathrm{\dots },v_{N-1}^m,v_N^m)\in V_h$ such that

(33) $$P{\delta }_h{𝐯}^{m+1/2}=A{𝐯}^{m+1/2}+{\text{S}}^m,1\le m\le M-1,$$

provided

(34) $${𝐯}^0={\text{v}}^0(x_i)foralli=i=0,1,\mathrm{\dots },N.$$

Then,

$$\Vert {𝐯}^m\Vert \le \frac{1}{{\mu }_1}\left({\mu }_2\Vert {𝐯}^0\Vert +\frac{\tau }{{\mu }_1}\sum _{l=0}^{m-1}\Vert {𝐒}^l\Vert \right),$$

where ${𝐒}^l={[S_0^l,S_1^l,\mathrm{\dots },S_N^l]}^{\text{T}}$ with $S_0^l=0$ and $S_N^l=0$ for all $l=0,1,\mathrm{\dots },M$.

Proof.

The negative definiteness of the operator $A$ implies that

.

Now, applying inner product on (33) with ${\text{v}}^{m+1/2}$ yields:

$(P{\delta }_{\tau }{\text{v}}^{m+1/2},{\text{v}}^{m+1/2})$ $=(A{\text{v}}^{m+1/2},{\text{v}}^{m+1/2})+({\text{S}}^m,{\text{v}}^{m+1/2})$

(35) $\le ({\text{S}}^m,{\text{v}}^{m+1/2}).$

Also,

$(P{\delta }_{\tau }{\text{v}}^{m+1/2},{\text{v}}^{m+1/2})$ $=(P\frac{1}{\tau }({\text{v}}^{m+1}-{\text{v}}^m),{\displaystyle \frac{1}{2}}({\text{v}}^{m+1}+{\text{v}}^m)$

$=\frac{1}{2\tau }\left({\Vert {\text{v}}^{m+1}\Vert }_P^2-{\Vert {\text{v}}^m\Vert }_P^2\right)$

(36) $=\frac{1}{2\tau }\left({\Vert {\text{v}}^{m+1}-\Vert {\text{v}}^m\parallel }_P\right)\left({\Vert {\text{v}}^{m+1}+\Vert {\text{v}}^m\parallel }_P\right)$

$\le ({\text{S}}^m,{\text{v}}^{m+1/2})\le \Vert {\text{S}}^m\Vert \Vert {\text{v}}^{m+1/2}\Vert$

$\le \frac{1}{{\mu }_1}\Vert {\text{S}}^m\Vert {\Vert {\text{v}}^{m+1/2}\Vert }_P$

(37) $\le \frac{1}{2{\mu }_1}\Vert {\text{S}}^m\Vert \left({\Vert {\text{v}}^{m+1}\Vert }_P+{\Vert {\text{v}}^m\Vert }_P\right).$

The inequality relating (4) and (37) reduces, for $0\le m\le M-1$, to

${\Vert {\text{v}}^{m+1}\Vert }_P\le \Vert {\text{v}}^m\Vert +\frac{\tau }{{\mu }_1}\Vert {\text{S}}^m\Vert$. Summing this for the first $m$ inequalities results:

$${\Vert {\text{v}}^m\Vert }_P\le {\Vert {\text{v}}^0\Vert }_P+\frac{\tau }{{\mu }_1}\sum _{l=0}^{m-1}\Vert {\text{S}}^l\Vert ,1\le m\le M-1.$$

Equivalence of the two norms concludes the proof. ∎

Lemma 12.

Leading to the above inner products and norms, the following results :

1. The operator ${\delta }_h^2$ is self adjoint on $V_h$.

2. ${|𝐮|}_1^2\le \frac{4}{h^2}$ for any $𝐮\in V_h$.

3. The operator $T_h=1+k{h}^2{\delta }_h^2$ is selft–adjoint on $V_h$, where $k$ is a given constant.

Proof.

1. Take any $𝐮,𝐯\in V_h$. Then, we must show that $⟨{\delta }_h^2𝐮,𝐯⟩=⟨𝐮,{\delta }_h^2𝐯⟩.$ We first note that $u_0{v}_1=v_0{u}_1$ and $u_N{v}_{N-1}=v_N{u}_{N-1}$, since the vectors $𝐮$ and $𝐯$ have zero boundary values; thereby, we have:

   $({\delta }_h^2𝐮,𝐯)$ $=h{\displaystyle \sum _{i=1}^{N-1}}{\delta }_h^2{u}_i{v}_i=h{\displaystyle \sum _{i=1}^{N-1}}\left(\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}\right)v_i$

   $=h{\displaystyle \sum _{i=1}^{N-1}}\frac{u_{i+1}{v}_i-2u_i{v}_i+u_{i-1}{v}_i}{h^2}=h{\displaystyle \sum _{i=1}^{N-1}}\frac{u_i{v}_{i+1}-2u_i{v}_i+u_i{v}_{i+1}}{h^2}$

   $=h{\displaystyle \sum _{i=1}^{N-1}}{u}_i\left(\frac{(v_{i+1}-2v_i+v_{i-1})}{h^2}\right)=(𝐮,{\delta }_h^2𝐯)$

   for all $𝐮,𝐯\in V_h$. That is, ${\delta }_h^2$ is self adjoint on $V_h$.

2. From Part (a), we have:

   ${|𝐮|}_1^2$ $=⟨{\delta }_h𝐮,{\delta }_h𝐮⟩=h{\displaystyle \sum _{i=1}^{N-1}}\left({\delta }_h{u}_{i-1/2}\right)\left({\delta }_h{u}_{i-1/2}\right)$

   $=h{\displaystyle \sum _{i=1}^{N-1}}{\left(\frac{u_i-u_{i-1}}{h}\right)}^2=h{\displaystyle \sum _{i=1}^{N-1}}\frac{(u_i^2-2u_i{u}_{i-1}+u_{i-1}^2)}{h^2}$

   $\le h{\displaystyle \sum _{i=1}^{N-1}}\frac{(u_i^2+(u_i^2+u_{i-1}^2)+u_{i-1}^2))}{h^2}=\frac{2}{h^2}\left(h{\displaystyle \sum _{i=1}^{N-1}}{u}_i^2+h{\displaystyle \sum _{i=1}^{N-2}}{u}_i^2\right)$

   $\le \frac{2}{h^2}\left(h{\displaystyle \sum _{i=1}^{N-1}}{u}_i^2+h{\displaystyle \sum _{i=1}^{N-1}}{u}_i^2\right)\le \frac{4}{h^2}{\Vert 𝐮\Vert }^2.$

3. Letting $T_h=1+k{h}^2{\delta }_h^2$ and using Part (a), we get:

   $(T_h𝐮,𝐯)$ $=(𝐮,𝐯)+b_2{h}^2({\delta }_h^2𝐮,𝐯)=(𝐮,𝐯)+b_2{h}^2(𝐮,{\delta }_h^2𝐯)$

   $=(𝐮,(1+b_2{h}^2{\delta }_h^2)𝐯)=(𝐮,T_h𝐯).$

   for all $𝐮,𝐯\in V_h$. Therefore, $T_h$ is a self–adjoint operator.

∎

4.1. Analysis of the second-order C–N scheme

This section gives the stability and convergence analysis of the C–N scheme presented in Section 3.4. First, Lemma 13 is presented along with its proof:

Lemma 13.

The matrices ${\widehat{A}}_{2,1}$ and ${\widehat{A}}_{2,1}^{\text{T}}$ are negative definite for $\frac{4}{3}\le \alpha \le 2$. Therefore, the corresponding operators ${\mathrm{\Delta }}_{+h,2,1}$ and ${\mathrm{\Delta }}_{-h,2,1}$ are also negative definite.

Proof.

The generating function of matrix ${\widehat{A}}_{2,1}$ is given by $G_{\alpha }={(1-e^{-ix})}^{\alpha }{({\beta }_0+{\beta }_3{e}^{-ix})}^{\alpha /2}{e}^{ix},\frac{4}{3}\le \alpha \le 2$.

$G_{\alpha }(x)$ $={\left(R_1{e}^{i{\theta }_1}\right)}^{\alpha }{\left(R_2{e}^{i{\theta }_2}\right)}^{\alpha /2}{e}^{ix}=R{e}^{i\left(\alpha {\theta }_1+\frac{\alpha }{2}{\theta }_2+x\right)},$

where ${\theta }_1=\frac{(\pi -x)}{2}$, ${\theta }_2=-{\tan }^{-1}\left(\frac{{\beta }_3\sin (x)}{{\beta }_0+{\beta }_3\cos (x)}\right)$,
and $R=R_1^{\alpha }{R}_2^{\alpha /2}={(2\sin (x/2))}^{\alpha }{\left({\beta }_0^2+2{\beta }_0{\beta }_3\cos (x)+{\beta }_3^2\right)}^{\alpha /4}$. The real part of $G_{\alpha }(z)$ is given below: $\mathrm{\Re }(G_{\alpha })(x)=R\cos \left(\alpha {\theta }_1+\frac{\alpha }{2}{\theta }_2+x\right)=R{H}_1(x,\alpha )$ for $\frac{4}{3}\le \alpha \le 2$, where $H_1(x,\alpha )=\cos \left(\alpha {\theta }_1+\frac{\alpha }{2}{\theta }_2+x\right)$. It must be shown that . Now, if and only of , since $R\ge 0$. Therefore, it will be proved that over the domain $[0,\pi ]\times [4/3,2]$. Let $Z(x,\alpha )=\alpha {\theta }_1+\frac{\alpha }{2}{\theta }_2+x$. For a fixed $\alpha \in [4/3,2]$, differentiating $Z$ with respect to $x$ gives

$$\frac{d}{dx}Z(x.\alpha )=-\frac{(1-\cos (x))(\alpha -1)(\alpha -2)(3\alpha -4)}{{\left[2(1-\alpha )+(2-\alpha )\cos (x)\right]}^2+{\left[(2-\alpha )\sin (x)\right]}^2}.$$

The foregoing derivative assumes positive values over the interval $(0,\pi )$ for an $\alpha \in (4/3,2)$ and thus, the function $Z$ is monotonically non-decreasing function over the interval $(0,\pi )$. Therefore, the maximum and minimum values of $Z$ are $Z_{\max }=Z(\pi ,\alpha )=\pi$ and $Z_{\min }=Z(0,\alpha )=\frac{\alpha \pi }{2}$, respectively. Therefore, $H_1(x,\alpha )$ is a non increasing function and thereby its maximum for . Now, by Lemma 8, we have .

Now, for any non–zero vector $v$, consider . Thus, the matrix ${\widehat{A}}_{2,1}$ is negative definite. Consequently, ${\widehat{A}}_{2,1}^{\text{T}}$, ${\mathrm{\Delta }}_{+h,2,1}$, and ${\mathrm{\Delta }}_{-h,2,1}$ are negative definite. This completes the proof. ∎

Remark 14.

Since the matrices ${\widehat{A}}_{2,1}$ and ${\widehat{A}}_{2,1}^{\text{T}}$ are negative definite, the matrix $c_1{\widehat{A}}_{2,1}+c_2{\widehat{A}}_{2,1}^{\text{T}}$ is also negative definite for any positive constants $c_1$ and $c_2$. Therefore, the operator ${\mathrm{\Delta }}_2=c_1{\mathrm{\Delta }}_{+h,2,1}+c_2{\mathrm{\Delta }}_{-h,2,1}$ is negative definite. Or, we may prove this using the linearity property of inner product: Let any $𝐮\in V_h$. Then, we have

Theorem 15.

The Crank–Nicholson scheme (27) with the approximation from the generating function $W_{2,1}(z)$ given by (17) for the space fractional diffusion equation is unconditionally stable for $\frac{4}{3}\le \alpha \le 2$.

Proof.

We have from Equation (27) that the iteration matrix of the C–N scheme, $M_2={(I-\widehat{B})}^{-1}(I+\widehat{B})$. To establish a stability criterion, we must have the spectral radius of $M_2$, . Now, for any $\lambda ({\widehat{B}}_{\alpha })\ne 1$, we have $\lambda (M_2)=\frac{1+\lambda ({\widehat{B}}_{\alpha })}{1-\lambda ({\widehat{B}}_{\alpha })}$. Then, if and only if if and only if . Also, we have ${\widehat{B}}_{\alpha }=\frac{\tau }{2}\left(K_1{\widehat{A}}_{2,1}+K_2{\widehat{A}}_{2,1}^{\text{T}}\right)$. Now, $\lambda ({\widehat{B}}_{\alpha })=\frac{\tau (K_1+K_2)}{2}\lambda ({\widehat{A}}_{2,1})$. Therefore, since $\tau ,K_1$, and $K_2$ are positive, if and only if . Now, from Lemma 13, we have . This completes the proof. ∎

Theorem 16.

The Crank-Nicholson finite difference scheme (24) with given initial and boundary conditions converges with order $2$ for $\frac{4}{3}\le \alpha \le 2$.

Proof.

Let $e_i^m=u_i^m-{\widehat{u}}_i^m$ be error at grid point $(x_i,t_m)$, where $u_i^m$ and ${\widehat{u}}_i^m$ denote the exact solution of the diffusion equation (1) and the corresponding approximate solution given by (24). Then, $e_0^m=0$ and $e_N^m=0$.

Also, let ${𝐞}^m=(e_0^m,e_1^m,\mathrm{\dots },e_{N-1}^m,e_N^m)$ and ${𝐑}^m=(R_0^m,R_1^m,\mathrm{\dots },R_{N-1}^m,R_N^m)$, where $R_i^m$ denotes the remainder term of (24) at $(x_i,t_m)$,  $0\le i\le N,\mathrm{\hspace{0.33em}0}\le m\le M-1$, $R_0^m=0$, and $R_N^m=0$.

Now,

${\Vert {𝐑}^m\Vert }^2=h{\sum }_{i=1}^{N-1}{\left(R_i^m\right)}^2=h{\sum }_{i=1}^{N-1}{|R_i^m|}^2\le h{\sum }_{i=1}^N{|R_i^m|}^2\le Nh{c}_1^2{({\tau }^2+h^2)}^2=(b-a)c_1^2{({\tau }^2+h^2)}^2$, where $c_1$ is a positive constant. Therefore, $\parallel {𝐑}^m\parallel \le \sqrt{(b-a)}{c}_1({\tau }^2+h^2)/]$.

Also, it is easy to see that, the error vector ${𝐞}^m$ governs the difference system:

(38) $${\delta }_{\tau }{𝐞}^{m+\frac{1}{2}}-{\mathrm{\Delta }}_{2,2}{𝐞}^{m+\frac{1}{2}}={𝐑}^m,{𝐞}^0=\mathrm{𝟎}.$$

In comparison with Equation (33), in Equation (38), $P$ is the identity operator and hence it is self–adjoint and its norm is equivalent to that of $𝐮$. So, ${\mu }_1=1={\mu }_2$. Also, from Remark 14, for $K_1,K_2\ge 0$, the operator ${\mathrm{\Delta }}_{2,2}=K_1{\mathrm{\Delta }}_{+h,2,1}+K_2{\mathrm{\Delta }}_{-h,2,2,1}$ is negative definite. Then, Theorem 11 views

$\Vert {𝐞}^m\Vert$ $\le \Vert {𝐞}^0\Vert +\tau {\displaystyle \sum _{l=0}^{m-1}}\Vert {𝐑}^l\Vert \le \tau {\displaystyle \sum _{l=0}^M}\Vert {𝐑}^l\Vert \le \tau M\sqrt{(b-a)}{c}_1({\tau }^2+h^2)$

$=T\sqrt{(b-a)}{c}_1({\tau }^2+h^2)=c_2({\tau }^2+h^2),$

where $c_2=T\sqrt{(b-a)}{c}_1$.

So, we complete the proof. ∎

4.2. Analysis of the third–order C–N quasi-compact scheme

In this section, the analysis of the proposed third order quasi–compact approximation is presented.

Lemma 17.

The QCD operator of order $3$ in (3.2) leads to the following for $\frac{4}{3}\le \alpha \le 2$:

1. for $r=1$.

2. The operator $P_h$ is self–adoint and $\frac{1}{3}{\Vert 𝐮\Vert }^2\le {\Vert 𝐮\Vert }_P^2\le {\Vert 𝐮\Vert }^2$ for $r=1$, where ${\Vert 𝐮\Vert }_P^2=(P_h𝐮,𝐮)$.

Proof.

(a) It is not hard to see that, the maximum of $a_2$ over the domain $[\frac{4}{3},2]$ is $\frac{3}{2}-\sqrt{\frac{11}{6}}$, occurring at $\alpha =\sqrt{\frac{24}{11}}$ and the minimum of $a_2$ over $[\frac{4}{3},2]$ is $\frac{1}{12}$, occurring at $\alpha =2$. Therefore, for $r=1.$

(b) Take any $𝐮,𝐯\in V_h$. Applying Part (c) of Lemma 12 with $k=a_2$ gives that the operator $P_h$ is self–adjoint.

Now, using Part (b) of Lemma 12, we have: $(P_h𝐮,𝐮)={\Vert 𝐮\Vert }^2-a_2{h}^2{|𝐮|}_1^2\ge {\Vert 𝐮\Vert }^2-\frac{4}{6}{\Vert 𝐮\Vert }^2=\frac{1}{3}{\Vert 𝐮\Vert }^2$. Hence, we complete the proof. ∎

Theorem 18.

The quasi compact Crank–Nicholson scheme (32) with the approximation from the generating function $W_{2,1}(z)$ is unconditionally stable for $\frac{4}{3}\le \alpha \le 2$.

Proof.

Consider the iteration matrix, $M_3$, of the C–N scheme set off in Equation (32) given by $M_3={({\widehat{P}}_{\alpha }-{\widehat{C}}_{\alpha })}^{-1}({\widehat{P}}_{\alpha }+{\widehat{C}}_{\alpha })={(I-{\widehat{P}}_{\alpha }^{-1}{\widehat{C}}_{\alpha })}^{-1}(I+{\widehat{P}}_{\alpha }^{-1}{\widehat{C}}_{\alpha })={(I-B_{\alpha ,3})}^{-1}(I+B_{\alpha ,3})$, where $B_{\alpha ,3}={\widehat{P}}_{\alpha }^{-1}{\widehat{C}}_{\alpha }$. Now, arguing analogously to Section 4.1, the spectral radius of matrix $M_3$, if and only if . The eigen-values of ${\widehat{P}}_{\alpha }$ are given by

$\lambda {({\widehat{P}}_{\alpha })}_m$ $=1-2a_2+2a_2\cos \left(\frac{m\pi }{N}\right)=1-2a_2(1-\cos \left(\frac{m\pi }{N}\right)$

$=1-4a_2{\sin }^2\left(\frac{m\pi }{2N}\right)=4a_2\left(\frac{1}{4a_2}-{\sin }^2\left(\frac{m\pi }{2N}\right)\right)$

$>0,m=1,2,\mathrm{\dots },N,$

since when $r=1$, for $\frac{4}{3}\le \alpha \le 2.$ Thus, $\lambda ({\widehat{P}}_{\alpha }^{\alpha })>0$; thereby we have if and only if .

Now, , since owing to Lemma 13. This results in giving . Therefore, the order $3$ Crank–Nicholson scheme is unconditionally stable for $4/3\le \alpha \le 2$. ∎

Theorem 19.

The quasi compact Crank–Nicholson finite difference scheme (31) with the given initial and boundary conditions converges with order $3$ to the exact solution of the diffusion problem (1) for $\frac{4}{3}\le \alpha \le 2$.

Proof.

In Part (b) of Lemma 17, we have proved that the operator $P_h$ is self–adjoint and $\frac{1}{3}{\Vert 𝐮\Vert }^2\le {\Vert 𝐮\Vert }_P^2\le {\Vert 𝐮\Vert }^2$. Also, it is not hard to see that the error vector ${𝐞}^m$ governs the difference system:

(39) $$P_h{\delta }_t{𝐞}^{m+\frac{1}{2}}-{\mathrm{\Delta }}_{2,2}{𝐞}^{m+\frac{1}{2}}={𝐑}^m,{𝐞}^0=\mathrm{𝟎},$$

where $𝐑=(R_0^m,R_1^m,\mathrm{\dots },R_{N-1}^m,R_N^m)$, where $R_i^m$ denotes the remainder term of (29) at $(x_i,t_m),0\le i\le N,0\le m\le M,R_0^m=0$ and $R_N^m=0$. Now, ${\Vert {𝐑}^m\Vert }^2=h{\sum }_{i=1}^{N-1}{(R_i)}^2\le h{\sum }_{i=1}^N{|R_i^m|}^2\le (b-a)c_3^2{({\tau }^2+h^3)}^2$. where $c_3$ is a positive constant. This implies that $\Vert {𝐑}^m\Vert \le \sqrt{(b-a)}{c}_3({\tau }^2+h^3)$.

Now, using Theorem 11 with ${\mu }_1=\frac{1}{\sqrt{3}},{\mu }_2=1$ and ${𝐒}^m={𝐑}^m$, we have:

$\Vert {𝐞}^m\Vert$ $\le \Vert {𝐞}^0\Vert +3\tau {\displaystyle \sum _{l=1}^{m-1}}\Vert {𝐑}^l\Vert \le \tau M\sqrt{b-a}{c}_3({\tau }^2+h^3)=T\sqrt{(b-a)}({\tau }^2+h^3)$

$=c_4({\tau }^2+h^3),$

where $c_4=T\sqrt{(b-a)}{c}_3$. So, we complete the proof. ∎

5. Numerical results

In this section, numerical examples are given to demonstrate the unconditional stability, convergence order, and accuracy of each scheme derived in Section 3. The following test example is considered:

Example 20.

Let $H(x,m,\alpha )=\frac{\mathrm{\Gamma }(m+1)}{\mathrm{\Gamma }(n+1-\alpha )}(x^{m-\alpha }+{(1-x)}^{m-\alpha })$ and $s_0(x)=x^5{(1-x)}^5$. The following example uses constant diffusion coefficients.

$K_1(x)$ $=1,K_2(x)=1,$

$f(x,t)$ $=-e^{-t}(s_0(x)+H(x,5,\alpha )-5H(x,6,\alpha )+10H(x,7,\alpha )$

$-10H(x,8,\alpha )+5H(x,9,\alpha )-H(x,10,\alpha ))$

$u(x,0)$ $=s_0(x),u(0,t)=0,u(1,t)=0.$

$Exact$ $Solutionu(x,t)=s_0(x)e^{-t}.$

Let, at a time final time $t=T$, the exact solution vector be defined by $U$ and a corresponding approximate solution vector be denoted by $\widehat{U}$. Then the maximum norm of error vector $\widehat{U}-U$ at grid size $h$ is given by $E_h={\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}={\max }_{1\le i\le n}|U_i-\widehat{U_i}\Vert$. The numerical convergence order $c$ is calculated by $c=\log (E_h/E_{h/2})/\log 2$. First, the second-order Crank-Nichoslon scheme is applied to Example 20 to calculate the errors and convergence orders for $\alpha =1.34,1.5$, and $1.9$. We choose $N=8,16,32,64,128,256,512,1024,2048=M$ with uniform sub-interval sizes $h=1/N$ and $\tau /M$. We then apply the new order $3$ quasi compact C-N-scheme presented in Section 3.5 to Example 20. The space domain is handled with $N$ sub-partitions and time domain is handled with $M=[N^{3/2}]+1$ sub-partitions. Table 1 and Table 2 demonstrate maximum absolute errors and convergence orders of these schemes.

$h=\tau$	$\alpha =1.34$		$\alpha =1.5$		$\alpha =1.9$
$N=M$	${\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}$	$c$	${\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}$	$c$	${\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}$	$c$
8	3.6245e-05	–	3.4278e-05	–	1.3912e-05	–
16	8.5169e-06	2.08	8.3193e-06	2.04	5.6605e-06	1.29
32	2.0773e-06	2.03	2.0709e-06	2.00	1.4222e-06	1.99
64	5.1597e-07	2.00	5.1907e-07	1.99	3.5503e-07	2.00
128	1.2880e-07	2.00	1.3012e-07	1.99	8.8753e-08	2.00
256	3.2192e-08	2.00	3.2587e-08	1.99	2.2192e-08	1.99
512	8.0477e-09	2.00	8.1546e-09	1.99	5.5488e-09	1.99
1024	2.0120e-09	1.99	2.0397e-09	1.99	1.3877e-09	1.99
2048	5.0323e-10	1.99	5.1009e-10	1.99	3.4964e-10	1.98

Table 1. Maximum errors and for order $2$ convergence of C-N scheme at $T=1$.

$h=\frac{1}{N}$	$\tau =\frac{1}{M}$	$\alpha =1.34$		$\alpha =1.5$		$\alpha =1.9$
$N$	M	${\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}$	$c$	${\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}$	$c$	${\Vert \widehat{U}-U\Vert }_{\mathrm{\infty }}$	$c$
8	23	3.3799e-06	–	1.3300e-06	–	1.3815e-06	–
16	65	1.1898e-07	4.82	2.1456e-07	2.63	3.0808e-08	5.48
32	182	5.3530e-09	4.47	3.1650e-08	2.76	1.9968e-09	3.94
64	513	2.8116e-10	4.25	4.2935e-09	2.88	4.6661e-10	2.09
128	1449	3.6421e-11	2.94	5.5870e-10	2.94	7.2068e-11	2.69
256	4097	6.5222e-12	2.48	7.1243e-11	2.97	9.8501e-12	2.87
512	11586	9.2923e-13	2.81	8.9930e-12	2.98	1.2345e-12	2.99

Table 2. Maximum errors and order $3$ convergence of C-N scheme at $T=1$.

Figure 1. Surface plot of the exact solution of Example 20 over $[0,1]\times [0,1].$

Both Table 1 and Table 2 confirm convergence orders, unconditional stability, and the accuracy of the second and third schemes, respectively. Furthermore, Fig. 1 exhibits the surface plot of the exact solution of the fractional diffusion equation in Example 20 over the domain $[0,1]\times [0,1]$.

6. Conclusion

In this paper, we present two new approximations for fractional derivatives, utilizing a recently developed unified explicit form. The first approximation achieves second-order accuracy, while the second approximation demonstrates third-order accuracy, derived from the former using a quasi-compact technique. These approximations were employed, together with the Crank-Nicholson method, to solve the space fractional diffusion equation. The unconditional stability and convergence of the resulting Crank-Nicholson schemes were established for fractional derivatives of order $\alpha$ in the interval $\frac{4}{3}\le \alpha \le 2$. Furthermore, numerical results confirm the unconditional stability and convergence orders.