Block splitting preconditioner for time-space fractional diffusion equations

1.   Introduction

In recent decades, the application of the fractional diffusion equation (FDE) in many fields has become more and more extensive ^[1,2,3]. FDE has become an indispensable tool for describing many of phenomena in mechanics and physics ^[4,5,6]. It has also attracted a lot of attention in biology ^[7], finance ^[8], image processing ^[9] and many other fields. In addition, FDE also has some very unique properties. For example, the spatial fractional diffusion equation can provide a sufficient and accurate description of the abnormal diffusion process, while the classic second-order diffusion equation often cannot accurately simulate this process. Therefore, more and more scholars have begun to study such important equations, and have obtained many excellent results ^{[10,11,12,13,14]}.

For the fractional diffusion equations, analytical solutions are usually not available, so numerical approximate solutions have become the main method. However, due to the non-local nature of fractional operators, the use of simple discretization, even if it is implicit, will lead to unconditional instability ^[15,16]. In addition, most FDE numerical methods tend to generate the full coefficient matrix, which requires the computational cost of ${\mathrm O}({{n}^{3 }})$ and storage of ${\mathrm O}({{n}^{2}})$ at each time step, where $n$ represents the number of spatial grids.

Recently, Meerschaet and Tadjeran ^[15,16] proposed that the fractional diffusion equation was discretized by using the implicit finite difference scheme of the Gr$\ddot{u}$nwald formula with displacement to overcome the difficulty of stability. Their approach has proven to be unconditional and stable. Later, Wang and Sircar ^[17] found that the coefficient matrix obtained by the Meerschaet-Tadjeran method has a Toeplitz-like structure, which can be expressed as the combination of diagonal matrix and Toeplitz matrix. Since the Toeplitz matrix has many good properties, such as it can be completely determined by $2n-1$ elements in line 1 and column 1, the storage requirements are greatly reduced from ${\mathrm O}({n}^{2})$ to ${\mathrm O}({n})$. In addition, as a special Toeplitz matrix, circulant matrix can be diagonalized by fast Fourier transform (FFT), while a Toeplitz matrix of $n$ order can be extended to a circular matrix of $2n$ order, so the matrix-vector multiplication for the Toeplitz-like can be obtained in $\text{O} (n\log n)$ operations by the FFT ^[18,19,20]. With this advantage, some scholars used the Conjugate Gradient Normal Residual (CGNR) method ^[21] to solve the linear system discretized by Meerschaet-Tadjeran method. Because its structure is similar to Toeplitz, the cost of each iteration of CGNR method is $\text{O} (n\log n)$. However, the numerical results show that the convergence effect of the CGNR method is good only when the diffusion coefficient function is small. To solve the above problems, Pang and Sun ^[22] proposed to use the multigrid method to solve the discretized FDE system obtained by the Meerschaet-Tadjeran method. This algorithm can keep the calculation cost of each iteration to $\text{O}(n\log n)$. The numerical results show that the multigrid method converges quickly even under the ill-condition cases. Although in very simple cases, the linear convergence of the multigrid method has not been proved in theoretical. The fast algorithm based on FFT and preprocessing technology has developed rapidly and constructing accelerated iteration of preprocessing has become a common method to solve these linear systems ^{[10,23,24,26,28,29]}.

2.   Discretization of time-space FDE

In this paper, we consider the time-space fractional diffusion equation as follows:

where $0 < \alpha < 1$, $1 < \beta < 2$ is the order of the fractional derivative, $f(x, t)$ is the source term, and the diffusion coefficient ${{e}_{1}}, {{ e}_{2}} > 0$, the definition of Caputo fractional derivative as follows:

$\Gamma (\cdot)$ is the gamma function, ${}_{a}D_{x}^{\beta }$ and ${}_{x}D_{b}^{\beta }$ are the left and right Riemann-Liouville fractional derivatives, respectively

Fix two positive integers $N$ and $M$, and divide the space domain $[0, L]$ and the time domain $[0, T]$ as follows:

Then use the shifted Gr$\ddot{u}$nwald approximations to discretize the left and right fractional derivatives in the space ^[15,16]

The alternate fractional binomial coefficient $g_{k}^{(\beta)}$ is given as follows:

and has the following properties ^[30]:

Then put the formula (2.2) into the Eq (2.1) to get the semi-discrete format in matrix-vector form

where $\mathbf{u}(t) = {{[{{u}_{1}}, {{u}_{1}}, \ldots, {{u}_{N-1}}]}}^{T}$, $_{0}^{C}\mathcal{D}_{t}^{\alpha}u(t) = {{[_{0}^{C}\mathcal{D}_{t}^{\alpha }{{u}_{1}}, _{0}^{C}\mathcal{D}_{t}^{\alpha }{{u}_{2}}, \ldots, _{0}^{C}\mathcal{D}_{t}^{\alpha }{{u}_{N-1}}]}^{T}}$, $\mathbf{f}(t) = {{[{{f}_{1}}, {{f}_{2}}, \ldots, {{f}_{N-1}}]}}^{T}$ with ${{f}_{i}} = f({{x}_{i}}, t) \ (0\le i\le N)$, $K = {{e}_{1}}{{T}_{\beta }}+{{e}_{2}}T_{\beta}^{T}$, Toeplitz matrix ${T}_{\beta }$ is given as follows:

On the time step, let $u_{i}^{j}\approx u({{x}_{i}}, {{t}_{j}})$ be the approximate solution. By using the $L2-{{1}_{\sigma}}$ formula ^[33], the temporal fractional derivative $_{0}^{C}\mathcal{D}_{t}^{\alpha} u(x, t)$ can be discretized as:

with $\sigma = 1-\frac{\alpha }{2}$, and for $j = 0$, $c_{0}^{(\alpha, \sigma)} = \frac{{{\tau }^{-\alpha }}}{\Gamma (2-\alpha)}a_{0}^{(\alpha, \sigma)}$, for $j\ge 1$,

with

Bring (2.7) into (2.5), the following discrete format can be obtained

with initial conditions $u_{i}^{0} = u_{0}\left(x_{i}\right)$ $(0 \leq i \leq N)$, where ${{\mathbf{u}}^{j+\sigma }} = \sigma {{\mathbf{u}}^{j+1}}+(1-\sigma){{\mathbf{u}}^{j}}$, ${{\mathbf{u}}^{j}} = {{\left[u_{1}^{j}, u_{2}^{j}, \ldots, u_{N-1}^{j} \right]}^{T}}$, ${{f}^{\ j+\sigma }} = {{\left[f_{1}^{\ j+\sigma }, f_{2}^{\ j+\sigma }, \ldots, f_{N-1}^{\ j+\sigma } \right]}^{T}}$ and $f_{i}^{\ j+\sigma } = f\left({{x}_{i}}, {{t}_{j+\sigma }} \right)$ $(0\le i\le N)$.

Let $\mathbf{0}$ and $I$ represent the zero matrix and the identity matrix, respectively. ${{A}_{0}} = {{h}^{\beta }}c_{0}^{(\alpha, \sigma)}I+\sigma K$, ${{\mathbf{b}}_{0}} = B{{\mathbf{u}}^{0}}+{{h}^{\beta }}{{\mathbf{f}}^{\ \sigma }},$

Let $v_{j}^{(\alpha, \sigma)} = \frac{\tau^{-\alpha}}{\Gamma(2-\alpha)}\left(a_{j}^{(\alpha, \sigma)}-b_{j}^{(\alpha, \sigma)}\right)$, then

Finally, put it into formula (2.8) to get

where $\mathbf{b} = {{\left[{{\mathbf{b}}_{1}}, {{\mathbf{b}}_{2}}, \ldots, {{\mathbf{b}}_{M-1}} \right]}^{T }}$,

For the above discrete linear system, we introduce the Kronecker product, and the Eq (2.9) can be expressed equivalently as ^[25]:

where

According to the above steps, Eq (2.1) is discretized a block lower triangular linear system (2.10).

3.   Two-Step Split (TSS) iteration method

In this section, we will be working on the block lower triangle Toeplitz linear system (2.10)

Let

then $\tilde{W} = {{T}_{1}}+{{T}_{2}}$ constitutes a split of the coefficient matrix $\tilde{W}$. The matrix $T_{1}$ is a block bi-diagonal matrix, with each block is a Toeplitz matrix. $T_{2}$ is a block lower triangular Toeplitz matrix, the form is the same as the matrix $\tilde{W}$. Next, we use the matrix splitting iteration method TSS to deal with the system (2.10).

(The TSS Iteration Method.) Given an initial guess ${{u}^{(0)}}\in {{\mathbb{C}}^{n}}$, for $k = 0, 1, 2, \ldots$, until the iteration sequence $\{{{u}^{(k)}}\}\in {{\mathbb{C}}^{n}}$ converges, compute the next iteration $\{{{u}^{(k+1)}}\}\in {{\mathbb{C}}^{n}}$ according to the following procedure

where $\gamma$ is a prescribed positive constant.

Remark 3.1. In the above TSS iterative method, the main amount of calculation comes from the product of matrix $(\gamma I+{{T}_{1} })$, $(\gamma I+{{T}_{2}})$ and vector. Since matrix $T_{1}$ is a block bi-diagonal matrix and each block is a Toeplitz matrix, matrix $T_{2}$ is the Toeplitz matrix of the lower triangle of the block, and each block is also a Toeplitz matrix, so we may use fast Fourier transform (FFT) to calculate them.

Next, we will discuss the convergence domain of the TSS iterative method and its optimal parameter.

By straightforward computations, the two half-iterates in Eq (3.3) can be integrated into a standard iteration form

where

Next, we will show that the TSS method is unconditionally convergent, and give the optimal parameters.

Theorem 3.1. Let $\tilde{W} = {T}_{1}+{T}_{2}$ form a matrix split of $\tilde{W}$, where matrix ${T}_{1}$, ${T}_{2}$ is defined by (3.1) and (3.2). Let $L(\gamma)$ be the iterative matrix of the TSS iterative method, then

and the spectral radius $\rho (L(\gamma))$ of $L(\gamma)$ satisfies

where

${{\lambda }_{\min }}$, ${{\lambda }_{\max }}$ denote the minimum and maximum eigenvalues of the matrix ${T}_{1}$, ${{\hat{\lambda }}_{\min }}$ and ${{\hat{\lambda }}_{\max }}$ represent the minimum and maximum eigenvalues of the matrix ${T}_{2}$, ${{t}_{1}}$ and ${{t}_{2}}$ represent an eigenvalue of matrices ${T}_{1}$ and ${{t}_{2}}$, respectively.

Proof. Since

Let

then the matrix $L(\gamma)$ is similar to $\hat{L}(\gamma)$, so

For ${{T}_{1}} = \tilde{B}\otimes {{e}_{1}}{{T}_{\beta }}$, ${T}_ {\beta }$ is a strictly diagonally dominant matrix, and the diagonal elements are greater than 0. So the eigenvalues of ${e}_{1}{T}_{\beta }$ are all greater than 0. In addition, since all the eigenvalues $\sigma$ of $\tilde{B}$ are positive, according to the properties of Kronecker product on eigenvalues, we can get $\lambda ({{T}_{1}}) = \lambda (\tilde{B}\otimes {{e}_{1}}{{T}_{\beta }}) > 0$, where $\lambda (.)$ indicates an eigenvalue.

In the same way, we can get $\lambda (\tilde{B}\otimes {{e}_{2}}T_{\beta }^{T}) > 0$. Because ${{T}_{2}} = {{h}^{\beta }}(\tilde{A}\otimes I)+\tilde{B}\otimes {{e}_{2}}T_ {\beta }^{T}$ is the lower triangular matrix of the block, and each block matrix on the diagonal is ${{h}^{\beta }}c_{0}^{(\alpha, \sigma)} I+\sigma {{e}_{2}}T_{\beta }^{T}$ is also a Toeplitz matrix with strictly diagonally dominant matrix, so the eigenvalue is greater than 0, $\lambda (T_{2}) > 0$ is established.

Since $\gamma > 0$, $t_{1} > 0$ and $t_{2} > 0$, so

That is $\varphi (\gamma) < 1$, the TSS iteration method converges unconditionally.

Because the convergence rate of TSS iteration is limited by $\varphi (\gamma)$, which depends not only on the spectral radius of matrix ${{T}_{1}}$ and ${{T}_{2}}$ but also on the parameter $\gamma$. The following lemma will give the optimal parameter ${{\gamma }^{*}}$ of the TSS iterative method.

Lemma 3.2^[26]. Let ${{\lambda }_{*}} = \sqrt{{{\lambda }_{\min }}{{\lambda }_{\max }}}$, ${{\hat{\lambda }}_{*}} = \sqrt{{{{\hat{\lambda }}}_{\min }}~~{{{\hat{\lambda }}}_{\max }}}$. Then the optimal parameter ${{\gamma }_{*}}$ that minimizes the function $\varphi (\gamma)$ is determined in the following three cases:

Case 1: If ${{\lambda }_{\text{*}}} = {{\hat{\lambda }}_{\text{*}}}$, then it holds that${{\gamma }_{*}} = {{\lambda }_{\text{*}}} = {{\hat{\lambda }}_{\text{*}}}$, and

Case 2: If ${{\lambda }_{\text{*}}} < {{\hat{\lambda }}_{\text{*}}}$, then

Case 2.1: for ${{\lambda }_{\min }}\ge {{\hat{\lambda }}_{\min }}$, or for ${{\hat{\lambda }}_{\max }} > {{\lambda }_{\max }}$, ${{\hat{\lambda }}_{\min }} > {{\lambda }_{\max }}$ and $\frac{{{\lambda }_{\max }}~~+{{{\hat{\lambda }}}_{\max }}}{{{\lambda }_{\min }}~~+{{{\hat{\lambda }}}_{\min }}}\ge \sqrt{\frac{{{\lambda }_{\max }}~~{{{\hat{\lambda }}}_{\max }}}{{{\lambda }_{\min }}~~{{{\hat{\lambda }}}_{\min }}}}$, it holds that ${{\gamma }_{*}} = {{\lambda }_{\text{*}}}$, and

Case 2.2: for ${{\lambda }_{\max }}\ge {{\hat{\lambda }}_{\max }}$, or for ${{\hat{\lambda }}_{\max }} > {{\lambda }_{\max }}$, ${{\hat{\lambda }}_{\min }} > {{\lambda }_{\max }}$ and $\frac{{{\lambda }_{\min }}~~+{{{\hat{\lambda }}}_{\min }}}{{{\lambda }_{\max }}~~+{{{\hat{\lambda }}}_{\max }}}\ge \sqrt{\frac{{{\lambda }_{\min }}~~{{{\hat{\lambda }}}_{\min }}}{{{\lambda }_{\max }}~~{{{\hat{\lambda }}}_{\max }}}}$, it holds that ${{\gamma }_{\text{*}}} = {{\hat{\lambda }}_{\text{*}}}$, and

Case 3: If ${{\hat{\lambda }}_{*}} < {{\lambda }_{*}}$, then

Case 3.1: for ${{\hat{\lambda }}_{\max }}\ge {{\lambda }_{\max }}$, or for ${{\lambda }_{\max }} > {{\hat{\lambda }}_{\max }}$, ${{\lambda }_{\max }} > {{\hat{\lambda }}_{\min }}$ and $\frac{{{\lambda }_{\min }}~~+{{{\hat{\lambda }}}_{\min }}}{{{\lambda }_{\max }}~~+{{{\hat{\lambda }}}_{\max }}}\ge \sqrt{\frac{{{\lambda }_{\min }}~~{{{\hat{\lambda }}}_{\min }}}{{{\lambda }_{\max }}~~{{{\hat{\lambda }}}_{\max }}}}$, it holds that ${{\gamma }_{\text{*}}} = {{\lambda }_{\text{*}}}$, and

Case 3.2: for ${{\hat{\lambda }}_{\min }}\ge {{\lambda }_{\min }}$, or for ${{\lambda }_{\min }} > {{\hat{\lambda }}_{\min }}$, ${{\lambda }_{\max }} > {{\hat{\lambda }}_{\max }}$ and $\frac{{{\lambda }_{\max }}~~+{{{\hat{\lambda }}}_{\max }}}{{{\lambda }_{\min }}~~+{{{\hat{\lambda }}}_{\min }}}\ge \sqrt{\frac{{{\lambda }_{\max }}~~{{{\hat{\lambda }}}_{\max }}}{{{\lambda }_{\min }}~~{{{\hat{\lambda }}}_{\min }}}}$, it holds that ${{\gamma }_{\text{*}}} = {{\hat{\lambda }}_{\text{*}}}$, and

4.   TSS preconditioner

From the TSS iterative method, we get the iterative matrix $M(\gamma) = \frac{1}{2\gamma }(\gamma I+{{T}_{1}})(\gamma I+{{T}_{2}})$ (see (3.5)) for the coefficient matrix $\tilde{W}$. By replacing the Toeplitz matrix ${{T}_{\beta }}$ in $M(\gamma)$ with the circulant matrix ${{S}_{\beta }}$ ^[34], we get the circulant preconditioner of the coefficient matrix $\tilde{W}$ of the linear system (2.10), denoted as

where

and the first columns of matrices ${{S}_{\beta }}$ and ${{S}_{\beta }}^{T}$ are respectively given by

Take the matrix $P(\gamma)$ as the preconditioner of the coefficient matrix $\tilde{W}$, to show its effectiveness, we need to consider the properties of the preconditioned matrix $P(\gamma)^{-1}\tilde{W}$. In the below discussion, let $\parallel \cdot \parallel$ be the infinite norm of the matrix.

Theorem 4.1. The preconditioner $P(\gamma) = \frac{1}{2\gamma }(\gamma I+{{S}_{1}})(\gamma I+{{S}_{2}})$ is reversible.

Proof. According to the properties of the binomial coefficient $g_{k}^{(\beta)}$ (see (2.4)), we have

By the Ger$\check{s}$gorin circle theorem^[38], we know that all eigenvalues of the circulant matrix ${{S}_{\beta }}$ and ${{S}_{\beta }}^{*}$ are within the open disc

so the circulant matrix ${{S}_{\beta }}$ and ${{S}_{\beta }}^{T}$ are strictly diagonally dominant. In addition, because all eigenvalues of $\tilde{B}$ are positive, by the properties of Kronecker product on eigenvalues, we can get $\lambda ({{S}_{1}}) > 0$, $\alpha I+{{S}_{1}}$ is reversible.

In the same way, $\lambda (\tilde{B}\otimes {{e}_{2}}T_{\beta }^{T}) > 0$. Because ${{S}_{2}} = {{h}^{\beta }}(\tilde{A}\otimes I)+\tilde{B}\otimes {{e}_{2}}S_ {\beta }^{T}$ is the lower triangular matrix of the block, and each block matrix on the diagonal is ${{h}^{\beta }}c_{0}^{(\alpha, \sigma)} I+\sigma {{e}_{2}}S_{\beta }^{T}$ is Toeplitz matrix with strictly diagonally dominant, so the eigenvalue is greater than 0, $\lambda (S_{2}) > 0$, $\alpha I+{{S}_{2}}$ is also reversible.

All in all, the preprocessing matrix $P(\gamma) = \frac{1}{2\gamma }(\gamma I+{{S}_{1}})(\gamma I+{{S}_{2} })$ is reversible.

Lemma 4.2^[35]. If the positive generating function $f$ of Toeplitz matrix ${{T}_{n}}\in {{\mathbb{C}}^{n\times n}}$ is in the Wiener class, ${{S}_{n}} \in {{\mathbb{C}}^{n\times n}}$ is a Strang's circulant preconditioner generated by ${T}_{n}$. Then, according to the equivalence of norm, for all $\varepsilon > 0$, there are positive integers ${N}^{'}$, ${M}^{'} > 0$, so that for all $n > {N}^{'}$, there are

with

Theorem 4.3. There exist matrices ${{E}_{L}}$, ${{F}_{L}}$, ${{E}_{R}}$ and ${{F}_{R}}$ such that

with

Proof. Since

and

According to the Lemma 4.2, let

with

Bring it into the above formula to get

with

and

Bring it into the above formula to get

with

and

Similarly

and

with

and

Bring it into the above formula to get

with

and

Theorem 4.4. There are matrices $\hat{E}$ and $\hat{F}$ such that $P{{(\gamma)}^{-1}}M(\gamma) = \hat{E}+\hat{ F}+I$, with

Proof.

with

and

Let $c = \max \{|c_{0}^{(\alpha, \sigma)}|, |c_{1}^{(\alpha, \sigma)}-c_{0}^{(\alpha, \sigma)}|, \ldots, |c_{M-2}^{(\alpha, \sigma)}-c_{M-3}^{(\alpha, \sigma)}|\}$, then

so

Theorem 4.5. There are matrices $E$ and $F$ such that $P{{(\gamma)}^{-1}}\tilde{W} = I+E(\gamma)+F(\gamma)$, with

Proof. Since

Let

then

5.   Numerical experiments

In this section, a numerical example used to show the performance of preconditioner ${{P}_{TSS}}(\gamma) = \frac{1}{2\gamma }(\gamma I+{S}_{1})(\gamma I+{S}_{2})$ (4.1) generated by the TSS iteration method. To illustrate the efficiency of ${{P}_{TSS}}(\gamma)$, the preconditioner ${{P}_{DTS}}(\gamma) = \frac{1 }{2\gamma }(\gamma I+{{h}^{\beta }}(\tilde{A}\otimes I))(\gamma I+\tilde{B}\otimes ({{e}_{1 }}{{S}_{\beta }}+{{e}_{2}}S_{\beta }^{T}))$ which generated by the diagonal and Toeplitz splitting (DTS) iteration method proposed by Bai ^[26] is tested. The Strang's circulant preconditioner ${{P}_{C}} = {{h}^{\beta }}(\tilde{A}\otimes I)+\tilde{B}\otimes ({{e} _{1}}{{S}_{\beta }}+{{e}_{2}}S_{\beta }^{T})$ which generated by directly replacing Toeplitz part of coefficient matrix with Strang's circulant matrix has been tested too. We choose the GMRES method to solve the Eq (2.1). The iteration terminated if the relative residual error satisfies $\frac{\parallel {{r}^{(k)}}{{\parallel }_{2}}}{\parallel {{r}^{(0)}}{{\parallel}_{2}}}\le {{10}^{-8}}$, where ${{r}^{(k)}}$ denotes the residual vector in the $kth$ iteration. The initial guess is chosen as the zero vector, and the experimental results include the CPU time and the number of iterations. All experiments are carried out via MATLAB 2020a on a Windows 10 (64 bit) PC.

Example 1^[25]. For Eq (2.1), consider the value of the fractional derivative $(\alpha, \beta) =$ $(0.1, 1.1)$, $(0.4, 1.7)$, $(0.7, 1.4)$, $(0.9, 1.9)$, the time domain is $[0, T] = [0, 1]$, and the space domain is $[0, L] = [0, 1]$, ${{e}_{1}} = 20, {{e}_{2}} = 0.02$, and source item is

where

the exact solution of the corresponding fractional diffusion equation is

Table 1 shows that compared with the preconditioner ${{P}_{DTS}}(\gamma)$, the preconditioner ${{P}_{TSS}}(\gamma)$ reduce the computational cost and number of iterations when we use the GMRES subspace iteration method to solve Example 1.

Figure 1 shows the eigenvalues distribution of the original coefficient matrix in Example 1 when ($(\alpha, \beta) = (0.7, 1.4)$, (N, M) = (33, 33). At this time, the eigenvalues are scattered in a large range. But after preprocessing, the eigenvalues gather obviously, as shown in Figure 2. Figure 2(a) on the left shows the result of the preconditioner generated by the TSS iterative method. On the right Figure 2(b) is the result of the preconditioner generated by the DTS iterative method. From the comparison results, the preconditioner generated by TSS makes the eigenvalues more concentrated.

6.   Conclusions

In this paper, we mainly researched the preconditioner of linear equations which are discretized from fractional diffusion equations. According to the two-Step Split (TSS) iteration method, we constructed a split iteration preconditioner that could reduce the computational cost and number of iterations when the Krylov subspace iteration method was used to solve this equation. Firstly, through the FDE discretization, we got the linear system $\tilde{W}\mathbf{u} = \mathbf{b}$ which contained the Kronecker product. Then, the coefficient matrix was split into two parts ${{T}_{1}}$, ${{T}_{ 2}}$, and introduced to a single parameter two-step split iteration method TSS. It was proved that the TSS iterative method was unconditionally convergent, the optimal parameter values were given, and the preconditioner $P(\gamma)$ of the coefficient matrix was constructed. Finally, to show the effectiveness of the preconditioner, theoretical analysis proved that preprocessed coefficient matrix could be expressed as the sum of an identity matrix, a low-rank matrix, and a small norm matrix. This showed that the preconditioner had a high degree of approximation to the coefficient matrix. It could be obtained by fast Fourier transform (FFT) in a small amount of calculation for the circulant matrix, calculating inversion, matrix-vector multiplication, etc. In the numerical experiment, a numerical example was used to demonstrate the effectiveness of the preconditioner proposed in this paper. From the analysis of the experimental results, the preconditioner we constructed could make the eigenvalue distribution of the coefficient matrix more concentrated than the preconditioner produced by the DTS iteration method. When using the Krylov subspace iteration method, such as the generalized minimal residual (GMRES) method, the number of iteration steps was reduced and the calculation speed was increased.

In addition, the proposed model is mainly for the constant-order fractional diffusion equations. Whether it can be extended to variable-order fractional diffusion problems ^{[42,43,44,45]} may also be considered in the future.

Acknowledgments

We would like to express our sincere thanks to the referees for insightful comments and invaluable suggestions that greatly improved the representation of this paper. The authors are partially supported by National Natural Science Foundation of China (11101071, 11271001).

Conflict of interest

The authors declare there is no conflicts of interest.