2 Grünwald-type approximations

In this section, we summarize the theory developed in [Reference Nasir and Nafa14] for constructing generating functions for SGD approximations of fractional derivatives with arbitrary orders and shifts.

Let $f(x)$ be defined in an infinite interval $(-\infty , a]$ or $[a, +\infty )$ . For a sequence of generic real weights $\{ w_{j}^{(\alpha )} \} $ with generating function

(2.1) $$ \begin{align}W(z) = \sum_{j=0}^{\infty} w_{j}^{(\alpha)} z^j , \end{align} $$

define the left and/or right GTD formula with shift r as

(2.2) $$ \begin{align} \Delta_{h\mp, r}^\alpha f(x) = \frac{1}{h^\alpha } \sum_{j=0}^\infty w_{j}^{(\alpha)} f(x \mp (j-r) h), \end{align} $$

where h is the subinterval size of a uniform partitioning of the domain of $f(x)$ into (infinite) subintervals. For simplicity of presentation, we do not distinguish the sequence of weights and its generating function.

We denote the generating function for an approximation with order p and shift $ r$ as $W_{p,r}(z) $ , and the approximation operator as $\Delta _{h\mp ,p,r}^\alpha $ . An equivalent characterization of the approximating generator $ W_{p,r}(z)$ for order $p\ge 1$ and shift r is given by the following theorem.

Theorem 2.2. Let $n =[\alpha ]+1, \; m$ be a nonnegative integer, $f(x)\in C^{m + n + 1}(\mathbb {R})$ and $D_x^k f(x) = d^k f(x)/dx^k \in L^1(\mathbb {R})$ for integers $k, \; 0 \le k \le m+n+1$ . Also, let the series of the generating function $W_{p,r}(z) $ given in (2.1) be absolutely convergent in a neighbourhood of the origin and

$$ \begin{align*} G_r(z) := \frac{1}{z^\alpha} W_{p,r}(e^{-z})e^{rz} \end{align*} $$

be analytic. Then, $W_{p,r}(z)$ approximates the left and right fractional differential operators $D_{x\mp }^{\alpha } $ with order p and shift r, $ 1\le p\le m $ , if and only if

(2.3) $$ \begin{align} G_r(z) = 1 + \sum_{l=p}^{\infty} a_{l} z^l = 1 + O(z^p) \end{align} $$

and

(2.4) $$ \begin{align} \Delta_{h\mp,p,r}^\alpha f(x) = D_{x\mp}^{\alpha} f(x) + h^p a_p D_{x\mp}^{\alpha + p} f(x) + h^{p+1} a_{p+1} D_{x\mp}^{\alpha + p+1} f(x) + \cdots + O(h^m), \end{align} $$

where $a_l \equiv a_l(r,\alpha )$ .

An immediate consequence of this is the following consistency condition for the approximation.

Using these results, Nasir and Nafa [Reference Nasir and Nafa14] constructed, by algebraic manipulations, generating functions that approximate the fractional differential operators $D_{x\mp }^\alpha $ with order p and shift r by considering fractional order powers of polynomials

$$ \begin{align*} W_{p,r} (z) = (b_0+b_1 z +\cdots + b_p z^p)^\alpha , \end{align*} $$

where the coefficients, $b_j\equiv b_j(\alpha ,r) $ , are to be determined from the condition (2.3).

For a second-order approximation ( $p = 2$ ), $W_{2,r}(z) = ( b_0+b_1 z +b_2 z^2 )^\alpha $ is obtained with

(2.5) $$ \begin{align} b_0 = \frac{3}{2} - \frac{r}{\alpha} , \quad b_1 = -2 + \frac{2r}{\alpha} \quad \text{ and } \quad b_2 = \frac{1}{2} - \frac{r}{\alpha}. \end{align} $$

We then have from (2.4), when $f(x)\in C^{ n+4}(\mathbb {R})$ ,

$$ \begin{align*} \Delta^\alpha_{h\mp,2,r}f(x) = D_{x\mp}^{\alpha}f(x) + h^2 a_2 D_{x\mp}^{\alpha+2}f(x) + h^3 a_3 D_{x\mp}^{\alpha+2}f(x) +O(h^4), \end{align*} $$

where the coefficients $a_2,\, a_3$ of the first two leading error terms are given by

$$ \begin{align*} a_2(r) = \frac{\alpha}{6}( -3 \lambda^{2} + 6 \lambda - 2) \quad \text{ and } \quad a_3(r) = \frac{\alpha (- 4 \lambda^{3} + 12 \lambda^{2} - 11 \lambda + 3)}{12}, \end{align*} $$

where $\lambda = {r}/{\alpha }$ .

Note that the second-order generating function with the coefficients in (2.5) reduces to the second-order Lubich generating function $ W_2(z)$ in (1.3) when there is no shift ( $r = 0$ ). We have from (1.3),

$$ \begin{align*} W_2(z) = \bigg( \sum_{k = 1}^{2} \frac{1}{k}(1-z)^k \bigg)^\alpha=\bigg( (1-z) + \frac{1}{2}(1-z)^2 \bigg)^\alpha=\bigg( \frac{3}{2} - 2z +\frac{1}{2}z^2 \bigg)^\alpha. \end{align*} $$

Hence, we may regard $W_{2,r}(z) $ as a shifted form of $W_2(z)$ . The coefficients $w_m^{(\alpha )} $ of the generating function $W_{2,r}(z) = \sum _{m = 0}^{\infty } w_m^{(\alpha )} z^m$ can be computed from the J.C.P. Miller formula,

$$ \begin{align*} w_0^{(\alpha)} & = b_0^\alpha, \quad w_1^{(\alpha)} = \alpha b_0^{\alpha-1} b_1, \\ w_m^{(\alpha)} & = \frac{1}{mb_0} [ (\alpha+1-m)w_{m-1}^{(\alpha)}b_1 + (2(\alpha+1)-m)w_{m-2}^{(\alpha)} b_2 ], \quad m = 2,3,\ldots. \end{align*} $$

Nasir and Nafa [Reference Nasir and Nafa13] used this second-order approximation to compute solutions for the space-fractional diffusion problems. They also constructed a third-order approximation using the QCA approach based on this second-order approximation.

3 Fourth-order approximations

In this section, we study the two fourth-order approximations for the fractional derivatives using the second-order shifted form presented in [Reference Gunarathna, Nasir and Daundasekara7].

The approximations are obtained by combining the WSGD and QCA approaches. We present the construction for the left fractional derivative, as analogous construction also applies to the right fractional derivative.

For two distinct integer shifts $r_1$ and $r_2$ , define an operator

$$ \begin{align*} B_{h-,(r_1,r_2)}^\alpha = \lambda_1 \Delta_{h-,r_1,2}^\alpha + \lambda_2 \Delta_{h-,r_2,2}^\alpha, \end{align*} $$

where the coefficients $\lambda _1$ and $\lambda _2$ are to be determined.

The generating function for $B_{h-,(r_1,r_2)}^\alpha $ is given by

$$ \begin{align*} W(z) = \lambda_1 W_{2,r_1}(z) + \lambda_2 W_{2,r_2}(z). \end{align*} $$

Expanding the WSGD generating function

$$ \begin{align*} G(z)&=\lambda_1 \frac{1}{z^{\alpha}} W_{2,r_1} (e^{-z}) + \lambda_2 \frac{1}{z^{\alpha}} W_{2,r_2} (e^{-z})\\ &=(\lambda_1 + \lambda_2) + (\lambda_1 a_2(r_1) + \lambda_2a_2(r_2)) z^2 + (\lambda_1 a_3(r_1) + \lambda_2a_3(r_2))z^3 + O(z^4), \end{align*} $$

and imposing conditions

(3.1a) $$ \begin{align} &\lambda_1 +\lambda_2 = 1, \end{align} $$

(3.1b) $$ \begin{align} &\lambda_1 a_3(r_1) +\lambda_2 a_3(r_2) = 0, \end{align} $$

we get

(3.2) $$ \begin{align} G(z) = 1 + {\tilde a}_2(r_1,r_2) z^2 + O(z^4), \end{align} $$

where ${\tilde a}_2(r_1,r_2) = \lambda _1 a_2(r_1) + \lambda _2a_2(r_2) $ .

Solving the system (3.1), we get expressions for $\lambda _1 $ and $\lambda _2$ :

$$ \begin{align*} \lambda_1(r_1,r_2) = \frac{-a_3(r_2)}{a_3(r_1) - a_3(r_2)}, \quad \lambda_2(r_1,r_2) = \frac{a_3(r_1)}{a_3(r_1) - a_3(r_2)}, \quad a_3(r_1) \ne a_3(r_2). \end{align*} $$

With the choices of shifts $(r_1,r_2) = (1,0)$ and $(1,-1)$ , these weights and coefficient ${\tilde a}_2(r_1,r_2) $ reduce to

(3.3) $$ \begin{align}\begin{aligned} \lambda_1(1,0) &= \frac{3 \alpha^{3}}{11 \alpha^{2} - 12 \alpha + 4 },\quad \lambda_2(1,0) = - \frac{(\alpha - 2) (\alpha - 1) (3 \alpha - 2)}{11 \alpha^{2} - 12 \alpha + 4}, \\[3pt]{\tilde a}_2(1,0) &=\frac{- 4 \alpha^{3} + 15 \alpha^{2} - 8 \alpha}{6(11 \alpha^{2} - 12 \alpha + 4) },\end{aligned} \end{align} $$

and

(3.4) $$ \begin{align}\begin{aligned} \lambda_1(1,-1) &= \frac{(\alpha + 1) (\alpha + 2) (3 \alpha + 2)}{2 (11 \alpha^{2} + 4)}, \, \lambda_{2}(1,-1) = \frac{-(\alpha - 2) (\alpha - 1) (3 \alpha - 2)}{2 (11 \alpha^{2} + 4)}, \\[3pt]{\tilde a}_2(1,-1) &= \frac{- 4 \alpha^{4} + 31 \alpha^2 - 12}{6\alpha (11 \alpha^{2} + 4)}, \end{aligned} \end{align} $$

respectively. In what follows, we refer to approximations (3.3) and (3.4) as QCA1 and QCA2, respectively.

From (3.2) and (2.4),

(3.5) $$ \begin{align} B_{h-,(r_1, r_2) }^\alpha f(x) &= D_{x-}^\alpha f(x) + h^2 {\tilde a}_2(r_1,r_2) D_{x-}^{2+\alpha} f(x) + O(h^{4})\nonumber\\ &=[I + h^2 {\tilde a}_2(r_1,r_2) D_x^2] D_{x-}^\alpha f(x) + O(h^{4}),\nonumber\\ &=P_x D_{x-}^\alpha f(x) + O(h^{4}),\end{align} $$

where $P_x = I + h^2 {\tilde a}_2(r_1,r_2) D_x^2, \, D^2_x = d^2/dx^2$ and I denotes the identity operator.

Using the discretization of the domain with uniform subintervals of size h, the second derivative $D_x^2$ is approximated by the second-order central difference approximation

$$ \begin{align*} \delta_h^2 f(x) = \frac{1}{h^2} (f(x-h) -2f(x)+f(x+h)) \end{align*} $$

so that

$$ \begin{align*} D_x^2 f(x) = \delta_h^2 f(x) + O(h^2) \end{align*} $$

and $P_x$ is approximated by

(3.6)

where $P_h = I + h^2 {\tilde a}_2 (r_1,r_2) \delta _h^2 $ .

From (3.5) and (3.6), with $(r_1,r_2) = (1,0) , (1,-1)$ , the fourth-order approximations for the operator $P_h D_{x-}^\alpha $ are given by

(3.7) $$ \begin{align} B_{h-,(r_1,r_2)}^\alpha f(x) = P_h D_{x-}^\alpha f(x) + O(h^4), \quad (r_1,r_2) = (1,0) , (1,-1). \end{align} $$

4 Approximation of fractional diffusion equation

For an application of our proposed higher-order approximation formula in (3.7), we consider the numerical approximation of the space fractional diffusion equation in the domain $[a,b]\times [0,T]$ :

(4.1) $$ \begin{align} \frac{\partial }{\partial t}u(x,t) = K_1 D_{x-}^{\alpha}u(x,t) +K_2 D_{x+}^{\alpha}u(x,t) + f(x,t), \quad x\in[a,b], \, t\in[0,T], \end{align} $$

with the initial and boundary conditions:

$$ \begin{align*} u(x,0)=s_0(x) , \quad x\in[a,b], \\ u(a,t)=\psi_1(t),\quad u(b,t)=\psi_2(t), \quad t\in[0,T], \end{align*} $$

where $u(x.t)$ is the unknown function to be determined; $K_1, K_2$ are nonnegative constant diffusion coefficients with $K_1 + K_2 \ne 0$ , that is, not both are simultaneously zero; $f(x,t)$ is a known source function.

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.

The space domain $[a,b]$ is partitioned by uniformly spaced $N+1$ points with sub-interval size $h = (b-a)/N$ . The time domain $[0,T]$ is discretized similarly by $M+1$ points with subinterval length $ \tau = T/M$ . These form a uniform mesh on the two-dimensional 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 $ . We also introduce the following notation:

$$ \begin{align*} u_i^m &= u(x_i,t^m),\quad t^{m+1/2} = \tfrac{1}{2}(t^{m+1}+t^m ),\quad f_i^{m+1/2} = f(x_i, t^{m+1/2}), \\ U^{m}&=[u_{0}^{m}, u_{1}^{m}, \ldots, u_{N}^{m}]^{T} \quad \text{and} \quad F^{m+1/2}=[f_{0}^{m+1/2},f_{1}^{m+1/2}, \ldots, f_{N}^{m+1/2}]^{T}. \end{align*} $$

Pre-multiplying (4.1) by $P_h$ , at time $t+\tau /2$ ,

$$ \begin{align*} P_{h}\frac{\partial u(x,t+\tau/2)}{\partial t} = K_1 P_h D_{x-}^{\alpha}u(x,t+\tau/2) + K_2 P_h D_{x+}^{\alpha}u(x,t+\tau/2) + P_h f(x,t+\tau/2). \end{align*} $$

Using the second-order approximations

$$ \begin{align*} \frac{\partial u(t+\tau/2)}{\partial t}&= \frac{1}{\tau}(u(x,t+ \tau)-u(x,t)) +O(\tau^2),\\ u(x,t+\tau/2)&= \frac{1}{2}(u(x,t+\tau)+u(x,t))+O(\tau^2), \end{align*} $$

the FDE at $(x_i, t^{m+1/2})$ gives the CN scheme:

(4.2) $$ \begin{align} P_h \frac{1}{\tau}({u_i^{m+1}-u_i^m}) = \frac{1}{2}L_{h}(u_i^{m+1}+u_i^m) + P_h f_i^{m+1/2} +O(\tau^2 +h^{4}), \end{align} $$

where $L_{h}=K_1 B_{h-,r_1,r_2} +K_2 B_{h+,r_1,r_2} $ .

Rearranging this,

$$ \begin{align*} \bigg(P_{h}-\frac{\tau}{2}L_{h}\bigg)u_{i}^{m+1}= \bigg(P_{h}+\frac{\tau}{2} L_{h}\bigg)u_{i}^{m} + \tau P_h f_i^{m+1/2} +O(\tau^3 +\tau h^{4}) \end{align*} $$

and the corresponding matrix form

(4.3) $$ \begin{align} (P_\alpha - DL_\alpha)U^{m+1} = (P_\alpha + DL_\alpha)U^{m} + \tau P_\alpha F^{m+1/2} + O(\tau^3+\tau h^4), \quad 1\le m \le M,\end{align} $$

where $D = (K1+K2)I$ ,

$$ \begin{align*} L_{\alpha}=\frac{\tau}{2}\bigg[\frac{K_1}{K1+K2}(\lambda_1A_{2,r_{1}}+ \lambda_2A_{2,r_{2}}) + \frac{K_2}{K1+K2}(\lambda_1A_{2,r_{1}}+ \lambda_2A_{2,r_{2}})^{T}\bigg], \end{align*} $$

and $P_{\alpha }=\text {Tri}[{\tilde a}_{2}, 1-2{\tilde a}_{2}, {\tilde a}_{2}]$ is a tridiagonal matrix of size $N+1$ corresponding to $P_{h}$ , and ${\tilde a}_{2} \equiv {\tilde a}_{2}(r_{1},r_{2})$ .

Now, let ${\hat U}^m = [{\hat u}_1^m, {\hat u} _2^m,\ldots , {\hat u}_{N-1}^m]^T $ be the solution of (4.3) after neglecting the error terms, where the entries ${\hat u}_i^m $ become the approximate values of the exact values $u_i^m$ . Also, let ${\hat P} _\alpha $ and ${\hat L}_\alpha $ be the reduced matrices from $P_\alpha $ and $ L_\alpha $ , respectively, by deleting the first and the last rows and columns, and ${\hat 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, (4.3) reduces to the ready-to-solve form

$$ \begin{align*} ({\hat P}_\alpha-D{\hat L}_\alpha){\hat U}^{m+1} = ({\hat P}_\alpha+D{\hat L} _\alpha){\hat U}^{m} + \tau {\hat P}_\alpha{\hat F}^{m+1/2} + \mathbf{\hat b} ^m, \quad m = 1,2,\ldots, M, \end{align*} $$

where the column vector $\mathbf {\hat b}^m = {\hat L}_0 (u^{m+1}_0 + u^{m}_0) + {\hat L}_N (u^{m+1}_N + u^{m}_N) $ and ${\hat L}_0 $ and ${\hat L}_N $ are respectively the first ( $ 0^{th}$ ) and last ( $N^{th}$ ) column vectors of the matrix $L_\alpha $ reduced as before.

4.1 DTTS iteration

To solve the linear system

$$ \begin{align*} (P_\alpha - DL_\alpha )w = F, \end{align*} $$

we re-write the matrix $P_\alpha - DL_\alpha $ as

$$ \begin{align*} P_\alpha - DL_\alpha = D(\sigma P_\alpha -L_\alpha) - (\sigma D - I) P_\alpha. \end{align*} $$

Then, the fixed-point iteration can be formed as in [Reference She16],

$$ \begin{align*} D(\sigma P_\alpha -L_\alpha)w_{k+1} &= (\sigma D - I) P_\alpha w_k + F. \end{align*} $$

Choose $\sigma $ such that the matrix $\sigma P_\alpha - L_\alpha $ is diagonally dominant.

5 Stability and convergence

In this section, we analyse the stability of the numerical schemes constructed in Section 3. We closely follow the analysis of Hao et al. [Reference Hao, Sun and Cao8].

Let $V_h = \{ \mathbf {v} = (v_0, v_1 , \ldots , v_N)| v_i \in \mathbb {R} , v_0 = v_N = 0 \} $ be the space of grid functions in the computational domain in space interval $[a, b]$ with N uniform subintervals of size h. Define a difference operator on the components of vector $\mathbf {v} \in V_h$ by

$$ \begin{align*} v_{i+1/2} =\frac{1}{2}( v_{i+1} + v_{i}) \quad \text{ and } \quad \delta_h v_{i-1/2} = \frac{1}{h}( v_{i} - v_{i-1}). \end{align*} $$

Also, consider a sequence $\mathbf {u}^m = (u_0^m, u_1^m,u_2^m, \ldots ,u_{N-1}^m, u_N^m) $ of vectors in $V_h$ . Let

$$ \begin{align*} u^{m+1/2} = \frac{1}{2}(u^{m+1} + u^m) \quad \text{and}\quad \delta_\tau u^{m + 1/2} = \frac{1}{\tau}( u^{m+1} - u^{m}).\end{align*} $$

For any $\mathbf {u},\mathbf {v} \in V_h$ , define the discrete inner products and norms as

$$ \begin{align*} (\mathbf{u},\mathbf{v})&= h\sum_{i=1}^{N-1} u_i v_i,\quad \|\mathbf{u}\|=\sqrt{(\mathbf{u},\mathbf{u})}, \\(\delta_{h} \mathbf{u}, \delta_{h} \mathbf{v})&= h\sum_{i=1}^{N-1} (\delta_{h} u_{i-1/2})(\delta_{h} v_{i-1/2}), \quad |\mathbf{u}|_1= \sqrt{({\delta_{h} \mathbf{u}, \delta_{h} \mathbf{u})}}. \end{align*} $$

Let $P $ be a self-adjoint difference operator on $V_h$ and the norm $\|\mathbf {u}\|_P^2 =(P\mathbf {u},\mathbf {u}) $ be equivalent to the norm $\|\mathbf {u}\|$ . That is, there exist positive constants $\mu _1 $ and $ \mu _2 $ such that

(5.1) $$ \begin{align} \mu_1 \|\mathbf{u}\| \le \|\mathbf{u}\|_P \le \mu_2 \|\mathbf{u}\|. \end{align} $$

Then, we have the following theorem.

Theorem 5.1. Let P be a self-adjoint difference operator on $V_h$ satisfying the inequality (5.1) and A be a negative-definite operator. Consider the discrete problem defined on $V_h$ .

Find $\mathbf {v}^m = (v_0^m, v_1^m,v_2^m, \ldots ,v_{N-1}^m, v_N^m) \in V_h $ that satisfies the system of difference equations with the given initial and boundary conditions:

(5.2) $$ \begin{align}\begin{aligned} & P \delta_\tau \mathbf{v}^{m+1/2} = A \mathbf{v}^{m+1/2} + \vec{S}^m, \\& \mathbf{v}^0 = \mathbf{v}^0(x_i) ,\quad 0\le i\le N. \end{aligned}\end{align} $$

Then,

$$ \begin{align*} \|\mathbf{v}^m\| \le \frac{1}{\mu_1} \bigg(\mu_2 \|\mathbf{v}^0\| + \frac{\tau}{\mu_1} \sum_{l=0}^{m-1} \|\vec{S}^l\| \bigg), \end{align*} $$

where $\vec {S}^l = [S_0^l, S_1^l, S_2^l, \ldots ,S_{N-1}^l, S_N^l]^T$ .

Proof. Since the operator A is negative definite, we have $(A \mathbf {v}^{m+1/2},\mathbf {v}^{m+1/2}) \le 0$ . Applying inner product on (5.2) with $\mathbf {v}^{m+1/2}$ gives

$$ \begin{align*} (P\delta_\tau \mathbf{v}^{m+1/2}, \mathbf{v}^{m+1/2}) = (A \mathbf{v}^{m+1/2}, \mathbf{v}^{m+1/2}) + (S^m,\mathbf{v}^{m+1/2})\le (S^m,\mathbf{v}^{m+1/2}). \end{align*} $$

Also,

(5.3) $$ \begin{align} (P\delta_\tau \mathbf{v}^{m+1/2}, \mathbf{v}^{m+1/2}) &= \bigg(P \frac{1}{\tau} (\mathbf{v}^{m+1} - \mathbf{v}^m) , \frac{1}{2} (\mathbf{v}^{m+1} + \mathbf{v}^m)\bigg)\nonumber\\ &=\frac{1}{2\tau}( (P_h v^{m+1} , v^{m+1}) - (P_h v^{m} , v^{m}) )\nonumber\\ &=\frac{1}{2\tau} ( \|\mathbf{v}^{m+1}\|_P^2 - \|\mathbf{v}^{m}\|_P^2 ) \nonumber\\ &=\frac{1}{2\tau} ( \|\mathbf{v}^{m+1}\|_P - \|\mathbf{v}^{m}\|_P )( \|\mathbf{v}^{m+1}\|_P + \|\mathbf{v}^{m}\|_P ). \end{align} $$

Again,

(5.4) $$ \begin{align} (P\delta_\tau \mathbf{v}^{m+1/2}, \mathbf{v}^{m+1/2}) & \le (\vec{S}^m, \mathbf{v}^{m+1/2}) \le \|\vec{S}^m\| \|\mathbf{v}^{m+1/2}\| \nonumber \\ &\le \frac{1}{\mu_1} \|\vec{S}^m\| \|\mathbf{v}^{m+1/2}\|_P \nonumber \\ &\le \frac{1}{2\mu_1} \|\vec{S}^m\| (\|\mathbf{v}^{m+1}\|_P +\|\mathbf{v}^{m}\|_P). \end{align} $$

The inequality relating (5.3) and (5.4) reduces, for $0\le m\le M-1$ , to $\|\mathbf {v}^{m+1}\|_P \le \|\mathbf {v}^{m}\|_P + ({\tau }/{\mu _1}) \|\vec {S}^m\|.$ Summing this for the first m inequalities,

$$ \begin{align*} \|\mathbf{v}^m\|_P \le \|\mathbf{v}^0\|_P + \frac{\tau}{\mu_1} \sum_{l=0}^{m-1 }\|\vec{S}^l\|, \quad 1 \le m \le M-1. \end{align*} $$

The equivalence of the two norms then gives

$$ \begin{align*} \|\mathbf{v}^m\| \le \frac{1}{\mu_1}\bigg( \mu_2 \|\mathbf{v}^0\| + \frac{\tau}{\mu_1} \sum_{l=0}^{m-1 }\|\vec{S}^l\| \bigg). \end{align*} $$

This completes the proof.

Proof. It is enough to prove that the spectral radius of the matrix $M=(I-B_\alpha )^{-1} (I+B_\alpha ) $ satisfies $\rho (M) \le 1$ , where $B_{\alpha }=P_{\alpha }^{-1}L_{\alpha }$ and $\rho (M)$ is the spectral radius of M. Note that for any eigenvalue $\lambda (B_{\alpha })(\neq 1) $ of $B_\alpha $ , $(1+\lambda (B_{\alpha }))/(1-\lambda (B_{\alpha }))$ is an eigenvalue of M. Thus, $ \lambda (B_{\alpha }) \le 0 $ if and only if $|1+\lambda (B_{\alpha })| \le | 1-\lambda (B_{\alpha })|$ , if and only if

$$ \begin{align*} \rho(M) \le \max_\lambda \bigg|\frac{1+\lambda(B_{\alpha})}{1-\lambda(B_{\alpha})} \bigg|\le 1.\\[-45pt] \end{align*} $$

Our analysis of stability follows the following lemmas on square matrices.

Lemma 5.3. For a complex square matrix A, let $H = (A+A^*)/2 $ be the Hermitian part of A. For any eigenvalue $\lambda (A)$ of A,

$$ \begin{align*} \lambda_{\min}(H) \le \Re(\lambda(A)) \le \lambda_{\max}(H), \end{align*} $$

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

We further require the following on Toeplitz matrices.

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

$$ \begin{align*} t_n = \frac{1}{2\pi} \int_{-\pi}^{\pi} G(x) e^{ \mathbf{i} nx}\,dx. \end{align*} $$

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

$$ \begin{align*} G_{\min} \le \lambda_{\min}(T) \le \lambda_{\max} (T) \le G_{\max}, \end{align*} $$

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

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

$$ \begin{align*} s_n = \frac{1}{2\pi}\int_{-\pi}^{\pi} G(x) e^{\mathbf{i}rx} e^{\mathbf{i}n x}\,dx = \frac{1}{2\pi}\int_{-\pi}^{\pi} G(x) e^{\mathbf{i}(n+r)x}\,dx = t_{n+r}. \end{align*} $$

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

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

$$ \begin{align*} \frac{1}{2\pi}\int_{-\pi}^{\pi} W_{2,r}(e^{-\mathbf{i} x}) e^{\mathbf{i} rx} e^{\mathbf{i} nx}\,dx = \frac{1}{2\pi}\sum_{k=0}^{\infty} \int_{-\pi}^{\pi} w_k e^{-\mathbf{i}k x} e^{\mathbf{i} (n+r) x}\,dx = w_{n+r}. \end{align*} $$

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. We further establish the following results on the generating function. For analytical simplicity, we write the generating function $W_{2,r}(z) $ in factored form as

$$ \begin{align*} W_{2,r}(z) = (1-z)^\alpha(b_{0,r} - b_{2,r} z)^\alpha , \end{align*} $$

where $b_{0,r} = 3/2-r/\alpha $ and $b_{2,r} = 1/2-r/\alpha $ .

Lemma 5.8. Let $W_{2,r}(z) = [(1-z)(b_{0,r} - b_{2,r} z)]^\alpha $ and $ G_r(x) = W_{2,r}(e^{-\mathbf {i} x}) e^{\mathbf {i} rx} $ . Then,

(5.5) $$ \begin{align} \Re(G_r(x)) = R_r(x) \cos (Z_r(x) ) , \quad 0\le x\le \pi, \end{align} $$

where

(5.6) $$ \begin{align} R_r(x) &= (2\sin(x/2) S_r(x))^\alpha \ge 0, \nonumber \\S_r(x) &= \sqrt{b_{0,r}^2 - 2b_{0,r} b_{2,r}\cos(x) +b_{2,r}^2} \;, \nonumber \\Z_r(x) &= \alpha \theta_r(x) + \alpha \frac{\pi-x}{2} + rx \nonumber \\\text{ and } \qquad \theta_r(x) &= \tan^{-1} \bigg(\frac{ b_{2,r} \sin(x) }{b_{0,r} - b_{2,r} \cos(x) }\bigg),\quad -\frac{\pi}{2} \le \theta_r(x) \le \frac{\pi}{2}. \end{align} $$

Proof. Note that

$$ \begin{align*}b_{0,r} - b_{2,r} e^{-\mathbf{i} x} = (b_{0,r} - b_{2,r} \cos(x)) + \mathbf{i} \beta \sin(x) = S_r(x) e^{\mathbf{i} \theta_r(x)}, \end{align*} $$

where

$$ \begin{align*}S_r(x)^2 = b_{0,r}^2 - 2b_{0,r} b_{2,r}\cos(x) +b_{2,r}^2 \quad \text{ and} \quad \tan(\theta_r(x)) = \frac{ b_{2,r} \sin(x) }{b_{0,r} - b_{2,r} \cos(x) }.\end{align*} $$

Now, for the factor $(1-e^{-\mathbf {i}x})$ in $G_r(x) $ (with $b_{0,r} = b_{2,r} = 1$ ), we have the radial function $S(x)$ given by $ S(x)^2 = 2 - 2\cos (x) = 4 \sin ^2(x/2) $ and the angular function $ \phi (x) $ given by

$$ \begin{align*} \tan(\phi(x)) = \frac{\sin(x)}{1-\cos(x)} = \cot(x/2) = \tan((\pi-x)/2),\end{align*} $$

which yields $\phi (x) = {\pi -x}/{2} $ .

Therefore,

$$ \begin{align*} G_r(x) &=(1-e^{-\mathbf{i}x})^\alpha (b_{0,r}-b_{2,r} e^{-\mathbf{i}x})^\alpha e^{\mathbf{i}r x} \\ &=(2\sin(x/2))^\alpha e^{\alpha \mathbf{i} (\pi-x)/2} S_r(x)^\alpha e^{\alpha \mathbf{i} \theta_r(x) } e^{\mathbf{i}r x} \\ &=R_r(x) e^{\mathbf{i}( \alpha \theta_r(x) + \alpha (\pi-x)/{2} + r x )}\\ &=R_r(x) e^{\mathbf{i} Z_r(x)}. \end{align*} $$

Hence, we have the generating function in polar form

$$ \begin{align*} G_{r}(x)=G_{r}(\alpha,x)=R_r(x) \cos (Z_r(x)) + \mathbf{i} R_r(x) \sin (Z_r(x)), \end{align*} $$

where $R_r(x) = (2\sin (x/2))^\alpha S_r^\alpha (x) $ . The real part of $G_r(x)$ is then given by (5.5). The proof is now complete.

Proof.

1. (a) It can be easily observed from (3.3) that

   $$ \begin{align*} \tfrac{1}{12}\leq a_{2}(1,0)\leq \tfrac{1}{6} < \tfrac{1}{5} \quad \text{and}\quad \tfrac{1}{12}\leq a_{2}(1,-1) < \tfrac{1}{5}.\end{align*} $$

2. (b) Defining the central difference operator $\delta _{h}^{2}$ by

   $$ \begin{align*}\delta_{h}^{2}v_{i}=\frac{1}{h^{2}}(v_{i+1}-2v_{i}+v_{i-1})\quad \text{with} \quad u_{-1}=u_{N+1}=0,\end{align*} $$

we have

$$ \begin{align*} (\delta_{h}^{2}\textbf{v}, \textbf{u}) &=\frac{1}{h^{2}} \sum_{i=0}^{N}(v_{i+1}u_{i}-2v_{i}u_{i}+v_{i-1}u_{i}) \\ &=\frac{1}{h^{2}} \sum_{i=0}^{N}(v_{i}u_{i-1}-2v_{i}u_{i}+v_{i}u_{i+1}) = (\textbf{v}, \delta_{h}^{2} \textbf{u}). \end{align*} $$

Also,

$$ \begin{align*} (\delta_{h}^{2}\textbf{v}, \textbf{v}) &= \frac{1}{h^{2}} \sum_{i=0}^{N}(v_{i+1}v_{i}-2v_{i}v_{i}+v_{i-1}v_{i})\\ & \le \frac{1}{h^{2}} \sum_{i=0}^{N}\bigg( \frac{1}{2} (|v_{i+1}|^2 +|v_{i}|^2) + 2|v_{i}|^2 + \frac{1}{2} ( |v_{i-1}|^2+|v_{i}|^2) \bigg) = \frac{4}{h^2} \|\textbf{v}\|^2. \end{align*} $$

We then have $(P_{h}\textbf {v}, \textbf {v})=((I+h^{2}a_{2}\delta _{h}^{2})\textbf {v}, \textbf {v})=(\textbf {v}, \textbf {v})+h^{2}a_{2}(\textbf {v}, \delta _{h}^{2}\textbf {v})=(\textbf {v},P_{h}\textbf {v}).$ Thus, $P_h$ is self-adjoint. Further,

$$ \begin{align*} |(P_h \mathbf{v},\mathbf{v}) - (\mathbf{v},\mathbf{v}) | = h^2 a_2 |(\delta_h^2 \mathbf{v},\mathbf{v})| \le h^2 \frac{1}{5} \frac{4}{h^2}(\mathbf{v},\mathbf{v}) = \frac{4}{5} \|\mathbf{v} \|^2. \end{align*} $$

Hence, $({1}/{5})\|\mathbf {v} \|^2 \le (P_h \mathbf {v},\mathbf {v}) \le ({9}/{5})\|\mathbf {v} \|^2 $ and thus, $P_h$ is positive definite. This completes the proof.

5.1 Stability analysis

We establish that CN-QCA1 and CN-QCA2 approximations are unconditionally stable for $1 \leq \alpha \leq 2$ .

The approximation matrix of CN-QCA1 is $ A_{2,1,0}=\lambda _{1,(1,0)}A_{2,1}+\lambda _{2,(1,0)}A_{2,0}$ with corresponding generating function given by

$$ \begin{align*} G_{c,(1,0)}(\alpha, x)=\lambda_{1} W_{2,1}(e^{-\text{i}x})e^{\text{i} x}+\lambda_{0}W_{2.0}(e^{-\text{i}x}), \end{align*} $$

where the convex weights $\lambda _{1(1,0)}$ and $\lambda _{2,(1,0)}$ are given in (3.3). The real part of this generating function is

$$ \begin{align*} \Re(G_{c,(1,0)}(\alpha,x))=(2\sin(x/2))^{\alpha}[\lambda_{1}S_{1}^{\alpha}( \alpha,x)\cos(Z_{1}(\alpha,x))+ \lambda_{0}S_{0}^{\alpha}(\alpha,x)\cos(Z_{0}(\alpha, x))]. \end{align*} $$

Let $H_{c,(1,0)}(\alpha , x)=\lambda _{1}S_{1}^{\alpha }(\alpha ,x)\cos (Z_{1}(\alpha , x))+\lambda _{0}S_{0}^{\alpha }(\alpha , x)\cos (Z_{0}(\alpha , x))$ . For $\alpha =1$ , we see that the convex weights have the values $ \lambda _{1}=1, \lambda _{0}=0$ . Also, for $r=1$ , since $b_{0}=-1/2$ and $ b_{2}=1/2$ , we have from (5.6),

$$ \begin{align*}\tan(\theta_{1}(\alpha, x)) =-\frac{ \sin(x)}{1+\sin(x)}=\tan(-x/2)\quad \text{resulting in} \quad \theta_{1}(\alpha,x)=-x/2.\end{align*} $$

Thus, $Z_{1}(1, x)=-{x}/{2}+({\pi -x})/{2}+x={\pi }/{2}$ and so, $ \cos (Z_{1}(1,x))=0$ and $H_{c,(1,0)}(1,x)=0$ for all $x \in [0, \pi ]$ . Therefore, we obtain $\Re (G_{c, (1,0)}(1,x))=0\:\text {for all} x\in [0,\pi ].$

However, for $\alpha =2$ , we again have $\lambda _{1}=1, \lambda _{0}=0$ and for $r=1$ , we have $b_{2}=0$ . Therefore, $\theta _{1}(2, x)=0$ . Now, $ Z_{1}(2,x)=2({\pi -x}0/{2}+x=\pi $ and $H_{c,(1,0)}(2,x)=-S_{}^{2}(2,x)<0$ for al $x\in [0, \pi ]$ . Hence,

$$ \begin{align*} \Re(G_{c,(1,0)}(2.x))=(2\sin(x/2))^{2}H_{c,(1,0)}(2,x)<0 \quad \text{for all }\: x\in[0,\pi]. \end{align*} $$

Showing that $H_{c(1,0)}(\alpha ,x)$ is decreasing in the region $1 < \alpha < 2$ for any x is mathematically cumbersome. Therefore, we use the contour plot of $H_{c,(1,0)}(\alpha ,x)$ shown in Figure 1. It follows from the figure that $H_{c,(1,0)}(\alpha ,x)<0$ for all $x \in [0, \pi ]$ . Since $ \sin ^{\alpha }(x/2)>0$ for all $0<x\leq \pi $ , we see that $ \Re (G_{c,(1,)}(\alpha ,x))$ is negative. As the real part of the generating function $G_{c,(1,0)}(\alpha .x)$ is differentiable, it remains negative in the region (see also the surface plot in Figure 1). Thus, QCA1 is unconditionally stable for any $\alpha $ lying in the interval $1\le \alpha \le 2$ .

Figure 1 Contour and surface plots of $H_{c,(1,0)}$ .

The generating function corresponding to the approximation matrix $A_{2,1,-1}=\lambda _{1}A_{2,1}+ \lambda _{-1}A_{2,-1}$ of CN-QCA2 is

$$ \begin{align*} G_{c,(1,-1)}(\alpha, x)=\lambda_{1} W_{2,1}(e^{-\text{i}x})e^{\text{i} x}+\lambda_{-1}W_{2,-1}(e^{-\text{i}x})e^{-\text{i}x}, \end{align*} $$

where the weights $\lambda _{1}$ and $\lambda _{-1}$ are given in (3.4).

The real part of the generating function is then

$$ \begin{align*} \Re(G_{c,(1,-1)}(\alpha,x))&=(2\sin(x/2))^{\alpha}[\lambda_{1}S_{1}^{\alpha}( \alpha,x)\cos(Z_{1}(\alpha,x)) \nonumber\\ &\quad +\lambda_{-1}S_{-1}^{\alpha}(\alpha,x)\cos(Z_{-1}(\alpha, x))]. \end{align*} $$

Let $H_{c,(1,-1)}(\alpha , x)=\lambda _{1}S_{1}^{\alpha }(\alpha ,x)\cos (Z_{1}(\alpha , x))+\lambda _{-1}S_{-1}^{\alpha }(\alpha , x)\cos (Z_{-1}(\alpha , x))$ . Then, for $\alpha =1$ , we have $\lambda _{1}=1$ and $\lambda _{-1}=0$ . Hence, the same analysis for CN-QCA1 is applied by using the contour plot of $H_{c,(1, -1)}(\alpha , x)$ over the same rectangular region displayed in Figure 2. Hence, we conclude that $ H_{c,(1,-1)}(\alpha ,x)$ is also nonpositive in the region for all $ 1\leq \alpha \leq 2$ and $0\leq x\leq \pi $ . It follows that CN-QCA2 is unconditionally stable for all $1\leq \alpha \leq 2$ .

Figure 2 Contour and surface plots of $H_{c,(1,-1)}$ .

5.2 Convergence

Here, we analyse the convergence of the CN-QCA1 and CN-QCA2 methods. First, observe that as the diffusion constants $K_{1}$ and $K_{2}$ are nonnegative, and the operators $\Delta _{h-,2,r}^{\alpha }$ and $ \Delta _{h+,2,r}^{\alpha }$ are negative definite [Reference Nasir and Nafa14], the operator $L_{h}$ is negative definite.

Theorem 5.10. The numerical solutions governed by the CN-QCA1 and CN-QCA2 schemes with the given boundary and initial conditions converge for $1\leq \alpha \leq 2$ .

Proof. Let $\textbf {e}^{m}=\textbf {U}^{m}-\hat {\textbf {U}}^{m}$ be the error vector of the solution, where $\textbf {u}$ and $\hat {\textbf {U}}^{m}$ are respectively the exact and approximate solution vectors of the diffusion problem (4.1). Then, from (4.2), we obtain the following for the error component $e_{i}^{m}$ of $\textbf {e}^{m}=[e_{0}^{m}, e_{1}^{m}, \ldots , e_{N}^{m}]^{T}$ :

$$ \begin{align*}P_{h}\frac{e_{i}^{m+1}-e_{i}^{m}}{\tau}=L_{h}\frac{e_{i}^{m+1}+e_{i}^{m}}{2}+O(\tau^{2}+h^{4}))\end{align*} $$

for all $i=0,1,\ldots , N$ . Denoting the residue term $O(\tau ^{2}+h^{4})$ by $S_{i}^{m}$ , and noting that $\delta _{\tau }e_{i}^{m+1/2}=({e_{i}^{m+1}-e_{i}^{m}})/{\tau }$ and $e^{m+1/2}=({e_{i}^{m+1}+e_{i}^{m}})/{2}$ ,

$$ \begin{align*}P_{h}\delta_{\tau}{\mathbf{e}}^{m+1/2}=L_{h}\mathbf{e}^{m+1/2}+\mathbf{S}^{m},\quad \mathbf{e}_{0}^{m}=0,\ \mathbf{e}_{N}^{m}=0,\ \mathbf{e}_{i}^{0}=0,\end{align*} $$

where $\mathbf {S}^{m}$ is the residue vector given by $[S_{0}^{m}, S_{1}^{m}, \ldots , S_{N}^{m}]^{T}$ for all $i, =0, 1, \ldots , N, m=0, 1, \ldots , M$ .

We know that, in view of part (b) of Lemma 5.9, $P_{h}$ is self-adjoint and, in fact, $L_{h}$ is negative definite. Now, Theorem 5.1 implies that

$$ \begin{align*} \|\textbf{e}^{m}\| \leq 25 \tau \sum_{l=0}^{m-1}\|\textbf{S}^{l}\|\leq 25 \tau \sum_{l=0}^{M-1}\|\textbf{S}^{l}\| \leq 25 \tau M\sqrt{h}\sqrt{N}c(\tau^{2}+h^{4})=K\tau \sqrt{h}(\tau^{2}+h^{4}), \end{align*} $$

where $K_\tau $ is a positive constant, which establishes the convergence as h and $\tau $ tend to zero. The proof is now complete.

6 Numerical results

To assess the convergence and performance of the CN-QCA1 and CN-QCA2 methods introduced in Section 5, we consider the following one-dimensional fractional diffusion test problem (4.1).

Test problem: Let

$$ \begin{align*} & G(x, m,\alpha)=\frac{\Gamma(m+1)}{ \Gamma(m+1-\alpha)} (x^{m-\alpha} + (1-x)^{m-\alpha}), \\ &s_0(x,m)=x^m(1-x)^m, \quad \text{and constant diffusion coefficients } K_1=K_2= 1, \\ & f(x,t)=-e^{-t}\bigg( s_0(x,m)+ \sum_{j=0}^m (-1)^j \binom{m}{j} G(x,m+j,\alpha) \bigg),\\ &u(x,0)=s_0(x,m),\,u(0,t) = 0 \quad \text{and} \quad u(1,t)=0. \end{align*} $$

The exact solution to this problem is $u(x,t) = s_0(x) e^{-t}$ .

The computation of approximate solutions for the above problem was performed with various computational domains of different grid sizes with $ N $ space discretization points and $ M= N^2 $ time discretization points. The numerical results for $ m = 4,6,2$ and 3, are displayed in Tables 1–4, respectively. Clearly, for $ m = 4 $ , $ m = 6$ and $1\leq \alpha \leq 2$ , both approximations have the order of convergence four, which confirms the analysis of the methods given in Section 5.

Table 1 Convergence of CN-QCA1 for the Test problem with $m=4$ .

We have computed the order of convergence for $m = 2$ and $ m = 3$ , and the fourth-order methods have regularity issues near the boundaries and the order of convergence is reduced to the degree m of the exact solution $s_0(x)$ . This phenomenon has been observed by many researchers [Reference Baeumer, Mihaly and Meerchaert1, Reference Baeumer, Mihály and Sankaranarayanan2, Reference She16] and the way around this is under investigation.

Table 2 Convergence of CN-QCA2 for the Test problem with $m=6$ .

Table 3 Convergence of CN-QCA1 for the Test problem with $m=2$ .

We have also implemented the DTTS iteration method discussed in Section 4.1 as described in [Reference She16]. The number of iterations for $\sigma = 1$ is two for all cases, for a tolerance of $Tol = 10^{-7}$ , and the solutions are the same as those obtained by the LU-decomposition.

Table 4 Convergence of CN-QCA2 for the Test problem with $m=3$ .