A two-grid mixed finite volume element method for nonlinear time fractional reaction-diffusion equations

1.   Introduction

In this paper, we consider the following nonlinear time fractional reaction-diffusion equations with the initial and Dirichlet boundary conditions

where $\Omega\subset \mathbb{R}^2$ is a bounded convex polygonal domain with the boundary $\partial\Omega$, $J = (0, T]$ with $0 < T < \infty$. Assume that the functions $u_0({\mathit{\boldsymbol{x}}})$, $g(u({\mathit{\boldsymbol{x}}}, t))$ and $f({\mathit{\boldsymbol{x}}}, t)$ are smooth enough, and there exists a constant $L > 0$ such that $|g(u)|\leq L|u|$. The diffusion coefficient matrix $\mathcal{A}({\mathit{\boldsymbol{x}}}) = (a_{ij}({\mathit{\boldsymbol{x}}}))_{2\times 2}$ is symmetric and uniformly positive definite, that is, there exist two constants $A_*, A^* > 0$ such that

Moreover, we should assume that $\mathcal{A}^{-1}({\mathit{\boldsymbol{x}}})$ satisfies the Lipschitz condition. In (1.1), the Caputo time fractional derivative $\frac{\partial^\alpha u({\mathit{\boldsymbol{x}}}, t)}{\partial t^\alpha}$ with order $\alpha\in (0, 1)$ is defined by

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

Fractional differential equations (FDEs) can be applied to simulate various natural phenomena in chemistry, physics and biology and so on ^[1,2,3], which have attracted great interest of more and more scholars. However, it is very difficult to obtain the exact solutions for a large number of FDEs due to the nonlocality of fractional integrals and derivatives and other reasons, such as complex nonlinear terms, initial or boundary conditions. Therefore, a lot of numerical algorithms have been proposed and applied to solve FDEs ^{[4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]}, including finite element (FE) methods, finite difference (FD) methods, finite volume/element (FV/FVE) methods, discontinuous Galerkin (DG) methods, spectral methods and so on. In this paper, we establish a two-grid algorithm to solve the nonlinear time fractional reaction-diffusion Eq (1.1).

The two-grid method is proposed and developed by Xu ^[24,25] to solve nonlinear elliptic partial differential equations based on FE methods. Because of the advantage of saving computing time, many scholars have extended and applied it to integer order partial differential equations. Dawson et al. ^[26] presented a two-grid mixed finite element (MFE) method for nonlinear parabolic equations which arises in flow through porous media, and gave the error analysis. Yan et al. ^[27] proposed a two-grid FVE method for the nonlinear Sobolev equations, and obtained optimal $H^1$-norm error estimate. Hou et al. ^[28] applied a two-grid expanded MFE method to solve semi-linear parabolic integro-differential equations, and gave the convergence analysis and some numerical results. Liu ^[29] presented a two-grid FVE method for semi-linear reaction-diffusion system of the solutes in the groundwater flow, and obtained the error estimates in $L^2$-norm and $H^1$-norm. In recent years, the two-grid method was also applied to solve fractional partial differential equations. Liu et al. ^[30] proposed a two-grid MFE algorithm for a nonlinear fourth-order reaction-diffusion model with the Caputo time fractional derivative, and obtained the unconditional stability and error estimates. Liu et al. ^[31] presented a two-grid FE algorithm for a time fractional Cable equation, in which the Riemann-Liouville fractional derivative was approximated by the second-order weighted and shifted Grünwald difference (WSGD) scheme. Li et al. ^[32] constructed a two-grid expanded MFE scheme to solve a semilinear time fractional reaction-diffusion equation, in which the Caputo fractional derivative was approximated by the $L1$-formula. Li et al. ^[33] proposed a two-grid FE method for a nonlinear time fractional diffusion equation, and gave some numerical results to confirm the theoretical results. Chen et al. ^[34] studied a two-grid modified method of characteristics scheme to solve nonlinear variable-order time fractional advection-diffusion equations, and obtained the optimal $L^2$-norm error estimates. Liu et al. ^[35] presented a two-grid FE fast algorithm to solve a nonlinear space-time fractional diffusion equation, and gave the stability and convergence analysis. From the current literatures, we find that there is no report about the two-grid fast algorithm based on the mixed finite volume element (MFVE) method ^{[36,37,38,39]} for solving the FDEs.

In this paper, we will construct a two-grid MFVE algorithm to solve the nonlinear time fractional reaction-diffusion equations. In temporal discretization, we select the classical $L1$-formula to approximate the Caputo time fractional derivative. In spatial discretization, we construct coarse and fine grids (containing primal and dual grids), and establish a two-grid MFVE scheme by introducing an auxiliary variable $\boldsymbol \lambda$ and using the transfer operator. The calculation process is divided into two steps: firstly, the coarse solution is computed iteratively by using the nonlinear MFVE scheme on the space coarse grid, then a linearized scheme is constructed by using the coarse solution, and finally solution on the space fine grid is obtained. In our theoretical analysis, we give the existence and uniqueness results of the fully discrete solutions for the two-grid MFVE scheme by applying the Browder fixed point theorem and the matrix theory, and obtain unconditional stability results and error estimates in $L^2(\Omega)$-norm for the variable $u$. Moreover, we derive the conditional stability results and error estimates in $(L^2(\Omega))^2$-norm and ${\mathit{\boldsymbol{H}}}({{{\rm{div}}}})$-norm for the variable $\boldsymbol\lambda$ by using a special analytical technique. Finally, we give some numerical results to verify the feasibility and effectiveness, and find that the proposed two-grid MFVE algorithm can greatly save the computing time.

The layout of this paper is as follows: By constructing coarse and fine grids (primal and dual) and introducing the transfer operator, a two-grid MFVE algorithm for the nonlinear time fractional reaction-diffusion equation is proposed in Section 2. Some properties of the transfer operator $\gamma_{\hbar}$ and the fractional Gronwall inequality are given, and the existence and uniqueness results are obtained in Section 3. In Sections 4 and 5, the stability and error estimates are derived in detailed. In Section 6, two numerical examples are given to verify the feasibility and effectiveness.

2.   Two-grid MFVE scheme

We shall use the standard Sobolev spaces $W^{m, p}(\Omega)$ with the norm $\|\cdot\|_{m, p}$. For $p = 2$, we define $H^m(\Omega) = W^{m, 2}(\Omega)$ with the norm $\|\cdot\|_m$, and $H^0(\Omega) = L^2(\Omega)$ with the inner product $(\cdot, \cdot)$ and the norm $\|\cdot\|$. We also use ${\mathit{\boldsymbol{H}}}({{{\rm{div}}}}, \Omega) = \{{\mathit{\boldsymbol{v}}} \in (L^2(\Omega))^2, {{{\rm{div}}}} {\mathit{\boldsymbol{v}}} \in L^2(\Omega)\}$ with the norm $\|\cdot\|_{{\mathit{\boldsymbol{H}}}({{{\rm{div}}}})}$. Furthermore, throughout this paper, the mark $C$ is a generic positive constant, which is independent of the mesh parameters.

In order to get the MFVE scheme, by introducing an auxiliary variable $\boldsymbol\lambda({\mathit{\boldsymbol{x}}}, t) = -\mathcal{A}({\mathit{\boldsymbol{x}}})\nabla u({\mathit{\boldsymbol{x}}}, t)$, we can rewrite the primal problem (1.1) as

Then, we can obtain the weak formulation of (2.1): Find $(\boldsymbol\lambda, u)\in {\mathit{\boldsymbol{V}}}\times W$ such that

where ${\mathit{\boldsymbol{V}}} = {\mathit{\boldsymbol{H}}}({{{\rm{div}}}}, \Omega)$ and $W = L^2(\Omega)$.

Now, we use $\mathscr{K}_h$ to denote a quasiuniform triangulation partition of the domain $\Omega$, that is $\mathscr{K}_h = \cup K_B$, where $K_B$ stands for the triangle with the barycenter $B$, referring to Figure 1. Let $h = \max\{h_{K_B}\}$, where $h_{K_B}$ is the diameter of the triangle $K_B$. Moreover, we should define the nodes of a triangular element to be its midpoints of three sides, and mark $P_1, P_2, ..., P_{M_S}$ as the inner nodes and $P_{M_S+1}, P_{M_S+2}, ..., P_M$ as the boundary nodes.

We select the lowest order Raviart-Thomas space ${\mathit{\boldsymbol{V}}}_h$ and piecewise constant function space $W_h$ as the trial function spaces for $\boldsymbol \lambda$ and $u$, respectively, where

Based on the primal partition $\mathscr{K}_h$, we construct the dual partition $\mathscr{K}^*_h$. Referring to Figure 1, the interior node $P_3$ belongs to the common side of two adjacent triangles $K_{B_1} = \Delta A_1A_2A_3$ and $K_{B_2} = \Delta A_1A_3A_5$, then we define the quadrilateral $A_1B_1A_3B_2$ to be the dual element for $P_3$. In general, for an interior node $P$, the dual element $K_P^*$ is the union of two triangles $K_L$ (with $\Delta A_1B_1A_3$) and $K_R$ (with $\Delta A_1A_3B_2$). For a boundary node such as $P_6$, the associated dual element is a triangle $K_I$ (with $\Delta A_5B_3 A_4$).

Integrating (2.1) on all the primal and dual partitions, respectively, we obtain

We define the transfer operator $\gamma_h: {\mathit{\boldsymbol{V}}}_h\rightarrow (L^2(\Omega))^2$ as follows

where $\chi^*_K$ is characteristic function of a set $K$. We use $\bar{{\mathit{\boldsymbol{Y}}}}_h = \gamma_h {{\mathit{\boldsymbol{V}}}}_h$ as the test function space, and rewrite (2.3) as

Similar to ^[37], making use of the operator $\gamma_h$ and the Green theorem, we have $(\nabla w_h, \gamma_h {\mathit{\boldsymbol{v}}}_h) = -({{{\rm{div}}}} {\mathit{\boldsymbol{v}}}_h, w_h), \forall {\mathit{\boldsymbol{v}}}_h\in {\mathit{\boldsymbol{V}}}_h, \forall w_h\in W_h.$ Then, we obtain the nonlinear semi-discrete MFVE scheme: For the selected appropriate $(\boldsymbol \lambda_h(0), u_h(0))$, find $(\boldsymbol \lambda_h(t), u_h(t))\in {\mathit{\boldsymbol{V}}}_h\times W_h$ such that

In order to approximate the Caputo time fractional derivative and give the fully discrete scheme, we should give the grid points $t_n = n\tau$ $(n = 0, 1, \cdots, N)$ in time interval $[0, T]$, where $N$ is a positive integer and $\tau = T/N$. We denote $\varphi^n = \varphi(\cdot, t_n)$ for a function $\varphi$. Following ^[4,5], we will approximate the fractional derivative $\frac{\partial^\alpha u({\mathit{\boldsymbol{x}}}, t)}{\partial t^\alpha}$ at $t = t_n$ by using the $L1$-formula as follows

where $b_k = (k+1)^{1-\alpha}-k^{1-\alpha}$, $b_0^n = (n-1)^{1-\alpha} -n^{1-\alpha}$, $b_n^n = 1$, $b_k^n = b_{n-k}-b_{n-k-1}$ $(0 < k < n)$. Setting $D_{\tau}^\alpha u^n = \frac{\tau ^{-\alpha}}{\Gamma(2-\alpha)}\sum\limits_{k = 0}^{n}b_k^n u^k$, we have $\frac{\partial^\alpha u({\mathit{\boldsymbol{x}}}, t_n)}{\partial t^\alpha} = D_{\tau}^\alpha u^n+R_t^n({\mathit{\boldsymbol{x}}})$. Following ^[4,5], we can get that if $u\in \mathcal{C}^2(\bar{J}, L^2(\Omega))$, then there exist a constant $C > 0$ independent of $\tau$ such that $\|R_t^n(\mathit{\boldsymbol{x}})\|\leq C\tau^{2-\alpha}$.

Let $\boldsymbol \lambda_h^n$ and $u_h^n$ be the numerical solutions of $\boldsymbol\lambda$ and $u$ at $t = t_n$, respectively. Then, we can obtain the nonlinear fully discrete MFVE scheme for the problem (1.1): For the properly selected $(\boldsymbol \lambda_h^0, u_h^0)$, find $(\boldsymbol \lambda_h^n, u_h^n)\in {{\mathit{\boldsymbol{V}}}}_h\times W_h$, $n = 1, 2, \cdots, N$, such that

For improving the nonlinear fully discrete MFVE scheme (2.7), we consider the following two-grid MFVE system based on the coarse grid $\mathscr{K}_H$ and the fine grid $\mathscr{K}_h$ with the corresponding dual grids $\mathscr{K}_H^*$ and $\mathscr{K}_h^*$, where $h\ll H$.

STEP I. On the coarse primal and dual grids $(\mathscr{K}_H$ and $\mathscr{K}_H^*)$, solve the following nonlinear system for $(\boldsymbol\lambda_H^n, u_H^n)\in {\mathit{\boldsymbol{V}}}_H \times W_H$, $n = 1, 2, \cdots, N$, such that

where $(\boldsymbol \lambda_H^0, u_H^0)\in {\mathit{\boldsymbol{V}}}_H \times W_H$ is defined in Section 5.

STEP II. On the fine primal and dual grids $(\mathscr{K}_h$ and $\mathscr{K}_h^*)$, solve the following linearized system for $(\hat{\boldsymbol\lambda}_h^n, \hat u_h^n)\in {\mathit{\boldsymbol{V}}}_h \times W_h$, $n = 1, 2, \cdots, N$, such that

where $(\hat{\boldsymbol\lambda}_h^0, \hat u_h^0)\in {\mathit{\boldsymbol{V}}}_h \times W_h$ is defined in Section 5.

Remark 2.1. In the actual numerical calculation of the two-grid systems (2.8) and (2.9), we can find a solution $(\boldsymbol \lambda_H^n, u_H^n)\in {\mathit{\boldsymbol{V}}}_H\times W_H$ on the coarse primal and dual grids $(\mathscr{K}_H$ and $\mathscr{K}_H^*)$ by calculating the nonlinear implicit system (2.8), then obtain the final solution $(\hat{\boldsymbol \lambda}_h^n, \hat u_h^n)\in {\mathit{\boldsymbol{V}}}_h\times W_h$ on the fine primal and dual grids $(\mathscr{K}_h$ and $\mathscr{K}_h^*)$ by calculating the linearized system (2.9). This calculation method will be more efficient than the standard nonlinear implicit system (2.7), and we will see this advantage from the numerical results.

3.   Existence and uniqueness

In the proof of existence and uniqueness and subsequent theoretical analysis, we need to use some important properties of transfer operator $\gamma_{\hbar}$ ($\hbar = H$ or $h$), which are as follows:

Lemma 3.1. ^[37] The transfer operator $\gamma_{\hbar}$ is bounded

Lemma 3.2. ^[38] The following symmetry relations holds

where $\bar{\mathcal{A}}^{-1}({\mathit{\boldsymbol{x}}}) = \mathcal{A}^{-1}(B)$, $\forall {\mathit{\boldsymbol{x}}} \in K_B$,

Lemma 3.3. ^[38] There exists three constants $\mu_1, \mu_2, \mu_3 > 0$ independent of $\hbar$ such that

Lemma 3.4. ^[38] There exists two constants $\mu_4, \mu_5 > 0$ independent of $\hbar$ such that

where $\|{\mathit{\boldsymbol{z}}}_{\hbar}\|_{1, \hbar}^2 = \|{\mathit{\boldsymbol{z}}}_{\hbar}\|^2+|{\mathit{\boldsymbol{z}}}|_{1, \hbar}^2$ and $|{\mathit{\boldsymbol{z}}}|_{1, \hbar}^2 = \sum\limits_{K\in \mathscr{K}_{\hbar}}(\|\nabla z_{\hbar}^1\|_{0, K}^2+\|\nabla z_{\hbar}^2\|_{0, K}^2)$, $\forall {\mathit{\boldsymbol{z}}}_{\hbar} = (z_{\hbar}^1, z_{\hbar}^2)\in {\mathit{\boldsymbol{V}}}_{\hbar}$.

Lemma 3.5. Let $\{\boldsymbol\lambda^n\}_{n = 0}^\infty$ be a function sequence on ${\mathit{\boldsymbol{V}}}_{\hbar}$, then we have

Lemma 3.6. ^[40] Let $\varphi^k\geq0$, $k = 0, 1, \ldots, N$, $\zeta > 0$ and $C_0\geq 1$ be two constants, which satisfy

Then, the following relation holds

Furthermore, the above result can be further written as

Lemma 3.7. ^[41] Let $\varphi^n$ be a function on $\Omega$, then

Lemma 3.8. ^[41] Let $\varphi^n, \varsigma^n\geq 0$, $n = 0, 1, \cdots,$ satisfy

where $\lambda_1, \lambda_2\geq 0$ are two constants independent of $\tau$. There exists a constant $\tau^* > 0$ such that, if $\tau \leq \tau^*$, then

where $E_\alpha(z) = \sum\limits_{k = 0}^\infty\frac{z^k}{\Gamma(1+k\alpha)}$ is the Mittag-Leffler function and $\lambda = \lambda_1+\frac{\lambda_2}{(2-2^{1-\alpha})}$.

Lemma 3.9. ^[42] [Browder fixed point theorem] Let $S$ be a finite dimensional space with the inner product $(\cdot, \cdot)_S$ and the norm $\|\cdot\|_S$, and the map $G: S\rightarrow S$ be continuous. Suppose the there exists $\mu > 0$ such that $(G(\xi), \xi)_S\geq 0$ for $\forall \xi \in S$ with $\|\xi\|_S = \mu$. Then, there exists $\xi^* \in S$ such that $G(\xi^*) = 0$ and $\|\xi^*\|_S\leq \mu$.

We first give the existence and uniqueness results for the nonlinear MFVE scheme (2.8) by using Lemma 3.9.

Theorem 3.1. Assume that $(\boldsymbol \lambda_H^i, u_H^i)$ $(i = 0, 1, \cdots, n-1)$ are given. There exists a constant $\tau_0 > 0$ such that, if $\tau < \tau_0$, then there exists a unique solution $(\boldsymbol \lambda_H^n, u_H^n)\in {\mathit{\boldsymbol{V}}}_H \times W_H$ for the nonlinear MFVE scheme (2.8) on the coarse primal and dual grids.

Proof. Let $G: {\mathit{\boldsymbol{V}}}_H \times W_H \rightarrow {\mathit{\boldsymbol{V}}}_H \times W_H$ be the map. For $\bar{\boldsymbol \lambda}_H, \bar{u}_H\in {\mathit{\boldsymbol{V}}}_H \times W_H$, we define $G(\bar{\boldsymbol \lambda}_H, \bar{u}_H)$ as follows:

The map $G$ is obviously continuous. Furthermore, setting ${\mathit{\boldsymbol{v}}}_H = \bar{\boldsymbol\lambda}_H, w_H = \bar{u}_H$ in (3.1), and applying Lemma 3.3, we have

Noting that $\sum\limits_{k = 0} ^{n-1}b_k^n = -1$, we have

Thus, there exists a constant $\tau_{0, 1} > 0$ such that, if $\tau < \tau_{0, 1}$, then

Because of the norm equivalence in finite dimensional normed linear space, there exists a constant $C_0 > 0$ such that $\|\bar{\boldsymbol\lambda}_H\|\geq C_0 \|\bar{\boldsymbol\lambda}_H\|_{{\mathit{\boldsymbol{H}}}({{{\rm{div}}}})}$. Thus, we have

where $\|(\bar{\boldsymbol\lambda}_H, \bar{u}_H)\|_{{\mathit{\boldsymbol{V}}}_H \times W_H}^2 = \|\bar{\boldsymbol\lambda}_H\|_{{\mathit{\boldsymbol{H}}}({{{\rm{div}}}})}^2+\|\bar{u}\|^2$ and $C_1 = \min\{\frac{1}{4\Gamma(2-\alpha)}, \mu_1 C_0^2 \tau^\alpha\}$. Let $\mu = \frac{1}{C_1}(1+\frac{\tau^\alpha}{2}\|f^n\|^2-\frac{1} {2\Gamma(2-\alpha)} \sum\limits_{k = 0}^{n-1}b_k^n\|u_H^k\|^2)$. Based on above analysis, we know that if $\|(\bar{\boldsymbol\lambda}_H, \bar{u}_H)\|_{{\mathit{\boldsymbol{V}}}_H \times W_H}^2 = \mu$, then $(G(\bar{\boldsymbol\lambda}_H, \bar{u}_H), (\bar{\boldsymbol\lambda}_H, \bar{u}_H))_{{\mathit{\boldsymbol{V}}}_H \times W_H}\geq 0$. Applying Lemma 3.9, we obtain that there exists $(\bar{\boldsymbol\lambda}_H^*, \bar{u}_H^*) \in {\mathit{\boldsymbol{V}}}_H \times W_H$ such that $G(\bar{\boldsymbol\lambda}_H^*, \bar{u}_H^*) = {\mathit{\boldsymbol{0}}}$. Then $(\boldsymbol\lambda_H^n, u_H^n) = (\bar{\boldsymbol\lambda}_H^*, \bar{u}_H^*)$ satisfies (2.8).

Next, we prove the uniqueness of the solution. Let $(\boldsymbol \Lambda_H^n, U_H^n)\in {\mathit{\boldsymbol{V}}}_H \times W_H$ be another solution of (2.8), and $(\boldsymbol \Lambda_H^0, U_H^0) = (\boldsymbol \lambda_H^0, u_H^0)$. Making use of (2.8), we obtain

where $p_H^n = u_H^n-U_H^n, {\mathit{\boldsymbol{q}}}_H^n = \boldsymbol\lambda_H^n-\boldsymbol\Lambda_H^n$. Choose $w_h = p_H^n, {\mathit{\boldsymbol{v}}}_h = {\mathit{\boldsymbol{q}}}_H^n$ in (3.6) to obtain

Applying Lemma 3.3, we have

There exits a constant $\tau_{0, 2} > 0$ such that, if $\tau \leq \tau_{0, 2}$, then $\|g\|_{1, \infty}\tau^\alpha \leq\frac{1} {2\Gamma(2-\alpha)}$, and

It follows that $\|p_H^n\| = 0$ and $\|{\mathit{\boldsymbol{q}}}_H^n\| = {\mathit{\boldsymbol{0}}}$. Setting $\tau_0 = \min\{\tau_{0, 1}, \tau_{0, 2}\}$, we have completed the proof of the theorem.

Now, we give the existence and uniqueness results for the linearized scheme (2.9).

Theorem 3.2. Assume that $(\hat{\boldsymbol \lambda}_h^i, \hat{u}_h^i)$ $(i = 0, 1, \cdots, n-1)$ are given. There exists a constant $\tau_1$ $(0 < \tau_1\leq\tau_0)$ such that, if $\tau < \tau_1$, then there exists a unique solution $(\hat{\boldsymbol \lambda}_h^n, \hat u_h^n)\in {\mathit{\boldsymbol{V}}}_h \times W_h$ for linearized scheme (2.9) on the fine primal and dual grids.

Proof. Let $\{\boldsymbol\phi_i\}_{i = 1}^{M_1}$ and $\{\varphi_j\}_{j = 1}^{M_2}$ be the basis functions of ${\mathit{\boldsymbol{V}}}_h$ and $W_h$, respectively. Then $(\hat{\boldsymbol\lambda}_h^n, \hat{u}_h^n)$ can be expressed as

Substituting (3.10) into (2.9), and taking ${\mathit{\boldsymbol{v}}}_h = \boldsymbol\phi_i$ $(i = 1, 2, \cdots, M_1)$ and $w_h = \varphi_j$ $(j = 1, 2, \cdots, M_2)$, we have

where

Noting that $B_1$ and $B_2$ are invertible, and applying the multiplication of partitioned matrices, we can get

Let $\psi(\tau) = \det(\frac{1}{\Gamma(2-\alpha)}B_1+\tau^\alpha B_3+\tau^\alpha C B_2^{-1}C^T)$, then $\psi(\tau)$ is a continuous function. Noting that $\psi(0) = \det(\frac{1} {\Gamma(2-\alpha)}B_1) > 0$. According to the property of continuous function, there exists a constant $\tau_1$ $(0 < \tau_1\leq\tau_0)$ such that, if $\tau < \tau_1$, then $\psi(\tau) > \frac{1}{2}\det(\frac{1} {\Gamma(2-\alpha)}B_1) > 0$. So the coefficient matrix of (3.11) is invertible, then there exists a unique solution for the linearized scheme (2.9).

4.   Stability analysis

We will give the stability results for the nonlinear MFVE scheme (2.8) and linearized MFVE scheme (2.9) on the coarse and fine grids, respectively.

Theorem 4.1. Let $(\boldsymbol \lambda_H^n, u_H^n)_{n = 0}^{N}\in {\mathit{\boldsymbol{V}}}_H \times W_H$ be the solution of system (2.8), then there exist a constant $C$ independent of $H$ and $\tau$ such that

Moreover, for a constant $c_0 > 0$, there exist a constant $\tau_2 > 0$ independent of $H$ and $\tau$ such that, if $H\leq c_0\tau\leq c_0\min\{\tau_2, \tau_0\}$ and $H\leq\hbar_0$, then

where $\tau_0$ is defined in Theorem 3.1, $\hbar_0 = \frac{\mu_1}{2\mu_3}$, $C > 0$ is a constant independent of $H$, $\tau$ and $c_0$.

Proof. Choosing $w_H = u_H^n$ and ${\mathit{\boldsymbol{v}}}_H = \boldsymbol\lambda_H^n$ in (2.8), we can get

Apply the Lemma 3.3 and Lemma 3.7 in (4.3) to obtain

Apply Lemma 3.8 in (4.4) to obtain

Now, making use of (2.8)$(b)$, we have

Choosing $w_H = D_{\tau}^\alpha u_H^n$ in (2.8)$(a)$ and ${\mathit{\boldsymbol{v}}}_H = \boldsymbol\lambda_H^n$ in (4.6), we have

Applying Lemma 3.5 in (4.7), and noting that $b_k^n < 0$ ($0\leq k < n$), we have

Apply Lemma 3.2 and Lemma 3.3 in (4.8) to obtain

Substituting (4.9) into (4.8), we have

Setting $\hbar_0 = \frac{\mu_1}{2\mu_3}$, when $H\leq\hbar_0$, we have $1-\frac{\mu_3}{\mu_1}H\geq \frac{1}{2}$ and

Applying Lemma 3.6 to have

Let $c_0 > 0$ be a constant. Selecting $H$ and $\tau$ to satisfy $H\leq c_0\tau$ in (4.12), we have

Noting that

then, there exists a constant $\tau_2 > 0$ such that, if $\tau < \min\{\tau_2, \tau_0\}$, where $\tau_0$ is defined in Theorem 3.1, then

Thus, we complete the proof of this theorem.

Theorem 4.2. Let $(\hat{\boldsymbol\lambda}_h^n, \hat u_h^n)_{n = 0}^{N}\in {\mathit{\boldsymbol{V}}}_h \times W_h$ be the solution of system (2.9), then there exist a constant $C > 0$ independent of $h$ and $\tau$ such that

Moreover, there exists a constant $C > 0$ independent of $h$, $\tau$ and $c_0$ such that, if $h\leq c_0\tau\leq c_0\min\{\tau_2, \tau_1\}$ and $h < \hbar_0$, then

where $\tau_1$ is defined in Theorem 3.2, $c_0, \hbar_0$ and $\tau_2$ are defined in Theorem 4.1.

Proof. Choosing $w_h = \hat u_h^n$ and ${\mathit{\boldsymbol{v}}}_h = \hat{\boldsymbol\lambda}_h^n$ in (2.9), we have

Apply Lemma 3.3 and Lemma 3.7 in (4.18) to obtain

Apply Lemma 3.8 and Theorem 4.1 to obtain

Now, making use of (2.9)$(b)$, we have

Choosing $w_H = D_{\tau}^\alpha \hat u_h^n$ in (2.9)$(a)$ and ${\mathit{\boldsymbol{v}}}_H = \hat{\boldsymbol\lambda}_h^n$ in (4.21), we have

Apply Lemma 3.5 to get

Noting that $b_k^n < 0$ $(0\leq k < n)$, we have

Making use of (4.1) and (4.20), we can apply the technique of (4.9) to obtain

Selecting $h$ to satisfy $h\leq \hbar_0$, where $\hbar_0 = \frac{\mu_1}{2\mu_3}$, we have $1-\frac{\mu_3} {\mu_1}h\geq \frac{1}{2}$ and

Applying the technique of (4.11)–(4.15), for the positive constant $c_0$ and $\tau_2$ which defined in Theorem 4.1, we can obtain that if $\tau < \min\{\tau_2, \tau_1\}$ and $h\leq c_0\tau$, then

We complete the proof of the stability.

Remark 4.1. (I) In Theorems 4.1 and 4.2, we also need to select $\tau$ to satisfy $\tau < \tau_0$ and $\tau < \tau_1$, respectively, because of the existence and uniqueness of the MFVE solutions in Theorems 3.1 and 3.2.

(II) From Theorems 4.1 and 4.2, we can see that $\|u_H^n\|$ and $\|\hat{u}_h^n\|$ are unconditionally stable, and $\|\boldsymbol\lambda_H^n\|$ and $\|\hat{\boldsymbol\lambda}_h^n\|$ are conditionally stable, because the bilinear $(\mathcal{A}^{-1}{\mathit{\boldsymbol{z}}}_\hbar, \gamma_q^*{\mathit{\boldsymbol{w}}}_\hbar)$ $(\forall{\mathit{\boldsymbol{z}}}_\hbar, {\mathit{\boldsymbol{w}}}_\hbar\in {\mathit{\boldsymbol{V}}}_\hbar, \hbar = H$ or $h)$ does not necessarily satisfy symmetry. When the coefficient $\mathcal{A}({\mathit{\boldsymbol{x}}})$ is a symmetry and positive definite constant matrix, making use of Lemma 3.2, we can see that $(\mathcal{A}^{-1}{\mathit{\boldsymbol{z}}}_\hbar, \gamma_\hbar^*{\mathit{\boldsymbol{w}}}_\hbar)$ is symmetry, Under this condition, we can obtain that $\|\boldsymbol\lambda_H^n\|$ and $\|\hat{\boldsymbol\lambda}_h^n\|$ are also unconditionally stable.

5.   Error estimates

In order to get the error estimates for two-grid MFVE systems (2.8) and (2.9), we should introduce the standard $L^2$-projection ^[44] $P_{\hbar} : W\rightarrow W_{\hbar}$, which satisfies

where $\hbar = H$ or $h$.

We introduce a generalized MFVE projection $(\tilde{\boldsymbol\lambda}_\hbar, \tilde{u}_\hbar): \bar{J}\rightarrow {\mathit{\boldsymbol{V}}}_\hbar\times W_\hbar$ such that, for $\hbar = H$ or $h$,

Then the above projection satisfies the following estimates.

Lemma 5.1. ^[43] There exists a constant $C > 0$ independent of $\hbar$ and $t$ such that, for $j = 0, 1$ and $\hbar = H$ or $h$

where ${{\mathit{\boldsymbol{H}}}}^1({{{\rm{div}}}}, \Omega) = \{{\mathit{\boldsymbol{v}}}\in (L^2(\Omega))^2: {{{\rm{div}}}} {\mathit{\boldsymbol{v}}}\in H^1(\Omega)\}$.

Lemma 5.2. ^[38] For $2 < q \leq \infty$ and $\hbar = H$ or $h$, the following $L^q$-estimate holds

Moreover, for $j = 0, 1$, the following superconvergence result holds

Now, let $\beta^n = u^n-\tilde{u}_H^n$, $\sigma^n = \tilde{u}_H^n -u_H^n$, $\boldsymbol\zeta^n = \boldsymbol\lambda^n-\tilde{\boldsymbol\lambda}_H^n$, $\boldsymbol\delta^n = \tilde{\boldsymbol\lambda}_H^n -\boldsymbol\lambda_H^n$, where $(\tilde{\boldsymbol\lambda}_H^n, \tilde{u}_H^n)\in {\mathit{\boldsymbol{V}}}_H\times W_H$ is the generalized MFVE projection of $(\boldsymbol\lambda, u)$, then we can obtain the error equations as follows

Theorem 5.1. Let $(\boldsymbol \lambda_H^n, u_H^n)\in {\mathit{\boldsymbol{V}}}_H \times W_H$ and $(\boldsymbol \lambda^n, u^n)\in {\mathit{\boldsymbol{V}}} \times W$ be the solutions of systems (2.8) and (2.2), respectively. Assume that $u, {{{\rm{div}}}}\boldsymbol\lambda\in \mathcal{C}^2(\bar{J}, H^1(\Omega))$, $\boldsymbol\lambda\in\mathcal{C}^2(\bar{J}, (H^1(\Omega))^2)$, and $(\boldsymbol \lambda_H^0, u_H^0) = (\tilde{\boldsymbol \lambda}_H^0, \tilde{u}_H^0)$, then there exists a constant $C > 0$ independent of $H$ and $\tau$ such that

Moreover, there exist a constant $C > 0$ independent of $H$, $\tau$ and $c_0$ such that, if $H\leq c_0\tau\leq c_0\min\{\tau_2, \tau_0\}$ and $H < \hbar_0$, then

where $\tau_0$ is defined in Theorem 3.1, $c_0, \hbar_0$ and $\tau_2$ are defined in Theorem 4.1.

Proof. Choosing ${\mathit{\boldsymbol{v}}}_H = \boldsymbol \delta^n$ and $w_H = \sigma^n$ in (5.4), we have

Making use of the Lagrange mean value theorem, $L^2$-projection $P_H$ and Lemma 5.2, we have

where $u_*^n$ is located between $u^n$ and $u_H^n$. And we can also obtain

Substituting (5.11) and (5.12) into (5.10), making use of Lemma 3.7, we obtain

Noting that $\sigma^0 = 0$, applying Lemma 3.8, we obtain

Making use of (5.4)$(b)$, we have

Choosing $w_H = D_{\tau}^\alpha \sigma^n$ in (5.4)$(a)$ and ${\mathit{\boldsymbol{v}}}_H = \boldsymbol\delta^n$ in (5.15), we have

For the term $-(g(u^n)-g(u^n_H), D_{\tau}^\alpha \sigma^n)$, similar to (5.11), we have

For the term $-(D_{\tau}^\alpha \beta^n, D_{\tau}^\alpha \sigma^n)$, similar to (5.12), we have

Substituting (5.17) and (5.18) into (5.16), applying Lemma 3.5, we obtain

Noting that $b_k^n < 0$ $(0\leq k < n)$, making use of (4.9) and (5.14), we have

Selecting $H$ to satisfy $H\leq \hbar_0$, where $\hbar_0 = \frac{\mu_1}{2\mu_3}$, we have $1-\frac{\mu_3} {\mu_1}H\geq \frac{1}{2}$ and

Applying the technique of (4.11)–(4.15), for the positive constant $c_0$ and $\tau_2$ which defined in Theorem 4.1, noting that $\boldsymbol\delta^0 = {\mathit{\boldsymbol{0}}}$, we can obtain that if $\tau < \min\{\tau_2, \tau_0\}$ and $H\leq c_0\tau$, then

Finally, we estimate $\|\boldsymbol\lambda^n-\boldsymbol\lambda_H^n\|_{{\mathit{\boldsymbol{H}}}({{{\rm{div}}}}, \Omega)}$. Choosing $w_H = {{{\rm{div}}}} \boldsymbol\delta^n$ in (5.4)$(a)$ and $v_H = \boldsymbol\delta^n$ in (5.15), we have

Noting that

similar to the proof of (5.17) and (5.18), we obtain

Making use of (5.14), we have

Then, apply Lemmas 5.1 and 5.2 to complete the proof.

Remark 5.1. For $2 < q\leq\infty$, making use of the inverse estimate and Lemma 5.2, we obtain

Moreover, when $q = 4$, we have

which will be applied to the following estimates for the linearized MFVE scheme (2.9).

Next, we give the error estimates for the linearized MFVE scheme (2.9). Let $\vartheta^n = u^n-\tilde{u}_h^n$, $\xi^n = \tilde{u}^n_h-\hat{u}_h^n$, $\boldsymbol\rho^n = \boldsymbol\lambda^n-\tilde{\boldsymbol\lambda}_h^n$, $\boldsymbol\theta^n = \tilde{\boldsymbol\lambda}^n_h-\hat{\boldsymbol\lambda}_h^n$, where $(\tilde{\boldsymbol\lambda}_h^n, \tilde{u}_h^n)\in {\mathit{\boldsymbol{V}}}_h\times W_h$ is the generalized MFVE projection of $(\boldsymbol\lambda, u)$, then we can obtain the error equations as follows

where $G = g(u^n)-g(u^n_H)-g'(u^n_H)(\hat u_h^n-u^n_H)$.

Theorem 5.2. Let $(\hat{\boldsymbol \lambda}_h^n, \hat u_h^n)\in {\mathit{\boldsymbol{V}}}_h \times W_h$ and $(\boldsymbol \lambda^n, u^n)\in {\mathit{\boldsymbol{V}}} \times W$ be the solutions of systems (2.9) and (2.2), respectively. Assume that $u\in \mathcal{C}^2(\bar{J}, W^{1, 4}(\Omega))$, $\boldsymbol\lambda\in\mathcal{C}^2(\bar{J}, (H^1(\Omega))^2)$, ${{{\rm{div}}}}\boldsymbol\lambda\in \mathcal{C}^2(\bar{J}, H^1(\Omega))$, and $(\hat{\boldsymbol \lambda}_h^0, \hat{u}_h^0) = (\tilde{\boldsymbol \lambda}_h^0, \tilde{u}_h^0)$, then there exists a constant $C > 0$ independent of $h$ and $\tau$ such that

Moreover, there exist a constant $C > 0$ independent of $h$, $\tau$ and $c_0$ such that, if $h\leq c_0\tau\leq c_0\min\{\tau_2, \tau_1\}$ and $h < \hbar_0$, then

where $\tau_1$ is defined in Theorem 3.2, $c_0, \hbar_0$ and $\tau_2$ are defined in Theorem 4.1.

Proof. Choosing ${\mathit{\boldsymbol{v}}}_h = \boldsymbol \theta^n$ and $w_h = \xi^n$ in (5.28), we have

Making use of the Taylor expansion for $g(u^n)$ on $u = u_H^n$, we have

where $u_\diamond^n$ is located between $u^n$ and $u_H^n$. Noting that

similar to (5.12) for $(D_{\tau}^\alpha \vartheta^n, \xi^n)$, applying Lemma 3.7, we obtain

Applying Remark 5.1, Lemma 5.1 and Lemma 3.8, noting that $\xi^0 = 0$, we obtain

Now, making use of (5.28)$(b)$, we have

Choosing $w_h = D_{\tau}^\alpha\xi^n$ in (5.28)$(a)$ and ${\mathit{\boldsymbol{v}}}_h = \boldsymbol \theta^n$ in (5.37), we have

For the term $-(G, D_{\tau}^\alpha \xi^n)$, similar to (5.34), we have

Similar to (5.12) for $(D_{\tau}^\alpha \vartheta^n, D_{\tau}^\alpha \xi^n)$, applying Lemma 3.5 in (5.38), we obtain

Noting that $b_k^n < 0$ $(0\leq k < n)$, and making use of (5.36), we get

Selecting $h$ to satisfy $h\leq \hbar_0$, where $\hbar_0 = \frac{\mu_1}{2\mu_3}$, we have $1-\frac{\mu_3}{\mu_1}h\geq \frac{1}{2}$ and

Applying the technique of (4.11)–(4.15), for the positive constant $c_0$ and $\tau_2$ which defined in Theorem 4.1, we can obtain that if $\tau < \min\{\tau_2, \tau_1\}$ and $h\leq c_0\tau$, then

Apply Lemma 5.1 to complete the proof of (5.29) and (5.30).

Then, we estimate $\|{{{\rm{div}}}}(\boldsymbol\lambda^n-\boldsymbol\lambda_h^n)\|$. Choosing $w_h = {{{\rm{div}}}} \boldsymbol \theta^n$ in (5.28) and ${\mathit{\boldsymbol{v}}}_h = \boldsymbol \theta^n$ in (5.37), we have

Apply Lemma 3.7 to obtain

Applying the technique of (5.24) and (5.25), we have

Finally, apply Lemma 5.1 and Remark 5.1 to complete the proof of (5.31).

Remark 5.2. (I) In Theorem 5.1, we should assume that $u, {{{\rm{div}}}}\boldsymbol\lambda\in \mathcal{C}^2(\bar{J}, H^1(\Omega))$, $\boldsymbol\lambda\in\mathcal{C}^2(\bar{J}, (H^1(\Omega))^2)$. We also need add the regularity $u\in \mathcal{C}^2(\bar{J}, W^{1, 4}(\Omega))$ in Theorem 5.2. Moreover, it should be pointed out that the solutions of FDEs usually show the initial weak singularity, some numerical methods ^{[45,46,47,48,49,50]} were proposed to deal with this problem.

(II) Similar to Remark 4.1, when the coefficient $\mathcal{A}({\mathit{\boldsymbol{x}}})$ is a symmetry and positive definite constant matrix, we can remove the conditions $H\leq c_0\tau$ and $h\leq c_0\tau$ in the analysis and results of Theorems 5.1 and 5.2, respectively.

6.   Numerical examples

In this section, we will give two examples with some numerical results to test the convergence rates and the influence of the fractional parameters. In (1.1), we choose $\Omega = (0, 1)^2$, $J = (0, T]$, $\mathcal{A}({\mathit{\boldsymbol{x}}})$ as the identity matrix, and the exact solution (similar as in ^[12,31])

where $\varpi$ is a parameter. Then, we can get the auxiliary variable

and the source function

Example 6.1. By choosing $T = 1$, $g(u) = \sin(u)$, and $\varpi = 2$, we carry out numerical simulation for some different fractional parameters $\alpha = 0.2, 0.4, 0.6.0.8$ and grid sizes. In Tables 1 and 2, we take $\tau = 1/5, 1/8, 1/10$, $h\thickapprox \sqrt{2}\tau^{2-\alpha}, H^2\thickapprox 2\tau^{2-\alpha}$ (in two-grid MFVE algorithm), and $h\thickapprox \sqrt{2}\tau^{2-\alpha}$ (in MFVE algorithm (2.7)), and obtain that the convergence rates in time direction are close to $2-\alpha$ for $u$ in $L^2(\Omega)$-norm and $\boldsymbol\lambda$ in $(L^2(\Omega))^2$ and ${\mathit{\boldsymbol{H}}}({{{\rm{div}}}}, \Omega)$-norms, which is consistent with the theoretical results in Theorems 5.1 and 5.2. For testing convergence rates in space direction, we fix the time step length $\tau = 1/100$, select the coarse and fine grid sizes to satisfy $h = H^2/\sqrt{2} = \sqrt{2}/4, \sqrt{2}/16, \sqrt{2}/25, \sqrt{2}/36$, and give numerical results and computing time for the two-grid MFVE algorithm in Table 3. At the same time, in Table 4, we give some numerical results for the MFVE algorithm (2.7) with grid sizes $h = \sqrt{2}/4, \sqrt{2}/16, \sqrt{2}/25, \sqrt{2}/36$. We can see that the convergence rates are close to $1$. Moreover, we choose the coarse and fine grid sizes to satisfy $h = H^2/\sqrt{2} = \sqrt{2}/4, \sqrt{2}/16, \sqrt{2}/25, \sqrt{2}/36$ and $h = \sqrt{2}\tau$, give numerical results for the two-grid MFVE algorithm in Table 5, and the corresponding numerical results for the MFVE algorithm (2.7) in Table 6. Then we obtain the same conclusions as that discussed in Tables 3 and 4. Furthermore, for the time parameter $t = 1$, we show the graphs of the exact solutions for $u$ and $\boldsymbol \lambda$ with $h = \sqrt{2}/32$ in Figures 2 and 4, respectively, also show the graphs of the numerical solutions based on the two-grid MFVE algorithm with $h = \sqrt{2}\tau = H^2/\sqrt{2} = \sqrt{2}/25$ in Figures 3 and 5. We find that the numerical solutions and the exact solutions have the same numerical behaviors.

Example 6.2. In this example, we take $T = 1$, $g(u) = u^3-u$, and $\varpi = 2+\alpha$, then obtain the exact solution $u({\mathit{\boldsymbol{x}}}, t) = t^{2+\alpha}\sin(2\pi x_1)\sin(2\pi x_2), {\mathit{\boldsymbol{x}}} = (x_1, x_2)\in [0, 1]^2, t\in [0, 1]$, the auxiliary variable $\boldsymbol \lambda({\mathit{\boldsymbol{x}}}, t) = -\nabla u({\mathit{\boldsymbol{x}}}, t)$. For some different fractional parameters $\alpha = 0.2, 0.4, 0.6.0.8$ and grid sizes, we conduct numerical experiments as in Example 6.1. For the two-grid MFVE algorithm and MFVE algorithm (2.7), we can see that the convergence rates in time direction are close to $2-\alpha$ (in Tables 7 and 8), and the convergence rates in space direction are close to $1$ (in Tables 9 and 10). Moreover, in Tables 11 and 12, we choose $h = \sqrt{2}\tau = H^2/\sqrt{2}$ (in two-grid MFVE algorithm) and $h = \sqrt{2}\tau$ (in MFVE algorithm), then obtain the same convergence rates as in Tables 9 and 10.

Base on the above the numerical results in Tables 1–6 and Figures 2–5 for Example 6.1 and Tables 7–12 for Example 6.2, we can know that the convergence rates are consistent with the theoretical results in Theorems 5.1 and 5.2. We also find that the two-grid MFVE algorithm can save the computing time compared with the MFVE algorithm while maintaining the same convergence rates. Finally, numerical results and the figures show that the proposed two-grid MFVE algorithm for the nonlinear time fractional reaction-diffusion equations is feasible and effective.

7.   Conclusions

In this paper, we construct the two-grid MFVE fast algorithm to solve the nonlinear time fractional reaction-diffusion equations with the Caputo time fractional derivative. We obtain the stability results and the optimal error estimates for $u$ (in $L^2(\Omega)$-norm) and $\boldsymbol \lambda$ (in $(L^2(\Omega))^2$-norm), and the sub-optimal error estimates for $\boldsymbol \lambda$ (in $H({{{\rm{div}}}}, \Omega)$-norm). Furthermore, we also give two numerical examples to verify that the proposed algorithm can greatly save the computing time. In future works, for the Caputo fractional derivative (1.2) with $\alpha\in (0, 1)$, we will try to use other approximation methods (such as $L1$-$2$, $L2$-$1_{\sigma}$, $L1$-$2$-$3$ formulas ^{[17,18,19,20]}) and the two-grid MFVE method to solve more fractional partial differential equations in scientific and engineering fields.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (11701299, 12161063, 11761053), the Natural Science Foundation of Inner Mongolia Autonomous Region (2020MS01003, 2021MS01018), the Prairie Talent Project of Inner Mongolia Autonomous Region, and the Postgraduate Scientific Research Innovation Support Project of Inner Mongolia University.

Conflict of interest

The authors declared that they have no conflicts of interest to this work.