Fast Crank-Nicolson compact difference scheme for the two-dimensional time-fractional mobile/immobile transport equation

1.   Introduction

In this paper, we focus on the efficient numerical scheme for solving the two-dimensional time-fractional mobile/immobile transport equation:

where $\Omega \subset \mathbb{R}^{2}$ is a bounded rectangle domain and $T$ is the positive final time. The corresponding initial and boundary conditions satisfy $u(\mathbf{x}, 0) = v(\mathbf{x}), \mathbf{x}\in \overline{\Omega}$ and $u(\mathbf{x}, t) = 0, (\mathbf{x}, t)\in\partial\Omega\times[0, T]$, respectively. Here, $\mu, \kappa_1$, and $\kappa_2$ are some fixed positive constants, $f$ and $v$ are two given smooth functions, $0 < \alpha < 1$, and $_{RL}D_{0, t}^{\alpha}$ is the Riemann-Liouville derivative defined by:

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

The fractional mobile/immobile transport Eq (1.1) can be viewed as the limiting equation that govern continuous time random walks with heavy tailed random waiting times ^[1]. The solution regularity of (1.1) can be derived by using the Laplace transform tools and operator approach, please refer to our earlier paper ^[2] for similar discussion. We do not intend to go into further discussion on this point, because this would detract from our focus in this paper.

Until now, the numerical solution of the Eq (1.1) has been partially investigated by some researchers ^[3,4]. By assuming the solution is sufficiently smooth, Liu et al. developed a fast compact finite difference scheme based on an equivalent form involving the Caputo derivative for solving the one-dimensional fractional mobile/immobile transport Eq ^[3]. Their fast solution technique is based on the fast Fourier transform to improve the computational performance in the time-marching system. Qiu et al. utilized the weighted and shifted Grünwald formula to obtain alternating direction implicit Galerkin finite element scheme for solving the distributed-order time-fractional mobile/immobile equation involving zero initial and Dirichlet boundary conditions ^[4]. We can see that the methods mentioned above are effective for solving smooth solution problems, but not for non-smooth solution ones, even though the non-smoothness of fractional models has gradually received considerable attention from researchers.

A variety of strategies can be found in ^[5,6,7,8] for dealing with the non-smooth solution case, just to name a few. One of the most commonly used is to employ the method of adding correction terms which is initiated in ^[9] and revisited in ^[10]. For the fractional model (1.1), there is a few literatures on such subject. In ^[11], Yin et al. derived the stable finite element scheme and considered the method of adding correction terms to overcome the initial singularity of the time fractional derivative. Inspired by this, and noting that the Crank-Nicolson scheme can be derived naturally by the modified L1 method ^[12], we shall therefore focus on the modified L1 method with correction terms to derive the efficient Crank-Nicolson scheme for solving the non-smooth problem in (1.1).

Fast algorithms for the fractional model (1.1) are also another focus of this paper. It is known that the nonlocal property of the fractional derivative always seriously impact the computational performance of the existing numerical schemes, especially for the high-dimensional problems. By employing the local radial basis functions, Nikan et al. proposed an efficient meshless approach for solving the two-dimensional time-fractional Klein-Kramers model ^[13]. In ^[14], the authors applied the fourth-order compact finite difference method to discrete the Poisson equation with Dirichlet boundary value condition, and solved the resulting discretized system efficiently with fast discrete sine transform (DST). This allows their algorithm to avoid computing matrix inversion directly. With the fast Fourier transform (FFT), the computational cost is reduced from $O(M_1^2 M_2^2)$ to $O(M_1M_2\log(M_1M_2))$ for two-dimensional problem. Here, $M_1$ and $M_2$ denote the total number grid points in $x$- and $y$-axis direction, respectively. So, for the model problem (1.1), we shall further employ the compact finite difference operator to discrete the Laplacian and then combining the DST technique to obtain the efficient fast Crank-Nicolson compact difference scheme. So far as we know, the use of DST technique with proper correction terms to efficiently solve the fractional model (1.1) has not been found in the existing literatures yet.

The contributions of this paper are listed as follows. First, we develop a fast Crank-Nicolson compact difference scheme by employing the modified L1 method and DST technique; see the scheme (4.1). Second, The non-smooth issue in the fractional model (1.1) is addressed by adding appropriate correction terms. And more importantly, the added correction terms have no impact on the stability of the original scheme; see the scheme (4.2). Third, we show theoretically and numerically that our scheme (4.2) solves the fractional model (1.1) rapidly and efficiently with the accuracy of $O(\tau^{2-\alpha}+h^4)$; see Theorems 3.1 and 3.2, and Tables 1-5. The organization of this paper is as follows. In Section 2, we derive the Crank-Nicolson compact difference scheme based on the modified L1 method. The stability and error estimates of the proposed scheme are proved rigorously in Section 3. In Section 4, we further apply the DST technology and the method of adding correction terms to deal with the non-smooth solution case. Numerical examples are demonstrated to conform the effectiveness of the scheme in Section 5. Finally, we give the conclusions of the paper in Section 6.

In what follows, the symbol $c$ (with or without subscript) is used to denote the constant, which may vary in different scenario but is independent of the temporal and spatial stepsizes.

2.   Derivation of the scheme

In this section, we provide the derivation of the numerical scheme for numerically solving (1.1).

Let $n_T$ be a positive integer. The temporal stepsize $\tau$ is given by $\tau = T/n_T$. We denote the grid point $t_n = n\tau$ for $n\geq 0$. If $g(t)\in C^2[0, T]$, we have the following modified L1 method for the approximation of Riemann-Liouville derivative at $t = t_{n+\frac{1}{2}}(: = (t_n+t_{n+1})/2)$:

where the error $|R^{n+\frac{1}{2}}|\leq c\tau^{2-\alpha}\max_{t\in[0, T]}|g''(t)|$ (cf. [12, Lemma 3.1]) and

Here, the weights $b_k, B_k,$ and $A_k$ are defined by

and

Letting $t = t_{n+1/2}$ in (1.1) yields

Notice that $\partial_{t} u(\mathbf{x}, t_{n+1/2}) = \delta_t u^{n+1/2} + O(\tau^2)$ where $\delta_t u^{n+1/2} = (u^{n+1}-u^n)/\tau$, we readily have the following form:

where the local truncation $R_\mathbf{x}^{n+1/2} = O(\tau^{2-\alpha})$ and $u^{n+1/2} = (u^{n+1}+u^n)/2$.

Next, we consider the fourth-order compact difference operator for spatial discretization.

To enable our numerical scheme and theoretical analysis to be easily extended to other high-dimensional problems, here we use partially the notations of the paper ^[15]. Denote the spatial stepsize $h_k = (x_{k}^R - x_k^L)/M_k$ where $M_k$ is a positive integer and the grid point $x_{k, j_k} = x_k^L + j_k h_k$ for $j_k = 0, 1, \cdots, M_k.$ The subscript $k(1\leq k\leq d)$ here represents the spatial direction at $k$th position. In this paper, we focus on the two-dimensional case: $d = 2$.

The discrete grids in space is given by $\overline{\Omega}_h = \{x_h = (x_{1, j_1}, x_{2, j_2}, \cdots, x_{d, j_d})|0\leq j_k\leq M_k, 1\leq k\leq d\}$. We further denote that $\Omega_h = \overline{\Omega}_h\cap \Omega$ and the boundary $\partial\Omega_h = \overline{\Omega}_h\cap \partial\Omega$. The space of grid function is denoted as $V_h = \{v| v = (v_h)_{\mathbf{x}_{h}} \; \rm{and}\; v_h = 0 \; \rm{for}\; \mathbf{x}_{h}\in\partial\Omega_h\}$. For simplicity, we introduce the following difference operator for the grid function $v_h = v(\mathbf{x}_h)$ with the index vector $h = (i_1, i_2, \cdots, i_d)$ at $k$th position:

with $\delta_k^2 v_{i_k} = (\delta_k v_{i_k+\frac{1}{2}} - \delta_k v_{i_k-\frac{1}{2}})/h_k$ and $\delta_{k} v_{i_k+\frac{1}{2}}: = (v_{i_{k}+1}-v_{i_{k}})/h_{k}.$ The compact difference operator is given by $\overline{\Delta}_{k} v_{i_{k}}: = \frac{\delta_{k}^{2}}{\mathcal{H}_{k}} v_{i_{k}}$. So we have the fourth-order spatial approximation of $\Delta v\left(\mathbf{x}_{h}\right)$ for $\mathbf{x}_{h}\in\Omega_h$ as follows: $\overline{\Delta}_h v_h: = \sum_k \overline{\Delta}_k v_h$.

Based on the above form (2.3), we apply the compact difference approximation in space to obtain

The local truncation error is given by $R_{\mathbf{x}t}^{n+\frac{1}{2}} = O(\tau^{2-\alpha} + h^4)$. Omitting the small term $R_{\mathbf{x}t}^{n+\frac{1}{2}}$, we obtain the following fully discrete Crank-Nicolson compact difference scheme for the model (1.1): Find the numerical solution $u_h^n$ of $u(\mathbf{x}_h, t_n)$ for $n\geq 1$, such that

The initial and boundary conditions are given by $u_h^{0} = v(\mathbf{x}_h)$ and $u(\mathbf{x}_{h})|_{\mathbf{x}_{h}\in\partial\Omega_h} = 0$, respectively.

Remark 2.1. When $\alpha\rightarrow 1$, the Crank-Nicolson compact difference scheme $(2.5)$ recovers the classical one

for solving the corresponding diffusion equation:

One the other hand, when $\alpha\rightarrow 0$, we have

from which we can numerically solve the classical diffusion equation with a reaction term:

Thus in this sense, we can say that the derived numerical scheme $(2.5)$ is compatible with integer-order one. Notice that the compatibility of the numerical scheme is important in solving nonlocal equations, one can refer to ^[16] for more details.

Remark 2.2. Formally, one can use the limiting properties of Riemann-Liouville derivative in the fractional model $(1.1)$ to naturally obtain the corresponding integer-order equation, however the underlying physical interpretation needs to be studied further, see the similar discussion in Section 4.2 of the paper ^[1].

3.   Stability and error estimates

In this section, we study the stability and error estimates for the Crank-Nicolson compact difference scheme (2.5).

For any grid function $v\in V_h$, the discrete $L^2$-norm is given by $\|v\| = \sqrt{(v, v)_h}$ with the discrete inner product $(u, v)_h = (\prod_{k = 1}^d h_k)\sum_{\mathbf{x}_h\in\Omega_h}u_h v_h$. The discrete $H^1$-seminorm and $H^1$-norm are denoted as $|v|_1 = \sqrt{\sum_{k = 1}^d\|\delta_{k} v_h\|^{2}}$ and $\|v\|_1 = \sqrt{\|v\|^{2}+|v|_1^{2}}$. In view of the embedding theorem, one can readily have the equivalence of $|v|_1$ and $\|v\|_1$ for any $v\in V_h$.

We shall need the boundedness of the discrete operator $\overline{D}_\tau^{\alpha}$, which is given by the following lemma.

Lemma 3.1. For any real-valued functions $g^n, n\geq 0$ defined on $\Omega$, the discrete operator $\overline{D}_\tau^{\alpha}$ defined by $(2.2)$ satisfies the following inequality:

Proof. One can refer to the Lemma 4.2 in ^[12] or Lemma 4.4 in ^[17] for similar discussion. Thus the proof is completed.

Now we are ready to present the stability of the scheme (2.5).

Theorem 3.1. The fully discrete Crank-Nicolson compact difference scheme $(2.5)$ is stable in the sense that

Proof. Taking the discrete inner product of (2.5) with $2\tau u_h^{n+1/2}$, one has

Since the matrix corresponding to $\mathcal{H}_{k}$ in $\overline{\Delta}_h$ is positive definite with the eigenvalues of the form: $\lambda_{k, j_k} = \frac{5}{6}+\frac{1}{6}\cos(\frac{j_k\pi}{M_k})$, one can obtain the boundedness of $\overline{\Delta}_h$ in discrete inner product, that is $\frac{3}{2}(\Delta_h u_h^{n+1/2}, u_h^{n+1/2})_h < (\overline{\Delta}_h u_h^{n+1/2}, u_h^{n+1/2})_h < (\Delta_h u_h^{n+1/2}, u_h^{n+1/2})_h$ with the notation $\Delta_h u_h^{n} = \sum_{k = 1}^d \delta_k^2 u_h^n$ (cf. [15, Theorem 3.1]). Notice that $(\Delta_h u_h^{n+1/2}, u_h^{n+1/2})_h = -|u_h^{n+1/2}|_1^2$, So discarding this nonpositive term and applying Lemma 3.1, we obtain

Denote $G^n = \kappa_1\|u_h^{n}\|^2 + \tau^{1-\alpha}\kappa_2\sum_{k = 1}^n b_{n-k}\|u_h^{k-1/2}\|^2$. Then the above inequality can be written as a more compact form:

Summing up $n$ from 1 to $m$ and replacing $m$ with $n$, we arrive at

By the Cauchy-Schwarz inequality, the third term on the right-hand side of the above inequality has the following estimates:

from which we have

It remains to estimate the $G^1$. To this end, we consider the case $n = 0$ for the scheme (2.5). Similarly, by taking the discrete inner product of (2.5) with $2\tau u_h^{1/2}$, we have

After expanding the above equation, we get

which leads to the inequality:

Applying the Cauchy-Schwarz inequality, we derive that

where the two positive constants $\epsilon_1$ and $\epsilon_2$ are chosen such that $\kappa_2 A_0\frac{\tau^{1-\alpha}}{\epsilon_1}\|u_h^{1/2}\|^2 + \frac{\tau^{1-\alpha}}{\epsilon_2}\|u_h^{1/2}\|^2\leq \kappa_2\tau^{1-\alpha}B_0\|u_h^{1/2}\|^2$. So, we can obtain

Thus, by (3.1), we have

This together the two inequalities $\sum_{k = 1}^n A_k\leq \frac{(1-\alpha)2^\alpha}{\Gamma(2-\alpha)}$ and $\tau^{\alpha}\sum_{k = 2}^{n+1}\frac{1}{b_{n+1-k}}\leq \frac{T^\alpha}{(1-\alpha)\Gamma(2-\alpha)}$ (cf. [17, Theorem 4.5]) yields the desired results.

The convergence result is presented by the following theorem.

Theorem 3.2. Suppose that $u\in C^2(0, T; C^6(\Omega)), f\in C(0, T; C^4(\Omega))$ and $v\in C(\Omega)$, then the numerical solution $u_h^n$ is convergent with respect to the discrete $L^2$-norm. Especially, the following error estimate holds:

where $n\geq 1$.

Proof. Denote the error $e_h^n = u(\mathbf{x}_h, t_n)-u_h^n$ for $\mathbf{x}_h\in\Omega_h$, then subtracting (2.5) from (2.4), we get the following error equation:

Apply the Theorem 3.1, we obtain

where the last inequality holds since $n\tau\leq cT$. Thus the proof is completed.

4.   Corrected scheme with fast solver

In order to avoid the direct calculations for the original scheme (2.5), we employ the DST in spatial direction to improve the computational performance. For the grid function $v_h$, the discrete sine transform of $v_h$ at the $k$th position is given by $v_{i_k} = \sum_{j_k = 1}^{M_k-1}\widehat{v_{j_k}}\sin(i_k j_k\pi/M_k)$. So, in view of the compact difference operator, we derive that

where $s_{j_{k}} = \cos(\frac{j_{k} \pi}{M_{k}})$ and $1\leq j_k\leq M_k-1$. One can refer to ^[14,15] for more details about the derivation. Therefore, the scheme (2.5) with fast solver has the following form:

from which we can obtain the numerical solution $u_h^n$ by the inverse DST of $\widehat{u}_\nu^{n}$ for $n\geq 1$. Here, the index set $\nu = \{(j_1, j_2, \cdots, j_d)|1\leq j_k\leq M_k-1, 1\leq k\leq d\}$. $\widehat{u}_\nu^{n}$ and $\widehat{f}^{n+\frac{1}{2}}$ are obtained from $u_h^{n}$ and $f^{n+\frac{1}{2}}$ by means of DST. From the scheme (4.1), we avoid directly computing the matrix inversion in the linear discrete system (2.5) at each temporal layer, thus greatly improving the computational efficiency, i.e., the computational cost will be $O(M_1M_2\log (M_1M_2))$ instead of $O(M_1^2 M_2^2)$.

Next, we further use the method of adding suitable correction terms to deal with non-smoothness of the solution. Consider the numerical approximation at $t = t_{n+\frac{1}{2}}$, we can modify the method (2.1) as

in which the starting weights $w_{n, k}^{(\alpha)}$ are chosen such that the above scheme is exact for some $g(t) = t^{\sigma_r}, (0\leq r\leq m)$, that is, they can be determined by the following linear system:

The corrected scheme for the first-order time derivative can be written as

Similarly, the starting weights $w_{n, k}^{(1)}$ are obtained by solving the linear system:

Thus, the fast Nicolson compact difference scheme (4.1) with correction terms has the following form:

It is worth noting that the numerical scheme (4.2) with suitable correction terms does not affect the stability of the original scheme (4.1), see ^[18] for similar discussion. Although calculating the starting weights raises some computational complexity, it improves the accuracy in the temporal direction for solution of equation (1.1) with low regularity.

5.   Numerical examples

In this part, we present two numerical examples to verify the theoretical result in Section 3 and the effectiveness of the scheme (4.2) in Section 4. We measure the $L^2$-norm error at $t = t_n$ by $e(n, h) = \|u(\mathbf{x}_h, t_n) - u_h^n\|$. The corresponding temporal and spatial convergent orders are calculated by $\log(e(n, h)/e(2n, h))$ and $\log(e(n, h)/e(n, h/2))$, respectively. The parameters $\mu$, $\kappa_1$ and $\kappa_2$ in equation (1.1) are all set to one. The computational domain are restricted to $\Omega = (0, 1)^2$. All the numerical results are obtained at $t = T$ with the final time $T = 1$. The numerical examples are tested in MATLAB software (R2020a) on an Apple OS platform with a quad-core 2.3 GHz processor and 8 GB of memory.

Example 5.1 (Smooth solution). Consider the following problem with zero Dirichlet boundary conditions:

where

The exact solution is $u = \sin(\pi x)\sin(\pi y)(1+t^\gamma)$. In order to test the correctness of the theoretical result in Section 3, we let $\gamma = 2.1$. By the fast Crank-Nicolson compact difference scheme $(4.1)$, the errors and convergence orders are obtained and demonstrated in Tables 1 and 2. From the data in these two tables, it can be observed that the temporal and spatial convergence orders of the numerical scheme in this example are consistent with the theoretical analysis.

Next, we compare the implemented CPU time applied by the fast Crank-Nicolson compact difference scheme $(4.1)$ with the original compact difference scheme $(2.5)$ in Table 3. As can be seen from the table, both schemes $(2.5)$ and $(4.1)$ achieve almost the same accuracy for the fixed temporal and spatial stepsizes and fractional order $\alpha$, however, it is clear that numerical scheme $(4.1)$ with DST requires less CPU execution time than $(2.5)$. And when the temporal stepsize $\tau$ is fixed, the computational advantage of numerical scheme $(4.1)$ is more obvious as the spatial stepsize $h$ keeps decreasing, which shows that numerical scheme with DST can effectively handle high-dimensional problems.

Example 5.2 (Non-smooth solution). Consider the following problem data:

with the zero initial value condition and zero boundary condition. Then the exact solution is given by $u = \sin(\pi x)\sin(\pi y)t^\gamma$.

In this example, we test the effectiveness of the fast Crank-Nicolson compact difference scheme with correction terms $(4.2)$ for dealing with the non-smooth solution case. To this end, we set $\gamma = 0.1$. So the first-order partial derivative of $u$ with respect to $t$ is unbounded at $t = 0$. The numerical results are presented in Tables 4 and 5. For Table 4, it can be noticed that when there is no correction term in $(4.2)$, the temporal accuracy of the numerical scheme is damaged as the regularity of the solution does not satisfy the smoothness requirement specified in the theoretical analysis. However, with the addition of the correction terms (i.e., $m = 1$ and $m = 3$, see the last four columns of Table 4), the temporal accuracy of the numerical scheme is somewhat maintained and therefore some improvement in the convergence orders can be observed. Similar results are also observed in Table 5. So, this indicates that the addition of suitable correction terms can surely improve the accuracy of the numerical scheme, thus confirming that the fast Crank-Nicolson compact difference scheme $(4.2)$ with correction terms is efficient for the non-smooth issue of the fractional model $(1.1)$.

6.   Conclusions

In this paper, we develop the fast Crank-Nicolson compact difference scheme for the time-fractional mobile/immobile equation (1.1). By using the DST technology and the method of adding correction terms, we improve the computational performance of the proposed scheme for solving the non-smooth solution problem in (1.1). The corresponding stability and error estimates are given rigorously. Numerical examples fully confirm the effectiveness of the proposed scheme.

Although our discussions only focus on two-dimensional problem, the theoretical analysis can be easily extended to other high-dimensional one. When the fractional order $\alpha$ goes to one or zero, we observe that the proposed scheme is compatible with the integer-order one. In the subsequent study, we shall further improve the computational effectiveness by combining fast algorithms in time, such as the fast convolution method ^[19], the sum-of-exponentials method ^[20] and so on. In addition, constructing efficient algorithms for high-dimensional fractional models with other boundary conditions is also in our scope of consideration.

Acknowledgments

The work is supported by the Guangxi Natural Science Foundation [grant numbers 2018GXNSFBA281020, 2018GXNSFAA138121] and the Doctoral Starting up Foundation of Guilin University of Technology [grant numbers GLUTQD2014045, GLUTQD2016044].

Conflict of interest

All authors declare no conflicts of interest in this paper.