A Fast Preconditioned Semi-Implicit Difference Scheme for Strongly Nonlinear Space-Fractional Diffusion Equations

Abstract

In this paper, we propose a semi-implicit difference scheme for solving one-dimensional nonlinear space-fractional diffusion equations. The method is first-order accurate in time and second-order accurate in space. It uses a fractional central difference formula and the backward Euler method to approximate its space and time derivatives, respectively. Stability and convergence properties of the proposed scheme are proved with the help of a discrete Grönwall inequality. Moreover, we extend the method to the solution of two-dimensional nonlinear models. A fast matrix-free implementation based on preconditioned Krylov subspace methods is presented for solving the discretized linear systems. The resulting fast preconditioned semi-implicit difference scheme reduces the memory requirement of conventional semi-implicit difference schemes from $\mathcal{O}\left(N_s^2\right)$ to $\mathcal{O}\left(N_s\right)$ and the computational complexity from $\mathcal{O}\left(N_s^3\right)$ to $\mathcal{O}(N_{s}log{N}_s)$ in each iterative step, where $N_s$ is the number of space grid points. Experiments with two numerical examples are shown to support the theoretical findings and to illustrate the efficiency of our proposed method.

1. Introduction

In this paper, we are interested in developing efficient numerical solutions of strongly nonlinear space-fractional diffusion equations (NSFDEs) of the following form:

$$\left\{\begin{array}{cc}\frac{\partial u(x,t)}{\partial t}=a\left(u\right)\frac{{\partial }^{\alpha }u(x,t)}{{\partial \left|x\right|}^{\alpha }}+f(x,t),\hfill & (x,t)\in (0,L_x)\times (0,T],\hfill \\ u(x,0)=\varphi \left(x\right),\hfill & x\in [0,L_x],\hfill \\ u(0,t)=\phi \left(t\right),u(L_x,t)=\psi \left(t\right),\hfill & t\in [0,T],\hfill \end{array}\right.$$

(1)

where $a\left(u\right),f(x,t),\varphi \left(x\right),\phi \left(t\right),\psi \left(t\right)$ are known functions. Function $u(x,t)$ may represent, for example, the concentration of a particle plume undergoing anomalous diffusion. The diffusion coefficient $a\left(u\right)$ is assumed to be positive, and the forcing function $f(x,t)$ represents the source or sink term [1,2]. In addition, we assume that Equation (1) has a smooth solution $u(x,t)$, and that there exist real positive numbers ${ϵ}_1,{ϵ}_2,ϵ,a_0,a_1$ and L satisfying the following:

$$\begin{array}{c}\hfill 0<a_0\le a(u(x,t)+{ϵ}_1)\le a_1,\end{array}$$

(2)

$$\begin{array}{c}\hfill |a(u(x,t)+{ϵ}_1)-a(u(x,t)+{ϵ}_2)|\le L|{ϵ}_1-{ϵ}_2|,\end{array}$$

(3)

where $|{ϵ}_1|\le ϵ,|{ϵ}_2|\le ϵ$, and $(x,t)\in [0,L_x]\times [0,T]$. Moreover, the symbol ${\displaystyle \frac{{\partial }^{\alpha }}{{\partial \left|x\right|}^{\alpha }}}$ denotes the Riesz fractional derivative operator [3,4,5] for $1<\alpha \le 2$, defined as follows:

$$\frac{{\partial }^{\alpha }u(x,t)}{{\partial \left|x\right|}^{\alpha }}=-\frac{1}{2cos\left(\frac{\alpha \pi }{2}\right)\Gamma (2-\alpha )}\frac{{\partial }^2}{\partial x^2}{\int }_{-\infty }^{\infty }{|x-\xi |}^{1-\alpha }u(\xi ,t)d\xi ,$$

(4)

where $\Gamma (·)$ is the Gamma function.

In the past few decades, a significant amount of research work has been devoted to the analysis and solution of fractional partial differential equations (fPDEs) especially for modeling anomalous diffusion [6,7], subdiffusion and and superdiffusion processes that can be observed in many real physical systems, such as heat transfer problems and wave analysis [8,9]. Space fractional equations are obtained by replacing the traditional space derivative with a fractional derivative in the integer-order diffusion equation and can describe heterogeneous substances with memory and genetic properties very effectively [10]. Closed-form analytic solutions are available only for limited classes of fPDEs. For solving more general equations, the only viable approach is to use discretization methods such as finite difference [11,12,13,14,15], finite element [16,17], finite volume [18,19] and spectral [20] methods.

The Riemann–Liouville fractional derivative is often used to model the dynamics of discrete systems with long-range interaction. Several studies have applied the weighted and shifted Grünwald difference (WSGD) formula to approximate the Riemann–Liouville derivative; see [14,21,22] for details. By using the WSGD formula for two-sided fractional derivatives combining the compact technique, Hao et al. presented a fourth-order difference approximation for the space fractional derivatives [23]. In addition, Çelik and Duman used a fractional centered difference operator to approximate the Riesz fractional derivative with second-order accuracy and to solve the Riesz SFDE (RSFDE) on a finite domain [24]. Later, Yang et al. [25] used the Diethelm method [15] to derive a new $\mathcal{O}(\Delta x^{3-\alpha })$ approximation for discretizing the Riesz fractional derivatives of order $\alpha \in (1,2)$, where $\Delta x$ is the step size. Then, a new finite difference scheme for solving the RSFDE was obtained by discretizing the first-order time derivative using the Crank–Nicolson (CN) method. In contrast with the classical diffusion operator $\Delta$, Riesz fractional derivative [3,5] is a special linear combination of left- and right-sided Riemann–Liouville fractional derivatives. By exploiting this property, Ding and Li [4] established a novel class of high-order numerical algorithms for Riesz derivatives through constructing new generating functions. In 2018, Lin et al. [26] studied the CN temporal discretization with various high-order spatial difference schemes for RSFDEs with variable diffusion coefficients and gave the unconditional stability and convergence analysis for temporally second-order and spatially jth-order ($j=1,2,3,4$) difference schemes for such equations with variable coefficients. In 2019, Lin et al. [27] studied the CN alternative direction implicit (CN-ADI) method for two-dimensional RSFDEs with nonseparable coefficients. They showed under mild assumptions the unconditional stability of the CN-ADI method in discrete ${\ell }^2$-norm and the consistency of cross perturbation terms arising from the CN-ADI method. We point the reader to [3,12,13,28] for other recent works on the numerical solution of fPDEs with the Riesz fractional derivative.

Due to their non-local nature, numerical approximations of space fractional operators generally lead to dense matrices. Conventional direct solution methods for linear systems based on variants of the Gaussian elimination algorithm may not be affordable to use because of their high computational complexity and large memory costs. The quest for fast algorithms for solving linear systems originated in the numerical discretization of fPDEs has become a topic of great interest in the recent literature on numerical methods for fPDEs. Lin et al. [29] proposed a splitting preconditioner for fast solution of Toeplitz-like linear systems arising from one- and two-dimensional time-dependent RSFDEs with variable diffusion coefficients [27]. Some other preconditioners are being applied also to solving the discretized (non-)linear systems from the solution of linear and non-linear problems; see, for example, [1,2,30,31] and the references therein.

In this paper, we study the efficient numerical solution of strongly NSFDEs because they are more suitable than linear models to describe some difficult physical processes, such as fractional distillation in nonlinear space, and when the diffusion coefficient is related to the concentration of a particle plume [31]. However, to the best of our knowledge, only a few studies on fast numerical schemes for nonlinear or semi-linear problems have been published in the literature; see, for example, refs. [32,33,34]. Our main contribution to the subject can be summarized as follows: first, for this class of strongly NSFDEs, we propose a semi-implicit difference scheme that can avoid solving the discretized nonlinear systems [35], and analyze its stability and convergence properties. Then, we extend the method to the solution of the two-dimensional formulation of Equation (1). In order to compute an efficient iterative solution of the discretized systems, a suitable preconditioning technique is developed and presented in the paper. The complexity and memory requirement of the resulting semi-implicit difference scheme are fairly lower, compared to the fully implicit difference scheme that requires to solve the discretized nonlinear systems, motivating our choice of the semi-implicit difference scheme for Equation (1) [36]. Moreover, unlike the traditional (semi-)implicit difference schemes, our proposed semi-implicit difference scheme can be implemented as a matrix-free method because it does not store any of the coefficient matrices involved in the whole solution process.

The rest of this article is organized as follows. In Section 2, we propose a semi-implicit difference scheme for Equation (1), and then we analyze its stability and convergence properties. In Section 3, we extend our method to the solution of the two-dimensional formulation of Equation (1). A fast matrix-free preconditioned iterative algorithm for solving the discretized linear system arising in the analysis of the two-dimensional case of Equation (1) is developed, and it is presented in Section 4. In Section 5, we report on numerical experiments to illustrate the efficiency of the proposed semi-implicit difference method. Some conclusions arising from this work are drawn in Section 6.

2. The Semi-Implicit Difference Scheme and Its Stability and Convergence Analysis

In this section, we propose a semi-implicit difference method to approximate the solution of Equation (1). The proposed method uses a fractional central difference formula to discretize the Riesz fractional derivative; it is stable and achieves first-order convergence in time and second-order convergence in space, respectively.

2.1. The Semi-Implicit Difference Scheme

We introduce the following notation. Let $M,N$ be positive integers, $h=L_x/N,x_i=ih,$$t_j=j\Delta t$ and $\Delta t=T/M$. We consider the sets of grid points ${\Omega }_h:=\left\{x_i|0\le i\le N\right\},{\Omega }_{\Delta t}:=\left\{t_j|0\le j\le M\right\},{\Omega }_{h\Delta t}:=\left\{(x_i,t_j)\right|0\le i\le N,0\le j\le M\}$, and let $v=\left\{v_i^j\right|0\le i\le N,0\le j\le M\}$ be a grid function defined on ${\Omega }_{h\Delta t}$. To simplify the notation, we write Equation (1) at the generic points $(x_i,t_j)$ as follows:

$$\frac{\partial u(x_i,t_j)}{\partial t}=a\left(u(x_i,t_j)\right)\frac{{\partial }^{\alpha }u(x_i,t_j)}{{\partial \left|x\right|}^{\alpha }}+f(x_i,t_j),1\le i\le N-1,1\le j\le M,$$

(5)

or, equivalently, as follows:

$$\frac{\partial u(x_i,t_j)}{\partial t}-a\left(u(x_i,t_{j-1})\right)\frac{{\partial }^{\alpha }u(x_i,t_j)}{{\partial \left|x\right|}^{\alpha }}=F^{\prime }(x_i,t_j),1\le i\le N-1,1\le j\le M,$$

(6)

where ${\displaystyle F^{\prime }(x_i,t_j)=f(x_i,t_j)+\left[a\left(u(x_i,t_j)\right)-a\left(u(x_i,t_{j-1})\right)\right]\frac{{\partial }^{\alpha }u(x_i,t_j)}{{\partial \left|x\right|}^{\alpha }}}$. Additionally, we introduce a grid function U defined as follows:

$$U_i^j:=u(x_i,t_j),f_i^j:=f(x_i,t_j),0\le i\le N,0\le j\le M.$$

(7)

By these notations, the partial derivative of $u(x,t)$ with respect to t evaluated at point $(x_i,t_j)$ is expressed as follows:

$$\frac{\partial u(x_i,t_j)}{\partial t}=\frac{u(x_i,t_j)-u(x_i,t_{j-1})}{\Delta t}+R_{i,j}^{\left(1\right)}\triangleq \frac{U_i^j-U_i^{j-1}}{\Delta t}+\mathcal{O}(\Delta t),$$

(8)

and the Riesz space fractional derivative as follows:

$$\frac{{\partial }^{\alpha }u(x_i,t_j)}{{\partial \left|x\right|}^{\alpha }}=-h^{-\alpha }{\Delta }_h^{\alpha }u(x_i,t_j)+R_{i,j}^{\left(2\right)}\triangleq {\delta }_h^{\alpha }{U}_i^j+\mathcal{O}\left(h^2\right),$$

(9)

where ${\displaystyle {\Delta }_h^{\alpha }{U}_i^j=\sum _{k=i-N}^i{g}_k^{\left(\alpha \right)}{U}_{i-k}^j}$ and ${\displaystyle g_k^{\left(\alpha \right)}=\frac{{(-1)}^k\Gamma (\alpha +1)}{\Gamma (\alpha /2-k+1)\Gamma (\alpha /2+k+1)}}$; refer to [24] for details. For brevity, we denote the following:

$$R_{ij}=-R_{i,j}^{\left(1\right)}-R_{i,j}^{\left(2\right)}+\left[a\left(u(x_i,t_j)\right)-a\left(u(x_i,t_{j-1})\right)\right]\frac{{\partial }^{\alpha }u(x_i,t_j)}{{\partial \left|x\right|}^{\alpha }}.$$

(10)

Finally, Equation (1) at grid points $(x_i,t_j)$ is written as follows:

$$\frac{U_i^j-U_i^{j-1}}{\Delta t}-a\left(U_i^{j-1}\right){\delta }_h^{\alpha }{U}_i^j=f_i^j+{\widehat{R}}_{ij},1\le i\le N-1,1\le j\le M,$$

(11)

where the terms $\left\{{\widehat{R}}_{ij}\right\}$ are small and satisfy the inequality as follows:

$$|{\widehat{R}}_{ij}|\le c_1(\Delta t+h^2),1\le i\le N-1,1\le j\le M.$$

(12)

Since the quantities $\left\{{\widehat{R}}_{ij}\right\}$ are small, we can omit them. By imposing the initial-boundary value conditions,

$$\begin{array}{c}\hfill U_i^0=\varphi \left(x_j\right),0\le i\le N,\end{array}$$

(13)

$$\begin{array}{c}\hfill U_0^j=\phi \left(t_j\right),U_N^j=\psi \left(t_j\right),1\le j\le M,\end{array}$$

(14)

the following semi-implicit difference scheme is obtained:

$$\left\{\begin{array}{cc}{\displaystyle \frac{u_i^j-u_i^{j-1}}{\Delta t}-a\left(u_i^{j-1}\right){\delta }_h^{\alpha }{u}_i^j=f_i^j,}\hfill & 1\le i\le N-1,1\le j\le M,\hfill \\ u_i^0=\varphi \left(x_j\right),\hfill & 0\le i\le N,\hfill \\ u_0^j=\phi \left(t_j\right),u_N^j=\psi \left(t_j\right),\hfill & 1\le j\le M,\hfill \end{array}\right.$$

(15)

that requires to solve the discretized linear system at each time step.

2.2. Stability and Convergence Analysis

In this subsection, we study the stability and convergence properties of the numerical scheme given by Equation (15). Let $u_i^j$ be the approximate solution of Equation (15), and introduce the error function as follows:

$$e_i^j=U_i^j-u_i^j,i=0,\cdots ,N,j=0,\cdots ,M,$$

defined as

$$\left\{\begin{array}{cc}{\displaystyle \frac{e_i^j-e_i^{j-1}}{\Delta t}-a\left(u_i^{j-1}\right){\delta }_h^{\alpha }{e}_i^j=\left[a\left(U_i^{j-1}\right)-a\left(u_i^{j-1}\right)\right]{\delta }_h^{\alpha }{U}_i^j+{\widehat{R}}_{ij},}\hfill & 1\le i\le N-1,1\le j\le M,\hfill \\ e_i^0=0,\hfill & 0\le i\le N,\hfill \\ e_0^j=0,e_N^j=0,\hfill & 1\le j\le M,\hfill \end{array}\right.$$

(16)

where $|{\widehat{R}}_{ij}|\le c_2(\Delta t+h^2)$. Denote $E_j={[e_1^j,e_2^j,\cdots ,e_{N-1}^j]}^{\top }$, and assume the following:

$$\parallel E^j{\parallel }_{\infty }=|e_{{\ell }_j}|=\underset{1\le \ell \le N-1}{max}\left|e_{\ell }^j\right|(0\le j\le M,1\le {\ell }_j\le N-1).$$

(17)

To prove the stability and assess the convergence order of the scheme (15), we need some preliminary theoretical results given in [24,35,37] that are recalled below.

Lemma 1.

Let

$$g_k^{\left(\alpha \right)}=\frac{{(-1)}^k\Gamma (\alpha +1)}{\Gamma (\alpha /2-k+1)\Gamma (\alpha /2+k+1)}$$

be the coefficients of the centered finite difference approximation (9) for $k=0,\pm 1,\pm 2,\dots$, and $\alpha >-1$. Then, we have the following:

(1) 

$g_0^{\left(\alpha \right)}\ge 0$ and $g_{-k}^{\left(\alpha \right)}=g_k^{\left(\alpha \right)}\le 0$ for all $\left|k\right|\ge 1$;

(2) 

${\displaystyle \sum _{-\infty }^{\infty }}{g}_k^{\left(\alpha \right)}=0$ and ${\displaystyle \sum _{k=i-N,k\ne 0}^i}{g}_k^{\left(\alpha \right)}\ge -g_0$.

Lemma 2.

(Discrete Grönwall inequality) Let $\left\{F^k\right|k\le 0\}$ be a sequence of non-negative numbers and let $\Delta t>0$ be such that the following is true:

$$F^{k+1}\le (1+c\Delta t)F^k+(\Delta t)g,k=0,1,2,\dots ,$$

(18)

where both c and g are non-negative scalars. Then, the following inequality holds:

$$F^k\le e^{ck\Delta t}\left(F^0+\frac{g}{c}\right)$$

(19)

with $k=0,1,\dots .$

At this stage, the following result on the numerical stability of the proposed one dimensional semi-implicit difference scheme can be proved:

Theorem 1.

Suppose that $\left\{u_i^j\right|0\le j\le M,0\le i\le N\}$ is the solution of the semi-implicit difference scheme (15), and define the following:

$$C_1=\frac{c_1}{L{M}_{\alpha }}exp\left(L{M}_{\alpha }T\right),M_{\alpha }=\left|\underset{\begin{array}{c}0\le x\le L,\\ 1\le t\le T\end{array}}{max}\frac{{\partial }^{\alpha }u(x,t)}{{\partial \left|x\right|}^{\alpha }}\right|.$$

(20)

Then, if $\Delta t$ and h satisfy the following condition

$$\Delta t\le \frac{ϵ}{2C_1},h\le \sqrt{\frac{ϵ}{2C_1}},$$

(21)

the semi-implicit difference scheme (15) is stable, and the following inequality holds:

$$\parallel e^k{\parallel }_{\infty }\le C_1(\Delta t+h^2),0\le k\le n.$$

(22)

Proof. 

According to Equation (15), we know that $\parallel e^k{\parallel }_{\infty }\le C_1(\Delta t+h^2)$ holds for $k=0$. Suppose that the inequality (22) holds for $0\le k\le l$, when $\Delta t$ and h satisfy the condition (21) and

$$\parallel e^k{\parallel }_{\infty }\le ϵ,0\le k\le l.$$

(23)

Therefore, the following holds:

$$a_0\le a\left(u_i^j\right)\le a_1,|a\left(U_i^j\right)-a\left(u_i^j\right)|\le L|e_i^j|,1\le i\le N-1,0\le j\le l.$$

(24)

From Equation (17), we have the following:

$$e_i^j-e_i^{j-1}=\Delta ta\left(u_i^{j-1}\right){\delta }_h^{\alpha }{e}_i^j+\Delta t\left[a\left(U_i^{j-1}\right)-a\left(u_i^{j-1}\right)\right]{\delta }_h^{\alpha }{U}_i^j+\Delta t{\widehat{R}}_{ij},$$

where ${\delta }_h^{\alpha }{e}_i^j=-h^{(-\alpha )}\sum _{k=i-N}^i{g}_k^{\left(\alpha \right)}{e}_{i-k}^j$. Thus we obtain the following:

$$e_i^j-e_i^{j-1}=-\frac{\Delta t}{h^{\alpha }}a\left(u_i^{j-1}\right)\left[\sum _{k=i-N,k\ne 0}^i{g}_k^{\left(\alpha \right)}{e}_{i-k}^j+g_0^{\left(\alpha \right)}{e}_i^j\right]+\Delta t\left[a\left(U_i^{j-1}\right)-a\left(u_i^{j-1}\right)\right]{\delta }_h^{\alpha }{U}_i^j+\Delta t{\widehat{R}}_{ij}.$$

So, we rewrite Equation (17) to obtain the following:

$$\left[1+\mu g_0^{\left(\alpha \right)}a\left(u_i^{j-1}\right)\right]e_i^j=\mu a\left(u_i^{j-1}\right)\sum _{k=i-N,k\ne 0}^i{g}_k^{\left(\alpha \right)}{e}_{i-k}^j+e_i^{j-1}+\Delta t\left[a\left(U_i^{j-1}\right)-a\left(u_i^{j-1}\right)\right]{\delta }_h^{\alpha }{U}_i^j+\Delta t{\tilde{R}}_{ij},$$

(25)

where $\mu =\Delta t/h^{\alpha }$. If we take infinite norms on both sides of Equation (25), and we set $i={\ell }_i$ obtained so that

$$\begin{array}{c}\hfill \left[1+\mu g_{0}a\left(u_i^{j-1}\right)\right]\parallel e^j{\parallel }_{\infty }\le \mu a\left(u_i^{j-1}\right)g_0^{\left(\alpha \right)}\parallel e^j{\parallel }_{\infty }+\parallel e^{j-1}{\parallel }_{\infty }+\Delta tL{M}_{\alpha }{\parallel e^{j-1}\parallel }_{\infty }+c_1\Delta t(\Delta t+h^2)\end{array}$$

for $1\le i\le N-1,1\le j\le l+1$. Thus, we have the following:

$$\parallel e^j{\parallel }_{\infty }=(1+\Delta tL{M}_{\alpha }){\parallel e^{j-1}\parallel }_{\infty }+c_1\Delta t(\Delta t+h^2),1\le j\le l+1.$$

By applying the discrete Grönwall inequality (cf. Lemma 2) to $F^k=e^{l-1}$, $c=L{M}_{\alpha }$ and $g=c_1(\Delta t+h^2)$, we obtain the inequality as follows:

$$\begin{array}{cc}\hfill \parallel e^{l+1}{\parallel }_{\infty }& =exp\left[L{M}_{\alpha }(l+1)\Delta t\right]·\left[\parallel e^0{\parallel }_{\infty }+\frac{c_1(\Delta t+h^2)}{L{M}_{\alpha }}\right]\hfill \\ \hfill & \le \frac{c_1}{L{M}_{\alpha }}exp\left(L{M}_{\alpha }T\right)(\Delta t+h^2)\hfill \\ \hfill & \triangleq C_1(\Delta t+h^2).\hfill \end{array}$$

The above proof shows that the semi-implicit difference scheme (15) is stable in the discrete $L_{\infty }$-norm.    □

3. Extensions to the Two-Dimensional Problem

In this section, we extend our method to the solution of two-dimensional problems. For this purpose, let $N_x,N_y$ be positive integers and define a spatial partition of nodes $h=L/N,x_i=i{h}_x(i=1,2,\cdots ,N_x),y_j=j{h}_y(j=1,2,\cdots ,N_y)$ with step sizes $h_x=\frac{L_x}{N_x}$, $h_y=\frac{L_y}{N_y}$ in the x-axis and y-axis, respectively, and a uniform time partition $t_k=k\Delta t$ of the interval $[0,T]$ with time step $\Delta t=T/M$. We consider the sets of grid points ${\Omega }_h:=\left\{(x_i,y_j)\right|0\le i\le N_x,0\le j\le N_y\},{\Omega }_{\Delta t}:=\left\{t_k\right|0\le k\le M\},{\Omega }_{h\Delta t}:=\left\{(x_i,y_j,t_k)\right|0\le i\le N_x,0\le j\le N_y,0\le k\le M\}$, and let $u=\left\{u_{i,j}^k\right|0\le i\le N_x,0\le j\le N_y,0\le k\le M\}$ be a grid function defined on ${\Omega }_{h\Delta t}$.

We consider the following two-dimensional problem:

$$\left\{\begin{array}{cc}{\displaystyle \frac{\partial u(x,y,t)}{\partial t}=a\left(u\right)\left[\frac{{\partial }^{\alpha }u(x,y,t)}{{\partial \left|x\right|}^{\alpha }}+\frac{{\partial }^{\alpha }u(x,y,t)}{{\partial \left|y\right|}^{\alpha }}\right]+f(x,y,t),}\hfill & (x,y,t)\in \Omega \times (0,T],\hfill \\ u(x,y,0)=\varphi (x,y),\hfill & (x,y)\in \Omega ,\hfill \\ u(x,y,t)=0,\hfill & (x,y,t)\in \partial \Omega \times (0,T].\hfill \end{array}\right.$$

(26)

where $\Omega =[0,L_x]\times [0,L_y]$. Upon evaluating Equation (26) at the points $(x_i,y_j,t_k)$, we obtain the following discrete equation:

$$\begin{array}{cc}\hfill \frac{\partial u(x_i,y_j,t_k)}{\partial t}& =a\left(u(x_i,y_j,t_k)\right)\left[\frac{\partial u(x_i,y_j,t_k)}{{\partial \left|x\right|}^{\alpha }}+\frac{\partial u(x_i,y_j,t_k)}{{\partial \left|y\right|}^{\alpha }}\right]+f(x_i,y_j,t_k),\hfill \\ \hfill & 1\le i\le N_x-1,1\le j\le N_y-1,1\le k\le M,\hfill \end{array}$$

(27)

which can be written equivalently as follows:

$$\begin{array}{cc}\hfill & \frac{\partial u(x_i,y_j,t_k)}{\partial t}-a\left(u(x_i,y_j,t_{k-1})\right)\left[\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|x\right|}^{\alpha }}+\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|y\right|}^{\alpha }}\right]=\hfill \\ \hfill & \left[a\left(u(x_i,y_j,t_k)\right)-a\left(u(x_i,y_j,t_{k-1})\right)\right]\left[\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|x\right|}^{\alpha }}+\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|y\right|}^{\alpha }}\right]+f(x_i,y_j,t_k).\hfill \end{array}$$

(28)

Let U be the grid function defined as follows:

$$U_{i,j}^k:=u(x_i,y_j,t_k),f_{i,j}^k:=f(x_i,y_j,t_k),0\le i\le N_x,0\le j\le N_y,0\le k\le M.$$

Thus, the partial derivative of t at the point $(x_i,y_j,t_k)$ is as follows:

$$\frac{\partial u(x_i,y_j,t_k)}{\partial t}=\frac{u(x_i,y_j,t_k)-u(x_i,y_j,t_{k-1})}{\Delta t}+R_{i,j,k}^{\left(1\right)}=\frac{u_{i,j}^k-u_{i,j}^{k-1}}{\Delta t}+\mathcal{O}(\Delta t),$$

and

$$\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|x\right|}^{\alpha }}=-h_x^{-\alpha }{\Delta }_{h_x}^{\alpha }u(x_i,y_j,t_k)+R_{i,j,k}^{\left(2\right)}\triangleq {\delta }_{h_x}^{\alpha }{U}_i^j+\mathcal{O}\left(h_x^2\right),$$

$$\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|y\right|}^{\alpha }}=-h_y^{-\alpha }{\Delta }_{h_y}^{\alpha }u(x_i,y_j,t_k)+R_{i,j,k}^{\left(3\right)}\triangleq {\delta }_{h_y}^{\alpha }{U}_i^j+\mathcal{O}\left(h_y^2\right),$$

where ${\displaystyle {\Delta }_{h_x}^{\alpha }{U}_{i,j}^k=\sum _{\ell =i-N_x}^i{g}_l^{\left(\alpha \right)}{U}_{i-\ell ,j}^k}$, ${\displaystyle {\Delta }_{h_y}^{\alpha }{U}_{i,j}^k=\sum _{\ell =j-N_y}^j{g}_{\ell }^{\left(\alpha \right)}{U}_{i,j-\ell }^k}$ and ${\displaystyle g_{\ell }^{\left(\alpha \right)}=}$$\frac{{(-1)}^{\ell }\Gamma (\alpha +1)}{\Gamma (\alpha /2-\ell +1)\Gamma (\alpha /2+\ell +1)}$.

We denote for simplicity the following:

$$\begin{array}{c}\hfill {\displaystyle R_{i,j,k}=-R_{i,j,k}^{\left(1\right)}-a\left(u(x_i,y_j,t_{k-1})\right)R_{i,j,k}^{\left(2\right)}-a\left(u(x_i,y_j,t_{k-1})\right)R_{i,j,k}^{\left(3\right)}+R^{\prime },}\end{array}$$

(29)

where ${\displaystyle R^{\prime }=\left[a\left(u(x_i,y_j,t_k)\right)-a\left(u(x_i,y_j,t_{k-1})\right)\right]\left[\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|x\right|}^{\alpha }}+\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{{\partial \left|y\right|}^{\alpha }}\right]}$, so that the equation can be rewritten as follows:

$$\frac{U_{i,j}^k-U_{i,j}^{k-1}}{\Delta t}-a\left(U_{i,j}^{k-1}\right)\left[{\delta }_{h_x}^{\alpha }{U}_{i,j}^k+{\delta }_{h_y}^{\alpha }{U}_{i,j}^k\right]=f_{i,j}^k+R_{i,j,k},$$

(30)

with $1\le i\le N_x-1,1\le j\le N_y-1,1\le k\le M$.

Observe that the terms $R_{i,j,k}$ are small. Omitting the quantities $R_{i,j,k}$ in the previous expression, and imposing the initial-boundary value conditions,

$$\begin{array}{c}\hfill U_{i,j}^0=\varphi (x_i,y_j),0\le i\le N_x,0\le j\le N_y;\\ \hfill U_{0,j}^k=U_{N_x,j}^k=\phi \left(t_k\right),U_{i,0}^k=U_{i,N_y}^k=\phi \left(t_k\right),0\le k\le M,\end{array}$$

(31)

the following semi-implicit difference scheme is obtained, where by $u_{i,j}^k$ we denote the approximate solution of Equation (31) at the points $(x_i,y_j,t_k)$:

$$\left\{\begin{array}{cc}\frac{u_{i,j}^k-u_{i,j}^{k-1}}{\Delta t}-a\left(u_{i,j}^{k-1}\right)\left[{\delta }_{h_x}^{\alpha }{u}_{i,j}^k+{\delta }_{h_y}^{\alpha }{u}_{i,j}^k\right]=f_{i,j}^k,\hfill & 1\le i\le N_x-1,1\le i\le N_y-1,1\le k\le M\hfill \\ u_{i,j}^0=\varphi (x_i,y_j),\hfill & 1\le i\le N_x-1,1\le i\le N_y-1,\hfill \\ u_{0,j}^k=u_{N_x,j}^k=\phi \left(t_k\right),u_{i,0}^k=u_{i,N_y}^k=\psi \left(t_k\right),\hfill & 1\le k\le M.\hfill \end{array}\right.$$

(32)

Below, we analyze stability and convergence properties of the two-dimensional semi-implicit difference scheme given by Equation (32). We define the following:

$$e_{i,j}^k=U_{i,j}^k-u_{i,j}^k,i=0,\cdots ,N_x,j=1,\cdots ,N_y,k=0,\cdots ,M,$$

the error satisfying the following equations:

$$\left\{\begin{array}{c}{\displaystyle \frac{e_{i,j}^k-e_{i,j}^{k-1}}{\Delta t}-a\left(u_{i,j}^{k-1}\right)\left[{\delta }_{h_x}^{\alpha }{e}_{i,j}^k+{\delta }_{h_y}^{\alpha }{e}_{i,j}^k\right]=\left[a\left(U_{i,j}^{k-1}\right)-a\left(u_{i,j}^{k-1}\right)\right]\left[{\delta }_{h_x}^{\alpha }{U}_{i,j}^k+{\delta }_{h_y}^{\alpha }{U}_{i,j}^k\right]+R_{i,j,k},}\hfill \\ 1\le i\le N_x-1,1\le j\le N_y-1,1\le k\le M,\hfill \\ e_{i,j}^0=0,0\le i\le N_x,0\le j\le N_y,\hfill \\ e_{0,j}^k=e_{N_x,j}^k=0,e_{i,0}^k=e_{i,N_y}^k=0,1\le k\le M.\hfill \end{array}\right.$$

(33)

By a similar proof to Theorem 1, the following result is established about the numerical stability of the semi-implicit difference scheme given by Equation (26).

Theorem 2.

Suppose that $\left\{u_{i,j}^k\right|0\le i\le N_x,0\le j\le N_y,0\le k\le M\}$ is the solution of the semi-implicit difference scheme (25), and define the following:

$$C_2=\frac{c_2}{L{M}_{\alpha }}exp\left(L{M}_{\alpha }T\right),M_{\alpha }=\left|\underset{\begin{array}{c}0\le x,y\le L,\\ 1\le t\le T\end{array}}{max}\left[\frac{{\partial }^{\alpha }u(x,t)}{{\partial \left|x\right|}^{\alpha }}+\frac{{\partial }^{\alpha }u(x,y,t)}{{\partial \left|y\right|}^{\alpha }}\right]\right|.$$

(34)

Then, when constants $\Delta t,h(=h_x=h_y)$ satisfy the following condition,

$$\Delta t\le \frac{ϵ}{2C_2},h\le \sqrt{\frac{ϵ}{2C_2}},$$

(35)

the semi-implicit difference scheme (32) is stable and as follows:

$$\parallel e^k{\parallel }_{\infty }\le C_2(\Delta t+h^2),0\le k\le M.$$

(36)

4. Fast Implementation of the Semi-Implicit Difference Scheme

In this section, we present a fast matrix-free implementation of the proposed one-dimensional and two-dimensional semi-implicit difference schemes that use the preconditioned biconjugate gradients stabilized (BiCGSTAB) method [38] to iteratively solve the nonsymmetric linear system with a Toeplitz-like structure arising from the discretization. Such an implementation enables us to reduce both the algorithmic complexity and the memory requirements for the solution significantly.

First, we analyze the properties of the coefficient matrices of the discretized linear systems.

4.1. One-Dimensional Coefficient Matrix

By rewriting Equation (15) in the following form,

$$\frac{u^j-u^{j-1}}{\Delta t}=\frac{1}{h^{\alpha }}{D}^{(j-1)}A{u}^j+\frac{1}{h^{\alpha }}{D}^{(j-1)}{v}^{j-1}+f^j,$$

(37)

where $v^{j-1}={[g_1^{\left(\alpha \right)}{u}_0^{j-1}+g_{-(N-1)}^{\left(\alpha \right)}{u}_N^j,g_2^{\left(\alpha \right)}{u}_0^{j-1}+g_{-(N-2)}^{\left(\alpha \right)}{u}_N^j,\cdots ,g_{N-1}^{\left(\alpha \right)}{u}_0^{j-1}+g_{-1}^{\left(\alpha \right)}{u}_N^j]}^{\top }$, we derive the following matrix form of the semi-implicit difference scheme:

$$\left(I-\frac{\Delta t}{h^{\alpha }}{D}^{(j-1)}A\right)u^j=u^{j-1}+\Delta t{F}^j,$$

(38)

where I is an identity matrix of suitable order, $F^j=\frac{1}{h^{\alpha }}{D}^{(j-1)}{v}^{j-1}+f^j$, $u^j={[u_1^j,u_2^j,\cdots ,u_N^j]}^{\top }$,

$$D^{(j-1)}=\left(\begin{array}{c}a\left(u_1^{j-1}\right)\\ & a\left(u_2^{j-1}\right)\\ & & \ddots \\ & & & a\left(u_{N-1}^{j-1}\right)\end{array}\right)\in {\mathbb{R}}^{(N-1)\times (N-1)},$$

is a diagonal matrix and

$$A=\left(\begin{array}{cccc}{g}_0^{\left(\alpha \right)}& g_{-1}^{\left(\alpha \right)}& \cdots & g_{-(N-2)}^{\left(\alpha \right)}\\ g_1^{\left(\alpha \right)}& g_0^{\left(\alpha \right)}& \cdots & g_{-(N-3)}^{\left(\alpha \right)}\\ ⋮& ⋮& & ⋮\\ g_{N-2}^{\left(\alpha \right)}& g_{N-3}^{\left(\alpha \right)}& \cdots & g_0^{\left(\alpha \right)}\end{array}\right)$$

is a symmetric negative definite Toeplitz matrix [3,35], which can be represented by a vector containing the elements of the first column, i.e., its storage requirement is of $\mathcal{O}\left(N\right)$.

4.2. Two-Dimensional Coefficient Matrix

The structure of the two-dimensional coefficient matrix is more complex than in the one-dimensional case. By writing Equation (26) in the following equivalent form

$$\frac{u^k-u^{k-1}}{\Delta t}={\tilde{D}}^{(k-1)}\left[W_x^{\left(\alpha \right)}\otimes I_y+I_x\otimes W_y^{\left(\alpha \right)}\right]u^k+f_{i,j}^k,$$

(39)

in which $u^k={[u_{11}^k,\cdots ,u_{N_x-1,1}^k,u_{1,2}^k,\cdots ,u_{N_x-1,2}^k,\cdots ,u_{1,N_y-1}^k,\cdots ,u_{N_x-1,N_y-1}^k]}^{\top }$, we can derive the discretized linear system as follows:

$$\left(I-\frac{\Delta t}{h^{\alpha }}{\tilde{D}}^{(k-1)}\tilde{A}\right)u^k=u^{k-1}+\Delta t{f}^k,$$

(40)

where $\tilde{A}=W_x^{\left(\alpha \right)}\otimes I_y+W_y^{\left(\alpha \right)}\otimes I_x$, $I_q(q\in \{x,y\})$ is also an identity matrix of suitable order and $f^k=[f_{11}^k,\cdots ,f_{N_x-1,1}^k,f_{12}^k,\cdots ,f_{N_x-1,2}^k,\cdots ,f_{1,N_y-1}^k,\cdots ,$$f_{N_x-1,N_y-1}^k{]}^{\top }$ and

$${\tilde{D}}^{(k-1)}=\left(\begin{array}{c}a\left(u_{11}^{k-1}\right)\\ & \ddots \\ & & a\left(u_{N_x-1,1}^{k-1}\right)\\ & & & \ddots \\ & & & & a\left(u_{1,N_y-1}^{k-1}\right)\\ & & & & & \ddots \\ & & & & & & a\left(u_{N_x-1,N_y-1}^{k-1}\right)\end{array}\right).$$

is a diagonal matrix. In order to simplify the calculation, we choose $h=h_x=h_y$ (i.e., $N_x=N_y=N$) and thus $W_x^{\alpha }=W_y^{\left(\alpha \right)}=W^{\left(\alpha \right)}=A$.

According to Lemma 2, $W^{\left(\alpha \right)}$ is a symmetric Toeplitz matrix of order $(N-1)$. We convert $W^{\left(\alpha \right)}$ into a $\tau$-matrix defined as follows:

$$\tau \left(W^{\left(\alpha \right)}\right)=W^{\left(\alpha \right)}-H_{N-1},$$

where $H_{N-1}$ is a Hankel matrix with constant entries along antidiagonals, namely $\left[g_2^{\left(\alpha \right)},g_3^{\left(\alpha \right)},\cdots ,g_{N-2}^{\left(\alpha \right)},0,0,0,g_{N-2}^{\left(\alpha \right)},\cdots ,g_3^{\left(\alpha \right)},g_2^{\left(\alpha \right)}\right].$

Notice that the $\tau$-matrix can be diagonalized by the sine transform matrix $S_{N-1}$ [39] as follows:

$$\tau \left(W^{\left(\alpha \right)}\right)=S_{N-1}\Lambda S_{N-1},$$

(41)

where $\Lambda$ is the diagonal matrix consisting of the eigenvalues of $\tau \left(W^{\left(\alpha \right)}\right)$. The entries of $S_{N-1}$ are given by the following:

$${\left[S_{N-1}\right]}_{i,j}=\sqrt{\frac{2}{N}}sin\left(\frac{\pi ij}{N}\right),1\le i,j\le N-1.$$

A simple calculation shows that $S_{N-1}$ is a symmetric orthogonal matrix, that is,

$$S_{N-1}^{-1}=S_{N-1}^{\top }=S_{N-1}.$$

Upon replacing $W^{\left(\alpha \right)}$ with $\tau \left(W^{\left(\alpha \right)}\right)$ in the expression of $\tilde{A}$, and using the properties of the Kronecker product, we derive the following:

$$\begin{array}{cc}\hfill P& =\tau \left(W^{\left(\alpha \right)}\right)\otimes I+I\otimes \tau \left(W^{\left(\alpha \right)}\right)\hfill \\ \hfill & =\left(S_{N-1}\otimes I\right)\left[(\Lambda \otimes I)+(I\otimes \Lambda )\right]\left(S_{N-1}\otimes I\right).\hfill \end{array}$$

(42)

Finally, we introduce matrices $\mathcal{A}$ and $\mathcal{M}$ defined as follows:

$$\mathcal{A}=I-\Delta t{\tilde{D}}^{(k-1)}\tilde{A}\mathrm{and}{P}_{\tau }=I-\Delta t{d}^{(k-1)}P,$$

(43)

where constant $d^{(k-1)}$ is defined as follows:

$$d^{(k-1)}:=\frac{a\left(u_{11}^{k-1}\right)+\cdots +a\left(u_{1,N-1}^{k-1}\right)+\cdots +a\left(u_{N-1,1}^{k-1}\right)+\cdots +a\left(u_{N-1,N-1}^{k-1}\right)}{{(N-1)}^2}>0,$$

(44)

due to the condition (2).

Under these assumptions, we can prove the following property of the preconditioner $P_{\tau }$ that guarantees the invertibility of $P_{\tau }$:

Proposition 1.

The preconditioner $P_{\tau }$ defined in (43) is symmetric positive definite.

Proof. 

From Equation (41) and Lemma 3.2 in [30], we know that all the eigenvalues of $\tau \left(W^{\left(\alpha \right)}\right)$ are negative as are those of matrix $W^{\left(\alpha \right)}$; that is, the matrix P defined in (42) is negative definite. Recalling the definition of the preconditioner $P_{\tau }$ given by Equation (43), and that $\Delta t,d^{(k-1)}>0$, we can conclude that the preconditioner $P_{\tau }$ is symmetric positive definite. □

Henceforth, for clarity, we rewrite Equation (40) as the following linear system:

$$\mathcal{A}u=b,$$

(45)

where the right-hand side vector is $b=u^{k-1}+\Delta t{f}^k$. Since the coefficient matrix $\mathcal{A}$ is non-symmetric, we can solve it by the preconditioned BiCGSTAB method described below: (Algorithm 1).

Algorithm 1 The preconditioned BiCGSTAB method.

1: $x_0$ is an initial guess; $r_0=b-\mathcal{A}{u}_0$; 2: ${\overline{r}}_0$ is an arbitrary vector, such that 3: $({\overline{r}}_0,r_0)\ne 0$, e.g., ${\overline{r}}_0=r_0$; ${\rho }_0=\alpha ={\omega }_0=1$; 4: $v_0=p_0=0$; 5: for$i=1,2,\dots ,$do 6:     ${\rho }_i=({\overline{r}}_0,r_{i-1});\beta =({\rho }_i/{\rho }_{i-1})(\alpha /{\omega }_{i-1})$; 7:     $p_i=r_{i-1}+\beta (p_{i-1}-{\omega }_{i-1}{v}_{i-1})$; 8:     Solve y from $P_{\tau }y=p_i$; 9:     $v_i=\mathcal{A}y$ 10:     $\alpha ={\rho }_i/({\overline{r}}_0,v_i)$; 11:     $s=r_{i-1}-\alpha v_i$; 12:     Solve z from $P_{\tau }z=s$; 13:     $t=\mathcal{A}z$; 14:     ${\omega }_i=(P_{\tau }^{-1}t,P_{\tau }^{-1}s)/(P_{\tau }^{-1}t,P_{\tau }^{-1}t)$; 15:     $u_i=u_{i-1}+\alpha y+{\omega }_{i}z$; 16:     if $u_i$ is accurate enough then quit; 17:     $r_i=s-{\omega }_{i}t$ 18: end for

At each iteration of the preconditioned BiCGSTAB algorithm, two matrix-vector products need to be computed, i.e., $\mathcal{A}v$ and $P_{\tau }^{-1}v$, where v represents an arbitrary vector. The matrix–vector product

$$\mathcal{A}v=(I-\Delta t{\tilde{D}}^{(k-1)}\tilde{A})v=v-\Delta t{\tilde{D}}^{(k-1)}\left(\tilde{A}v\right),$$

(46)

is calculated as follows:

$$\begin{array}{cc}\hfill \tilde{A}v& =(W^{\left(\alpha \right)}\otimes I+I\otimes W^{\left(\alpha \right)})v\hfill \\ \hfill & =(W^{\left(\alpha \right)}\otimes I)v+(I\otimes W^{\left(\alpha \right)})v\hfill \\ \hfill & =W^{\left(\alpha \right)}·V+V·W^{\left(\alpha \right)},\hfill \end{array}$$

(47)

where $V=[v_1,v_2,\cdots ,v_{N-1}],v_i\in {\mathbb{R}}^{N-1}$ and $v=vec\left(V\right)$. Obviously, this problem is transformed into a matrix–vector operation with matrix $W^{\left(\alpha \right)}$ and vector v. On the one hand, the calculation $W^{\left(\alpha \right)}·V$ is expressed as $W^{\left(\alpha \right)}·[v_1,v_2,\cdots ,v_{N-1}]$ one by one and costs $\mathcal{O}(nlogn)$ arithmetic operations, where $n={(N-1)}^2$. On the other hand, the operation $V·W^{\left(\alpha \right)}$ is computed as $W^{\left(\alpha \right)}·V^{\top }$, and thus, it can be carried out at the same with $\mathcal{O}(nlogn)$ complexity. Similarly, since matrix P can be diagonalized, the preconditioner $P_{\tau }^{-1}$ can be diagonalized too, so the inverse of $P_{\tau }$ exists. The preconditioning operation $P_{\tau }^{-1}v$ is as follows:

$$\begin{array}{cc}\hfill P_{\tau }^{-1}v& ={\left[I-\Delta t{d}^{(k-1)}P\right]}^{-1}v\hfill \\ \hfill & =(S_{N-1}\otimes I){\left[I-\Delta t{d}^{(k-1)}(\Lambda \otimes I+I\otimes \Lambda )\right]}^{-1}(S_{N-1}\otimes I)v.\hfill \end{array}$$

In conclusion, each iteration of the preconditioned BiCGSTAB iterative algorithm to solve Equation (45) requires $\mathcal{O}(nlogn)$ arithmetic operations in total and can be implemented as a matrix-free manner.

It is worth noting that, although the preconditioner $P_{\tau }^{-1}$ is well defined according to Proposition 1, due to the presence of nonlinear diffusion coefficients and the nonsymmetry of the discretized linear systems (40), it is not simple to derive theoretically the eigenvalue distribution of the preconditioned matrix [33]. However, in the next section, we provide numerical evidence of the favorable clustering of the eigenvalues of the preconditioned matrices for some problems.

5. Numerical Experiments

In this section, two numerical experiments are presented to illustrate the efficiency of the proposed preconditioned semi-implicit difference scheme, and to support our theoretical findings. All experiments are performed in MATLAB R2019b on Intel(R) Core(TM) i5-10210U CPU @ 1.60∼2.11 GHz and 8.00 GB of RAM. On occasion, one can reduce non-homogeneous Dirichlet boundary conditions to homogeneous Dirichlet boundary conditions by using the following transformation:

$$u(x,t)=W(x,t)+V(x,t),$$

(48)

where ${\displaystyle V(x,t)=\varphi \left(t\right)+\frac{[\Phi (t)-\phi (t\left)\right]x}{L_x}}$ and $W(x,t)$ is an unknown function such that $W(0,t)=0,W(L_x,t)=0$.

Example 1.

For the first problem, we consider Equation (1) with $L_x=T=1$, $a\left(u\right)=u^2$, $\varphi \left(x\right)=x^2{(1-x)}^2$, $\phi \left(t\right)=\psi \left(t\right)=0$, and the source term $f(x,t)$ is computed as follows:

$$\begin{array}{cc}\hfill f(x,t)=& \alpha {(1+t)}^{\alpha -1}{x}^2{(1-x)}^2+\frac{a\left(u\right){(1+t)}^{\alpha }}{2cos(\pi \alpha /2)}\{\frac{\Gamma \left(3\right)}{\Gamma (3-\alpha )}\left[x^{2-\alpha }+{(1-x)}^{2-\alpha }\right]-\hfill \\ \hfill & \frac{2\Gamma \left(4\right)}{\Gamma (4-\alpha )}\left[x^{3-\alpha }+{(1-x)}^{3-\alpha }\right]+\frac{\Gamma \left(5\right)}{\Gamma (5-\alpha )}\left[x^{4-\alpha }+{(1-x)}^{4-\alpha }\right]\}.\hfill \end{array}$$

This problem can be viewed as a modification of [1,24], and its exact solution is $u(x,t)={(1+t)}^{\alpha }{x}^2{(1-x)}^2$. Denoting by $u_j^n$ the numerical solution, we define the $L_{\infty }$-norm errors and the $L_2$-norm errors as follows:

$$\begin{array}{c}{E}_{\infty }(\Delta t,h)=\underset{1\le j\le M}{max}\parallel u(ih,j\Delta t)-u_i^j{\parallel }_{\infty },E_2(\Delta t,h)=\underset{1\le j\le M}{max}\parallel u(ih,j\Delta t)-u_i^j\parallel ,\hfill \end{array}$$

respectively, and the space and time convergence orders as follows:

$$\begin{array}{c}rat{e}_1^t={log}_2\left(\frac{E_{\infty }(2\Delta t,h)}{E_{\infty }(\Delta t,h)}\right),rat{e}_1^x={log}_2\left(\frac{E_{\infty }(4h^2,2h)}{E_{\infty }(h^2,h)}\right),\hfill \\ rat{e}_2^t={log}_2\left(\frac{E_2(2\Delta t,h)}{E_2(\Delta t,h)}\right),rat{e}_2^x={log}_2\left(\frac{E_2(4h^2,2h)}{E_2(h^2,h)}\right).\hfill \end{array}$$

The data reported in

Table 1. The errors and temporal convergence orders when $N=512$.

$\mathit{\alpha }$	M	${\mathit{E}}_{\mathit{\infty }}(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_1^{\mathit{t}}$	${\mathit{E}}_2(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_2^{\mathit{t}}$
	8	1.6125e − 03	–	8.6746e − 04	–
	16	7.9547e − 04	1.0194	4.2894e − 04	1.0160
1.20	32	3.9507e − 04	1.0097	2.1329e − 04	1.0080
	64	1.9688e − 04	1.0048	1.0635e − 04	1.0040
	128	9.8288e − 05	1.0022	5.3109e − 05	1.0018
	8	4.7362e − 03	–	2.6365e − 03	–
	16	2.3256e − 03	1.0261	1.2984e − 03	1.0219
1.50	32	1.1525e − 03	1.0128	6.4439e − 04	1.0107
	64	5.7374e − 04	1.0063	3.2102e − 04	1.0053
	128	2.8628e − 04	1.0030	1.6023e − 04	1.0025
	8	1.2993e − 02	–	7.4538e − 03	–
	16	6.3112e − 03	1.0417	3.6319e − 03	1.0373
1.90	32	3.1120e − 03	1.0201	1.7935e − 03	1.0179
	64	1.5455e − 03	1.0098	8.9131e − 04	1.0088
	128	7.7025e − 04	1.0047	4.4437e − 04	1.0042

and

Table 2. The errors and spatial convergence orders when $M=N^2$.

$\mathit{\alpha }$	N	${\mathit{E}}_{\mathit{\infty }}(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_1^{\mathit{x}}$	${\mathit{E}}_2(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_2^{\mathit{x}}$
	4	1.1988e − 03	–	6.3875e − 04	–
	8	2.9694e − 04	2.0133	1.5674e − 04	2.0269
1.20	16	7.4054e − 05	2.0035	3.9116e − 05	2.0025
	32	1.8502e − 05	2.0009	9.7746e − 06	2.0006
	64	4.6248e − 06	2.0002	2.4434e − 06	2.0001
	4	3.5560e − 03	–	1.9379e − 03	–
	8	8.7005e − 04	2.0489	4.7737e − 04	2.0213
1.50	16	2.1633e − 04	2.0120	1.1889e − 04	2.0055
	32	5.4010e − 05	2.0030	2.9695e − 05	2.0013
	64	1.3498e − 05	2.0007	7.4219e − 06	2.0004
	4	1.0923e − 02	–	6.1965e − 03	–
	8	2.6106e − 03	2.0649	1.5174e − 03	1.9285
1.90	16	6.4524e − 04	2.0165	3.7556e − 04	1.9873
	32	1.6085e − 04	2.0041	9.3655e − 05	1.9954
	64	4.0184e − 05	2.0010	2.3399e − 05	1.9992

refer to the convergence results obtained by the fast preconditioned iterative algorithm and by a direct method (i.e., the backslash “∖” operator in MATLAB) for solving the non-symmetric discretized linear system, respectively. Because of the small size of the linear systems presented in Table 2, the cost of direct solvers is affordable [16]. Table 1 presents the maximum errors and the corresponding convergence orders in time of the semi-implicit difference scheme when $N=512$, for $\alpha =1.20,1.50,1.90$, respectively. The numerical results show a convergence order approximately equal to 1 in time for the semi-implicit difference scheme. Meanwhile, the experiments reported in Table 2 confirm that the semi-implicit difference scheme is convergent with second-order accuracy in space. These results are consistent with the theoretical analysis presented in Section 2.

Example 2.

For the second problem, we consider the two-dimensional nonlinear Riesz space-fractional diffusion problem, which is modified from [34] as follows:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial u(x,y,t)}{\partial t}=a\left(u\right)\left[\frac{{\partial u(x,y,t)}^{\alpha }}{{\partial \left|x\right|}^{\alpha }}+\frac{{\partial u(x,y,t)}^{\alpha }}{{\partial \left|y\right|}^{\alpha }}\right]+f(x,y,t),(x,y)\in \Omega ,0<t\le 1}\hfill \\ u(0,y,t)=u(1,y,t)=u(x,0,t)=u(x,1,t)=0,0\le t\le 1\hfill \\ u(x,y,0)=x^3{(1-x)}^3{y}^3{(1-y)}^3,(x,y)\in \Omega ,\hfill \end{array}\right.$$

(49)

where $\Omega =(0,1)\times (0,1)$, $a\left(u\right)=\left|u\right|+10$, and

$$\begin{array}{cc}\hfill f(x,y,t)=& -u(x,y,t)+\frac{a\left(u\right)}{2cos\left(\frac{\alpha \pi }{2}\right)}{e}^{-t}\{y^3{(1-y)}^3[g_1\left(x\right)-3g_2\left(x\right)+3g_3\left(x\right)-g_4\left(x\right)]\hfill \\ \hfill & +x^3{(1-x)}^3[g_1\left(y\right)-3g_2\left(y\right)+3g_3\left(y\right)-g_4\left(y\right)]\}\hfill \end{array}$$

(50)

where

$$\begin{array}{c}\hfill g_i\left(x\right)=\frac{\Gamma (i+3)}{\Gamma (i+3-\alpha )}\left[x^{i+2-\alpha }+{(1-x)}^{i+2-\alpha }\right],i=1,2,3,4.\end{array}$$

(51)

The exact solution for this problem is $u(x,y,t)=e^{-t}{x}^3{(1-x)}^3{y}^3{(1-y)}^3$.

The data reported in

Table 3. The errors and temporal convergence orders when $N=32$.

$\mathit{\alpha }$	M	${\mathit{E}}_{\mathit{\infty }}(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_1^{\mathit{t}}$	${\mathit{E}}_2(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_2^{\mathit{t}}$
1.20	4	1.8489e − 05	–	6.3044e − 06	–
	8	9.4453e − 06	0.9690	3.2208e − 06	0.9689
	16	4.7734e − 06	0.9846	1.6277e − 06	0.9846
	32	2.3998e − 06	0.9921	8.1828e − 07	0.9922
1.50	4	1.8488e − 05	–	6.3044e − 06	–
	8	9.4452e − 06	0.9689	3.2207e − 06	0.9690
	16	4.7732e − 06	0.9846	1.6276e − 06	0.9846
	32	2.3996e − 06	0.9922	8.1822e − 07	0.9922
1.90	4	1.8489e − 05	–	6.3045e − 06	–
	8	9.4456e − 06	0.9690	3.2208e − 06	0.9690
	16	4.7737e − 06	0.9845	1.6278e − 06	0.9845
	32	2.4001e − 06	0.9920	8.1838e − 07	0.9921

and

Table 4. The errors and spatial convergence orders when $M=N^2$.

$\mathit{\alpha }$	N	${\mathit{E}}_{\mathit{\infty }}(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_1^{\mathit{x}}$	${\mathit{E}}_2(\mathbf{\Delta }\mathit{t},\mathit{h})$	${\mathit{rate}}_2^{\mathit{x}}$
1.20	4	4.7734e − 06	–	1.6181e − 06	–
	8	1.2036e − 06	1.9877	4.1036e − 07	1.9793
	16	3.0226e − 07	1.9943	1.0299e − 07	1.9943
	32	7.6388e − 08	1.9844	2.6021e − 08	1.9855
1.50	4	4.7732e − 06	–	1.6180e − 06	–
	8	1.2034e − 06	1.9878	4.1030e − 07	1.9795
	16	3.0205e − 07	1.9943	1.0298e − 07	1.9943
	32	7.6178e − 08	1.9873	2.5954e − 08	1.9883
1.90	4	4.7738e − 06	–	1.6182e − 06	–
	8	1.2039e − 06	1.9874	4.1046e − 07	1.9791
	16	3.0264e − 07	1.9920	1.0314e − 07	1.9926
	32	7.6769e − 08	1.9790	2.6120e − 08	1.9814

are convergence results obtained with the preconditioned iterative algorithm. Table 3 presents the maximum errors and the corresponding convergence orders in time of the semi-implicit difference scheme when $N=32$, for $\alpha =1.20,1.50,1.90$, respectively. The numerical results exhibit a convergence order approximately equal to 1 in time for the semi-implicit scheme. Meanwhile, Table 4 shows that the semi-implicit scheme is convergent with second-order accuracy in space. These results are consistent with the theoretical analysis presented in Section 3.

Table 5. Total number of iterations and CPU time (in seconds) when $M=N^2$.

$\mathit{\alpha }$	N	Direct	I		${\mathit{P}}_{\mathit{\tau }}$
		CPU	Iter	CPU	Iter	CPU
	16	0.026	15.0	0.041	4.0	0.021
	32	1.635	24.0	0.350	4.0	0.094
1.20	64	130.982	37.0	2.354	5.0	0.532
	128	23,599.059	50.4	41.207	5.0	5.248
	256	out of memory	59.9	844.715	5.0	38.040
	16	0.024	18.0	0.044	4.0	0.023
	32	1.614	33.0	0.479	4.0	0.094
1.60	64	117.905	59.4	3.764	4.0	0.430
	128	24,804.137	95.1	78.685	4.8	4.962
	256	out of memory	154.3	706.824	5.0	38.935
	16	0.020	19.0	0.050	3.0	0.016
	32	1.562	36.0	0.522	4.0	0.096
1.70	64	751.896	64.0	4.065	4.0	0.428
	128	32,105.020	109.4	89.902	4.0	4.185
	256	out of memory	188.2	858.853	4.0	30.393
	16	0.021	20.0	0.051	3.0	0.017
	32	1.689	43.7	0.636	3.0	0.071
1.90	64	263.742	75.3	4.785	4.0	0.426
	128	22,645.044	144.7	387.991	4.0	4.211
	256	out of memory	273.0	5633.458	4.0	30.324

shows comparative experiments between a direct method and two fast iterative algorithms for the solution of the discretized linear system at each step of the implicit difference scheme for $\alpha =1.20,1.60,1.70,1.90$. It presents the total number of iterations (Iter) and the solution time (CPU in seconds) when $M=N^2$. Obviously, we can see from the table that the two iterative algorithms are significantly faster than a direct algorithm, motivating our quest for fast solvers. As seen from Figure 1, the proposed preconditioner $P_{\tau }$ can be indeed efficient to accelerate the convergence of the preconditioned BiCGSTAB method in terms of the clustering eigenvalues, i.e., most eigenvalues of the preconditioned matrices increasingly gather around 1 when we take $\alpha \to 2$.

Since it is proved that the convergence order of the proposed semi-implicit different scheme is $\mathcal{O}(\Delta t+h^2)$, the Richardson extrapolation can be used to improve the temporal convergence order. For the sake of clarity, we do the following:

(1)

We apply the numerical scheme (24) with time step $\Delta t$, and we compute the solution $U_{ij}^n(h_x,h_y,\Delta t)$. According to the error analysis of the method, the following holds:

$$u(x_i,y_j,t_n)-U_{i,j}^n(h_x,h_y,\Delta t)=\mathcal{O}(h_x^2+h_y^2)+c\Delta t+\mathcal{O}(\Delta t^2),$$

(52)

(2)

Another run with the scheme (24) is carried out using $\frac{\Delta t}{2}$, and a new solution $U_{ij}^{2n}\left(h_x,h_y,\frac{\Delta t}{2}\right)$ is computed on the fine temporal grid. According to the error analysis, we have the following:

$$u(x_i,y_j,t_n)-U_{i,j}^{2n}\left(h_x,h_y,\frac{\Delta t}{2}\right)=\mathcal{O}(h_x^2+h_y^2)+c\frac{\Delta t}{2}+\mathcal{O}\left({\left(\frac{\Delta t}{2}\right)}^2\right),$$

(53)

(3)

We compute $\left(5.6\right)\times 2-\left(5.5\right)$ and define ${\tilde{U}}_{i,j}^n=2U_{i,j}^{2n}\left(h_x,h_y,\frac{\Delta t}{2}\right)-U_{i,j}^n(h_x,h_y,\frac{\Delta t}{2})$:

$$u(x_i,y_j,t_n)-{\tilde{U}}_{i,j}^n=\mathcal{O}(h_x^2+h_y^2+\Delta t^2).$$

(54)

This means that Richardson extrapolation indeed can improve the temporal convergence order of the proposed semi-implicit difference scheme from 1 to 2. However, it is always sensitive to the parameter selection used for solving the model problem and thus, it is numerically unstable (we omit the specific numerical results here).

6. Concluding Remarks

In this paper, a fast preconditioned semi-implicit difference scheme for a one-dimensional strongly NSFDE is first presented, which is also extended to solve a two-dimensional NSFDE. In addition, the stability and convergence properties of the new numerical scheme are proved. At the same time, a fast preconditioned iterative algorithm is designed to solve the discretized linear systems on each time step efficiently and rapidly without necessarily storing the coefficient matrix of the linear system at each time level. In the future work, we will design the fast numerical methods and explore their stability and convergence analyses for solving a class of strongly NSFDEs. Moreover, due to the nonlocality of the Riesz fractional derivative, developing the fast matrix-free iterative methods accelerated by suitable preconditioners (e.g., the recent parallel-in-time preconditioners [36,40]) for solving the discretized linear systems should be also very meaningful.