Approximate inverse preconditioners for linear systems arising from spatial balanced fractional diffusion equations

1.   Introduction

Fractional diffusion equations (FDEs) have been utilized to model anomalous diffusion phenomena in the real world, see for instance ^{[1,2,4,18,19,22,23,25]}. One of the main features of the fractional differential operator is nonlocality. It brings big challenge for finding the numerical solution of FDEs, as the coefficient matrix of the discretized FDEs is typically dense, which requires $O(N^3)$ of computational cost and $O(N^2)$ of memory storage if a direct solution method is employed, where $N$ is the number of unknowns. However, by making use of the Toeplitz-like structure of the coefficient matrices, many efficient algorithms have been developed, see for instance ^{[6,9,12,13,14,15,16,21]}.

In this paper, we consider the following initial-boundary value problem of spatial balanced FDEs ^[7,24]:

where $d_{\pm}(x, t) > 0$ are diffusion coefficients, $f(x, t)$ is the source term, and $\alpha, \beta$ are the fractional orders satisfying $\dfrac{1}{2} < \alpha, \beta < 1$. Here ${}_{{x_L}}^{{{\rm{RL}}}}{{\rm{D}}}_{x}^{{\gamma}}$ and ${}_{{x}}^{{{\rm{RL}}}}{{\rm{D}}}_{x_R}^{{\gamma}}$ denote the left-sided and right-sided Riemann-Liouville fractional derivatives for $0 < \gamma < 1$, respectively, and are defined by ^[17]

where $\Gamma(\cdot)$ is the Gamma function, while ${}_{{x_L\!}}^{{{\rm{C}}}}{{\rm{D}}}_{x}^{{\gamma}}$ and ${}_{{x}}^{{{\rm{C}}}}{{\rm{D}}}_{x_R}^{{\gamma}}$ denote the left-sided and right-sided Caputo fractional derivatives for $0 < \gamma < 1$, respectively, and are defined by ^[17]

The fractional differential operator in (1.1) is called the balanced central fractional derivative, which was studied in ^[24]. One advantage of such fractional differential operator is that its variational formulation has a symmetric bilinear form, which can greatly benefit theoretical investigation.

Recently, a finite volume approximation for the spatial balanced FDEs (1.1) is proposed in ^[7]. By applying a standard first-order difference scheme for the time derivative and a finite volume discretization scheme for the spatial balanced fractional differential operator, a series of systems of linear equations are generated, whose coefficient matrices share the form of the sum of a tridiagonal matrix and two Toeplitz-times-diagonal-times-Toeplitz matrices. One attractive feature of these coefficient matrices is that they are symmetric positive definite, so that the linear systems can be solved by CG method, in which the three-term recurrence can significantly reduce the computational and storage cost. However, due to the ill-conditioning of the coefficient matrices, CG method, when applied to solve the resulted linear systems, usually converges very slow. Therefore, preconditioners should be applied to improve the computational efficiency. In ^[7], the authors proposed two preconditioners: circulant preconditioner for the constant diffusion coefficient case and banded preconditioner for the variational diffusion coefficient case.

In this paper, we consider the approximate inverse preconditioners for the resulting linear systems arising from the finite volume discretization of the spatial balanced FDEs (1.1). Our preconditioner is based on the symmetric approximate inverse strategy studied in ^[11] and the Sherman-Morrison-Woodburg formula. Rigorous analysis shows that the preconditioned matrix can be written as the sum of the identity matrix, a small norm matrix, and a low-rank matrix. Therefore, the quick convergence of the CG method for solving the preconditioned linear systems is expected. Numerical examples, for both one-dimensional and two-dimensional cases, are given to demonstrate the robustness and effectiveness of the proposed preconditioner. We remark that our preconditioner can also be applied to another class of conservative balanced FDEs which was studied in ^[8,10].

The rest of the paper is organized as follows. In Section 2, we present the discretized linear systems of the balanced FDEs. Our new preconditioners are given in Section 3, and their properties are investigated in detail in Section 4. In Section 5, we carry out the numerical experiments to demonstrate the performance of the proposed preconditioners. Finally, we give the concluding remarks in Section 6.

2.   Discretization of the spatial balanced FDEs

Let $\Delta t = T/M_t$ be the time step where $M_t$ is a given positive integer. We define a temporal partition $t_{j} = j\Delta t$ for $j = 0, 1, 2, \dots, M_t$. The first-order time derivative in (1.1) can be discretized by the standard backward Euler scheme, and we obtain the following semidiscrete form:

for $j = 1, 2, \cdots, M_t$. Let $\Delta x = (x_R-x_L)/(N+1)$ be the size of the spatial grid where $N$ is a positive integer. We define a spatial partition $x_{i} = x_L+i\Delta x$ for $i = 0, 1, 2, \dots, N+1$, and denote by $x_{i-\frac{1}{2}} = \dfrac{x_{i-1}+x_{i}}{2}$ the midpoint of the interval $[x_{i-1}, x_{i}]$. Integrating both sides of (2.1) over $[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}]$ gives

Let $\mathcal{S}_{\Delta x}(x_L, x_R)$ be the space of continuous and piecewise-linear functions with respect to the spatial partition, and define the nodal basis functions $\phi_k(x)$ as

for $k = 1, 2, \dots, N$, and

The approximate solution $u_{\Delta x}(x, t_j)\in \mathcal{S}_{\Delta x}(x_L, x_R)$ can be expressed as

Therefore, the corresponding finite volume scheme leads to

which can be further approximated and we obtain ^[7]

where $\eta_{\alpha} = \dfrac{\Delta t}{\Gamma(2-\alpha)^2 \Delta x^{2\alpha}}$, $\eta_{\beta} = \dfrac{\Delta t}{\Gamma(2-\beta)^2 \Delta x^{2\beta}}$, and

for $\gamma = \alpha, \beta$, with

The initial value and boundary condition are

Collecting all components $i$ into a single matrix system, we obtain

where

are diagonal matrices, $M = \dfrac{1}{8}{\mathrm{tridiag}}(1, 6, 1)$, $u^{(j)} = \left[u_{1}^{(j)}, u_{2}^{(j)}, \ldots, u_{N}^{(j)}\right]^{ { \intercal }}$, $f^{(j)} = \left[f_{1}^{(j)}, f_{2}^{(j)}, \ldots, \right.$ $\left.f_{N}^{(j)}\right]^{ { \intercal }}$ with

and

with

The matrices $\tilde{G}_{\alpha}$, $\tilde{G}_{\beta}$ are $N$-by-$(N+1)$ Toeplitz matrices defined by

We remark that the entries of $\tilde{G}_{\alpha}$ and $\tilde{G}_{\beta}$ are independent on the time partition $t_j$.

We denote the coefficient matrix by $A^{(j)}$, that is,

As $M$ and $\tilde{D}_{+}^{(j)}, \tilde{D}_{-}^{(j)}$ are symmetric positive definite, we can see that the coefficient matrix $A^{(j)}$ is symmetric positive definite too.

For the entries of $\tilde{G}_{\alpha}$ and $\tilde{G}_{\beta}$, we have the following result, which is directly obtained from Lemma 2.3 in ^[15].

Lemma 2.1. Let $g_{k}^{(\gamma)}$ be defined by (2.5) with $\dfrac{1}{2} < \gamma < 1$. Then

${ \rm(1)}$ $g_{0}^{(\gamma)} > 0$, $g_{1}^{(\gamma)} < g_{2}^{(\gamma)} < \cdots < g_{k}^{(\gamma)} < \cdots < 0$;

${ \rm(2)}$ $\left|g_{k}^{(\gamma)}\right| < \dfrac{c_{\gamma}}{(k-1)^{\gamma+1}}$ for $k = 2, 3, \dots$, where $c_{\gamma} = \gamma(1-\gamma)$;

${ \rm(3)}$ $\lim\limits_{k\to\infty} \left|g_{k}^{(\gamma)}\right| = 0$, $\sum\limits_{k = 0}^{\infty} g_{k}^{(\gamma)} = 0$, and $\sum\limits_{k = 0}^{p} g_{k}^{(\gamma)} > 0$ for $p\geq 1$.

3.   The approximate inverse preconditioner

As the coefficient matrix is symmetric positive definite, we can apply CG method to solve the linear systems (2.6). In order to improve the performance and reliability of the CG method, preconditioning is necessarily to employed. In this section, we develop the approximate inverse preconditioner for the linear system (2.6).

For the convenience of analysis, we omit the superscript $"(j)"$ for the $j$th time step, that is, we denote the coefficient matrix (2.7) as

where $\tilde{D}_{+} = \tilde{D}_{+}^{(j)}$, $\tilde{D}_{-} = \tilde{D}_{-}^{(j)}$.

Let $g_\alpha$ and $G_\alpha$ be the first column and the last $N$ columns of $\tilde{G}_\alpha$, respectively, that is, $\tilde{G}_\alpha = \big[ g_\alpha, G_\alpha\big]$ where

Analogously, we denote by $g_\beta$ and $G_\beta$ the last column and the first $N$ columns of $\tilde{G}_\beta$, respectively, that is, $\tilde{G}_\beta = \left[G_\beta, g_\beta\right]$ where

Then we have

where

Therefore, $A$ can be written as

where

and

According to Sherman-Morrison-Woodburg Theorem, we have

In the following, we consider the approximation of $\hat{A}^{-1}$. To this end, we first define a matrix $\tilde{A}$ with

where $D_{+} = {\mathrm{diag}}\left(\left\{d_{+}\left(x_{k}\right)\right\}_{k = 1}^{N}\right)$ and $D_{-} = {\mathrm{diag}}\left(\left\{d_{-}\left(x_{k}\right)\right\}_{k = 1}^{N}\right)$.

Assume $d_{+}(x), d_{-}(x)\in C[x_L, x_R]$, then it is easy to see that for any $\epsilon > 0$, when the spatial gird $\Delta x$ is small enough, we have

Meanwhile, we learn from Lemma 2.1 that there exists $c_{\alpha}, c_{\beta} > 0$ such that $\|G_{\alpha}\|_{2}\leq c_{\alpha}, \|G_{\beta}\|_{2}\leq c_{\beta}$. Let

Then it holds that

This indicates that $\tilde{A}$ can be a good approximation of $\hat{A}$ when $\Delta x$ is very small. However, it is not easy to compute the inverse of $\tilde{A}$.

Motivated by the symmetric approximate inverse preconditioner proposed in ^[11], we consider the following approximations

where $e_i$ is the $i$-th column of the identity matrix $I$ and

which is symmetric positive definite. That is, we approximate the $i$-th column of $\tilde{A}^{-1/2}$ by the $i$-th column of $K_{i}^{-1/2}$. Then we propose the following symmetric approximate inverse preconditioner

Although $K_i$ is symmetric positive definite, it is not easy to compute the inverse of its square root. Hence we further approximate $K_i$ by circulant matrices whose inverse square roots can be easily obtained by FFT.

Let $C_{M}$, $C_{\alpha}$ and $C_{\beta}$ be Strang's circulant approximations ^[5] of $M$, $G_{\alpha}$ and $G_{\beta}$, respectively. Then we obtain the following preconditioner

where $C_{i} = C_{M}+\eta_{\alpha}d_{+}(x_{i}) C_{\alpha}C_{\alpha}^{ { \intercal }} +\eta_{\beta} d_{-}(x_{i}) C_{\beta}C_{\beta}^{ { \intercal }}$ is circulant matrix.

In order to make the preconditioner more practical, similar to the idea in ^[11,13], we utilize the interpolation technique. Let $\left\{\tilde{x}_{k}\right\}_{k = 1}^{\ell}$ be $\ell$ distinct interpolation points in $[x_L, x_R]$ with a small integer $\ell$ ($\ell\ll N$), and denote by $\lambda\triangleq\left(\lambda_{M}, \lambda_{\alpha}, \lambda_{\beta}\right)$, where $\lambda_{M}$, $\lambda_{\alpha}$, $\lambda_{\beta}$ are certain positive real numbers. Then we define the function

Let

be the piecewise-linear interpolation for $g_{\lambda}(x)$ based on the $\ell$ points $\left\{(\tilde{x}_k, g_{\lambda}(\tilde{x}_k))\right\}_{k = 1}^{\ell}$. Define

Then $\tilde{C}_{k}$ is symmetric positive definite and can be diagonalized by FFT, that is,

where $F$ is the Fourier matrix and $\tilde{\Lambda}_k$ is a diagonal matrix whose diagonal entries are the eigenvalues of $\tilde{C}_{k}$. By making use of the interpolation technique, we approximate $C_{i}^{-1/2}$ by

By substituting these approximations into $P_{2}^{-1}$, we obtain the practical preconditioner

where $\Phi_{k} = {\mathrm{diag}}\left(\phi_{k}(x_1), \phi_{k}(x_2), \dots, \phi_{k}(x_N)\right)$. Now applying $P_{3}^{-1}$ to any vector requires about $O(\ell N\log N)$ operations, which is acceptable for a moderate number $\ell$.

Finally, by substituting $\hat{A}^{-1}$ in the Sherman-Morrison-Woodburg formula (3.5) with $P_{3}^{-1}$, we obtain the preconditioner

We remark that both $P_3^{-1}$ and $P_4^{-1}$ can be taken as preconditioners. It is clear that implementing $P_{4}^{-1}$ requires more computational work than $P_{3}^{-1}$. However, we note that $P_4^{-1}$ is a rank-2 update of $P_3^{-1}$, and $g_\alpha$, $g_\beta$ are independent on $t_j$, which indicates that, during the implementation of the preconditioner $P_4^{-1}$ to the CG method, we can precompute $P_{3}^{-1} U$ ahead, as well as the inner product $U^{ { \intercal }}P_{3}^{-1} U$ and the inverse of the 2-by-2 matrix $I+U^{ { \intercal }}P_{3}^{-1} U$. Therefore, at each CG iteration with preconditioner $P_4^{-1}$, besides the matrix-vector product with $P_{3}^{-1}$, only two inner-products and two vector updates are needed. Therefore, it is expected that $P_4^{-1}$ may have better performance for the one dimensional problems.

4.   Analysis of the preconditioner

Since $P_4^{-1}$ is obtained by substituting $\hat{A}^{-1}$ with $P_{3}^{-1}$ in the expression of $A^{-1}$, the approximation property of $P_4^{-1}$ to $A^{-1}$ is dependent on how close $P_3^{-1}$ to $\hat{A}^{-1}$ will be. Therefore, in this section, we study the difference between $P_3^{-1}$ and $\hat{A}^{-1}$.

We first introduce the off-diagonal decay property ^[20], which is crucial for studying the circulant approximation of the Toeplitz matrix.

Definition 4.1. Let $A = \left[a_{i, j}\right]_{i, j\in\mathcal{I}}$ be a matrix, where the index set is $\mathcal{I} = \mathbb{Z}, \mathbb{N}$ or $\{1, 2, \dots, N\}$. We say $A$ belongs to the class $\mathcal{L}_{s}$ if

holds for $s > 1$ and some constant $c > 0$, and say $A$ belongs to the class $\mathcal{E}_{r}$ if

holds for $r > 0$ and some constant $c > 0$.

The following results hold for the off-diagonal decay matrix class $\mathcal{L}_{s}$ and $\mathcal{E}_{r}$ ^[3,11,15,20].

Lemma 4.1. Let $A = \left[a_{i, j}\right]_{i, j\in\mathcal{I}}$ be a nonsingular matrix, where the index set is $\mathcal{I} = \mathbb{Z}, \mathbb{N}$ or $\{1, 2, \dots, N\}$.

${ \rm(1)}$ If $A\in\mathcal{L}_{s}$ for some $s > 1$, then $A^{-1} \in \mathcal{L}_{s}$.

${ \rm(2)}$ If $A\in\mathcal{L}_{1+s_1}$ and $B\in\mathcal{L}_{1+s_2}$ are finite matrices, then $AB\in\mathcal{L}_{1+s}$, where $s = \min\{s_1, s_2\}$.

${ \rm(3)}$ If $A\in\mathcal{E}_{r_1}$ and $B\in\mathcal{E}_{r_2}$ are finite matrices, then $AB\in\mathcal{E}_{r}$ for some constant $0 < r < \min\{r_1, r_2\}$.

${ \rm(4)}$ Let $A$ be a banded finite matrix and Hermitian positive definite, and let $f$ be an analytic function on $[\lambda_{\min}(A), \lambda_{\max}(A)]$ and $f(\lambda)$ is real for real $\lambda$, where $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the minimal and maximal eigenvalues of $A$, respectively. Then $f(A)$ has the off-diagonal exponential decay property (4.2). In particular, let $f(x) = x^{-1}, x^{-1/2}, \, \mathit{{\text{and}}}\, x^{1/2}$, respectively, then $A^{-1}$, $A^{-1/2}$, and $A^{1/2}$ have the off-diagonal exponential decay property (4.2).

Assume that $d_{+}(x), d_{-}(x)\in C[x_L, x_R]$, it follows from Lemma 2.1 and Lemma 4.1 that we have

for $\frac{1}{2} < \alpha, \beta < 1$. Therefore, it holds that $\tilde{A}, \ \tilde{A}^{-1} \in \mathcal{L}_{1+\min\{\alpha, \beta\}}.$

4.1.   Approximation property of $P_{1}^{-1}$ to $\tilde{A}^{-1}$

Given an integer $m$, let $G_{\alpha, m}$ and $G_{\beta, m}$ be the $(2m+1)$-banded approximations of $G_{\alpha}$ and $G_{\beta}$, respectively, that is,

In the following discussion, we let

which can be regarded as some kind of banded approximation of the $K_i$ defined in (3.8). We remark that, in the actual applications, we still employ $K_i$ in (3.8) to construct the preconditioners.

Define

As $\tilde{A}_{m}$ and $K_{i}$ are banded matrices and symmetric positive definite, it follows from Lemma 4.1 that $\tilde{A}_{m}^{-1}$, $K_{i}^{-1}$ and $K_{i}^{-1/2}$ have off diagonal exponential decay property. Meanwhile, we have

and hence

For a given $\epsilon > 0$, it follows from Theorem 12 in ^[7] that there exists an integer $m_1 > 0$ such that for all $m\geq m_1$ we have

Therefore, it holds that

and

Now, we turn to estimate the upper bound of $\|P_{1}^{-1}-\tilde{A}_{m}^{-1}\|_{2}$. Since both $P_1^{-1}$ and $\tilde{A}_{m}^{-1}$ are symmetric, we have

For the first item, we have the following estimation.

Lemma 4.2. Let $K_{j}$ be defined by (4.3). Assume that $0 < d_{\min} \leq d_{+}(x), \, d_{-}(x) \leq d_{\max}$. Then for a given $\epsilon > 0$, there exist an integer $N_1 > 0$ and two positive constants $c_3, c_4$, such that

where $\Delta_{N_1}d_{+, j} = \max\limits_{j-N_1\leq k \leq j+N_1} |d_{+, k}-d_{+, j}|$, $\Delta_{N_1}d_{-, j} = \max\limits_{j-N_1\leq k \leq j+N_1} |d_{-, k}-d_{-, j}|$.

Proof. We have

As it was shown in Theorem 2.2 of ^[3], for a given $\epsilon > 0$, we can find a polynomial $p_{k}(t) = \sum_{\ell = 0}^{k} a_\ell t^\ell$ of degree $k$ such that

where $\big[\, \cdot\, \big]_{l, r}$ denotes the $(l, r)$-th entry of a matrix. Then we can write (4.7) as

where

Note that $p_{k}(K_i)$ is $(2km+1)$-banded matrix and $K_{i}^{-1/2}$ has the off diagonal exponential decay property, we know that $K_{i}^{-1/2}-p_{k}(K_i)$ also has the off diagonal exponential decay property. According to Lemma 4.1, we learn that $H_{ij}^{(1)}$ has the off diagonal exponential decay property. Hence, analogous to the proof of Lemma 3.8 in ^[11], we can find constants $\hat{c} > 0$ and $\hat{r} > 0$ such that

On the other hand, we denote the $j$-th column of $H_{ij}^{(1)}$ by $[h_{1, j}, h_{2, j}, \cdots, h_{N, j}]^{ { \intercal }}$, and let $\tilde{K}_{i} \triangleq K_{i}^{-1 / 2}-p_{k}\left(K_{i}\right)$. Observe that

Then based on (4.8) and the fact that $K_{j}^{-1/2}$ has off diagonal exponential property, we can show that there exists a constant $\hat{c}_1$ such that

Denote

From (4.10) and (4.11), we can show that there exists a constant $c_1 > 0$, such that

Analogously, we can show that there exists a constant $c_2 > 0$ such that

Now we turn to $H_{ij}^{(2)}$. It follows from the proof of Lemma 3.11 in ^[11] that $p_k(K_i)-p_k(K_j)$ has the form

where

It is easy to see that both $\tilde{H}_{ij}^{(2)}$ and $\hat{H}_{ij}^{(2)}$ have the off diagonal exponential decay property. As

there exists an integer $N_1$ such that

where $\Delta_{N_1}d_{+, j} = \max\limits_{j-N_1\leq k \leq j+N_1} |d_{+, k}-d_{+, j}|$ and $\Delta_{N1}d_{-, j} = \max\limits_{j-N_1\leq k \leq j+N_1} |d_{-, k}-d_{-, j}|$.

Let $c_3 = \tilde{c}$ and $c_4 = c_1+c_2+4d_{\max}$. It follows from (4.12), (4.13) and (4.15) that

□

Now we consider the second item in (4.6).

Lemma 4.3. Let $K_j$ and $\tilde{A}_{m}$ be defined by (4.3) and (4.4), respectively. Assume that $0 < d_{\min} \leq d_{+}(x), d_{-}(x) \leq d_{\max}$. Then for a given $\epsilon > 0$, there exists an integer $N_2 > 0$ and constants $c_{5}, c_{6} > 0$ such that

Proof. Let

then we have

On the one hand, observe that

we have

where

and

As $K_j^{-1}$ has off diagonal exponential decay property, similar to the derivation of (4.15), we can show that, given a $\epsilon > 0$, there exists a integer $\tilde{N}_1 > 0$ and a constant $\tilde{c}_{1} > 0$, such that

and

Therefore, there exists positive constants $\tilde{c}_2$ and $\tilde{c}_3$, such that

On the other hand, observe that

we have

Since $G_{\alpha}$ and $G_{\beta}$ have off diagonal decay property, then as shown in the proof of Lemma 4.7 in ^[15], given a $\epsilon > 0$, there exists a integer $\tilde{N}_2 > 0$ and a constant $\tilde{c}_{4} > 0$ such that

and

where $\Delta_{\tilde{N}_2}d_{+, k} = \max\limits_{k-\tilde{N}_2\leq l \leq k+\tilde{N}_2} |d_{+, l}-d_{+, k}|$ and $\Delta_{\tilde{N}_2}d_{-, k} = \max\limits_{k-\tilde{N}_2\leq l \leq k+\tilde{N}_2} |d_{-, l}-d_{-, k}|$. Therefore, there exists positive constants $\tilde{c}_5$ and $\tilde{c}_6$ such that

Now, let

It follows from (4.18) and (4.19) that, given a $\epsilon > 0$, there exists a integer $N_2 > 0$ and constants $c_{5}, c_{6} > 0$ such that

The proof is completed. □

In combination with (4.5), Lemma 4.2 and Lemma 4.3, we obtain the following theorem.

Theorem 4.1. Assume $0\leq d_{\min} \leq d_+(x), d_-(x)\leq d_{\max}$. Then given an $\epsilon > 0$, there exist an integer $N_3 > 0$ and constants $c_{7}, c_{8} > 0$ such that

Remark 1. If $d_{+}(x), \, d_{-}(x)\in C[x_L, x_R]$, then $\max\limits_{1\leq k\leq N}\max\{\Delta_{N_3}d_{+, k}, \Delta_{N_3}d_{-, k}\}$ can be very small for sufficiently large $N$, which implies that $P_{1}^{-1}$ then will be a good approximation to $\tilde{A}^{-1}$.

4.2.   Approximation of $P_{2}^{-1}$ to $P_{1}^{-1}$

Now we consider the difference between $P_{2}^{-1}$ to $P_{1}^{-1}$. As shown in Section 3.3 of ^[11], we know that, for a given $\epsilon > 0$, there exists a polynomial $p_{k}(t) = \sum_{\ell = 0}^{k} a_\ell t^\ell$ of degree $k$ such that

Define

and

Then we have

Lemma 4.4. Let $\tilde{P}_1$ and $\tilde{P}_2$ be defined by (4.22) and (4.23), respectively. Then we have

Proof. Observe that

It is easy to see that

For the matrix $G_{\alpha}$ and its Strang's circulant approximation $C_{\alpha}$, they have the following form

where $K_0, K_1\in\mathbb{R}^{m \times m}$. Therefore, we have

Similarly, we can show that

Therefore, we can see from (4.25), (4.26), (4.27) and (4.28) that $K_{i}-C_{i}$ is of the form

where $"+"$ denotes a $m$-by-$m$ block matrix.

Analogous to the proof of Lemma 3.11 in ^[11], it follows from

and

that ${\mathrm{rank}}(\tilde{P}_2-\tilde{P}_1)\leq 8km$. □

On the other hand, for $P_1-\tilde{P}_1$ and $P_2-\tilde{P}_2$, we can directly apply the proof of Lemma 3.12 in ^[11] to obtain the following lemma.

Lemma 4.5. Let $\tilde{P}_1$ and $\tilde{P}_2$ be defined by (4.22) and (4.23), respectively. Then given a $\epsilon > 0$, there exist constants $\tilde{c}_7 > 0$ and $\tilde{c}_8 > 0$ such that

By Lemma 4.4 and Lemma 4.5, we obtain the following theorem.

Theorem 4.2. Let $P_{1}^{-1}$ and $P_{2}^{-1}$ be defined by (3.9) and (3.10), respectively. Then given a $\epsilon > 0$, there exists constants $\tilde{c}_9 > 0$ and $k > 0$ such that

where $E = (P_2^{-1}-\tilde{P}_2)+(\tilde{P}_1-P_1^{-1})$ is a small norm matrix with $\|E\|_{2} < \tilde{c}_9\epsilon$, and $M = \tilde{P}_2-\tilde{P}_1$ is a low rank matrix with ${\mathrm{rank}}(M)\leq 8km$.

4.3.   Approximation of $P_{3}^{-1}$ to $P_{2}^{-1}$

We can directly apply the analysis of section 3.4 in ^[11] to reach the following conclusion.

Lemma 4.6. Let $P_{2}^{-1}$ and $P_{3}^{-1}$ be defined by (3.10) and (3.11), respectively. Then, for a given $\epsilon > 0$, there exist an integer $\ell_{0} > 0$ and a constant $c_{9} > 0$ such that

holds when the number of interpolation points in $P_{3}^{-1}$ satisfies $\ell\geq \ell_{0}$.

Summarizing our analysis, we have the following theorem.

Theorem 4.3. Let $P_{3}^{-1}$ and $A$ be defined by (3.11) and (3.2), respectively. Then we have

where $E$ is a low rank matrix and $S$ is a small norm matrix.

Proof. Observe that

By (3.2), (3.7), Theorems 4.1, 4.2 and Lemma 4.6, we can show that

where $E_1$ and $E_2$ are two low rank matrices, and $S_1$ and $S_2$ are two small norm matrices. Therefore, it is easy to see from (4.29) that

where $E$ is a low rank matrix and $S$ is a small norm matrix. □

5.   Numerical experiments

In this section, we carry out numerical experiments to show the performance of the proposed preconditioners $P_3$ and $P_4$. We use the preconditioned conjugate gradient (PCG) method to solve the symmetric positive definite linear systems (2.6). In all experiments, we set the stopping criterion as

where $r_{k}$ is the residual vector after $k$ iterations and $b$ is the right-hand side.

For the sake of comparison, we also test two other types of preconditioners which are proposed in ^[7]. One is the Strang's circulant preconditioner $C_s$ defined as

where $s(M)$, $s(G_{\alpha})$ and $s(G_{\beta})$ are the strang-circulant approximations of the Toeplitz matrices $M$, $G_{\alpha}$ and $G_{\beta}$, and $\tilde{d}_{+}$ and $\tilde{d}_{-}$ are the mean values of the diagonals of $\tilde{D}_{+}$ and $\tilde{D}_{-}$, respectively. Another is the banded preconditioner $B(\ell)$ defined as

where $\tilde{B}_{\alpha}$ and $\tilde{B}_{\beta}$ are the banded approximations of $\tilde{G}_{\alpha}$ and $\tilde{G}_{\beta}$ which are of the form

respectively.

Example 1. We consider balanced fractional diffusion equation (1.1) with $(x_L, x_R)$ $= (0, 1)$ and $T = 1$. The diffusion coefficients are given by

and $f(x, t) = 1000$.

In the numerical testing, we set the number of temporal girds with $M_t = N/2$. The results for different values of $\alpha$ and $\beta$ are reported in Table 1. We do not report the numerical results of the unpreconditioned CG method as it converges very slow.

In the table, we use $"P_3(\ell)"$ and $"P_4(\ell)"$ to denote our proposed preconditioners with $\ell$ being the number of interpolation points, "Iter" to denote the average number of iteration steps required to solve the discrete fractional diffusion equation (2.6) at each time step, and "CPU" to denote the total CPU time in seconds for solve the whole discretized system.

From Table 1, we can see that the banded preconditioner $B(\ell)$ and our proposed preconditioners $P_3(\ell)$ and $P_4(\ell)$ perform much better than the circulant preconditioner $C_s$ in terms of both iteration number and elapsed time. The banded preconditioner $B(\ell)$ has better performance when $\alpha = \beta = 0.9$. In this case, the off-diagonals of the coefficient matrix decay to zero very quickly. For other cases, $P_3(\ell)$ and $P_4(\ell)$ perform much better.

We also see that the performance of $P_3(\ell)$ and $P_4(\ell)$ are very robust in the sense that the average iteration numbers are almost the same for different values of $\alpha$ and $\beta$, as well as for the different mesh size.

Example 2. Consider the two dimensional balanced fractional diffusion equation

where $\frac{1}{2} < \alpha_1, \beta_1, \alpha_2, \beta_2 < 1$, and $d_{\pm}(x, y, t) > 0$, $e_{\pm}(x, y, t) > 0$ for $(x, y, t)\in\Omega \times (0, T]$.

As was shown in ^[7], the finite volume discretization leads to

where

Here $M_{x}, M_{y}$ are tridiagonal matrices and $D_{\pm}^{(j)}$, $E_{\pm}^{(j)}$ are block diagonal matrices

where $N_{x}$ and $N_{y}$ denote the number of grid points in the $x$-direction and $y$-direction, respectively.

Similar to the construction of $P_{3}$ of (3.11), we can define an approximate inverse preconditioner for the two dimensional problems. For convenience of expression, we omit the superscript in $A^{(j)}$. To construct the preconditioner, we define the following circulant matrices

Then we choose the interpolation points $\{(\tilde{x}_i, \tilde{y}_j\}$, $i = 1, 2, \ldots, \ell_x$, $j = 1, 2, \ldots, \ell_y$, and we approximate $\tilde{C}_{i, j}^{-1/2}$ by

where $F$ refers to the two dimensional discrete Fourier transform matrix of size $N_{x}N_{y}$, and $\tilde{\Lambda}_{k, l}$ is the diagonal matrix whose diagonals are the eigenvalues of $\tilde{C}_{k, l}$ for $1\leq k\leq \ell_x$ and $1\leq l\leq \ell_y$. Then we obtain the resulting preconditioner

where $\Phi_{k, l} = {\rm{diag}}\left(\phi_{k, l}(x_1, y_1), \phi_{k, l}(x_1, y_2), \cdots, \phi_{k, l}(x_{N_x}, y_{N_y})\right)$ for $1\leq k\leq \ell_x$ and $1\leq l\leq \ell_y$.

In the tests, we set $\Omega = (0, 1)\times(0, 1)$, $T = 1$,

and $f(x, y, t) = 1000$. We also set $N_{x} = N_{y} = N$, $M = N/2$ and $\ell_{x} = \ell_{y} = \ell$. The maximum iteration number is set to be 500. As it is expensive to compute the second item of the Sherman-Morrison-Woodburg formula 3.5 for the two dimensional problem, we only test the preconditioner $P_3(\ell)$.

The results are reported in Table 2. We can see that, for the two dimensional problem, although our proposed preconditioners $P_3(3)$ and $P_3(4)$ take more iteration steps to converge than the banded preconditioners $B(3)$ and $B(4)$, it is much cheaper to implement our new preconditioners than banded preconditioner and circulant preconditioner in terms of CPU time.

6.   Concluding remarks

In this paper, we consider the preconditioners for the linear systems resulted from the discretization of the balanced fractional diffusion equations. The coefficient matrices are symmetric positive definite and can be written as the sum of a tridiagonal matrix and two Toeplitz-times-diagonal-times-Toeplitz matrices. We investigate the approximate inverse preconditioners and show that the spectra of the preconditioned matrices are clustered around 1. Numerical experiments have shown that the conjugate gradient method, when applied to solve the preconditioned systems, converges very quickly. Besides, we extend our preconditioning technique to the case of two dimensional problems.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is supported by NSFC grant 11771148 and Science and Technology Commission of Shanghai Municipality grant 22DZ2229014.

Conflict of interest

All authors declare no conflicts of interest in this paper.