An efficient second-order convergent scheme for one-side space fractional diffusion equations with variable coefficients

TL;DR

This paper introduces a second-order finite difference scheme for one-sided space fractional diffusion equations with variable coefficients, achieving higher accuracy and efficiency through stability analysis and a Toeplitz preconditioner.

Contribution

The paper develops a second-order convergent scheme with an efficient Toeplitz preconditioner for variable coefficient fractional diffusion equations, improving upon existing methods.

Findings

Unconditionally stable with second-order convergence in time and space.

Preconditioned linear systems have bounded condition numbers, ensuring fast convergence.

Numerical results confirm the scheme's accuracy and efficiency.

Abstract

In this paper, a second order finite difference scheme is investigated for time-dependent one-side space fractional diffusion equations with variable coefficients. The existing schemes for the equation with variable coefficients have temporal convergence rate no better than second order and spatial convergence rate no better than first order, theoretically. In the presented scheme, the Crank-Nicolson temporal discretization and a second-order weighted-and-shifted Gr\"unwald-Letnikov spatial discretization are employed. Theoretically, the unconditional stability and the second-order convergence in time and space of the proposed scheme are established under some conditions on the diffusion coefficients. Moreover, a Toeplitz preconditioner is proposed for linear systems arising from the proposed scheme. The condition number of the preconditioned matrix is proven to be bounded by a constant independent of the discretization step-sizes so that the Krylov subspace solver for the preconditioned linear systems converges linearly. Numerical results are reported to show the convergence rate and the efficiency of the proposed scheme.

Tables

Table 1. Table 1: Numerical results for Example 1 when τ = 2 − 10 𝜏 superscript 2 10 \tau=2^{-10} .

		PGMRES-T				PLU
$\alpha$	$h$	$\mathrm{iter}$	CPU	$E(h,\tau )$	$Rat{e}_h$	CPU	$E(h,\tau )$	$Rat{e}_h$
1.2	$2^{-8}$	2.1	3.12s	3.47e-5	–	48.40s	3.49e-5	–
	$2^{-9}$	2.2	6.27s	8.42e-6	2.04	121.34s	8.60e-6	2.02
	$2^{-10}$	2.3	12.07s	1.87e-6	2.17	327.90s	1.99e-6	2.11
1.5	$2^{-8}$	3.4	4.12s	3.18e-5	–	48.79s	3.18e-5	–
	$2^{-9}$	3.6	7.89s	7.82e-6	2.03	123.67s	7.81e-6	2.02
	$2^{-10}$	3.9	15.31s	1.66e-6	2.23	350.71s	1.81e-6	2.11
1.8	$2^{-8}$	4.5	5.11s	2.51e-5	–	53.24s	2.50e-5	–
	$2^{-9}$	4.7	9.51s	6.23e-6	2.01	138.52s	6.12e-6	2.03
	$2^{-10}$	4.9	17.37s	1.57e-6	1.99	350.04s	1.40e-6	2.13

Table 2. Table 2: Numerical results for Example 1 when τ = 2 − 11 𝜏 superscript 2 11 \tau=2^{-11} .

		PGMRES-T				PLU
$\alpha$	$\tau$	$\mathrm{iter}$	CPU	$E(h,\tau )$	$Rat{e}_{\tau }$	CPU	$E(h,\tau )$	$Rat{e}_{\tau }$
1.2	$2^{-7}$	10.6	4.81s	2.01e-5	–	120.74s	2.01e-5	–
	$2^{-8}$	6.7	5.63s	4.80e-6	2.07	252.02s	4.74e-6	2.08
	$2^{-9}$	3.8	8.20s	1.16e-6	2.05	493.55s	9.68e-7	2.29
1.5	$2^{-7}$	12.4	4.47s	2.17e-5	–	124.21s	2.17e-5	–
	$2^{-8}$	9.5	6.30s	5.19e-6	2.06	252.72s	5.19e-6	2.06
	$2^{-9}$	6.5	9.27s	1.08e-6	2.27	486.69s	1.12e-6	2.22
1.8	$2^{-7}$	12.6	4.55s	2.34e-5	–	119.27s	2.34e-5	–
	$2^{-8}$	10.1	6.47s	5.73e-6	2.03	238.38s	5.68e-6	2.04
	$2^{-9}$	7.4	9.67s	1.45e-6	1.98	478.66s	1.28e-6	2.15

Table 3. Table 3: Numerical results for Example 2 when τ = 2 − 7 𝜏 superscript 2 7 \tau=2^{-7} .

		PGMRES-T				PLU
$(\alpha ,\beta )$	$h$	$\mathrm{iter}$	CPU	$E(h,\tau )$	$Rat{e}_h$	CPU	$E(h,\tau )$	$Rat{e}_h$
(1.01,1.09)	$2^{-4}$	4.0	6.63s	3.38e-3	–	35.73s	3.38e-3	–
	$2^{-5}$	4.1	9.80s	8.14e-4	2.05	234.97s	8.14e-4	2.05
	$2^{-6}$	4.1	16.19s	1.90e-4	2.10	5017.20s	1.90e-4	2.10
(1.5,1.3)	$2^{-4}$	4.4	7.43s	3.00e-3	–	34.20s	3.00e-3	–
	$2^{-5}$	4.6	10.62s	7.34e-4	2.03	219.68s	7.34e-4	2.03
	$2^{-6}$	5.8	20.95s	1.72e-4	2.09	4875.58s	1.72e-4	2.09
(1.5,1.6)	$2^{-4}$	4.9	7.59s	3.00e-3	–	34.00s	3.00e-3	–
	$2^{-5}$	5.8	12.71s	7.41e-4	2.03	220.50s	7.41e-4	2.03
	$2^{-6}$	6.8	23.99s	1.73e-4	2.09	4769.31s	1.73e-4	2.10
(1.5,1.9)	$2^{-4}$	5.8	8.75s	3.00e-3	–	33.81s	3.00e-3	–
	$2^{-5}$	6.7	14.47s	7.33e-4	2.03	219.61s	7.33e-4	2.03
	$2^{-6}$	7.4	25.85s	1.71e-4	2.10	4746.60s	1.71e-4	2.10
(1.2,1.2)	$2^{-4}$	4.1	6.63s	3.30e-3	–	33.50s	3.30e-3	–
	$2^{-5}$	4.2	9.71s	8.03e-4	2.03	217.33s	8.03e-4	2.03
	$2^{-6}$	4.6	17.13s	1.89e-4	2.09	4653.91s	1.89e-4	2.09
(1.5,1.5)	$2^{-4}$	4.5	7.13s	3.00e-3	–	33.71s	3.00e-3	–
	$2^{-5}$	5.4	12.09s	7.40e-4	2.03	218.36s	7.40e-4	2.03
	$2^{-6}$	6.5	23.10s	1.73e-4	2.10	4663.71s	1.73e-4	2.10
(1.8,1.8)	$2^{-4}$	6.2	9.22s	2.30e-3	–	33.73s	2.30e-3	–
	$2^{-5}$	7.2	15.40s	5.71e-3	2.04	218.98s	5.71e-4	2.04
	$2^{-6}$	8.3	28.54s	1.31e-4	2.12	4841.52s	1.31e-4	2.12

Table 4. Table 4: Numerical results for Example 2 when h = 2 − 6 ℎ superscript 2 6 h=2^{-6} .

		PGMRES-T				PLU
$(\alpha ,\beta )$	$\tau$	$\mathrm{iter}$	CPU	$E(h,\tau )$	$Rat{e}_h$	CPU	$E(h,\tau )$	$Rat{e}_h$
(1.01,1.09)	$2^{-3}$	219.1	62.13s	6.77e-3	–	334.49s	6.77e-3	–
	$2^{-4}$	82.9	43.48s	1.61e-3	2.07	636.48s	1.61e-3	2.07
	$2^{-5}$	33.1	30.14s	3.48e-4	2.21	1249.74s	3.48e-4	2.21
(1.5,1.3)	$2^{-3}$	14.5	3.50s	6.80e-3	–	338.66s	6.80e-3	–
	$2^{-4}$	12.7	5.71s	1.60e-3	2.07	612.05s	1.60e-3	2.07
	$2^{-5}$	10.2	9.00s	3.50e-4	2.21	1226.04s	3.50e-4	2.21
(1.5,1.6)	$2^{-3}$	12.3	3.01s	6.80e-3	–	331.95s	6.80e-3	–
	$2^{-4}$	11.3	5.16s	1.60e-3	2.07	644.72s	1.60e-3	2.07
	$2^{-5}$	10.1	8.96s	3.45e-4	2.23	1220.89s	3.45e-4	2.23
(1.5,1.9)	$2^{-3}$	11.0	2.73s	6.80e-3	–	336.93s	6.80e-3	–
	$2^{-4}$	10.3	4.77s	1.60e-3	2.07	638.47s	1.60e-3	2.07
	$2^{-5}$	9.5	8.42s	3.46e-4	2.23	1217.05s	3.46e-4	2.23
(1.2,1.2)	$2^{-3}$	15.4	3.63s	6.80e-3	–	338.42s	6.80e-3	–
	$2^{-4}$	12.7	5.70s	1.60e-3	2.07	629.27s	1.60e-3	2.07
	$2^{-5}$	10.2	9.01s	3.47e-4	2.22	1238.65s	3.47e-4	2.22
(1.5,1.5)	$2^{-3}$	13.4	3.23s	6.80e-3	–	334.47s	6.80e-3	–
	$2^{-4}$	12.1	5.47s	1.60e-3	2.07	646.24s	1.60e-3	2.07
	$2^{-5}$	10.5	9.19s	3.46e-4	2.22	1223.25s	3.46e-4	2.22
(1.8,1.8)	$2^{-3}$	12.3	3.00s	6.80e-3	–	335.93s	6.80e-3	–
	$2^{-4}$	11.5	5.23s	1.60e-3	2.06	625.35s	1.60e-3	2.06
	$2^{-5}$	10.6	9.30s	3.58e-4	2.20	1221.17s	3.58e-4	2.20

Equations

$$\begin{array}[]{l}\displaystyle{\frac{\partial u(x,t)}{\partial t}=d(x)\,{{}_{x_{L}}}D_{x}^{\alpha}u(x,t)+f(x,t)},\qquad x\in(x_{L},x_{R}),~{}t\in(0,T],\\ \vspace{2mm}u(x_{L},t)=u(x_{R},t)=0,\qquad t\in[0,T],\\ u(x,0)=\varphi(x),\qquad x\in[x_{L},x_{R}],\end{array}$$

$$\begin{array}[]{l}\vspace{2mm}\displaystyle{\,{{}_{x_{L}}}D_{x}^{\alpha}u(x,t)=\dfrac{1}{\Gamma(2-\alpha)}\dfrac{\partial^{2}}{\partial x^{2}}\int_{x_{L}}^{x}\dfrac{u(\xi,t)}{(x-\xi)^{\alpha-1}}\,d\xi},\end{array}$$

$$\begin{array}[]{l}\vspace{2mm}\displaystyle{{{}_{x_{L}}}D_{x}^{\alpha}u(x_{i})=\frac{1}{h^{\alpha}}\sum_{k=0}^{i}w_{k}^{(\alpha)}u(x_{i-k+1})+{\cal O}(h^{2})},\end{array}$$

$$w_0^{(\alpha )}=\frac{\alpha }{2}{g}_0^{(\alpha )},w_k^{(\alpha )}=\frac{\alpha }{2}{g}_k^{(\alpha )}+\frac{2-\alpha }{2}{g}_{k-1}^{(\alpha )}\text{mbox}fork\ge 1,$$

$$\displaystyle g_{0}^{(\alpha)}=1,\quad g_{k}^{(\alpha)}=\Big{(}1-\frac{\alpha+1}{k}\Big{)}g_{k-1}^{(\alpha)},~{}\mbox{ for }~{}k=1,2,\ldots.$$

$$\frac{u_i^n-u_i^{n-1}}{\tau }=\frac{1}{2h^{\alpha }}{d}_i\sum _{k=0}^i{w}_k^{(\alpha )}\left(u_{i-k+1}^{n-1}+u_{i-k+1}^n\right)+f_i^{n-\frac{1}{2}}+R_i^{n-\frac{1}{2}},1\le i\le M,$$

$$\displaystyle D=\mathrm{diag}(d_{1},d_{2},...,d_{M}),\quad G_{\alpha}=\left[\begin{array}[]{ccccc}w_{1}^{(\alpha)}&w_{0}^{(\alpha)}&0&\cdots&0\\ w_{2}^{(\alpha)}&w_{1}^{(\alpha)}&w_{0}^{(\alpha)}&\ddots&\vdots\\ \vdots&w_{2}^{(\alpha)}&w_{1}^{(\alpha)}&\ddots&0\\ \vdots&\ddots&\ddots&\ddots&w_{0}^{(\alpha)}\\ w_{M}^{(\alpha)}&\cdots&\cdots&w_{2}^{(\alpha)}&w_{1}^{(\alpha)}\\ \end{array}\right],$$

$$\frac{1}{\tau }\left(u^n-u^{n-1}\right)=\frac{1}{2h^{\alpha }}D{G}_{\alpha }\left(u^{n-1}+u^n\right)+f^{n-\frac{1}{2}},n=1,2,\dots ,N.$$

$${\lambda }_n{w}^{T}w\le w^{T}Hw\le {\lambda }_1{w}^{T}w.$$

$${\Vert u^n\Vert }_{D^{-1}}^2\le \exp (2T){\Vert \phi \Vert }_{D^{-1}}^2+[\exp (2T)-1]\underset{1\le k\le n}{\max }{\Vert f^{k-\frac{1}{2}}\Vert }_{D^{-1}}^2,n=1,2,\dots ,N,$$

$$\frac{1}{\tau }h{\left(u^{n-1}+u^n\right)}^T{D}^{-1}\left(u^n-u^{n-1}\right)=$$

$$+h{\left(u^{n-1}+u^n\right)}^T{D}^{-1}{f}^{n-\frac{1}{2}}.$$

$$\frac{1}{2h^{\alpha }}h{\left(u^{n-1}+u^n\right)}^T{G}_{\alpha }\left(u^{n-1}+u^n\right)$$

$$=\frac{1}{2h^{\alpha }}h{\left(u^{n-1}+u^n\right)}^T\mathcal{H}(G_{\alpha })\left(u^{n-1}+u^n\right)\le 0.$$

$$h(u^n{)}^T{D}^{-1}{u}^n-h(u^{n-1}{)}^T{D}^{-1}{u}^{n-1}$$

$$\le \tau h(u^n{)}^T{D}^{-1}{f}^{n-\frac{1}{2}}+\tau h(u^{n-1}{)}^T{D}^{-1}{f}^{n-\frac{1}{2}}.$$

$${\Vert u^n\Vert }_{D^{-1}}^2\le {\Vert u^{n-1}\Vert }_{D^{-1}}^2+\frac{\tau }{2}{\Vert u^n\Vert }_{D^{-1}}^2+\frac{\tau }{2}{\Vert u^{n-1}\Vert }_{D^{-1}}^2+\tau {\Vert f^{n-\frac{1}{2}}\Vert }_{D^{-1}}^2,$$

$${\Vert u^n\Vert }_{D^{-1}}^2\le \frac{2+\tau }{2-\tau }{\Vert u^{n-1}\Vert }_{D^{-1}}^2+\frac{2\tau }{2-\tau }{\Vert f^{n-\frac{1}{2}}\Vert }_{D^{-1}}^2.$$

$${\Vert u^n\Vert }_{D^{-1}}^2\le$$

$$+\frac{2\tau }{2-\tau }$$

$$\displaystyle\Big{(}\frac{2+\tau}{2-\tau}\Big{)}^{n}=\Big{(}1+\frac{2\tau}{2-\tau}\Big{)}^{n}\leq(1+2\tau)^{n}\leq\lim\limits_{N\rightarrow+\infty}\left(1+\frac{2T}{N}\right)^{N}=\exp(2T),$$

$$\displaystyle\frac{2\tau}{2-\tau}\sum_{k=1}^{n}\Big{(}\frac{2+\tau}{2-\tau}\Big{)}^{k-1}=\left(\frac{2+\tau}{2-\tau}\right)^{n}-1\leq\exp(2T)-1.$$

$$\Vert e^n\Vert \le c_2({\tau }^2+h^2),$$

$$\frac{1}{\tau }\left(e^n-e^{n-1}\right)=\frac{1}{2h^{\alpha }}D{G}_{\alpha }\left(e^{n-1}+e^n\right)+R^{n-\frac{1}{2}},1\le n\le N,$$

$$e_0^n=e_{M+1}^n=0,1\le n\le N,e_i^0=0,0\le i\le M+1.$$

$${\Vert e^n\Vert }_{D^{-1}}^2\le [\exp (2T)-1]\underset{1\le k\le n}{\max }{\Vert R^{k-\frac{1}{2}}\Vert }_{D^{-1}}^2,n=1,2,\dots ,N.$$

$$\Vert e^n{\Vert }^2\le [c_2({\tau }^2+h^2){]}^2,n=1,2,\dots ,N.$$

$${\mathbf{u}}^T{G}_{\alpha }\mathbf{v}=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }\sum _{k=1}^M{u}_k{e}^{-\mathbf{i}kx}\sum _{k=1}^M{v}_k{e}^{\mathbf{i}kx}g(\alpha ,x)dx.$$

$${\varsigma }_{\alpha }\triangleq \underset{x}{\min }\frac{\Re [-g(\alpha ,x)]}{|g(\alpha ,x)|}=|\cos (\frac{\alpha }{2}\pi )|,$$

$${\mathbf{u}}^T\tilde{D}G\tilde{D}\mathbf{u}\le 2\underset{i}{\max }\{|{\tilde{d}}_i{|}^2\}{\mathbf{u}}^{T}G\mathbf{u},$$

$$\left(\kappa -\frac{\sqrt{2}({\kappa }_{max}-{\kappa }_{min})}{{\varsigma }_{\alpha }}\right){\mathbf{u}}^{T}G\mathbf{u}$$

$$\le {\mathbf{u}}^T(-D{G}_{\alpha }-G_{\alpha }^{T}D)\mathbf{u}\le \left(\kappa +\frac{\sqrt{2}({\kappa }_{max}-{\kappa }_{min})}{{\varsigma }_{\alpha }}\right){\mathbf{u}}^{T}G\mathbf{u},$$

$${\mathbf{u}}^T(-D{G}_{\alpha }-G_{\alpha }^{T}D)\mathbf{u}=\kappa {\mathbf{u}}^{T}G\mathbf{u}+{\mathbf{u}}^T\left(-\tilde{D}{G}_{\alpha }-G_{\alpha }^T\tilde{D}\right)\mathbf{u}.$$

Full text

††journal: BIT

1 Introduction

In the paper, we study an efficient numerical method for solving the one-side space fractional diffusion equation (OSFDE) with variable coefficients. To begin with, we firstly present the one-dimensional OSFDE (the two dimension case will be discussed in Section 3) tadjeran2006205 ; sousa2011numerical ; sousaelic :

[TABLE]

where $d(x)$, which satisfies $0<{d_{-}}\leq d(x)\leq{d_{+}}<\infty$, is a strictly positive known function, $\varphi$ and $f$ are both known functions, $u$ is unknown to be solved, ${{}_{x_{L}}}D_{x}^{\alpha}u(x,t)$ is the Riemann-Liouville (RL) fractional derivative of order $\alpha\in(1,2)$ defined as (podlubny1999 ; samko1993fractional )

[TABLE]

with $\Gamma(\cdot)$ denoting the gamma function.

Due to the nonlocal dependence, fractional derivatives model many challenging phenomena more accurately than integer-order derivatives do, which has therefore attracted lots of interests in recent years. As an illustration of this fact, some applications of fractional calculus and anomalous diffusion have been discussed in the books baleanu2010new ; klages2008anomalous ; klafter2012fractional ; mainardi2010fractional ; sabatier2007advances ; ortigueira2011fractional and the references therein. It is well-known that closed-form analytic solutions of fractional diffusion equations are usually not available especially in the existence of variable coefficients. Moreover, because of the nonlocality of the fractional derivative and the existence of the variable coefficients, the discretization of the OSFDE tends to generate dense matrix with high displacement rank 111see mng2004 , for which the discrete linear systems related to the variable-coefficients OSFDE are time-consuming to directly solve. Therefore, studying reliable discretization schemes and the corresponding fast iterative solvers for the OSFDE becomes an urgent topic.

There have been many schemes applicable to or solely developed for OSFDEs; see, e.g., tadjeran2006205 ; meerschaert2004finite ; tadjeran2007813 ; tianwy2015 ; chendeng2014 ; sousa2011numerical ; sousaelic ; qu2014note ; vong2016high ; zhaozsunwcao2015 ; leislhyc2016 . In zhaozsunwcao2015 ; leislhyc2016 ; vong2016high , numerical schemes with spatial fourth-order convergence for a space fractional diffusion equation is developed by applying the technique of compact operators, which is however only available for constant-coefficient case. Another spatially fourth-order accurate scheme is studied in chendeng2014 by implementing weighted-and-shifted Lubich difference operators whose convergent property is established only for constant diffusion coefficients. Some second-order numerical schemes are proposed in tadjeran2006205 ; tadjeran2007813 ; sousa2011numerical ; sousaelic for solving OSFDEs with variable coefficients, which however does not provide convergence proof. In qu2014note , the stability and convergence of the second-order numerical scheme for variable coefficient equations are established for $\alpha\in(\alpha_{0},2)$, where $\alpha_{0}\approx 1.5546$ is a solution of the equation $3^{3-\gamma}-4\times 2^{3-\gamma}+6=0$. In linstbcvg2017 , a series of numerical schemes for Riesz space fractional diffusion equation have been proven to be convergent and stable. Nevertheless, the proof technique used in linstbcvg2017 heavily depends on the symmetry of discretization matrix of Riesz fractional derivative, which is not applicable to the OSFDE that involves the non-symmetric one-sided fractional derivatives weighted by variable coefficients.

In this paper, we propose a second-order scheme for the one- and two-dimensional OSFDEs weighted by variable coefficients. The Crank-Nicolson method and a second-order weighted-and-shifted Grünwald-Letnikov difference (WSGD, see tianwy2015 ) operator are employed to discretize temporal and spatial derivatives, respectively. For the one-dimensional OSFDE, the proposed scheme is proven to be unconditionally stable and second-order convergent in time and space without additional assumption on the diffusion coefficient. It has been shown in vonglyu2017 that the symmetric part of the discretization matrix of $d(x){{}_{x_{L}}}D_{x}^{\alpha}$ is negative semi-definite under some conditions on $d(x)$. We extend this one-dimension result to two-dimension case under additional assumptions on the diffusion coefficients 222see the assumptions in Lemma 2, based on which the proposed scheme for the two-dimension OSFDE is proven to be unconditionally stable and second-order convergent in time and space.

As mentioned above, the direct solver for the linear systems arising from variable-coefficients OSFDE requires too much computational time. Fortunately, the discretization matrix has Toeplitz-like structure due to which its matrix-vector multiplication can be fast computed via fast Fourier transforms (FFTs). Because of the fast matrix-vector multiplication, fast iterative solvers for the linear systems can be possibly developed. However, the discretization matrix of the OSFDE is ill-conditioned when $\tau/h^{\alpha}$ is large, where $\tau$ and $h$ represent the temporal and spatial step-sizes, respectively. Thus, a Toeplitz preconditioner is proposed to reduce the condition number of the one- and two-dimensional discretization matrices. Theoretically, we show that the condition number of the preconditioned matrix is uniformly bounded by a constant independent of $\tau$ and $h$ under certain conditions on the diffusion coefficients 333see the assumptions in Theorems 2.4, 3.3 so that the Krylov subspace method for the preconditioned linear systems converges linearly no matter the unpreconditioned matrix is ill-conditioned or not.

This paper is organized as follows. In Section 2, we propose a second-order scheme and its corresponding Toeplitz preconditioner for the one-dimensional

OSFDE, analyze the unconditional stability and convergence of the proposed scheme, estimate the condition number of the preconditioned matrix. In Section 3, we extend the scheme and the preconditioner to two-dimensional case. In Section 4, numerical results are reported to show the efficiency and accuracy of the proposed scheme.

2 Stability and Convergence of Discrete One-Dimensional OSFDE and Its Preconditioning

We need some notations to describe the discretization for (1). Let $h=(x_{R}-x_{L})/(M+1)$ and $\tau=T/N$ be the space and time step sizes, respectively, where $M$ and $N$ are given positive integers. And denote $x_{i}=x_{L}+ih$ for $i=0,1,\ldots,M+1$, $t_{n}=n\tau$ for $n=0,1,\ldots,N$. Throughout this paper, the discretization on RL fractional derivative is based on the following second-order WSGD formula tianwy2015 :

[TABLE]

which is under the smooth assumptions $u$, ${}_{-\infty}D_{x}^{\alpha+2}u$ and Fourier transform of

${}_{-\infty}D_{x}^{\alpha+2}u$ belong to $L^{1}(\mathbb{R})$. The coefficients $w_{k}^{(\alpha)}$ were defined by (tianwy2015 )

[TABLE]

where $g_{k}^{(\alpha)}$ are the coefficients of the power series of $(1-z)^{\alpha}$, and they can be obtained recursively as

[TABLE]

Next, we introduce the finite difference scheme for solving (1). Let $u^{n}_{i}$ be the numerical approximation of $u(x_{i},t_{n})$, and denote $d_{i}=d(x_{i})$, $\varphi_{i}=\varphi(x_{i})$, $f^{n-\frac{1}{2}}_{i}=f(x_{i},t_{n-\frac{1}{2}})$, where $t_{n-\frac{1}{2}}=(t_{n-1}+t_{n})/2$ and $n=1,2,\ldots,N$. Then, applying the Crank-Nicolson technique and approximation (3) to the time derivative and the space fractional derivatives of (1) respectively, we get

[TABLE]

where $|R_{i}^{n-\frac{1}{2}}|\leq c_{1}(\tau^{2}+h^{2})$ for a positive constant $c_{1}$; see, e.g., tianwy2015 .

Denote $u^{n}=[u_{1}^{n},u_{2}^{n},\ldots,u_{M}^{n}]^{T}$, $f^{n-\frac{1}{2}}=[f_{1}^{n-\frac{1}{2}},f_{2}^{n-\frac{1}{2}},\ldots,f_{M}^{n-\frac{1}{2}}]^{T}$, and

[TABLE]

where $\{w_{k}^{(\alpha)}\}_{k=0}^{M}$ are the coefficients given in (4).

Omitting the small term $R_{i}^{n+\frac{1}{2}}$ in (5), then equation (1) can be solved numerically by the following finite difference scheme in matrix form

[TABLE]

2.1 Stability and convergence

Some General Notations:

$\bullet$ $\mathbb{C}^{m\times n}$ ($\mathbb{R}^{m\times n}$, respectively) denotes the set of all $m\times n$ complex (real, respectively) matrices.

$\bullet$ ${\cal H}(X)$ denotes the symmetric part of a square matrix $X$.

Lemma 1

(tianwy2015 ) The matrix $G_{\alpha}+G_{\alpha}^{T}$ is negative definite.

Lemma 2

(laub2005matrix ) Let symmetric matrix $H\in\mathbb{R}^{n\times n}$ with eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}$. Then for all $w\in\mathbb{R}^{n\times 1}$,

[TABLE]

Now we show the stability and convergence of scheme (12) by energy method.

Theorem 2.1

The finite difference scheme (12) is unconditionally stable and its solution satisfies the following estimate

[TABLE]

where $\|\cdot\|_{D^{-1}}$ is the norm induced by the inner product, $\langle v_{1},v_{2}\rangle_{D^{-1}}:=hv_{1}^{\rm T}D^{-1}v_{2}$.

Proof

Some steps of this proof are similar to those of Theorem 3.8 in vonglyu2017 . Multiplying

$h\left(u^{n-1}+u^{n}\right)^{T}D^{-1}$ on the both sides of (12), we get

[TABLE]

Notice that $w^{T}G_{\alpha}w=w^{T}{\cal H}(G_{\alpha})w$ for any real vector $w$. Therefore, by Lemma 1, the first term on the right hand side of (13) can be estimated as

[TABLE]

As a result

[TABLE]

Applying Cauchy-Schwarz inequality on the right hand side of (2.1), we get

[TABLE]

which is equivalent to

[TABLE]

Iterating (15) for $n$ times, we obtain

[TABLE]

For the small $\tau~{}(\tau\leq 1)$, we have

[TABLE]

and

[TABLE]

The result follows from (16)–(18).

Theorem 2.2

Let $u(x_{i},t_{n})$ be the exact solution of (1) and $u_{i}^{n}$ be the solution of finite difference scheme (12). Denote $e_{i}^{n}=u(x_{i},t_{n})-u_{i}^{n}$, $0\leq i\leq M+1$, $0\leq n\leq N$. Then there exists a positive constant $c_{2}$ such that

[TABLE]

where $e^{n}=[e_{1}^{n},e_{2}^{n},\ldots,e_{M}^{n}]^{T}$ and $\|\cdot\|$ denotes the discrete $L^{2}$ norm, i.e. $\|v\|=\sqrt{hv^{T}v}$,

Proof

Denote $R^{n-\frac{1}{2}}=[R_{1}^{n-\frac{1}{2}},R_{2}^{n-\frac{1}{2}},\ldots,R_{M}^{n-\frac{1}{2}}]^{T}$. We can easily show that $e^{n}$ and $e_{i}^{n}$ satisfy the following error equations

[TABLE]

By Theorem 2.1, we have

[TABLE]

As $D^{-1}$ is a positive diagonal matrix, utilizing Lemma 2, we get

[TABLE]

2.2 An estimate on the field of values of $DG_{\alpha}+G_{\alpha}^{T}D$

In this subsection, we focus on estimating the field of values of $DG_{\alpha}+G_{\alpha}^{T}D$, the results of which will be further applied to the analysis of one-dimensional preconditioning and the extension to two-dimensional OSFDE. First, we denote $g(\alpha,x)$ as the generating function mng2004 of the Toeplitz matrix $G_{\alpha}$. The next two lemmas describe some properties concerning $g(\alpha,x)$, which will be useful to obtaining the desired estimation.

Lemma 3

(mng2004 ) Let ${\bf u}=[u_{1},u_{2},\ldots,u_{M}]^{T},{\bf v}=[v_{1},v_{2},\ldots,v_{M}]^{T}\in\mathbb{R}^{M\times 1}$. Then we have

[TABLE]

Lemma 4

(vonglyu2017 ) It holds that

[TABLE]

where $\Re[g(\alpha,x)]$ denotes the real part of $g(\alpha,x)$.

The following lemma provides a novel bound to the field of values of ${\tilde{D}}G{\tilde{D}}$, where $G=-G_{\alpha}-G_{\alpha}^{T}$ and ${\tilde{D}}$ is a diagonal matrix satisfying some properties.

Lemma 5

*(vonglyu2017 ) Denote $G=-G_{\alpha}-G_{\alpha}^{T}$. Suppose that

${\tilde{D}}=\mathrm{diag}(\tilde{d}(x_{1}),\tilde{d}(x_{2}),...,\tilde{d}(x_{M}))$ for some function $\tilde{d}(x)$ defined on $(x_{L},x_{R})$. For any real vector ${\bf u}=[u_{1},u_{2},\ldots,u_{M}]^{T}$, we have*

[TABLE]

if ${\tilde{d}}(x)$ is convex and ${\tilde{d}}(x)\geq 0$, or ${\tilde{d}}(x)$ is concave and ${\tilde{d}}(x)\leq 0$.

Assuming $0\leq\kappa_{min}\leq d(x)\leq\kappa_{max}<\infty$. The following theorem reveals some inclusion relations between numerical ranges of $G$ and $-DG_{\alpha}-G_{\alpha}^{\rm T}D$, which acts an important role in the analysis of the proposed preconditioner.

Theorem 2.3

For any ${\bf u}=[u_{1},u_{2},\ldots,u_{M}]^{T}$, we have

[TABLE]

where $\kappa=\kappa_{max}$ when $d(x)$ is concave, and $\kappa=\kappa_{min}$ when $d(x)$ is convex.

Proof

Denote ${\tilde{D}}=D-\kappa I$, then $DG_{\alpha}+G_{\alpha}^{T}D=\kappa(G_{\alpha}+G_{\alpha}^{T})+{\tilde{D}}G_{\alpha}+G_{\alpha}^{T}{\tilde{D}}$. And, for any ${\bf u}=[u_{1},u_{2},\ldots,u_{M}]^{T}$, we have

[TABLE]

Denote $u(x)=\sum_{k=1}^{M}{u}_{k}e^{{\bf i}kx}$ and $v(x)=\sum_{k=1}^{M}{({\tilde{D}}u)}_{k}e^{{\bf i}kx}$, it follows by Lemmas 3 and 5 that

[TABLE]

Using Lemma 3 again and applying Cauchy-Schwarz inequality, Lemma 4, (21) and (2.2), we get

[TABLE]

Thus, the desired result can be obtained just by utilizing the following inequality

[TABLE]

2.3 Toeplitz preconditioner for the discrete one-dimensional fractional diffusion equation

To solve (12) is equivalent to recursively solve the following linear systems

[TABLE]

where ${\bf A}={I}_{M}-\eta{D}{G}_{\alpha}$, $\eta=\tau/(2h^{\alpha})$, ${I}_{k}$ denotes $k\times k$ identity matrix, ${b}^{n}=({I}_{M}+\eta{D}{G}_{\alpha}){u}^{n-1}+\tau{\bf f}^{n-\frac{1}{2}}$. As explained in the introduction section, a good preconditioner is required for the linear systems in (23).

For any $m\times m$ diagonal matrix, ${\bf C}=\mathrm{diag}(c_{1},c_{2},...,c_{m})$, denote ${\tt mean}({\bf C})=\frac{1}{m}\sum\limits_{i=1}^{m}c_{i}$. In this subsection, we propose a Toeplitz preconditioner for the linear systems in (23) such that

[TABLE]

where $\bar{d}={\tt mean}({D})$. In the following, we discuss a computationally effective representation of ${\bf P}^{-1}$, which allows fast matrix-vector multiplication of ${\bf P}^{-1}$.

Let ${\bf v}=(v_{1},v_{2},...,v_{M})^{\rm T}$ and $\tilde{{\bf v}}=(\tilde{v}_{1},\tilde{v}_{2},...,\tilde{v}_{M})^{\rm T}$ be solutions of following linear systems

[TABLE]

According to the Gohberg-Semencul-type formula gohbergolshevsky , ${\bf P}^{-1}$ can be expressed as follows

[TABLE]

where ${\bf S}_{1}$, ${\bf S}_{2}$ are skew-circulant matrices with ${\bf v}$, $\bar{{\bf v}}=(-\tilde{v}_{M},\tilde{v}_{1},...,\tilde{v}_{M-1})^{\rm T}$ as their first columns, respectively; ${\bf C}_{1}$, ${\bf C}_{2}$ are circulant matrices with

$\hat{{\bf v}}=(\tilde{v}_{M},\tilde{v}_{1},...,\tilde{v}_{M-1})^{\rm T}$, ${\bf v}$ as their first columns, respectively. From (25), we see that $v_{1}$ is the first diagonal entry of ${\bf P}^{-1}$. From Lemma 1, we see that ${\bf P}+{\bf P}^{\rm T}$ is positive definite. Thus,

[TABLE]

which means (26) is applicable. Moreover, the Toeplitz linear systems in (25) can be efficiently solved by the super fast direct solver proposed in hoog1987 .

For ${\bf C}\in\mathbb{C}^{m\times n}$, denote by $\Sigma({\bf C})$, the set of singular values of ${\bf C}$. Also denote $\Sigma^{2}({\bf C})=\{x^{2}|x\in\Sigma({\bf C})\}$. For any matrix ${\bf C}\in\mathbb{C}^{m\times m}$, denote by $\sigma({\bf C})$, the spectrum of ${\bf C}$. For a number $\lambda$, denote by $\Re(\lambda)$, the real part of $\lambda$.

For any invertible matrix ${\bf C}\in\mathbb{C}^{m\times m}$, define its condition number as

[TABLE]

Lemma 6

(see (tianwy2015, , Lemma 2.7))* For any ${\bf C}\in\mathbb{C}^{m\times m}$, it holds*

[TABLE]

As a preconditioner, the invertibility is essential.

Proposition 1

${\bf P}$* is invertible for any $\alpha\in(1,2)$.*

Proof

By lemme 1, it is easy to see that ${\bf P}+{\bf P}^{\rm T}$ is positive definite and thus has positive eigenvalues. From Lemma 6, we see that $\{\Re(\lambda)|\lambda\in\sigma({\bf P})\}\subset(0,+\infty)$. Therefore, ${\bf P}$ is invertible.

For any Hermitian matrices ${\bf H}_{1},{\bf H}_{2}\in\mathbb{C}^{m\times m}$, denote ${\bf H}_{1}\prec{\bf H}_{2}$ or ${\bf H}_{2}\succ{\bf H}_{1}$ if ${\bf H}_{2}-{\bf H}_{1}$ is Hermitian positive definite. Especially, we denote ${\bf O}\prec{\bf H}_{1}$ or ${\bf H}_{1}\succ{\bf O}$, when ${\bf H}_{1}$ itself is Hermitian positive definite. Also, we use ${\bf H}_{1}\preceq{\bf H}_{2}$ or ${\bf H}_{2}\succeq{\bf H}_{1}$ to denote a Hermitian positive semidefinite ${\bf H}_{2}-{\bf H}_{1}$ and use ${\bf O}\preceq{\bf H}_{1}$ or ${\bf H}_{1}\succeq{\bf O}$ to denote a Hermitian positive semidefinite ${\bf H}_{1}$.

Next, we are to estimate the condition number of the preconditioned matrix ${\bf A}{\bf P}^{-1}$.

Proposition 2

For positive numbers $\xi_{i}$, $\zeta_{i}$ ($1\leq i\leq m$), it obviously holds that

[TABLE]

Theorem 2.4

Assume

(i)

for any $x\in(x_{L},x_{R})$, $d(x)\in[\kappa_{\min},\kappa_{\max}]$ for positive constants $\kappa_{\min}$ and $\kappa_{\max}$,

(ii)

$\kappa_{\max}-\nu_{\alpha}>0$, with $\nu_{\alpha}=\sqrt{2}(\kappa_{\max}-\kappa_{\min})/\varsigma_{\alpha}$,

(iii)

$d(x)$* is concave.*

Then, for any $N\geq 1$, any $M\geq 1$, $\Sigma^{2}({\bf A}{\bf P}^{-1})\subset[\check{s},\hat{s}]$ and thus

[TABLE]

where $\check{s}$ and $\hat{s}$ are positive constants independent of $\tau$, $h$ and given by

[TABLE]

Proof

By straightforward calculation,

[TABLE]

By Theorem 2.3, we see that

[TABLE]

For any non-zero vector ${\bf y}\in\mathbb{R}^{M\times 1}$, denote ${\bf z}={\bf P}^{-1}{\bf y}$. Then, it holds

[TABLE]

By (27),

[TABLE]

By Proposition (2) and (2.3),

[TABLE]

During the proof above, there is no constraint on $M$ and $N$. Thus, for any $N\geq 1$, any $M\geq 1$, $\Sigma^{2}({\bf A}{\bf P}^{-1})\subset[\check{s},\hat{s}]$ and $\sup\limits_{N,M\geq 1}\rm cond({\bf A}{\bf P}^{-1})\leq\sqrt{\hat{s}/\check{s}}$.

Similar to proof of Theorem 2.4, one can prove following theorem.

Theorem 2.5

Assume

(i)

for any $x\in(x_{L},x_{R})$, $d(x)\in[\kappa_{\min},\kappa_{\max}]$ for positive constants $\kappa_{\min}$ and $\kappa_{\max}$,

(ii)

$\kappa_{\min}-\nu_{\alpha}>0$, with $\nu_{\alpha}=\sqrt{2}(\kappa_{\max}-\kappa_{\min})/\varsigma_{\alpha}$,

(iii)

$d(x)$* is convex.*

Then, for any $N\geq 1$, any $M\geq 1$, $\Sigma^{2}({\bf A}{\bf P}^{-1})\subset[\check{s},\hat{s}]$ and thus

[TABLE]

where $\check{s}$ and $\hat{s}$ are positive constants independent of $\tau$, $h$ and given by

[TABLE]

Remark 1

Theorems 2.4–2.5 show that $\rm cond({\bf A}{\bf P}^{-1})$ has an upper bound independent of $\tau$ and $h$ under certain assumptions on the coefficient function $d$. Thus, Krylov subspace method for such preconditioned linear systems converges linearly and independently on the discretization step-sizes.

3 Extension to Two-dimensional OSFDE

In this section, we study the following two-dimensional OSFDE tadjeran2007813 :

[TABLE]

where $\alpha,\beta\in(1,2)$, $\Omega=(x_{L},x_{R})\times(y_{L},y_{R})$ and $d(x,y)$, $e(x,y)$ are nonnegative functions, ${{}_{x_{L}}}D_{x}^{\alpha}u(x,y,t)$ denotes the $\alpha$-order RL derivative with respect to $x$ direction defined as

[TABLE]

${{}_{y_{L}}}D_{y}^{\beta}u(x,y,t)$ can be defined in a similar way.

To state a finite difference scheme for (29), we need more notations. Let $\tau=T/N$, $h_{1}=\frac{x_{R}-x_{L}}{M_{1}+1}$, $h_{2}=\frac{y_{R}-y_{L}}{M_{2}+1}$, where $M_{1}$, $M_{2}$ and $N$ are some positive integers. For $i=0,1,\ldots,M_{1}+1$, $j=0,1,\ldots,M_{2}+1$ and $n=0,1,\ldots,N$, denote $x_{i}=ih_{1}$, $y_{j}=jh_{2}$ and $t_{n}=n\tau$. Denote $t_{n-\frac{1}{2}}=\frac{t_{n}+t_{n-1}}{2}$ for $n=1,2,...,N$. Let $\bar{\Omega}_{h}=\{(x_{i},y_{j})|0\leq i\leq M_{1}+1,0\leq j\leq M_{2}+1\}$, ${\Omega}_{h}={\bar{\Omega}_{h}}\cap{\Omega}$, $\partial{{\Omega}_{h}}={\bar{\Omega}_{h}}\cap{\partial{\Omega}}$. Furthermore, denote $d_{i,j}=d(x_{i},y_{j})$, $e_{i,j}=e(x_{i},y_{j})$, $f_{i,j}^{n-\frac{1}{2}}=f(x_{i},y_{j},t_{n-\frac{1}{2}})$, and $\varphi_{i,j}=\varphi(x_{i},y_{j})$, and let $u_{i,j}^{n}$ be the numerical approximation of $u(x_{i},y_{j},t_{n})$. Then, in a similar way with the one-dimensional case, we can derive the Crank-Nicolson scheme for the two-dimensional problem (29) as the following

[TABLE]

where ${\hat{R}}_{i,j}^{n-\frac{1}{2}}\leq c_{3}(\tau^{2}+h_{1}^{2}+h_{2}^{2})$ for a positive constant $c_{3}$.

Take

[TABLE]

Omitting the small term ${\hat{R}}_{i,j}^{n-\frac{1}{2}}$ in (30), the finite difference scheme in matrix form for (29) can be given as:

[TABLE]

where $I$ is the identity matrix, the symbol ‘$\otimes$’ denotes the Kronecker product, and $G_{\beta}$ has the similar definition to $G_{\alpha}$.

3.1 Stability and Convergence of the Two-dimensional Problem

To discuss the stability and convergence of scheme (3), we denote

[TABLE]

and introduce a set:

[TABLE]

Now we present the stability of the scheme (3).

Theorem 3.1

For any $Q\in{\mathfrak{D}}$, the finite difference scheme (3) is unconditionally stable and its solution satisfies the following estimate

[TABLE]

where $\|\cdot\|_{Q}$ is defined as $\|v\|_{Q}^{2}:=hv^{T}Qv$.

Proof

Multiplying $h\left(u^{n-1}+u^{n}\right)^{T}Q$ on the both sides of (3), we get

[TABLE]

Since ${\cal H}(QA)$ is negative semi-definite, we have

[TABLE]

Then it follows

[TABLE]

The rest of the proof is similar to that in Theorem 2.1.

With Theorem 3.1, the convergence of scheme (3) can be directly obtained:

Theorem 3.2

Let $u(x_{i},y_{j},t_{n})$ be the exact solution of (29) and smooth enough, $u_{i,j}^{n}$ be the solution of finite difference scheme (3). Denote $e_{i,j}^{n}=u(x_{i},y_{j},t_{n})-u_{i,j}^{n}$, $0\leq i\leq M_{1}+1$, $0\leq j\leq M_{2}+1$, $0\leq n\leq N$. For any $Q\in{\mathfrak{D}}$, there exists a positive constant $c_{4}$ such that

[TABLE]

The remaining and important thing is to give the feature of the set ${\mathfrak{D}}$. However, it seems difficult to depict all the elements of ${\mathfrak{D}}$. In the following Corollaries 1 and 2, we show that there are some matrices belong to ${\mathfrak{D}}$ when the variable coefficients $d(x,y),e(x,y)$ satisfy some certain conditions, this ensures that ${\mathfrak{D}}$ is not an empty set which is necessary for the stability and convergence. We discuss the existence of those matrices in two cases:

$\bullet$ Case 1 When $d(x,y),e(x,y)$ are separable respect to $x$ and $y$.

In this case, we denote $d(x,y)={\tilde{d}}(x){\hat{d}}(y)$ and $e(x,y)={\tilde{e}}(x){\hat{e}}(y)$, and take

[TABLE]

Then $D={\hat{D}}\otimes{\tilde{D}}$, $E={\hat{E}}\otimes{\tilde{E}}$.

Corollary 1

If ${\tilde{d}}_{-}\leq{\tilde{d}}(x)\leq{\tilde{d}}_{+}$ and ${\hat{e}}_{-}\leq{\hat{e}}(y)\leq{\hat{e}}_{+}$ for some positive constants ${\tilde{d}}_{-},{\tilde{d}}_{+},{\hat{e}}_{-}$ and ${\hat{e}}_{+}$, then ${\hat{E}}^{-1}\otimes{\tilde{D}}^{-1}\in{\mathfrak{D}}$.

Proof

We have $A=\frac{1}{2h_{1}^{\alpha}}({\hat{D}}\otimes{\tilde{D}}G_{\alpha})+\frac{1}{2h_{2}^{\beta}}({\hat{E}}G_{\beta}\otimes{\tilde{E}})$, then

[TABLE]

which is negative semi-definite. Thus ${\hat{E}}^{-1}\otimes{\tilde{D}}^{-1}\in{\mathfrak{D}}$.

$\bullet$ Case 2 When $d(x,y)$ and $e(x,y)$ are non-separable.

As in Lemma 4, we denote

[TABLE]

where $g(\beta,x)$ is the generating function of matrix $G_{\beta}$.

Corollary 2

(i)

Assume that

[TABLE]

[TABLE]

Then, $I\in{\mathfrak{D}}$ if the following conditions are fulfilled

[TABLE]

where $\kappa^{d}(y)=\kappa_{max}^{d}(y)$ when $d(x,y)$ is a concave function of $x$, $\kappa^{d}(y)=\kappa_{min}^{d}(y)$ when $d(x,y)$ is a convex function of $x$, $\kappa^{e}(x)=\kappa_{max}^{e}(x)$ when $e(x,y)$ is a concave function of $y$, $\kappa^{e}(x)=\kappa_{min}^{e}(x)$ when $e(x,y)$ is a convex function of $y$;

(ii)

Assume that

[TABLE]

Then $D^{-1}\in{\mathfrak{D}}$ if the following condition is fulfilled

[TABLE]

where $\kappa^{e^{\prime}}(x)=\kappa_{max}^{e^{\prime}}(x)$ when $\frac{e(x,y)}{d(x,y)}$ is a concave function of $y$, and $\kappa^{e^{\prime}}(x)=\kappa_{min}^{e^{\prime}}(x)$ when $\frac{e(x,y)}{d(x,y)}$ is a convex function of $y$;

(iii)

Assume that

[TABLE]

Then $E^{-1}\in{\mathfrak{D}}$ if the following condition is fulfilled

[TABLE]

where $\kappa^{d^{\prime}}(y)=\kappa_{max}^{d^{\prime}}(y)$ when $\frac{d(x,y)}{e(x,y)}$ is a concave function of $x$, and $\kappa^{d^{\prime}}(y)=\kappa_{min}^{d^{\prime}}(y)$ when $\frac{d(x,y)}{e(x,y)}$ is a convex function of $x$.

Proof

We firstly prove ${\bf(i)}$. Denote

[TABLE]

Take ${\tilde{D}}=D-K^{d}\otimes I$ and ${\tilde{E}}=E-I\otimes K^{e}$. Then

[TABLE]

For any ${\bf u}=[u_{1,1},u_{2,1},\ldots,u_{M_{1},1},u_{1,2},\ldots,u_{M_{1},2},\ldots,u_{1,M_{2}},\ldots,u_{M_{1},M_{2}}]^{T}$, we have

[TABLE]

Referring to the proof of Theorem 2.3, it is easy to obtain

[TABLE]

where

[TABLE]

[TABLE]

So

[TABLE]

Which implies that ${\cal H}(A)$ is negative semi-definite if the conditions in (32)-(33) hold. Hence $I\in{\mathfrak{D}}$.

Similarly, one can show ${\bf(ii)}$ and ${\bf(iii)}$.

3.2 The Two-Dimensional Toeplitz Preconditioner

In this subsection, we extend the Toeplitz preconditioner to two-dimensional case. To solve (3) is equivalent to solve the following $N$ linear systems

[TABLE]

where $I_{k}$ denotes the $k\times k$ identity, ${\bf A}={I}_{\hat{M}}+{D}{B}_{x}+{E}{B}_{y}$, ${B}_{x}=-\eta_{x}({I}_{M_{2}}\otimes{G}_{\alpha})$, ${B}_{y}=-\eta_{y}({G}_{\beta}\otimes{I}_{M_{1}})$ $\hat{M}=M_{1}M_{2}$, $\eta_{x}=\tau/(2h_{1}^{\alpha})$, $\eta_{y}=\tau/(2h_{2}^{\beta})$, ${\bf b}^{n}=({I}_{\hat{M}}-{D}{B}_{x}-{E}{B}_{y}){u}^{n-1}+\tau{\bf f}^{n-\frac{1}{2}}$. Our two-level Toeplitz preconditioner for preconditioning (34) is defined as follows

[TABLE]

where $\bar{d}={\tt mean}({D})$, $\bar{e}={\tt mean}({E})$. The preconditioned Krylov subspace method with preconditioner ${\bf P}$ is employed to solve the linear systems in (34). Hence, in each iteration, it requires to compute some matrix-vector multiplications like ${\bf P}^{-1}{\bf z}$ for some randomly given ${\bf z}$, i.e., it requires to solve the linear system of the form

[TABLE]

Next, we introduce a multigrid method to solve (36).

For the choices of coarse-gird matrices, interpolation and restriction, we refer to the geometric grid coarsening, piecewise linear interpolation and its transpose. For the choice of pre-smoothing iteration, we refer to the block Jacobi iteration, i.e.,

[TABLE]

where ${\bf T}_{x}={I}_{\hat{M}}+\bar{d}{B}_{x}$ is the block diagonal part of ${\bf P}$, ${\bf x}^{k}$ is an initial guess of ${\bf x}$ in (36). Since ${\bf T}_{x}$ is a block diagonal matrix with identical Toeplitz blocks, its inversion, ${\bf T}_{x}^{-1}$ can be computed efficiently with the help of Gohberg-Semencul-type formula as discussed in Section 2. For the choice of post-smoother, we refer to the block Jacobi iteration for the permuted linear system, i.e,

[TABLE]

where ${\bf T}_{y}={I}_{\hat{M}}+\bar{e}{B}_{y}$, ${\bf x}^{k}$ is an initial guess of ${\bf x}$ in (36). One can easily find a $x$-$y$ ordering permutation matrix ${\bf Q}$ such that

[TABLE]

Thus, ${\bf T}_{y}^{-1}={\bf Q}^{\rm T}({I}_{\hat{M}}-\bar{e}\eta_{y}{I}_{M_{1}}\otimes{G}_{\beta})^{-1}{\bf Q}$, which means the implementation of (38) still requires to compute an inversion of a block diagonal matrix with identical Toeplitz blocks. Therefore, (38) can still be fast implemented using the Gohberg-Semencul-type formula. Similar to proof of Proposition 1, one can prove the following proposition.

Proposition 3

${\bf P}$* defined in (35) is invertible for any $\alpha\in(1,2)$.*

Theorem 3.3

Let $d(x,y)\equiv\nu_{1}a(x,y)$ and $e(x,y)\equiv\nu_{2}a(x,y)$ for any $(x,y)\in\Omega$ with nonnegative constants $\nu_{1}$ and $\nu_{2}$. Assume

(i)

$a(x,y)\in[\check{a},\hat{a}]$* with $\check{a}>0$ for any $(x,y)\in\Omega$,*

(ii)

For any $x\in(x_{L},x_{R})$, $a(x,\cdot)$ is convex or concave on $y\in(y_{L},y_{R})$; and for any $y\in(y_{L},y_{R})$, $a(\cdot,y)$ is convex or concave on $x\in(x_{L},x_{R})$.

(iii)

$\check{c}_{1}=\inf\limits_{y\in(y_{L},y_{R})}[\tilde{\mathcal{M}}_{1}(y)-\sqrt{2}(\hat{\mathcal{M}}_{1}(y)-\check{\mathcal{M}}_{1}(y))/\varsigma_{\alpha}]>0$* with $\hat{\mathcal{M}}_{1}(y)=:\sup\limits_{x\in(x_{L},x_{R})}a(x,y)$ and $\check{\mathcal{M}}_{1}(y)=:\inf\limits_{x\in(x_{L},x_{R})}a(x,y)$,*

[TABLE]

$\check{c}_{2}=\inf\limits_{x\in(x_{L},x_{R})}[\tilde{\mathcal{M}}_{2}(x)-\sqrt{2}(\hat{\mathcal{M}}_{2}(x)-\check{\mathcal{M}}_{2}(x))/\varsigma_{\beta}]>0$* with $\hat{\mathcal{M}}_{2}(x)=:\sup\limits_{y\in(y_{L},y_{R})}a(x,y)$ and $\check{\mathcal{M}}_{2}(x)=:\inf\limits_{y\in(y_{L},y_{R})}a(x,y)$,*

[TABLE]

Then, for any positive integers, $N$, $M_{1}$ and $M_{2}$, it holds $\Sigma^{2}({\bf A}{\bf P}^{-1})\subset[\check{s},\hat{s}]$ and thus

[TABLE]

where $\check{s}$, $\hat{s}$ are positive constants independent of $\tau$, $h_{1}$ and $h_{2}$:

[TABLE]

Proof

Denote

[TABLE]

with $a_{i,j}=a(x_{i},y_{j})$. Also, denote $\bar{a}={\tt mean}({D}_{a})$. By straightforward calculation,

[TABLE]

where ${W}=\nu_{1}{B}_{x}+\nu_{2}{B}_{y}$. Rewrite ${D}_{a}=\mathrm{diag}({D}_{a,1},{D}_{a,2},...,{D}_{a,M_{2}})$ with

${D}_{a,i}=\mathrm{diag}(a_{1,i},a_{2,i},...,a_{M_{1},i})$. Then, it is easy to see that ${B}_{x}^{\rm T}{D}_{a}+{D}_{a}{B}_{x}=\mathrm{diag}({H}_{1},{H}_{2},...,{H}_{M_{2}})$ with ${H}_{i}=-\eta_{x}({D}_{a,i}{G}_{\alpha}+{G}_{\alpha}^{\rm T}{D}_{a,i})$. Denote $l_{1}(y)=\tilde{\mathcal{M}}_{1}(y)-\sqrt{2}(\hat{\mathcal{M}}_{1}(y)-\check{\mathcal{M}}_{1}(y))/\varsigma_{\alpha}$ and $s_{1}(y)=\tilde{\mathcal{M}}_{1}(y)+\sqrt{2}(\hat{\mathcal{M}}_{1}(y)-\check{\mathcal{M}}_{1}(y))/\varsigma_{\alpha}$ . Then, applying Theorem 2.3 to ${\bf(i)}$–${\bf(iii)}$, we have

[TABLE]

Therefore,

[TABLE]

where the first ‘$\prec$’ is obvious. Recall the permutation matrix defined in (39). Denote $\tilde{B}_{y}:={\bf Q}{B}_{y}{\bf Q}^{\rm T}=-\eta_{y}{I}_{M_{1}}\otimes{G}_{\beta}$, $\tilde{D}_{a}=\mathrm{diag}(\tilde{D}_{a,1},\tilde{D}_{a,2},...,\tilde{D}_{a,M_{1}})$ with $\tilde{D}_{a,i}=\mathrm{diag}(a_{i,1},a_{i,2},...,a_{i,M_{2}})$. Then, it is easy to check that

[TABLE]

Similarly to proof of (42), applying Theorem 2.3 to ${\bf(i)}$, ${\bf(ii)}$ and ${\bf(iii)}$ yields

[TABLE]

Moreover, it is easy to see that

[TABLE]

By (42)–(44),

[TABLE]

For any non-zero vector ${y}\in\mathbb{R}^{M\times 1}$, denote ${z}={\bf P}^{-1}{y}$. Then, it holds

[TABLE]

By (45),

[TABLE]

By Proposition (2) and (3.2),

[TABLE]

During the proof above, there is no constraint on $M$ and $N$. Thus, for any $N\geq 1$, any $M\geq 1$, $\Sigma^{2}({\bf A}{\bf P}^{-1})\subset[\check{s},\hat{s}]$ and $\sup\limits_{N,M\geq 1}\rm cond({\bf A}{\bf P}^{-1})\leq\sqrt{\hat{s}/\check{s}}$.

4 Numerical experiments

In this section, we test several examples to support theoretical results of Theorems 2.2, 3.2 and to show the efficiency of the Toeplitz preconditioner. We employ generalized minimal residual (PGMRES) method with the Toeplitz preconditioner to solve (12) and (3). We denote PGMRES method with Toeplitz preconditioner by PGMRES-T. The stopping criterion for PGMRES-T is set as $\frac{||{\bf r}_{k}||_{2}}{||{\bf r}_{0}||_{2}}\leq$1e-7, where ${\bf r}_{k}$ denotes the residual vector at $k$-th iteration. Also, to illustrate the efficiency of PGMRES-T, we compare it with the direct solver, pivoted LU factorization (PLU). All numerical experiments are performed via MATLAB R2015a on a PC with the configuration: Intel(R) Core(TM) i7-4720 CPU 2.60 GHz and 8 GB RAM.

Recall that $h$ is the spatial step-size for one-dimensional discretization. We also set $h_{1}=h_{2}=h$ in two-dimensional discretization for the related experiments in this section. Define the error as

[TABLE]

Then, the spatial and temporal convergence rates are measured as follows

[TABLE]

Denote by CPU, the running time by unit seconds. Denote by ’iter’, the average of iteration numbers of PGMRES method for the $N$ linear systems in (12) or (3).

Example 1

Consider a one-dimensional OSFDE with $[x_{L},x_{R}]=[0,1]$, $T=1$ and

[TABLE]

The explicit expression of exact solution for the example is $u(x,t)=2^{6}x^{3}(1-x)^{3}t^{3}$.

We employ both PGMRES-T and PLU to solve the linear systems (23) arising from Example 1, the results of which are listed in Tables 1–2.

From Tables 1–2, we see that the CPU cost of PGMRES-T is much less than that of PLU solver while the error, $E(h,\tau)$ of the two solvers are almost the same, which demonstrate the efficiency of the Toeplitz preconditioner. Also, as $\tau$ or $h$ changes in Tables 1–2, the iteration number of PGMRES-T varies slightly, which shows a linear convergence of PGMRES-T. Moreover, the temporal convergence rate, $Rate_{h}$ and the spatial convergence rate,$Rate_{\tau}$ from Tables 1–(2) are always close to 2, which supports the theoretical result of Theorem 2.2.

Example 2

Consider a two-dimensional OSFDE with $[x_{L},x_{R}]=[y_{L},y_{R}]=[0,2]$, $T=1$ and

[TABLE]

The explicit expression of exact solution for the example is $u(x,y,t)=x^{4}(2-x)^{4}y^{4}(2-y)^{4}t^{3}$.

We employ both PGMRES-T and PLU to solve the linear systems (34) arising from Example 2, the results of which are listed in Tables 3–4.

Tables 3–4 shows that the CPU cost of PGMRES-T is much less than that of PLU, which demonstrates the efficiency of the Toeplitz preconditioner in two-dimensional case. Again, the iteration number of PGMRES-T changes slightly as $\tau$ or $h$ changes in Tables 3–4, which shows a linear convergence of PGMRES-T. Moreover, the convergence rates, $Rate_{h}$ and $Rate_{\tau}$ shown in Tables 3–4 are close to 2, which coincides with the theoretical result of Theorem 3.2.

5 Concluding remarks

We study second-order schemes for time-dependent one- and two-dimensional OSFDEs with variable diffusion coefficients, in which implicit Crank-Nicolson scheme and WSGD formula are employed to discretize the temporal and the spatial derivatives, respectively. Theoretically, we have established the unconditional stability and second-order convergence for the one-dimensional scheme without additional assumption, and for the two-dimensional scheme with certain assumptions on diffusion coefficients presented in Corollaries 1–2. To accelerate the solution process, Toeplitz preconditioners have been proposed for both one- and two-dimensional schemes. The condition numbers of the preconditioned matrices have been proven to be bounded by a constant independent of discretization step-sizes under certain assumptions on the diffusion coefficients presented in Theorems 2.4, 2.5, 3.3. Numerical results reported have shown the second-order convergence rate of the proposed schemes and the efficiency of the proposed preconditioners.