1 Introduction
Difference approximations for fractional derivatives have gained popularity as shown in [Reference Lubich10, Reference Nasir and Nafa13, Reference Tian, Zhou and Deng18] and the references therein. The left
$(-)$
and right
$(+)$
Riemann–Liouville (R-L) fractional derivatives are defined, as in [Reference Podlubny15], by
$$ \begin{align*} D_{x\mp}^\alpha f(x) = \frac{(\pm 1)^n} {\Gamma(n-\alpha)} \frac{d^n}{d x^n} \int_{I_{x\mp}} (\pm (x-\xi))^{n-\alpha -1} f(\xi) \, d\xi, \end{align*} $$
where
$\Gamma $
denotes the Gamma function,
$I_{x-} = (-\infty ,x]$
,
$I_{x+} = [x,\infty )$
and the integer n is the ceiling of
$\alpha $
satisfying
$n-1 \le \alpha < n $
.
The corresponding Grünwald difference (GD) approximations for
$D_{x\mp }^\alpha f(x)$
are defined by
$$ \begin{align} \delta_{x\mp,h,r}^\alpha f(x) = \frac{1}{h^\alpha} \sum_{k=0}^{\infty} g_k^{(\alpha)} f(x \mp k h), \end{align} $$
where
$$ \begin{align*}g_k^{(\alpha)} = (-1)^k \frac{ \Gamma(\alpha+1)}{k!\Gamma{(\alpha+1-k)}}, \quad k\ge 0\end{align*} $$
are the coefficients of the power series expansion of the generating function
${W_1(z) = (1-z)^\alpha }$
and h is the subinterval size of the equi-spaced discretization
$I_{h,x\mp } = \{ x_k = x \mp kh \vert k=0,1,\cdots \} $
of the interval domain
$ I_{x\mp } $
.
The GD approximation is of first-order accuracy, but when applied to the fractional diffusion equation (FDE), displays unstable solutions even for implicit Euler and Crank–Nicolson (CN) numerical schemes, although these methods are unconditionally stable for classical integer-order diffusion equations. To overcome this difficulty, a shifted Grünwald difference approximation (SGD) was introduced by Meerschaert and Tadjeran [Reference Meerschaert and Tadjeran11], for the fractional derivative R-L, and is given by
$$ \begin{align} \delta_{x\mp,h,r}^\alpha f(x) = \frac{1}{h^\alpha} \sum_{k=0}^{\infty} g_k^{(\alpha)} f(x \mp (k-r) h). \end{align} $$
The SGD approximation in (1.2), with shift
$r=1$
, restores the stability of the methods mentioned above with the same first-order accuracy.
For higher order approximations, Meerchaert et al. [Reference Tadjeran, Meerschaert and Scheffler17] used an extrapolation improvement for space discretization to obtain second-order accuracy for solving the FDE. Nasir et al. [Reference Nasir, Gunawardana and Abeyrathna12] obtained second-order accuracy through a noninteger shift
$r=\alpha /2$
in (1.2), displaying super convergence. Convex combinations of various shifts in (1.2) were used to obtain some second-order approximations [Reference Tian, Zhou and Deng18]. These are called the weighted and shifted Grünwald difference (WSGD) approximations. A third-order approximation through WSGD was not successful as it fails to give desired stability for
$1< \alpha \le 2 $
[Reference Tian, Zhou and Deng18]. However, it is shown in [Reference Li and Deng9] that the unconditional stability for the third-order method is possible for a restricted range of the fractional order
$\alpha $
. A fourth-order approximation was obtained by Hao [Reference Hao, Sun and Cao8] through a quasi-compact difference approximation (QCA) on a WSGD approximation.
All the above higher-order approximations were obtained directly from the first-order Grünwald formula (1.1) only with the Grünwald generating function
$W_{1}(z)=(1-z)^{\alpha }$
with various shifts.
Earlier, Lubich [Reference Lubich10] derived some higher-order Grünwald type difference (GTD) approximations in terms of some generating functions of the form, for order p,
$$ \begin{align} W_p(z) = \bigg( \sum_{k = 1}^{p} \frac{1}{k}(1-z)^k \bigg)^\alpha , \quad 1\le p\le 6. \end{align} $$
In these cases, the coefficients
$g_k^{(\alpha )}$
in (1.1) are replaced by the coefficients of the power series expansions of
$W_p(z)$
to obtain order p approximations without shift. In fact, these higher-order generating functions were presented in [Reference Lubich10] to be used in a convolution quadrature formula for fractional integrals where
$ \alpha <0 $
.
Applying these generators for fractional derivatives, however, suffers from the same stability difficulties in solving the FDE. Shifted formulae with these generating functions reduce the approximation order to one and hence become uninteresting. Ding and Li [Reference Ding and Li4, Reference Ding and Li5] used a second-order Lubich-type approximation with shift
$r=1$
for the left and right fractional derivatives to approximate Riesz derivatives.
Nasir and Nafa [Reference Nasir and Nafa14] constructed a new set of Lubich-type generating functions for approximation orders that will retain the order with any shift r. These generating functions are also in the form of fractional-order power of a polynomial with coefficients dependent on
$\alpha $
and shift r. For the latest review of recent difference approximations for fractional derivatives, we refer to [Reference Cai and Li3].
Among these generating functions, the second-order generating function
$ W_{2,r}(z) $
is found to be reliable for
$1 \le \alpha \le 2 $
for stability and consistency. A theoretical justification and numerical verification for this approximation are given in [Reference Nasir and Nafa14]. Using this second-order generating function, Nasir and Nafa [Reference Nasir and Nafa13] obtained a third-order approximation using the QCA approach that turned out to be unconditionally stable for the CN scheme for the fractional diffusion equation. An explicit formulation for the generating function of the SGD for any order p and any arbitrary shift r is given by Guranthna et al. [Reference Gunarathna, Nasir and Daundasekera6].
Gunarathna et al. [Reference Gunarathna, Nasir and Daundasekara7] constructed order-4 approximations using both WSGD and QCA approaches and the second-order Lubich-type generating function
$W_{2,r}(z) $
. In this paper, we analyse the stability and convergence of these fourth-order methods for the FDE with constant diffusion coefficients combined with the CN scheme, and prove that they are unconditionally stable.
The paper is organized as follows. In Section 2, the previous work of constructing higher-order generating functions are summarized. In Section 3, the new set of fourth-order approximations is derived. In Section 4, the derived approximations are applied to space-fractional diffusion problems with the CN numerical scheme. Stability and convergence of the schemes are analysed in Section 5. Numerical results are presented in Section 6 and conclusions are drawn in Section 7.
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
$$ \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
$$ \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.
Definition 2.1. Let
$W(z)$
and
$\Delta _{h\mp ,r}f(x)$
be given as in (2.1) and (2.2), respectively.
-
(a)
$W(z)$
is said to approximate the left and right fractional derivatives
$D_{x\mp }^{\alpha } f(x) $
at x with shift r in the sense of Grünwald if
$$ \begin{align*} D_{x\mp}^{\alpha} f(x) = \lim_{h \rightarrow 0} \Delta_{h\mp , r}^\alpha f(x). \end{align*} $$
-
(b) The GTD in (2.2) is consistent of order
$p\ge 1$
with
$D_{x\mp }^{\alpha } f(x) $
if
$$ \begin{align*} \Delta_{h\mp,r}^\alpha f(x) - D_{x\mp}^{\alpha} f(x) = O(h^p). \end{align*} $$
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
$$ \begin{align} G_r(z) = 1 + \sum_{l=p}^{\infty} a_{l} z^l = 1 + O(z^p) \end{align} $$
and
$$ \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.
Corollary 2.3. If a generating function
$W_{p,r}(z)$
approximates the fractional differential operator
$D_{x\mp }^\alpha $
, then
$W_{p,r}(1) = 0 $
.
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
$$ \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
$$ \begin{align} &\lambda_1 +\lambda_2 = 1, \end{align} $$
$$ \begin{align} &\lambda_1 a_3(r_1) +\lambda_2 a_3(r_2) = 0, \end{align} $$
we get
$$ \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
$$ \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
$$ \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.
$$ \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
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
$$ \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]$
:
$$ \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:
$$ \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
$$ \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
$$ \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:
$$ \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,
$$ \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,
$$ \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.
Theorem 5.2. Let
$u \in L^1(\mathbb {R})$
. Then, the Crank–Nicolson scheme (4.3) using the approximation of the generating function
$W_{2,1}(z)$
for the SFD is unconditionally stable if and only if the matrix
$P_\alpha ^{-1} L_{\alpha }$
is negative definite.
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.
Lemma 5.6. If
$G(x)$
is the generating function for a Toeplitz matrix
$T = [t_{i-j}]$
, then
$G(x) e^{\mathbf {i}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^{\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$
.
Lemma 5.7. The generating functions of the matrices
$A_{2,r} $
and
$ A_{2,r}^T $
corresponding to the approximation operators
$\Delta _{h-,r,2} $
and
$\Delta _{h+,r,2} $
of the second-order approximation with shift r are given by
$W_{2,r}(e^{-\mathbf {i} x}) e^{\mathbf {i} rx} $
and its conjugate
$ W_{2,r}(e^{\mathbf {i} x}) e^{-\mathbf {i} rx} $
, respectively.
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,
$$ \begin{align} \Re(G_r(x)) = R_r(x) \cos (Z_r(x) ) , \quad 0\le x\le \pi, \end{align} $$
where
$$ \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.
Lemma 5.9. For the QCA1 and QCA2 operators of order 4 presented in Section 3, we have the following for
$1 \le \alpha \leq 2$
:
-
(a)
$\tfrac {1}{12} \leq a_{2}(r_{1}, r_{2})\leq \tfrac {1}{5}$
for
$(r_{1}, r_{2})=(1,0)$
and
$ (1,-1)$
; -
(b) the operator
$P_{h}(r_{1}, r_{2})$
is positive definite, self-adjoint and satisfies
$\tfrac {1}{5}\|\mathbf {v}\| \leq (P_h \mathbf {v},\mathbf {v}) \leq \tfrac {9}{5}\|\mathbf {v}\|$
.
Proof.
-
(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*} $$
-
(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$
.
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$
.
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.
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.
Convergence of CN-QCA2 for the Test problem with
$m=6$
.
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.
Convergence of CN-QCA2 for the Test problem with
$m=3$
.
7 Conclusion
Two fourth-order approximations for the Riemann–Liouville fractional derivatives are constructed using a second-order shifted Lubich-type approximation introduced by Nasir and Nafa [Reference Nasir and Nafa13]. These approximations are used with the Crank–Nicolson scheme to approximate the solution of a one-dimensional fractional diffusion problem. The results found for the test problem show that both schemes converge and the order of convergence is equal to four. Furthermore, the analysis given in this paper for the construction of higher-order methods has the mathematical merits of being involving, sound and interesting. Also, it sheds light on the use of this class of methods.
Acknowledgements
We are grateful to the editor and the referees for the valuable and constructive comments they have made. Their comments have improved the presentation of the paper.
Open access
