A Mixed Finite Element Method for the Multi-Term Time-Fractional Reaction–Diffusion Equations

Abstract

In this work, a fully discrete mixed finite element (MFE) scheme is designed to solve the multi-term time-fractional reaction–diffusion equations with variable coefficients by using the well-known $L1$ formula and the Raviart–Thomas MFE space. The existence and uniqueness of the discrete solution is proved by using the matrix theory, and the unconditional stability is also discussed in detail. By introducing the mixed elliptic projection, the error estimates for the unknown variable u in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm and for the auxiliary variable $\mathit{\lambda }$ in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms are obtained. Finally, three numerical examples are given to demonstrate the theoretical results.

1. Introduction

Fractional calculus and fractional partial differential equations (FPDEs) have been confirmed to be very important tools in describing some anomalous phenomena and processes with memory and nonlocal properties [1,2,3,4,5,6]. Moreover, some underlying and complex processes can be described more appropriately by multi-term FPDEs [7,8,9], as they contains multiple fractional derivative or calculus terms. In recent years, many numerical methods have been increasingly used by scholars to solve multi-term FPDEs. Liu et al. [10] constructed some finite difference (FD) schemes to solve the multi-term time-fractional wave-diffusion equations by using two fractional predictor–corrector methods. Dehghan et al. [11] devised two high-order numerical schemes to solve the multi-term time-fractional diffusion-wave equations by using the compact FD method and Galerkin spectral technique. Ren and Sun [12] established an efficient compact FD scheme for the multi-term time-fractional diffusion-wave equation by using the $L1$ formula. Zheng et al. [13] proposed a high-order space–time spectral method for the multi-term time-fractional diffusion equations by using the Legendre polynomials in the temporal direction and the Fourier-like basis functions in the spatial direction. Du and Sun [14] constructed an FD scheme for multi-term time-fractional mixed diffusion and wave equations by using the $L2-1_{\sigma }$ formula. Hendy and Zaky [15] proposed a spectral method for a coupled system of nonlinear multi-term time–space fractional diffusion equations by using the $L1$ formula on a time-graded mesh. Liu et al. [16] developed an ADI Legendre spectral method for solving a multi-term time-fractional Oldroyd-B fluid-type diffusion equation. Wei and Wang [17] constructed a higher-order numerical scheme for the multi-term variable-order time-fractional diffusion equation by using the local discontinuous Galerkin method. She et al. [18] considered a spectral method for solving the multi-term time-fractional diffusion problem by using a modified $L1$ formula.

Meanwhile, many scholars selected the finite element (FE) method for solving the multi-term FPDEs and have achieved excellent results. Jin et al. [19] developed an FE method for a multi-term time-fractional diffusion equation and considered the case of smooth and nonsmooth initial data. Li et al. [20] proposed an FE method to solve a higher-dimensional multi-term fractional diffusion equation on nonuniform time meshes. Zhou et al. [21] developed a weak Galerkin FE method for solving multi-term time-fractional diffusion equations by using a convolution quadrature formula. Bu et al. [22] proposed a space–time FE method for solving the multi-term time–space fractional diffusion equation based on the suitable graded time mesh. Feng et al. [23] proposed an FE method for a multi-term time-fractional mixed subdiffusion and diffusion-wave equation on the convex domain by using mixed L-type schemes. Meng and Stynes [24] considered an $L1$ FE method for a multi-term time-fractional initial-boundary value problem on the temporal graded mesh. Yin et al. [25] constructed a class of efficient time-stepping FE schemes for multi-term time-fractional reaction–diffusion-wave equations by using the shifted convolution quadrature method. Huang et al. [26] proposed an $\alpha$-robust FE method for a multi-term time-fractional diffusion problem on a graded mesh by using the $L1$ formula. Liu et al. [27] proposed an FE method for solving a multi-term variable-order time-fractional diffusion equation and developed an efficient parallel-in-time algorithm to reduce the computational costs.

In this work, we will construct a fully discrete mixed finite element (MFE) scheme for the following multi-term time-fractional reaction–diffusion (TFRD) equations with variable coefficients:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}P\left(D_t\right)u(\mathit{x},t)-\mathrm{div}(\mathcal{A}\left(\mathit{x}\right)\nabla u(\mathit{x},t))+p\left(\mathit{x}\right)u(\mathit{x},t)=f(\mathit{x},t),\hfill & (\mathit{x},t)\in \Omega \times J,\hfill \\ u(\mathit{x},t)=0,\hfill & (\mathit{x},t)\in \partial \Omega \times \overline{J},\hfill \\ u(\mathit{x},0)=u_0\left(\mathit{x}\right),\hfill & \mathit{x}\in \overline{\Omega },\hfill \end{array}\right.\end{array}$$

(1)

where $\Omega \subset R^2$ is a convex and bounded polygon region with boundary $\partial \Omega$, $J=(0,T]$ with $0<T<\infty$. We assume that the source function $f(\mathit{x},t)$, initial data $u_0\left(\mathit{x}\right)$, and non-negative coefficient $p\left(\mathit{x}\right)$ are smooth enough. Specifically, for the symmetric diffusion coefficient matrix $\mathcal{A}\left(\mathit{x}\right)$, we should assume that there exist two constants $A_0$, $A_1>0$ such that

$$\begin{array}{c}\hfill A_0{\mathit{z}}^T\mathit{z}\le {\mathit{z}}^T\mathcal{A}\left(\mathit{x}\right)\mathit{z}\le A_1{\mathit{z}}^T\mathit{z},\forall \mathit{z}\in R^2,\forall \mathit{x}\in \overline{\Omega }.\end{array}$$

Moreover, the multi-term time-fractional derivative $P\left(D_t\right)u(\mathit{x},t)$ is defined by

$$\begin{array}{c}\hfill P\left(D_t\right)u(\mathit{x},t)=\sum _{i=1}^m{b}_i{D}_t^{{\alpha }_i}u(\mathit{x},t),0<{\alpha }_m<{\alpha }_{m-1}<\cdots <{\alpha }_1<1,\end{array}$$

where $b_i$   $(i=1$, 2, ⋯, $m)$ are the positive real numbers and $D_t^{{\alpha }_i}u$ is the Caputo time-fractional derivative as follows:

$$\begin{array}{cc}\hfill D_t^{{\alpha }_i}u(\mathit{x},t)=& \frac{1}{\Gamma (1-{\alpha }_i)}{\int }_0^t\frac{\partial u(x,s)}{\partial s}\frac{1}{{(t-s)}^{{\alpha }_i}}\mathrm{d}s,\hfill \end{array}$$

where $\Gamma (·)$ denotes the $\Gamma$-function.

It should be noted that the MFE method, as an important numerical calculation method, has been widely used to solve FPDEs [28,29,30,31,32], and some scholars have also used this method to solve the multi-term FPDEs [33,34,35]. In [33], Shi et al. proposed an $H^1$ -Galerkin mixed finite element (MFE) method for the multi-term time-fractional diffusion equations and gave a superconvergence result. In [34], Li et al. proposed an MFE method for the multi-term time-fractional diffusion and diffusion-wave equations by using an MFE space contained in ${\left(L^2(\Omega )\right)}^d\times H_0^1(\Omega )$, where $d=2,3$. In [35], Cao et al. constructed a nonconforming MFE scheme for the multi-term time-fractional mixed diffusion and diffusion-wave equations. Motivated by the above excellent works, we will construct a fully discrete MFE scheme for the multi-term TFRD equation (1) by using the Raviart–Thomas MFE space and the $L1$ formula, analyze the existence, uniqueness, and unconditional stability in detail, and give error estimates for u (in discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and auxiliary variable $\mathit{\lambda }$ (in discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and discrete $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms). Finally, we give numerical experiments to demonstrate the efficiency of the proposed method.

The remainder of this paper is arranged as follows. In Section 2, we construct a fully discrete MFE scheme for the multi-term TFRD equations by using the Raviart–Thomas MFE space and the $L1$ -formula. In Section 3, we give a fractional Grönwall inequality and analyze the existence and uniqueness of the discrete solution. We derive the unconditional stability results and a priori error estimates in detail in Section 4 and Section 5, respectively. Finally, three numerical examples are given to verify the theoretical results.

2. Mixed Finite Element Method

We introduce the flux $\mathit{\lambda }(\mathit{x},t)=-\mathcal{A}\left(\mathit{x}\right)\nabla u(\mathit{x},t)$ as the auxiliary variable. Then, the original problem (1) can be rewritten as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\left(a\right)P\left(D_t\right)u(\mathit{x},t)+\mathrm{div}\mathit{\lambda }(\mathit{x},t)+p\left(\mathit{x}\right)u(\mathit{x},t)=f(\mathit{x},t),\hfill & (\mathit{x},t)\in \Omega \times J,\hfill \\ \left(b\right){\mathcal{A}}^{-1}\left(x\right)\mathit{\lambda }(\mathit{x},t)+\nabla u(\mathit{x},t)=0,\hfill & (\mathit{x},t)\in \Omega \times J,\hfill \\ \left(c\right)u(\mathit{x},t)=0,\hfill & (\mathit{x},t)\in \partial \Omega \times \overline{J},\hfill \\ \left(d\right)u(\mathit{x},0)=u_0\left(\mathit{x}\right),\hfill & \mathit{x}\in \overline{\Omega }.\hfill \end{array}\right.\end{array}$$

(2)

Let $V=L^2(\Omega )$ and $\mathit{W}=\mathit{H}(\mathrm{div},\Omega )=\{\mathit{w}\in {\left(L^2(\Omega )\right)}^2:\mathrm{div}\mathit{w}\in L^2(\Omega )\}.$ Then, we obtain the mixed variational formulation of (2): find $(u,\mathit{\lambda })\in V\times \mathit{W}$ such that

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\left(a\right)(P\left(D_t\right)u,v)+(\mathrm{div}\mathit{\lambda },v)+(pu,v)=(f,v),\hfill & \forall v\in V,\hfill \\ \left(b\right)({\mathcal{A}}^{-1}\mathit{\lambda },\mathit{w})-(u,\mathrm{div}\mathit{w})=0,\hfill & \forall \mathit{w}\in \mathit{W},\hfill \\ \left(c\right)u(\mathit{x},0)=u_0\left(\mathit{x}\right),\hfill & \forall \mathit{x}\in \overline{\Omega }.\hfill \end{array}\right.\end{array}$$

(3)

Let $K_h$ be a quasi-uniform triangulation of the domain $\Omega$, $h_T$ be the diameter of the triangle $T\in K_h$ and denote $h=max\left\{h_T\right\}$. We select the Raviart–Thomas MFE space $V_h\times {\mathit{W}}_h\subset V\times \mathit{W}$, that is,

$$\begin{array}{cc}\hfill & V_h\left(K\right)=\{v_h\in V:v_h{|}_T\in P_r\left(T\right),\forall T\in K_h\},\hfill \\ \hfill & {\mathit{W}}_h\left(K\right)=\{{\mathit{w}}_h\in \mathit{W}:{\mathit{w}}_h{|}_T\in {\left(P_r\left(T\right)\right)}^2\oplus \left(\mathit{x}{P}_r\left(T\right)\right),\forall T\in K_h\},\hfill \end{array}$$

where the notation ⊕ indicates a direct sum, $\mathit{x}{P}_r\left(T\right)=(x_1{P}_r\left(T\right),x_2{P}_r\left(T\right)),\mathit{x}=(x_1,x_2)$ and $r\ge 0$ is a given integer.

Let $\tau =T/N$ and $t_n=n\tau$ for $n=0$, 1, 2, ⋯, N, where N is a positive integer. For the parameters ${\alpha }_i$ and $i=1$, 2, ⋯, m, the Caputo time-fraction derivative $D_t^{{\alpha }_i}u(\mathit{x},t)$ at $t=t_n$ is approximated by using the well-known $L1$ formula [36,37] as follows:

$$\begin{array}{cc}\hfill D_t^{{\alpha }_i}u(\mathit{x},t_n)=& \frac{1}{\Gamma (1-{\alpha }_i)}{\int }_0^{t_n}\frac{\partial u(x,s)}{\partial s}\frac{1}{{(t_n-s)}^{{\alpha }_i}}ds\hfill \\ \hfill =& \frac{1}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}\frac{u^{k+1}-u^k}{\tau }[{(t_n-t_k)}^{1-{\alpha }_i}-{(t_n-t_{k+1})}^{1-{\alpha }_i}]+Q_{{\alpha }_i}^n\left(x\right)\hfill \\ \hfill =& \frac{1}{\Gamma (2-{\alpha }_i)}[d_{{\alpha }_i,1}^n{u}^n+\sum _{k=1}^{n-1}(d_{{\alpha }_i,k+1}^n-d_{{\alpha }_i,k}^n)u^{n-k}-d_{{\alpha }_i,n}^n{u}^0]+Q_{{\alpha }_i}^n\left(x\right)\hfill \\ \hfill =& \frac{1}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^n{\tilde{d}}_{{\alpha }_i,k}^n{u}^k+Q_{{\alpha }_i}^n\left(x\right),\hfill \end{array}$$

(4)

where $d_{{\alpha }_i,k}^n=\frac{{(t_n-t_{n-k})}^{1-{\alpha }_i}-{(t_n-t_{n-k+1})}^{1-{\alpha }_i}}{\tau }$, ${\tilde{d}}_{{\alpha }_i,0}^n=-d_{{\alpha }_i,n}^n$, ${\tilde{d}}_{{\alpha }_i,n}^n=-d_{{\alpha }_i,1}^n$, and ${\tilde{d}}_{{\alpha }_i,k}^n=d_{{\alpha }_i,n-k+1}^n-d_{{\alpha }_i,n-k}^n(0<k\le n-1)$. Setting $D_N^{{\alpha }_i}{u}^n=\frac{1}{\Gamma (2-{\alpha }_i)}{\displaystyle \sum _{k=0}^n}{\tilde{d}}_{{\alpha }_i,k}^n{u}^k$, we have $D_t^{{\alpha }_i}u(\mathit{x},t_n)=D_N^{{\alpha }_i}{u}^n+Q_{{\alpha }_i}^n\left(x\right)$, where $Q_{{\alpha }_i}^n\left(x\right)$ is the truncation error.

Based on the above definitions, and setting $u_h^n$ and ${\mathit{\lambda }}_h^n$ to be the discrete solutions of u and $\mathit{\lambda }$ at $t=t_n$, respectively, then we can design a fully discrete MFE scheme for the original problem (1): find $(u_h^n,{\lambda }_h^n)\in V_h\times {\mathit{W}}_h$ such that

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\left(a\right)({\displaystyle \sum _{i=1}^m}{b}_i{D}_N^{{\alpha }_i}{u}_h^n,v_h)+(\mathrm{div}{\mathit{\lambda }}_h^n,v_h)+(p{u}_h^n,v_h)=(f^n,v_h),\hfill & \forall v_h\in V_h,\hfill \\ \left(b\right)({\mathcal{A}}^{-1}{\mathit{\lambda }}_h^n,{\mathit{w}}_h)-(u_h^n,\mathrm{div}{\mathit{w}}_h)=0,\hfill & \forall {\mathit{w}}_h\in {\mathit{W}}_h,\hfill \end{array}\right.\end{array}$$

(5)

where $(u_h^0,{\mathit{\lambda }}_h^0)\in V_h\times {\mathit{W}}_h$ satisfies

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\left(a\right)({\mathcal{A}}^{-1}{\mathit{\lambda }}_h^0,{\mathit{w}}_h)-(u_h^0,\mathrm{div}{\mathit{w}}_h)=0,\hfill & \forall {\mathit{w}}_h\in {\mathit{W}}_h,\hfill \\ \left(b\right)(\mathrm{div}{\mathit{\lambda }}_h^0,v_h)=(\mathrm{div}{\mathit{\lambda }}_0,v_h),\hfill & \forall v_h\in V_h,\hfill \end{array}\right.\end{array}$$

(6)

where ${\mathit{\lambda }}_0\left(\mathit{x}\right)=-\mathcal{A}\left(\mathit{x}\right)\nabla u_0\left(\mathit{x}\right)$.

Remark 1.

(I) In the MFE scheme (5)–(6), we particularly emphasize the calculation of initial values $(u_h^0,{\mathit{\lambda }}_h^0)$, as this calculation will be used in stability and convergence analyses. Moreover, from the mixed elliptic projection $R_h$ defined in Section 5, we can see that $(u_h^0,{\mathit{\lambda }}_h^0)=(R_h{u}_0,R_h{\mathit{\lambda }}_0)$.

(II) Compared with the standard FE methods, it is well known that the MFE method can not only reduce the smoothness requirement of the finite element space, but also simultaneously calculate multiple physical quantities. These advantages are very important and popular in practical applications.

3. Existence and Uniqueness

In this section, we shall prove the existence and uniqueness for the MFE scheme (5)–(6). We first give some lemmas, which are important in subsequent theoretical analysis.

Lemma 1

([38]). There exist two positive constants ${\mu }_0$ and ${\mu }_1$ such that

$$\begin{array}{c}\hfill {\mu }_0{\parallel \mathit{w}\parallel }^2\le {\parallel \mathit{w}\parallel }_{{\mathcal{A}}^{-1}}^2\le {\mu }_1{\parallel \mathit{w}\parallel }^2,\mathit{where}{\parallel \mathit{w}\parallel }_{{\mathcal{A}}^{-1}}^2=({\mathcal{A}}^{-1}\mathit{w},\mathit{w}),\forall \mathit{w}\in \mathit{W}.\end{array}$$

Lemma 2

([28]). Let ${\left\{{\mathit{z}}^n\right\}}_{n=0}^{\infty }$ be a sequence on ${\mathit{W}}_h$. Then, the following identity holds:

$$\begin{array}{cc}\hfill \sum _{k=0}^n{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{z}}^k,{\mathit{z}}^n)=& \frac{1}{2}[{\tilde{d}}_{{\alpha }_i,n}^n({\mathcal{A}}^{-1}{\mathit{z}}^n,{\mathit{z}}^n)+\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{z}}^k,{\mathit{z}}^k)\hfill \\ \hfill & +\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}({\mathit{z}}^n-{\mathit{z}}^k),{\mathit{z}}^n-{\mathit{z}}^k].\hfill \end{array}$$

Lemma 3. 

Let $\{{\phi }^k$ : $0\le k\le N\}$ be a non-negative sequence, $\{{\mathit{\xi }}^k$ : $0\le k\le N\}$ be a nondecreasing positive sequence, and $C_0\ge 1$ be a constant, which satisfy

$$\begin{array}{c}\hfill \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n{\phi }^n\le -C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n{\phi }^k+{\mathit{\xi }}^n,1\le n\le N.\end{array}$$

(7)

Then, we have

$$\begin{array}{c}\hfill {\phi }^n\le C_0^n({\phi }^0+\frac{1}{{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,n}^n}{\mathit{\xi }}^n),1\le n\le N.\end{array}$$

(8)

Further, we can further write (8) as

$$\begin{array}{c}\hfill {\phi }^n\le C_0^n({\phi }^0+\sum _{i=1}^m\frac{\Gamma (1-{\alpha }_i)t_n^{{\alpha }_i}}{b_i}{\mathit{\xi }}^n),1\le n\le N.\end{array}$$

(9)

Proof. 

When $n=1$ in (7), we have

$$\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,1}^1{\phi }^1\le -C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,0}^1{\phi }^0+{\mathit{\xi }}^1.$$

(10)

Noting that ${\tilde{d}}_{{\alpha }_i,0}^n=-d_{{\alpha }_i,n}^n,{\tilde{d}}_{{\alpha }_i,n}^n=d_{{\alpha }_i,1}^n$, we have

$$\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,1}^1{\phi }^1\le C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,1}^1{\phi }^0+{\mathit{\xi }}^1.$$

(11)

Then, we can obtain

$${\phi }^1\le C_0({\phi }^0+\frac{1}{{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,1}^1}{\mathit{\xi }}^1).$$

(12)

It means that the conclusion (8) is valid for the case of $n=1$. Assume that (8) is valid for $n=1,2,\cdots ,r$. We now need to prove that it also holds for $n=r+1$. Selecting $n=r+1$ in (7), we obtain

$$\begin{array}{cc}\hfill & \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,j+1}^{j+1}{\phi }^{j+1}\hfill \\ \hfill \le & -C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^j{\tilde{d}}_{{\alpha }_i,k}^{j+1}{\phi }^k+{\mathit{\xi }}^{j+1}\hfill \\ \hfill =& C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=1}^j(d_{{\alpha }_i,j-k+1}^{j+1}-d_{{\alpha }_i,j-k+2}^{j+1}){\phi }^k+C_0\sum _{i=1}^m{b}_i{d}_{{\alpha }_i,j+1}^{j+1}{\phi }^0+{\mathit{\xi }}^{j+1}\hfill \\ \hfill =& C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{j-1}(d_{{\alpha }_i,k+1}^{j+1}-d_{{\alpha }_i,k+2}^{j+1}){\phi }^{j-k}+C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,j+1}^{j+1}{\phi }^0+{\mathit{\xi }}^{j+1}.\hfill \end{array}$$

(13)

Noting that $0<d_{{\alpha }_i,k+1}^{k+1}<d_{{\alpha }_i,k}^k$ and $0<d_{{\alpha }_i,k+1}^n<d_{{\alpha }_i,k}^n,(k=0,1\cdots j)$, we have

$$\begin{array}{cc}\hfill & \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,j+1}^{j+1}{\phi }^{j+1}\hfill \\ \hfill \le & C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{j-1}(d_{{\alpha }_i,k+1}^{j+1}-d_{{\alpha }_i,k+2}^{j+1})[C_0^{j-k}({\phi }^0+\frac{1}{{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,j-k}^{j-k}}{\mathit{\xi }}^{j-k})]\hfill \\ \hfill & +C_0\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\mathrm{d}}_{{\alpha }_i,j+1}^{j+1}({\phi }^0+\frac{1}{{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,j+1}^{j+1}}{\mathit{\xi }}^{j+1})\hfill \\ \hfill \le & C_0^{j+1}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}[\sum _{k=0}^{j-1}(d_{{\alpha }_i,k+1}^{j+1}-d_{{\alpha }_i,k+2}^{j+1})+d_{{\alpha }_i,j+1}^{j+1}]({\phi }^0+\frac{1}{{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,j+1}^{j+1}}{\mathit{\xi }}^{j+1})\hfill \\ \hfill \le & C_0^{j+1}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,1}^{j+1}({\phi }^0+\frac{1}{{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{d}_{{\alpha }_i,j+1}^{j+1}}{\mathit{\xi }}^{j+1}).\hfill \end{array}$$

(14)

Therefore, using the mathematical induction method, we can complete the proof of (8). From [37], we know $n^{{\alpha }_i}(n^{1-{\alpha }_i}-{(n-1)}^{1-{\alpha }_i})\ge 1-{\alpha }_i$, and then

$$\begin{array}{cc}\hfill d_{{\alpha }_i,n}^n& =\frac{{\left(n\tau \right)}^{1-{\alpha }_i}-{(n\tau -\tau )}^{1-{\alpha }_i}}{\tau }=\frac{(n^{1-{\alpha }_i}-{(n-1)}^{1-{\alpha }_i})}{{\tau }^{{\alpha }_i}}\ge \frac{(1-{\alpha }_i)}{{\tau }^{{\alpha }_i}{n}^{{\alpha }_i}}.\hfill \end{array}$$

(15)

Thus, making use of (8) and (15), we can complete the proof of (9). □

Next, we give the existence and uniqueness results for the MFE scheme (5).

Theorem 1. 

The MFE scheme (5) has a unique solution.

Proof. 

Let $V_h=\mathrm{span}\{{\varphi }_1,{\varphi }_2,\cdots ,{\varphi }_{M_1}\}$ and ${\mathit{W}}_h=\mathrm{span}\{{\mathit{\psi }}_1,{\mathit{\psi }}_2,\cdots ,{\mathit{\psi }}_{M_2}\}$. Then, $u_h^n$ and ${\mathit{\lambda }}_h^n$ can be written as

$$\begin{array}{c}\hfill u_h^n=\sum _{i=1}^{M_1}{\tilde{u}}_i^n{\varphi }_i,{\mathit{\lambda }}_h^n=\sum _{j=1}^{M_2}{\tilde{s}}_j^n{\mathit{\psi }}_j.\end{array}$$

(16)

Substituting (16) into (5) and selecting $v_h={\varphi }_i$   $(i=1,2,\cdots ,M_1)$ and ${\mathit{w}}_h={\mathit{\psi }}_j$   $(j=1,2,\cdots ,M_2)$, we have

$$\left[\begin{array}{cc}{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n{B}_1+B_3& D^T\\ -D& B_2\end{array}\right]\left[\begin{array}{c}{U}^n\\ L^n\end{array}\right]=\left[\begin{array}{c}{F}^n-{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{\displaystyle \sum _{k=0}^{n-1}}{\tilde{d}}_{{\alpha }_i,k}^n{B}_1{U}^k\\ 0\end{array}\right],$$

(17)

where

$$\begin{array}{cc}{U}^n={({\tilde{u}}_1^n,{\tilde{u}}_2^n,\cdots ,{\tilde{u}}_{M_1}^n)}^T,\hfill & L^n={({\tilde{s}}_1^n,{\tilde{s}}_2^n,\cdots ,{\tilde{s}}_{M_2}^n)}^T,\hfill \\ B_1={\left(({\varphi }_i,{\varphi }_j)\right)}_{M_1\times M_1},\hfill & B_2={\left(({\mathcal{A}}^{-1}{\mathit{\psi }}_i,{\mathit{\psi }}_j)\right)}_{M_2\times M_2},\hfill \\ B_3={\left((p{\varphi }_i,{\varphi }_j)\right)}_{M_1\times M_1},\hfill & D={\left((\mathrm{div}{\mathit{\psi }}_i,{\varphi }_j)\right)}_{M_2\times M_1},\hfill \\ F^n={\left((f^n,{\varphi }_i)\right)}_{M_1\times 1},\hfill \end{array}$$

Noting that $B_1$ and $B_2$ are symmetric positive definite matrices and $B_3$ is a symmetric semi-positive matrix, we have

$$\left[\begin{array}{cc}E& -D^T{B}_2^{-1}\\ 0& E\end{array}\right]\left[\begin{array}{cc}{\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n{B}_1+B_3& D^T\\ -D& B_2\end{array}\right]=\left[\begin{array}{cc}G& 0\\ -D& B_2\end{array}\right].$$

(18)

where $G={\displaystyle \sum _{i=1}^m}\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n{B}_1+B_3+D^T{B}_2^{-1}D$. It is easy to see that G is invertible, so the coefficient matrix of linear Equation (17) is invertible. This means that the MFE scheme (5) has a unique solution. □

Remark 2. 

For Lemma 3, when $C_0=1$, a similar conclusion can be seen from the proof of Theorem 3.1 in [20]. When $C_0>1$, some special applications can be seen from [39]. It should be noted that this lemma can be considered a fractional Grönwall inequality without any other conditions for its existence, which will play a crucial role in the subsequent proof process of stability and convergence analyses.

4. Stability Analysis

In this section, we will discuss the unconditional stability for the MFE scheme (5)–(6).

Theorem 2. 

Let ${(u_h^n,{\mathit{\lambda }}_h^n)}_{n=1}^N$ be the solutions of the MFE scheme (5). Then, there exists a constant $C>0$ independent of h and N such that

$$\begin{array}{cc}\hfill & \parallel u_h^n\parallel \le \parallel u_h^0\parallel +\sum _{i=1}^m\frac{\Gamma (1-{\alpha }_i)t_n^{{\alpha }_i}}{b_i}\underset{t\in [0,T]}{sup}\parallel f\left(t\right)\parallel \triangleq U_h^{\diamond },\hfill \\ \hfill & \parallel {\mathit{\lambda }}_h^n\parallel \le C\left(\parallel {\mathit{\lambda }}_h^0\parallel +{\left(\sum _{i=1}^m\frac{\Gamma (1-{\alpha }_i)t_n^{{\alpha }_i}}{b_i}\right)}^{1/2}(\underset{t\in [0,T]}{sup}\parallel f\left(t\right)\parallel +{\parallel p\parallel }_{\infty }{U}_h^{\diamond })\right).\hfill \end{array}$$

Proof. 

Taking $v_h=u_h^n$ and $w_h={\mathit{\lambda }}_h^n$ in (5), we have

$$(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n,u_h^n)+({\mathcal{A}}^{-1}{\mathit{\lambda }}_h^n,{\mathit{\lambda }}_h^n)+(p{u}_h^n,u_h^n)=(f^n,u_h^n).$$

(19)

Using Lemma 1 and the definition of $D_N^{{\alpha }_i}{u}_h^n$, we have

$$\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n(u_h^n,u_h^n)+{\mu }_0{\parallel {\mathit{\lambda }}_h^n\parallel }^2\le \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n(u_h^k,u_h^n)+(f^n,u_h^n).$$

(20)

Applying the Cauchy–Schwarz inequality yields

$$\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n\parallel u_h^n\parallel \le \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n\parallel u_h^k\parallel +\parallel f^n\parallel .$$

(21)

Using Lemma 3, we obtain

$$\parallel u_h^n\parallel \le \parallel u_h^0\parallel +\sum _{i=1}^m\frac{\Gamma (1-{\alpha }_i)t_n^{{\alpha }_i}}{b_i}\underset{t\in [0,T]}{sup}\parallel f\left(t\right)\parallel \triangleq U_h^{\diamond }.$$

(22)

Next, using (5) $\left(b\right)$ and (6) $\left(a\right)$, we have

$$(A^{-1}\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\mathit{\lambda }}_h^n,{\mathit{w}}_h)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n,\mathrm{div}{\mathit{w}}_h)=0,\forall {\mathit{w}}_h\in {\mathit{W}}_h.$$

(23)

Choosing $v_h={\displaystyle \sum _{i=1}^m}{b}_i{D}_N^{{\alpha }_i}{u}_h^n$ and ${\mathit{w}}_h={\mathit{\lambda }}_h^n$ in (5) $\left(a\right)$ and (23), respectively, we obtain

$${\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n\parallel }^2+(A^{-1}\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\mathit{\lambda }}_h^n,{\mathit{\lambda }}_h^n)+(p{u}_h^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n)=(f^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n).$$

(24)

Using Lemma 2 in (24) yields

$$\begin{array}{cc}\hfill & {\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n\parallel }^2+\frac{1}{2}\sum _{i=0}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}[{\tilde{d}}_{{\alpha }_i,n}^n(A^{-1}{\mathit{\lambda }}_h^n,{\mathit{\lambda }}_h^n)+\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n(A^{-1}{\mathit{\lambda }}_h^k,{\mathit{\lambda }}_h^k)\hfill \\ \hfill & -\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n(A^{-1}({\mathit{\lambda }}_h^n-{\mathit{\lambda }}_h^k),{\mathit{\lambda }}_h^n-{\mathit{\lambda }}_h^k)]=(f^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n)-(p{u}_h^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n).\hfill \end{array}$$

(25)

Because of ${\tilde{d}}_{{\alpha }_i,k}^n<0,0<k\le n-1$, we have

$$\begin{array}{cc}\hfill & {\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n\parallel }^2+\frac{1}{2}\sum _{i=0}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n(A^{-1}{\mathit{\lambda }}_h^n,{\mathit{\lambda }}_h^n)\hfill \\ \hfill \le & -\frac{1}{2}\sum _{i=0}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n(A^{-1}{\mathit{\lambda }}_h^k,{\mathit{\lambda }}_h^k)\hfill \\ \hfill & +(f^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n)-(p{u}_h^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n).\hfill \end{array}$$

(26)

Apply the Cauchy–Schwarz inequality and the Young inequality in (26) to obtain

$$\begin{array}{cc}\hfill & {\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n\parallel }^2+\frac{1}{2}\sum _{i=0}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n(A^{-1}{\mathit{\lambda }}_h^n,{\mathit{\lambda }}_h^n)\hfill \\ \hfill \le & -\frac{1}{2}\sum _{i=0}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n(A^{-1}{\mathit{\lambda }}_h^k,{\mathit{\lambda }}_h^k)\hfill \\ \hfill & +\frac{1}{2}{\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{u}_h^n\parallel }^2+{\parallel f^n\parallel }^2+{\parallel p\parallel }_{\infty }^2{\parallel u_h^n\parallel }^2.\hfill \end{array}$$

(27)

Using Lemma 3 in (27), we obtain

$$\begin{array}{c}\hfill \parallel {\mathit{\lambda }}_h^n\parallel \le C(\parallel {\mathit{\lambda }}_h^0\parallel +{\left(\sum _{i=1}^m\frac{\Gamma (1-{\alpha }_i)t_n^{{\alpha }_i}}{b_i}\right)}^{1/2}(\underset{t\in [0,T]}{sup}\parallel f\left(t\right)\parallel +{\parallel p\parallel }_{\infty }{U}_h^{\diamond })).\end{array}$$

(28)

Thus, we complete the proof. □

5. Convergence Analysis

In this section, we will present the convergence results. For this purpose, we first introduce the mixed elliptic projection $(R_{h}u,R_h\mathit{\lambda })\in V_h\times {\mathit{W}}_h$ defined by

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\left(a\right)({\mathcal{A}}^{-1}(\mathit{\lambda }-R_h\mathit{\lambda }),{\mathit{w}}_h)-(u-R_{h}u,\mathrm{div}{\mathit{w}}_h)=0,\hfill & \forall w_h\in {\mathit{W}}_h,\hfill \\ \left(b\right)(\mathrm{div}(\mathit{\lambda }-R_h\mathit{\lambda }),v_h)=0,\hfill & \forall v_h\in V_h.\hfill \end{array}\right.\end{array}$$

(29)

Then, the above projection satisfies the classical estimates as follows.

Lemma 4

([40,41]). There exists a constant $C>0$ independent of h and N such that

$$\begin{array}{cc}\hfill & \parallel \mathit{\lambda }-R_h\mathit{\lambda }\parallel \le C{h}^{r+1}{\parallel \mathit{\lambda }\parallel }_{r+1},\mathit{for}\mathit{\lambda }\in {\left(H^{r+1}(\Omega )\right)}^2,\hfill \\ \hfill & \parallel \mathrm{div}(\mathit{\lambda }-R_h\mathit{\lambda })\parallel \le C{h}^{r+1}{\parallel \mathrm{div}\mathit{\lambda }\parallel }_{r+1},\mathit{for}\mathrm{div}\mathit{\lambda }\in H^{r+1}(\Omega ),\hfill \\ \hfill & \parallel u-R_{h}u\parallel \le C{h}^{r+1}({\parallel u\parallel }_{r+1}+{\parallel \mathit{\lambda }\parallel }_{r+1}),\mathit{for}u\in H^{r+1}(\Omega ),\mathit{\lambda }\in {\left(H^{r+1}(\Omega )\right)}^2.\hfill \end{array}$$

For the truncation error $Q_{{\alpha }_i}^n$   $(i=1,2,\cdots ,m)$ of the $L1$ formula, from [36,37], we give the following estimates.

Lemma 5. 

Let $u\in {\mathcal{C}}^2(\overline{J},L^2(\Omega ))$. Then, we have

$$\begin{array}{cc}\hfill & \parallel Q_{{\alpha }_i}^n\parallel \le C{N}^{-(2-{\alpha }_i)},i=1,2,\cdots ,m,\hfill \\ \hfill & \parallel \sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n\parallel \le C{N}^{-(2-{\alpha }_1)},\hfill \end{array}$$

where $C>0$ is a constant independent of h and N.

Now, we write the errors $u\left(t_n\right)-u_h^n=u\left(t_n\right)-R_{h}u\left(t_n\right)+R_{h}u\left(t_n\right)-u_h^n={\rho }^n+{\theta }^n$ and $\mathit{\lambda }\left(t_n\right)-{\mathit{\lambda }}_h^n=\mathit{\lambda }\left(t_n\right)-R_h\mathit{\lambda }\left(t_n\right)+R_h\mathit{\lambda }\left(t_n\right)-{\mathit{\lambda }}_h^n={\mathit{\xi }}^n+{\mathit{\eta }}^n$. From (3) and (5), making use of the mixed elliptic projection $R_h$, we have the following error equations:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\left(a\right)({\displaystyle \sum _{i=1}^m}{b}_i{D}_N^{{\alpha }_i}({\theta }^n+{\rho }^n),v_h)+(\mathrm{div}{\mathit{\eta }}^n,v_h)+(p({\theta }^n+{\rho }^n),v_h)\hfill \\ =-({\displaystyle \sum _{i=1}^m}{b}_i{Q}_{{\alpha }_i}^n,v_h),\hfill & \forall v_h\in V_h,\hfill \\ \left(b\right)({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{w}}_h)-({\theta }^n,\mathrm{div}{\mathit{w}}_h)=0,\hfill & \forall {\mathit{w}}_h\in {\mathit{W}}_h.\hfill \end{array}\right.\end{array}$$

(30)

Noting that $(u_h^0,{\mathit{\lambda }}_h^0)=(R_h{u}_0,R_h{\mathit{\lambda }}_0)$, we have ${\theta }^0=0$ and ${\mathit{\eta }}^0=\mathbf{0}$. We next give the convergence results for the MFE scheme (5)–(6).

Theorem 3. 

Let $(u^n,{\mathit{\lambda }}^n)\in V\times \mathit{W}$ and $(u_h^n,{\mathit{\lambda }}_h^n)\in V_h\times {\mathit{W}}_h$ be the solutions of (3) and (5), respectively. Assume that $u,\mathrm{div}\mathit{\lambda }\in {\mathcal{C}}^2(\overline{J},H^{r+1}(\Omega ))$, $\mathit{\lambda }\in {\mathcal{C}}^2(\overline{J},{\left(H^{r+1}(\Omega )\right)}^2)$. Then, we have

$$\begin{array}{cc}\hfill & \underset{1\le n\le N}{max}\parallel u\left(t_n\right)-u_h^n\parallel +\underset{1\le n\le N}{max}\parallel \mathit{\lambda }\left(t_n\right)-{\mathit{\lambda }}_h^n\parallel \le C(h^{r+1}+N^{-(2-{\alpha }_1)}),\hfill \\ \hfill & \underset{1\le n\le N}{max}{\parallel \mathit{\lambda }\left(t_n\right)-{\mathit{\lambda }}_h^n\parallel }_{\mathit{H}(\mathrm{div},\Omega )}\le C(1+N^{\frac{{\alpha }_m}{2}})(h^{r+1}+N^{-(2-{\alpha }_1)}),\hfill \end{array}$$

where $C>0$ is a constant independent of h and N.

Proof. 

Taking $v_h={\theta }^n$ and ${\mathit{w}}_h={\mathit{\eta }}^n$ in (30), we can obtain

$$\begin{array}{cc}\hfill & (\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n,{\theta }^n)+({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{\eta }}^n)+(p{\theta }^n,{\theta }^n)\hfill \\ \hfill =& -(\sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n,{\theta }^n)-(p{\rho }^n,{\theta }^n)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n,{\theta }^n).\hfill \end{array}$$

(31)

Noting that $p\left(\mathit{x}\right)\ge 0$, using the Lemma 1 and the definition of $D_N^{{\alpha }_i}{u}_h^n$, we have

$$\begin{array}{cc}\hfill & \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n({\theta }^n,{\theta }^n)+{\mu }_0{\parallel {\mathit{\eta }}^n\parallel }^2\hfill \\ \hfill \le & -\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\theta }^k,{\theta }^n)-(\sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n,{\theta }^n)-(p{\rho }^n,{\theta }^n)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n,{\theta }^n).\hfill \end{array}$$

(32)

Applying the Cauchy–Schwarz inequality, we obtain

$$\begin{array}{cc}\hfill & \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n{\parallel {\theta }^n\parallel }^2+{\mu }_0{\parallel {\mathit{\eta }}^n\parallel }^2\hfill \\ \hfill \le & -\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n\parallel {\theta }^k\parallel \parallel {\theta }^n\parallel \hfill \\ \hfill & +(\parallel \sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n\parallel +{\parallel p\parallel }_{\infty }\parallel \rho \parallel +\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n\parallel )\parallel {\theta }^n\parallel ,\hfill \end{array}$$

(33)

and then

$$\begin{array}{cc}\hfill & \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n{\parallel {\theta }^n\parallel }^2\hfill \\ \hfill \le & -\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n\parallel {\theta }^k\parallel +(\parallel \sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n\parallel +{\parallel p\parallel }_{\infty }\parallel \rho \parallel +\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n\parallel ).\hfill \end{array}$$

(34)

Using Lemmas 4 and 5, we obtain

$$\begin{array}{cc}\hfill \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n\parallel {\theta }^n\parallel & \le -\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n\parallel {\theta }^k\parallel +C(N^{-(2-{\alpha }_1)}+h^{r+1}).\hfill \end{array}$$

(35)

Noting that ${\theta }^0=0$ and using Lemma 3, we obtain

$$\parallel {\theta }^n\parallel \le C(h^{r+1}+N^{-(2-{\alpha }_1)}).$$

(36)

Now, from (30) $\left(b\right)$, we obtain

$$({\mathcal{A}}^{-1}\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\mathit{\eta }}^n,{\mathit{w}}_h)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n,\mathrm{div}{\mathit{w}}_h)=0,\forall {\mathit{w}}_h\in {\mathit{W}}_h.$$

(37)

Choosing $v_h={\displaystyle \sum _{i=1}^m}{b}_i{D}_N^{{\alpha }_i}{\theta }^n$ and ${\mathit{w}}_h={\mathit{\eta }}^n$ in (30) $\left(a\right)$ and (37), respectively, we can obtain

$$\begin{array}{cc}\hfill & {\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n\parallel }^2+({\mathcal{A}}^{-1}\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\mathit{\eta }}^n,{\mathit{\eta }}^n)\hfill \\ \hfill =& -(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n)-(p({\rho }^n+{\theta }^n),\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n)-(\sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n).\hfill \end{array}$$

(38)

Using Lemma 2, we have

$$\begin{array}{cc}\hfill & \parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n{\parallel }^2+\frac{1}{2}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}[{\tilde{d}}_{{\alpha }_i,n}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{\eta }}^n)-\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n\left({\mathcal{A}}^{-1}({\mathit{\eta }}^n-{\mathit{\eta }}^k,{\mathit{\eta }}^n-{\mathit{\eta }}^k)\right]\hfill \\ \hfill =& -\frac{1}{2}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^k,{\mathit{\eta }}^k)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n)\hfill \\ \hfill & -(\sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n)-(p({\rho }^n+{\theta }^n),\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n).\hfill \end{array}$$

(39)

Noting that ${\tilde{d}}_{{\alpha }_i,k}^n<0,0<k\le n-1$ and using Lemma 1, we obtain

$$\begin{array}{cc}\hfill & {\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n\parallel }^2+\frac{1}{2}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{\eta }}^n)\hfill \\ \hfill \le & -\frac{1}{2}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^k,{\mathit{\eta }}^k)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n)\hfill \\ \hfill & -(\sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n,\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n)-(p({\rho }^n+{\theta }^n),\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n).\hfill \end{array}$$

(40)

Applying the Cauchy–Schwarz and the Young inequality in (40) yields

$$\begin{array}{cc}\hfill & {\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\theta }^n\parallel }^2+\frac{1}{2}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{\eta }}^n)\hfill \\ \hfill \le & -\frac{1}{2}\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^k,{\mathit{\eta }}^k)+\frac{1}{2}{\parallel \sum _{i=1}^m{b}_i{D}_N^{\alpha }{\theta }^n\parallel }^2\hfill \\ \hfill & +2(\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n\parallel +{\parallel \sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n\parallel }^2+{\parallel p\parallel }_{\infty }^2({\parallel {\rho }^n\parallel }^2+{\parallel {\theta }^n\parallel }^2)).\hfill \end{array}$$

(41)

Using Lemmas 4 and 5, we obtain

$$\begin{array}{cc}\hfill & \sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{\tilde{d}}_{{\alpha }_i,n}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{\eta }}^n)\hfill \\ \hfill \le & -\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^k,{\mathit{\eta }}^k)+C(N^{-2(2-{\alpha }_1)}+h^{2r+2}).\hfill \end{array}$$

(42)

Noting that ${\mathit{\eta }}^0=\mathbf{0}$ and using Lemma 3, we obtain

$$\parallel {\mathit{\eta }}^n\parallel \le C(h^{r+1}+N^{-(2-{\alpha }_1)}).$$

(43)

We now estimate $\parallel {\mathit{\lambda }}^n-{\mathit{\lambda }}_h^n{\parallel }_{\mathit{H}(\mathrm{div},\Omega )}$. Taking $v_h={\displaystyle \sum _{i=1}^m}{b}_i{D}_N^{{\alpha }_i}{\theta }^n$ and ${\mathit{w}}_h={\mathit{\eta }}^n$ in (30) $\left(a\right)$ and (37), respectively, we have

$$\begin{array}{cc}\hfill \parallel \mathrm{div}{\mathit{\eta }}^n{\parallel }^2=& -({\mathcal{A}}^{-1}\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\mathit{\eta }}^n,{\mathit{\eta }}^n)-(\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n,\mathrm{div}{\mathit{\eta }}^n)\hfill \\ \hfill & -(p({\rho }^n+{\theta }^n),\mathrm{div}{\mathit{\eta }}^n)-(\sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n,\mathrm{div}{\mathit{\eta }}^n).\hfill \end{array}$$

(44)

For the term $-({\mathcal{A}}^{-1}{\displaystyle \sum _{i=1}^m}{b}_i{D}_N^{{\alpha }_i}{\mathit{\eta }}^n,{\mathit{\eta }}^n)$, noting that ${\displaystyle \sum _{k=0}^{n-1}}(-{\tilde{d}}_{{\alpha }_i,k}^n)=T^{-{\alpha }_i}{N}^{{\alpha }_i}$, we obtain

$$\begin{array}{cc}\hfill -({\mathcal{A}}^{-1}\sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\mathit{\eta }}^n,{\mathit{\eta }}^n)=& -\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\left(\sum _{k=0}^{n-1}{\tilde{d}}_{{\alpha }_i,k}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^k,{\mathit{\eta }}^n)+{\tilde{d}}_{{\alpha }_i,n}^n({\mathcal{A}}^{-1}{\mathit{\eta }}^n,{\mathit{\eta }}^n)\right)\hfill \\ \hfill \le & {\mu }_1\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}\sum _{k=0}^{n-1}(-{\tilde{d}}_{{\alpha }_i,k}^n)\parallel {\mathit{\eta }}^k\parallel \parallel {\mathit{\eta }}^n\parallel \hfill \\ \hfill \le & C\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{T}^{-{\alpha }_i}{N}^{{\alpha }_i}{(N^{-(2-{\alpha }_1)}+h^{r+1})}^2.\hfill \end{array}$$

(45)

Then, it holds from (45) that

$$\begin{array}{cc}\hfill \parallel \mathrm{div}{\mathit{\eta }}^n{\parallel }^2=& 2(\parallel \sum _{i=1}^m{b}_i{D}_N^{{\alpha }_i}{\rho }^n\parallel +{\parallel \sum _{i=1}^m{b}_i{Q}_{{\alpha }_i}^n\parallel }^2+{\parallel p\parallel }_{\infty }^2({\parallel {\rho }^n\parallel }^2+{\parallel {\theta }^n\parallel }^2))\hfill \\ \hfill & +C\sum _{i=1}^m\frac{b_i}{\Gamma (2-{\alpha }_i)}{T}^{-{\alpha }_i}{N}^{{\alpha }_i}{(N^{-(2-{\alpha }_1)}+h^{r+1})}^2+\frac{1}{2}{\parallel \mathrm{div}{\mathit{\eta }}^n\parallel }^2.\hfill \end{array}$$

(46)

Using Lemmas 4 and 5, we have

$$\parallel \mathrm{div}{\mathit{\eta }}^n\parallel \le C(1+N^{\frac{{\alpha }_m}{2}})(h^{r+1}+N^{-(2-{\alpha }_1)}).$$

(47)

Then, we finish the proof. □

Remark 3.

(I) For variables u and λ, we define the discrete norms of the errors as follows:

$$\begin{array}{cc}\hfill & \parallel u-u_h{\parallel }_{{\widehat{L}}^{\infty }\left(L^2(\Omega )\right)}=\underset{1\le n\le N}{max}\parallel u\left(t_n\right)-u_h^n\parallel ,\hfill \\ \hfill & \parallel \mathit{\lambda }-{\mathit{\lambda }}_h{\parallel }_{{\widehat{L}}^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)}=\underset{1\le n\le N}{max}\parallel \mathit{\lambda }\left(t_n\right)-{\mathit{\lambda }}_h^n\parallel ,\hfill \\ \hfill & \parallel \mathit{\lambda }-{\mathit{\lambda }}_h{\parallel }_{{\widehat{L}}^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)}=\underset{1\le n\le N}{max}{\parallel \mathit{\lambda }\left(t_n\right)-{\mathit{\lambda }}_h^n\parallel }_{\mathit{H}(\mathrm{div},\Omega )}.\hfill \end{array}$$

From Theorem 3, we obtain the optimal a priori error estimate results for u in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm and λ in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ norm and obtain the suboptimal error estimate for λ in the discrete $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norm. In the actual calculation in the next section, we achieve the optimal convergence rates for variables u and λ based on the above discrete norms.

(II) It should be pointed out that the solutions of many FPDEs have an initial layer at $t=0$ (see [42,43]). To overcome this difficulty, some scholars have adopted nonuniform mesh methods and achieved excellent results [24,26,42,44,45,46]. Moreover, it is noted that $\{{\xi }^k$ : $0\le k\le N\}$ in Lemma 3 is required to be a nondecreasing positive sequence, so the error estimates for the MFE scheme (5)–(6) with the temporal nonuniform method should adopt some other techniques. It is gratifying that the numerical results in Example 3 show that the MFE scheme (5)–(6) with the temporal graded mesh is feasible and effective.

6. Numerical Examples

In this section, we given three test examples to verify the effectiveness and convergence accuracy of the proposed MFE scheme (5)–(6) and adopt the lowest-order Raviart–Thomas MFE space for variables u and $\mathit{\lambda }$ in the numerical experiments.

Example 1. 

Consider the following two-term TFRD equation:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{D}_t^{{\alpha }_1}u(\mathit{x},t)+D_t^{{\alpha }_2}u(\mathit{x},t)-\Delta u(\mathit{x},t)+p\left(\mathit{x}\right)u(\mathit{x},t)=f(\mathit{x},t),\hfill & (\mathit{x},t)\in \Omega \times J,\hfill \\ u(\mathit{x},t)=0,\hfill & (\mathit{x},t)\in \partial \Omega \times \overline{J},\hfill \\ u(\mathit{x},0)=u_0\left(\mathit{x}\right),\hfill & \mathit{x}\in \overline{\Omega },\hfill \end{array}\right.\end{array}$$

(48)

where $J=(0,1]$, $\Omega ={(0,1)}^2$, $p\left(\mathit{x}\right)=1+x_1^2+x_2^2$, $\mathit{x}=(x_1,x_2)\in \Omega$, $u(\mathit{x},0)=0$, and the source function f is taken by

$$\begin{array}{cc}\hfill f(\mathit{x},t)=& \left(\frac{\Gamma (3+{\alpha }_1+{\alpha }_2)}{\Gamma (3+{\alpha }_2)}{t}^{2+{\alpha }_2}+\frac{\Gamma (3+{\alpha }_1+{\alpha }_2)}{\Gamma (3+{\alpha }_1)}{t}^{2+{\alpha }_1}+(2{\pi }^2+p\left(\mathit{x}\right))t^{2+{\alpha }_1+{\alpha }_2}\right)\hfill \\ \hfill & \times sin\left(\pi x_1\right)sin\left(\pi x_2\right).\hfill \end{array}$$

And we can find the analytical solutions for variables u and λ as follows:

$$\begin{array}{cc}\hfill & u(\mathit{x},t)=t^{2+{\alpha }_1+{\alpha }_2}sin\left(\pi x_1\right)sin\left(\pi x_2\right),\hfill \\ \hfill & \mathit{\lambda }(\mathit{x},t)=-\pi t^{2+{\alpha }_1+{\alpha }_2}(cos\left(\pi x_1\right)sin\left(\pi x_2\right),sin\left(\pi x_1\right)cos\left(\pi x_2\right)).\hfill \end{array}$$

In the numerical simulation, we select fractional parameters ${\alpha }_1=0.9$, $0.7$, $0.5$ and ${\alpha }_2$ = 0.1, 0.4 in Equation (48) and know that among these different fractional parameters, the convergence rates are only related to the largest fractional parameter ${\alpha }_1$ from Theorem 3. By taking $N=5$, 8, 10, 16 and the corresponding $h=\sqrt{2}/N^{2-{\alpha }_1}$, we give the error results and convergence rates in

Table 1. Numerical results with $h\approx \sqrt{2}/N^{2-{\alpha }_1}$ and ${\alpha }_1=0.9$ in Example 1.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	5	$8.6930\times {10}^{-2}$	−	$3.3630\times {10}^{-1}$	−	$1.7545\times {10}^{+0}$	−
	8	$5.2424\times {10}^{-2}$	1.0760	$2.0229\times {10}^{-1}$	1.0814	$1.0569\times {10}^{+0}$	1.0784
	10	$4.0389\times {10}^{-2}$	1.1687	$1.5578\times {10}^{-1}$	1.1709	$8.1390\times {10}^{-1}$	1.1708
	16	$2.3908\times {10}^{-2}$	1.1157	$9.2174\times {10}^{-2}$	1.1165	$4.8141\times {10}^{-1}$	1.1172
0.4	5	$8.7218\times {10}^{-2}$	−	$3.3776\times {10}^{-1}$	−	$1.7631\times {10}^{+0}$	−
	8	$5.2624\times {10}^{-2}$	1.0750	$2.0332\times {10}^{-1}$	1.0799	$1.0619\times {10}^{+0}$	1.0786
	10	$4.0555\times {10}^{-2}$	1.1674	$1.5663\times {10}^{-1}$	1.1692	$8.1786\times {10}^{-1}$	1.1704
	16	$2.4011\times {10}^{-2}$	1.1153	$9.2701\times {10}^{-2}$	1.1160	$4.8371\times {10}^{-1}$	1.1175

,

Table 2. Numerical results with $h\approx \sqrt{2}/N^{2-{\alpha }_1}$ and ${\alpha }_1=0.7$ in Example 1.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	5	$6.5262\times {10}^{-2}$	−	$2.5184\times {10}^{-1}$	−	$1.3146\times {10}^{+0}$	−
	8	$3.4916\times {10}^{-2}$	1.3308	$1.3451\times {10}^{-1}$	1.3344	$7.0291\times {10}^{-1}$	1.3321
	10	$2.6204\times {10}^{-2}$	1.2863	$1.0092\times {10}^{-1}$	1.2875	$5.2740\times {10}^{-1}$	1.2873
	16	$1.4175\times {10}^{-2}$	1.3073	$5.4581\times {10}^{-2}$	1.3077	$2.8520\times {10}^{-1}$	1.3080
0.4	5	$6.5379\times {10}^{-2}$	−	$2.5244\times {10}^{-1}$	−	$1.3184\times {10}^{+0}$	−
	8	$3.4998\times {10}^{-2}$	1.3296	$1.3493\times {10}^{-1}$	1.3328	$7.0499\times {10}^{-1}$	1.3318
	10	$2.6269\times {10}^{-2}$	1.2858	$1.0125\times {10}^{-1}$	1.2869	$5.2895\times {10}^{-1}$	1.2875
	16	$1.4212\times {10}^{-2}$	1.3070	$5.4771\times {10}^{-2}$	1.3073	$2.8602\times {10}^{-1}$	1.3081

and

Table 3. Numerical results with $h\approx \sqrt{2}/N^{2-{\alpha }_1}$ and ${\alpha }_1=0.5$ in Example 1.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	5	$4.7518\times {10}^{-2}$	−	$1.8311\times {10}^{-1}$	−	$9.5625\times {10}^{-1}$	−
	8	$2.2765\times {10}^{-2}$	1.5657	$8.7633\times {10}^{-2}$	1.5679	$4.5800\times {10}^{-1}$	1.5663
	10	$1.6367\times {10}^{-2}$	1.4786	$6.2996\times {10}^{-2}$	1.4792	$3.2925\times {10}^{-1}$	1.4791
	16	$8.1859\times {10}^{-3}$	1.4742	$3.1504\times {10}^{-2}$	1.4743	$1.6465\times {10}^{-1}$	1.4744
0.4	5	$4.7568\times {10}^{-2}$	−	$1.8337\times {10}^{-1}$	−	$9.5799\times {10}^{-1}$	−
	8	$2.2802\times {10}^{-2}$	1.5645	$8.7824\times {10}^{-2}$	1.5662	$4.5893\times {10}^{-1}$	1.5658
	10	$1.6396\times {10}^{-2}$	1.4781	$6.3145\times {10}^{-2}$	1.4785	$3.2993\times {10}^{-1}$	1.4790
	16	$8.2016\times {10}^{-3}$	1.4738	$3.1585\times {10}^{-2}$	1.4739	$1.6499\times {10}^{-1}$	1.4744

for the MFE scheme (5)–(6), which show that the convergence rates in the temporal direction for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms) are close to $2-{\alpha }_1$. Moreover, in order to test convergence rates in the spatial direction, by fixing $N=100$ and taking $h=\sqrt{2}/4$, $\sqrt{2}/8$, $\sqrt{2}/16$, $\sqrt{2}/32$, we give the error results and convergence rates in

Table 4. Numerical results with $\tau =T/N=1/100$ and ${\alpha }_1=0.9$ in Example 1.

${\mathit{\alpha }}_2$	h	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	$\sqrt{2}/4$	$1.2927\times {10}^{-1}$	−	$5.0254\times {10}^{-1}$	−	$2.5970\times {10}^{+0}$	−
	$\sqrt{2}/8$	$6.5234\times {10}^{-2}$	0.9867	$2.5171\times {10}^{-1}$	0.9975	$1.3118\times {10}^{+0}$	0.9853
	$\sqrt{2}/16$	$3.2696\times {10}^{-2}$	0.9965	$1.2589\times {10}^{-1}$	0.9996	$6.5762\times {10}^{-1}$	0.9963
	$\sqrt{2}/32$	$1.6361\times {10}^{-2}$	0.9989	$6.2964\times {10}^{-2}$	0.9996	$3.2908\times {10}^{-1}$	0.9988
0.4	$\sqrt{2}/4$	$1.2926\times {10}^{-1}$	−	$5.0244\times {10}^{-1}$	−	$2.5971\times {10}^{+0}$	−
	$\sqrt{2}/8$	$6.5231\times {10}^{-2}$	0.9866	$2.5169\times {10}^{-1}$	0.9973	$1.3119\times {10}^{+0}$	0.9853
	$\sqrt{2}/16$	$3.2696\times {10}^{-2}$	0.9964	$1.2589\times {10}^{-1}$	0.9995	$6.5766\times {10}^{-1}$	0.9962
	$\sqrt{2}/32$	$1.6363\times {10}^{-2}$	0.9987	$6.2973\times {10}^{-2}$	0.9994	$3.2913\times {10}^{-1}$	0.9987

,

Table 5. Numerical results with $\tau =T/N=1/100$ and ${\alpha }_1=0.7$ in Example 1.

${\mathit{\alpha }}_2$	h	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	$\sqrt{2}/4$	$1.2929\times {10}^{-1}$	−	$5.0266\times {10}^{-1}$	−	$2.5969\times {10}^{+0}$	−
	$\sqrt{2}/8$	$6.5241\times {10}^{-2}$	0.9867	$2.5174\times {10}^{-1}$	0.9976	$1.3118\times {10}^{+0}$	0.9853
	$\sqrt{2}/16$	$3.2698\times {10}^{-2}$	0.9966	$1.2590\times {10}^{-1}$	0.9997	$6.5757\times {10}^{-1}$	0.9963
	$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2954\times {10}^{-2}$	0.9999	$3.2900\times {10}^{-1}$	0.9991
0.4	$\sqrt{2}/4$	$1.2928\times {10}^{-1}$	−	$5.0258\times {10}^{-1}$	−	$2.5970\times {10}^{+0}$	−
	$\sqrt{2}/8$	$6.5239\times {10}^{-2}$	0.9867	$2.5173\times {10}^{-1}$	0.9975	$1.3118\times {10}^{+0}$	0.9853
	$\sqrt{2}/16$	$3.2697\times {10}^{-2}$	0.9966	$1.2590\times {10}^{-1}$	0.9996	$6.5758\times {10}^{-1}$	0.9963
	$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2954\times {10}^{-2}$	0.9999	$3.2900\times {10}^{-1}$	0.9990

and

Table 6. Numerical results with $\tau =T/N=1/100$ and ${\alpha }_1=0.5$ in Example 1.

${\mathit{\alpha }}_2$	h	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	$\sqrt{2}/4$	$1.2930\times {10}^{-1}$	−	$5.0273\times {10}^{-1}$	−	$2.5969\times {10}^{+0}$	−
	$\sqrt{2}/8$	$6.5243\times {10}^{-2}$	0.9868	$2.5176\times {10}^{-1}$	0.9977	$1.3118\times {10}^{+0}$	0.9852
	$\sqrt{2}/16$	$3.2699\times {10}^{-2}$	0.9966	$1.2591\times {10}^{-1}$	0.9997	$6.5756\times {10}^{-1}$	0.9963
	$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2955\times {10}^{-2}$	0.9999	$3.2899\times {10}^{-1}$	0.9991
0.4	$\sqrt{2}/4$	$1.2929\times {10}^{-1}$	−	$5.0266\times {10}^{-1}$	−	$2.5969\times {10}^{+0}$	−
	$\sqrt{2}/8$	$6.5242\times {10}^{-2}$	0.9867	$2.5175\times {10}^{-1}$	0.9976	$1.3118\times {10}^{+0}$	0.9853
	$\sqrt{2}/16$	$3.2698\times {10}^{-2}$	0.9966	$1.2590\times {10}^{-1}$	0.9997	$6.5757\times {10}^{-1}$	0.9963
	$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2955\times {10}^{-2}$	0.9999	$3.2899\times {10}^{-1}$	0.9991

, which show that the convergence rates in the spatial direction for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms) are close to 1.

Example 2. 

Consider the following three-term TFRD equation:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{D}_t^{{\alpha }_1}u(\mathit{x},t)+D_t^{{\alpha }_2}u(\mathit{x},t)+D_t^{{\alpha }_3}u(\mathit{x},t)-\Delta u(\mathit{x},t)+p\left(\mathit{x}\right)u(\mathit{x},t)=f(\mathit{x},t),\hfill & (\mathit{x},t)\in \Omega \times J,\hfill \\ u (\mathit{x},t)=0,\hfill & (\mathit{x},t)\in \partial \Omega \times \overline{J},\hfill \\ u(\mathit{x},0)=u_0\left(\mathit{x}\right),\hfill & \mathit{x}\in \overline{\Omega },\hfill \end{array}\right.\end{array}$$

(49)

where the spatial domain Ω, temporal domain J, coefficient $p\left(\mathit{x}\right)$, and initial data $u(\mathit{x},0)$ are as in Example 1 and the source function f is taken by

$$\begin{array}{cc}\hfill f(\mathit{x},t)=& (\frac{\Gamma (3+{\alpha }_1+{\alpha }_2+{\alpha }_3)}{\Gamma (3+{\alpha }_2+{\alpha }_3)}{t}^{2+{\alpha }_2+{\alpha }_3}+\frac{\Gamma (3+{\alpha }_1+{\alpha }_2+{\alpha }_3)}{\Gamma (3+{\alpha }_1+{\alpha }_3)}{t}^{2+{\alpha }_1+{\alpha }_3}\hfill \\ \hfill & +\frac{\Gamma (3+{\alpha }_1+{\alpha }_2+{\alpha }_3)}{\Gamma (3+{\alpha }_1+{\alpha }_2)}{t}^{2+{\alpha }_1+{\alpha }_2}+(2{\pi }^2+p\left(\mathit{x}\right))t^{2+{\alpha }_1+{\alpha }_2+{\alpha }_3})sin\left(\pi x_1\right)sin\left(\pi x_2\right).\hfill \end{array}$$

And we can also find the analytical solutions for variables u and λ as follows:

$$\begin{array}{cc}\hfill & u(\mathit{x},t)=t^{2+{\alpha }_1+{\alpha }_2+{\alpha }_3}sin\left(\pi x_1\right)sin\left(\pi x_2\right),\hfill \\ \hfill & \mathit{\lambda }(\mathit{x},t)=-\pi t^{2+{\alpha }_1+{\alpha }_2+{\alpha }_3}(cos\left(\pi x_1\right)sin\left(\pi x_2\right),sin\left(\pi x_1\right)cos\left(\pi x_2\right)).\hfill \end{array}$$

In this example, since the Equation (49) contains three Caputo time-fractional derivative terms, we specifically take the fractional parameters ${\alpha }_1=0.9$, $0.7$, $0.5$ and $({\alpha }_2,{\alpha }_3)=(0.4,0.2),(0.3,0.1)$. From Theorem 3, we also point out that the convergence rates are only related to the maximum fractional parameter ${\alpha }_1$. In

Table 7. Numerical results with $h\approx \sqrt{2}/N^{2-{\alpha }_1}$ and ${\alpha }_1=0.9$ in Example 2.

${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.4	0.2	5	$8.7400\times {10}^{-2}$	−	$3.3868\times {10}^{-1}$	−	$1.7681\times {10}^{+0}$	−
		8	$5.2741\times {10}^{-2}$	1.0747	$2.0392\times {10}^{-1}$	1.0795	$1.0648\times {10}^{+0}$	1.0790
		10	$4.0650\times {10}^{-2}$	1.1670	$1.5711\times {10}^{-1}$	1.1686	$8.2004\times {10}^{-1}$	1.1704
		16	$2.4067\times {10}^{-2}$	1.1152	$9.2988\times {10}^{-2}$	1.1159	$4.8494\times {10}^{-1}$	1.1177
0.3	0.1	5	$8.7138\times {10}^{-2}$	−	$3.3735\times {10}^{-1}$	−	$1.7608\times {10}^{+0}$	−
		8	$5.2567\times {10}^{-2}$	1.0753	$2.0303\times {10}^{-1}$	1.0804	$1.0605\times {10}^{+0}$	1.0788
		10	$4.0507\times {10}^{-2}$	1.1679	$1.5638\times {10}^{-1}$	1.1698	$8.1673\times {10}^{-1}$	1.1706
		16	$2.3980\times {10}^{-2}$	1.1154	$9.2546\times {10}^{-2}$	1.1162	$4.8304\times {10}^{-1}$	1.1175

,

Table 8. Numerical results with $h\approx \sqrt{2}/N^{2-{\alpha }_1}$ and ${\alpha }_1=0.7$ in Example 2.

${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.4	0.2	5	$6.5465\times {10}^{-2}$	−	$2.5288\times {10}^{-1}$	−	$1.3208\times {10}^{+0}$	−
		8	$3.5052\times {10}^{-2}$	1.3291	$1.3521\times {10}^{-1}$	1.3321	$7.0629\times {10}^{-1}$	1.3319
		10	$2.6310\times {10}^{-2}$	1.2857	$1.0146\times {10}^{-1}$	1.2868	$5.2989\times {10}^{-1}$	1.2877
		16	$1.4234\times {10}^{-2}$	1.3070	$5.4886\times {10}^{-2}$	1.3073	$2.8651\times {10}^{-1}$	1.3083
0.3	0.1	5	$6.5343\times {10}^{-2}$	−	$2.5225\times {10}^{-1}$	−	$1.3173\times {10}^{+0}$	−
		8	$3.4972\times {10}^{-2}$	1.3300	$1.3479\times {10}^{-1}$	1.3334	$7.0434\times {10}^{-1}$	1.3321
		10	$2.6247\times {10}^{-2}$	1.2860	$1.0114\times {10}^{-1}$	1.2872	$5.2845\times {10}^{-1}$	1.2876
		16	$1.4200\times {10}^{-2}$	1.3071	$5.4707\times {10}^{-2}$	1.3075	$2.8575\times {10}^{-1}$	1.3082

and

Table 9. Numerical results with $h\approx \sqrt{2}/N^{2-{\alpha }_1}$ and ${\alpha }_1=0.5$ in Example 2.

${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.4	0.2	5	$4.7614\times {10}^{-2}$	−	$1.8360\times {10}^{-1}$	−	$9.5930\times {10}^{-1}$	−
		8	$2.2831\times {10}^{-2}$	1.5638	$8.7972\times {10}^{-2}$	1.5654	$4.5961\times {10}^{-1}$	1.5656
		10	$1.6417\times {10}^{-2}$	1.4779	$6.3255\times {10}^{-2}$	1.4782	$3.3041\times {10}^{-1}$	1.4790
		16	$8.2126\times {10}^{-3}$	1.4738	$3.1641\times {10}^{-2}$	1.4738	$1.6523\times {10}^{-1}$	1.4745
0.3	0.1	5	$4.7550\times {10}^{-2}$	−	$1.8327\times {10}^{-1}$	−	$9.5743\times {10}^{-1}$	−
		8	$2.2788\times {10}^{-2}$	1.5650	$8.7752\times {10}^{-2}$	1.5669	$4.5859\times {10}^{-1}$	1.5661
		10	$1.6385\times {10}^{-2}$	1.4783	$6.3087\times {10}^{-2}$	1.4788	$3.2967\times {10}^{-1}$	1.4791
		16	$8.1952\times {10}^{-3}$	1.4740	$3.1552\times {10}^{-2}$	1.4742	$1.6486\times {10}^{-1}$	1.4745

, for different $N=5$, 8, 10, 16, we give the error results and convergence rates for the MFE scheme (5)–(6), where the spatial grid sizes are also taken as $h=\sqrt{2}/N^{2-{\alpha }_1}$. We can also see that the convergence rates in the temporal direction for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms) are close to $2-{\alpha }_1$. Furthermore, in

Table 10. Numerical results with $\tau =T/N=1/100$ and ${\alpha }_1=0.9$ in Example 2.

${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	h	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.4	0.2	$\sqrt{2}/4$	$1.2924\times {10}^{-1}$	−	$5.0230\times {10}^{-1}$	−	$2.5973\times {10}^{+0}$	−
		$\sqrt{2}/8$	$6.5229\times {10}^{-2}$	0.9865	$2.5168\times {10}^{-1}$	0.9970	$1.3119\times {10}^{+0}$	0.9853
		$\sqrt{2}/16$	$3.2696\times {10}^{-2}$	0.9964	$1.2589\times {10}^{-1}$	0.9994	$6.5768\times {10}^{-1}$	0.9962
		$\sqrt{2}/32$	$1.6364\times {10}^{-2}$	0.9986	$6.2979\times {10}^{-2}$	0.9992	$3.2916\times {10}^{-1}$	0.9986
0.3	0.1	$\sqrt{2}/4$	$1.2925\times {10}^{-1}$	−	$5.0236\times {10}^{-1}$	−	$2.5972\times {10}^{+0}$	−
		$\sqrt{2}/8$	$6.5231\times {10}^{-2}$	0.9865	$2.5169\times {10}^{-1}$	0.9971	$1.3119\times {10}^{+0}$	0.9853
		$\sqrt{2}/16$	$3.2696\times {10}^{-2}$	0.9964	$1.2589\times {10}^{-1}$	0.9995	$6.5765\times {10}^{-1}$	0.9962
		$\sqrt{2}/32$	$1.6362\times {10}^{-2}$	0.9988	$6.2971\times {10}^{-2}$	0.9994	$3.2912\times {10}^{-1}$	0.9987

,

Table 11. Numerical results with $\tau =T/N=1/100$ and ${\alpha }_1=0.7$ in Example 2.

${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	h	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.4	0.2	$\sqrt{2}/4$	$1.2926\times {10}^{-1}$	−	$5.0244\times {10}^{-1}$	−	$2.5971\times {10}^{+0}$	−
		$\sqrt{2}/8$	$6.5237\times {10}^{-2}$	0.9866	$2.5172\times {10}^{-1}$	0.9972	$1.3118\times {10}^{+0}$	0.9853
		$\sqrt{2}/16$	$3.2697\times {10}^{-2}$	0.9965	$1.2590\times {10}^{-1}$	0.9996	$6.5758\times {10}^{-1}$	0.9963
		$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2954\times {10}^{-2}$	0.9999	$3.2901\times {10}^{-1}$	0.9990
0.3	0.1	$\sqrt{2}/4$	$1.2927\times {10}^{-1}$	−	$5.0249\times {10}^{-1}$	−	$2.5970\times {10}^{+0}$	−
		$\sqrt{2}/8$	$6.5238\times {10}^{-2}$	0.9866	$2.5172\times {10}^{-1}$	0.9973	$1.3118\times {10}^{+0}$	0.9853
		$\sqrt{2}/16$	$3.2697\times {10}^{-2}$	0.9965	$1.2590\times {10}^{-1}$	0.9996	$6.5758\times {10}^{-1}$	0.9963
		$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2954\times {10}^{-2}$	0.9999	$3.2900\times {10}^{-1}$	0.9991

and

Table 12. Numerical results with $\tau =T/N=1/100$ and ${\alpha }_1=0.5$ in Example 2.

${\mathit{\alpha }}_2$	${\mathit{\alpha }}_3$	h	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.4	0.2	$\sqrt{2}/4$	$1.2927\times {10}^{-1}$	−	$5.0252\times {10}^{-1}$	−	$2.5970\times {10}^{+0}$	−
		$\sqrt{2}/8$	$6.5240\times {10}^{-2}$	0.9866	$2.5173\times {10}^{-1}$	0.9973	$1.3118\times {10}^{+0}$	0.9853
		$\sqrt{2}/16$	$3.2698\times {10}^{-2}$	0.9965	$1.2590\times {10}^{-1}$	0.9996	$6.5757\times {10}^{-1}$	0.9963
		$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2955\times {10}^{-2}$	0.9999	$3.2899\times {10}^{-1}$	0.9991
0.3	0.1	$\sqrt{2}/4$	$1.2928\times {10}^{-1}$	−	$5.0257\times {10}^{-1}$	−	$2.5970\times {10}^{+0}$	−
		$\sqrt{2}/8$	$6.5241\times {10}^{-2}$	0.9866	$2.5174\times {10}^{-1}$	0.9974	$1.3118\times {10}^{+0}$	0.9853
		$\sqrt{2}/16$	$3.2698\times {10}^{-2}$	0.9966	$1.2590\times {10}^{-1}$	0.9996	$6.5757\times {10}^{-1}$	0.9963
		$\sqrt{2}/32$	$1.6359\times {10}^{-2}$	0.9991	$6.2955\times {10}^{-2}$	0.9999	$3.2899\times {10}^{-1}$	0.9991

, we also fix $N=100$ and take $h=\sqrt{2}/4$, $\sqrt{2}/8$, $\sqrt{2}/16$, $\sqrt{2}/32$, give the error results and convergence rates, and see that the convergence rates in the spatial direction for u and $\mathit{\lambda }$ in the above corresponding discrete norms are also close to 1.

Based on the numerical results in Table 1, Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8, Table 9, Table 10, Table 11 and Table 12 obtained from the above two test examples, we can see that the convergence rates in the spatial and temporal directions for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ norm) are in agreement with the theoretical results in Theorem 3, and those for $\mathit{\lambda }$ (in the discrete $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norm) are higher than the theoretical result. These results fully demonstrate that the proposed MFE method for the multi-term TFRD equations is effective.

Example 3. 

Consider the two-term TFRD equation in Example 1 with weak regularity solutions near the initial time $t=0$, where the source function f is taken by

$$\begin{array}{cc}\hfill f(\mathit{x},t)=& (\frac{2}{\Gamma (3-{\alpha }_1)}{t}^{2-{\alpha }_1}+\Gamma (1+{\alpha }_1)+\frac{\Gamma (2+{\alpha }_2)}{\Gamma (2+{\alpha }_2-{\alpha }_1)}{t}^{1+{\alpha }_2-{\alpha }_1}\hfill \\ \hfill & +\frac{2}{\Gamma (3-{\alpha }_2)}{t}^{2-{\alpha }_2}+\frac{\Gamma (1+{\alpha }_1)}{\Gamma (1+{\alpha }_1-{\alpha }_2)}{t}^{{\alpha }_1-{\alpha }_2}+\Gamma (2+{\alpha }_2)t\hfill \\ \hfill & +(2{\pi }^2+p\left(\mathit{x}\right))(t^2+t^{{\alpha }_1}+t^{1+{\alpha }_2}))sin\left(\pi x_1\right)sin\left(\pi x_2\right).\hfill \end{array}$$

And we can also find the analytical solutions for variables u and λ as follows:

$$\begin{array}{cc}\hfill & u(\mathit{x},t)=(t^2+t^{{\alpha }_1}+t^{1+{\alpha }_2})sin\left(\pi x_1\right)sin\left(\pi x_2\right),\hfill \\ \hfill & \mathit{\lambda }(\mathit{x},t)=-\pi (t^2+t^{{\alpha }_1}+t^{1+{\alpha }_2})(cos\left(\pi x_1\right)sin\left(\pi x_2\right),sin\left(\pi x_1\right)cos\left(\pi x_2\right)).\hfill \end{array}$$

In this example, we will select the graded mesh to discretize the interval $[0,T]$ and set $t_n=T{(n/N)}^{\gamma }$, for $n=0$, 1, 2, ⋯, N, where constant $\gamma \ge 1$ is the temporal graded mesh parameter. The ideal optimal error estimates for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms) should be $O(N^{-min\{\gamma {\alpha }_1,2-{\alpha }_1\}}+h)$. Here, we will mainly test the convergence rates in the temporal direction with the graded mesh parameter $\gamma =1$ and $(2-{\alpha }_1)/{\alpha }_1$. We first conduct numerical experiments with $\gamma =1$. Then, the optimal convergence rate in the temporal direction is ${\alpha }_1$. For fractional parameters ${\alpha }_1=0.9$, $0.7$, $0.5$ and ${\alpha }_2=0.1$, $0.4$, we take the time mesh parameter $N=20$, 40, 80, 160 and special spatial grid parameters: (i) when ${\alpha }_1=0.9$, take $h\approx 2\sqrt{2}/N^{{\alpha }_1}$ ; (ii) when ${\alpha }_1=0.7$, take $h\approx \sqrt{2}/N^{{\alpha }_1}$ ; (iii) when ${\alpha }_1=0.5$, take $h\approx \sqrt{2}/\left(2N^{{\alpha }_1}\right)$. Then, we give the numerical results in

Table 13. Numerical results with ${\alpha }_1=0.9$ and graded mesh parameter $\gamma =1$ in Example 3.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	20	$1.9572\times {10}^{-1}$	−	$7.5520\times {10}^{-1}$	−	$3.9354\times {10}^{+0}$	−
	40	$1.1208\times {10}^{-1}$	0.8042	$4.3165\times {10}^{-1}$	0.8070	$2.2539\times {10}^{+0}$	0.8041
	80	$6.0396\times {10}^{-2}$	0.8920	$2.3245\times {10}^{-1}$	0.8930	$1.2146\times {10}^{+0}$	0.8920
	160	$3.2722\times {10}^{-2}$	0.8842	$1.2591\times {10}^{-1}$	0.8845	$6.5806\times {10}^{-1}$	0.8842
0.4	20	$1.9571\times {10}^{-1}$	−	$7.5516\times {10}^{-1}$	−	$3.9354\times {10}^{+0}$	−
	40	$1.1208\times {10}^{-1}$	0.8042	$4.3164\times {10}^{-1}$	0.8069	$2.2540\times {10}^{+0}$	0.8041
	80	$6.0396\times {10}^{-2}$	0.8920	$2.3244\times {10}^{-1}$	0.8929	$1.2146\times {10}^{+0}$	0.8920
	160	$3.2722\times {10}^{-2}$	0.8842	$1.2591\times {10}^{-1}$	0.8845	$6.5806\times {10}^{-1}$	0.8842

,

Table 14. Numerical results with ${\alpha }_1=0.7$ and graded mesh parameter $\gamma =1$ in Example 3.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	20	$1.9573\times {10}^{-1}$	−	$7.5526\times {10}^{-1}$	−	$3.9353\times {10}^{+0}$	−
	40	$1.2068\times {10}^{-1}$	0.6976	$4.6487\times {10}^{-1}$	0.7002	$2.4269\times {10}^{+0}$	0.6974
	80	$7.4765\times {10}^{-2}$	0.6908	$2.8779\times {10}^{-1}$	0.6918	$1.5035\times {10}^{+0}$	0.6907
	160	$4.4872\times {10}^{-2}$	0.7365	$1.7268\times {10}^{-1}$	0.7369	$9.0241\times {10}^{-1}$	0.7365
0.4	20	$1.9572\times {10}^{-1}$	−	$7.5523\times {10}^{-1}$	−	$3.9354\times {10}^{+0}$	−
	40	$1.2068\times {10}^{-1}$	0.6976	$4.6486\times {10}^{-1}$	0.7001	$2.4269\times {10}^{+0}$	0.6974
	80	$7.4765\times {10}^{-2}$	0.6908	$2.8779\times {10}^{-1}$	0.6918	$1.5035\times {10}^{+0}$	0.6907
	160	$4.4872\times {10}^{-2}$	0.7365	$1.7268\times {10}^{-1}$	0.7369	$9.0241\times {10}^{-1}$	0.7365

and

Table 15. Numerical results with ${\alpha }_1=0.5$ and graded mesh parameter $\gamma =1$ in Example 3.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	20	$1.7410\times {10}^{-1}$	−	$6.7141\times {10}^{-1}$	−	$3.5006\times {10}^{+0}$	−
	40	$1.2069\times {10}^{-1}$	0.5286	$4.6487\times {10}^{-1}$	0.5303	$2.4269\times {10}^{+0}$	0.5285
	80	$8.7212\times {10}^{-2}$	0.4687	$3.3576\times {10}^{-1}$	0.4694	$1.7538\times {10}^{+0}$	0.4686
	160	$6.2812\times {10}^{-2}$	0.4735	$2.4175\times {10}^{-1}$	0.4739	$1.2632\times {10}^{+0}$	0.4735
0.4	20	$1.7410\times {10}^{-1}$	−	$6.7139\times {10}^{-1}$	−	$3.5006\times {10}^{+0}$	−
	40	$1.2069\times {10}^{-1}$	0.5286	$4.6487\times {10}^{-1}$	0.5303	$2.4269\times {10}^{+0}$	0.5285
	80	$8.7212\times {10}^{-2}$	0.4687	$3.3576\times {10}^{-1}$	0.4694	$1.7538\times {10}^{+0}$	0.4686
	160	$6.2812\times {10}^{-2}$	0.4735	$2.4175\times {10}^{-1}$	0.4739	$1.2632\times {10}^{+0}$	0.4735

, which show that the convergence rates in the temporal direction for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms) are close to ${\alpha }_1$.

Next, we conduct numerical experiments with $\gamma =(2-{\alpha }_1)/{\alpha }_1$. Then, the optimal convergence rate is $2-{\alpha }_1$. We take the time mesh parameter $N=5$, 8, 10, 16 and the spatial grid parameter $h=\sqrt{2}/N^{2-{\alpha }_1}$. Then, we give the numerical results in

Table 16. Numerical results with ${\alpha }_1=0.9$ and graded mesh parameter $\gamma =\frac{2-{\alpha }_1}{{\alpha }_1}$ in Example 3.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	5	$2.6024\times {10}^{-1}$	−	$1.0062\times {10}^{+0}$	−	$5.2342\times {10}^{+0}$	−
	8	$1.5673\times {10}^{-1}$	1.0789	$6.0412\times {10}^{-1}$	1.0855	$3.1526\times {10}^{+0}$	1.0787
	10	$1.2067\times {10}^{-1}$	1.1718	$4.6479\times {10}^{-1}$	1.1750	$2.4273\times {10}^{+0}$	1.1718
	16	$7.1368\times {10}^{-2}$	1.1175	$2.7470\times {10}^{-1}$	1.1189	$1.4355\times {10}^{+0}$	1.1176
0.4	5	$2.6021\times {10}^{-1}$	−	$1.0060\times {10}^{+0}$	−	$5.2351\times {10}^{+0}$	−
	8	$1.5673\times {10}^{-1}$	1.0787	$6.0410\times {10}^{-1}$	1.0852	$3.1531\times {10}^{+0}$	1.0787
	10	$1.2067\times {10}^{-1}$	1.1716	$4.6479\times {10}^{-1}$	1.1748	$2.4276\times {10}^{+0}$	1.1718
	16	$7.1372\times {10}^{-2}$	1.1174	$2.7472\times {10}^{-1}$	1.1188	$1.4357\times {10}^{+0}$	1.1176

,

Table 17. Numerical results with ${\alpha }_1=0.7$ and graded mesh parameter $\gamma =\frac{2-{\alpha }_1}{{\alpha }_1}$ in Example 3.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	5	$1.9568\times {10}^{-1}$	−	$7.5499\times {10}^{-1}$	−	$3.9360\times {10}^{+0}$	−
	8	$1.0462\times {10}^{-1}$	1.3323	$4.0284\times {10}^{-1}$	1.3365	$2.1043\times {10}^{+0}$	1.3323
	10	$7.8497\times {10}^{-2}$	1.2872	$3.0217\times {10}^{-1}$	1.2887	$1.5789\times {10}^{+0}$	1.2873
	16	$4.2450\times {10}^{-2}$	1.3079	$1.6336\times {10}^{-1}$	1.3086	$8.5380\times {10}^{-1}$	1.3081
0.4	5	$1.9566\times {10}^{-1}$	−	$7.5492\times {10}^{-1}$	−	$3.9370\times {10}^{+0}$	−
	8	$1.0462\times {10}^{-1}$	1.3320	$4.0289\times {10}^{-1}$	1.3361	$2.1049\times {10}^{+0}$	1.3323
	10	$7.8506\times {10}^{-2}$	1.2870	$3.0221\times {10}^{-1}$	1.2885	$1.5793\times {10}^{+0}$	1.2874
	16	$4.2457\times {10}^{-2}$	1.3078	$1.6339\times {10}^{-1}$	1.3084	$8.5400\times {10}^{-1}$	1.3081

and

Table 18. Numerical results with ${\alpha }_1=0.5$ and graded mesh parameter $\gamma =\frac{2-{\alpha }_1}{{\alpha }_1}$ in Example 3.

${\mathit{\alpha }}_2$	N	$\mathit{u}–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\mathit{L}}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left({\left({\mathit{L}}^2\right)}^2\right)$	Rates	$\mathit{\lambda }–{\widehat{\mathit{L}}}^{\mathit{\infty }}\left(\mathit{H}\left(\mathbf{div}\right)\right)$	Rates
0.1	5	$1.4253\times {10}^{-1}$	−	$5.4923\times {10}^{-1}$	−	$2.8673\times {10}^{+0}$	−
	8	$6.8272\times {10}^{-2}$	1.5661	$2.6278\times {10}^{-1}$	1.5685	$1.3733\times {10}^{+0}$	1.5663
	10	$4.9083\times {10}^{-2}$	1.4788	$1.8890\times {10}^{-1}$	1.4794	$9.8727\times {10}^{-1}$	1.4790
	16	$2.4548\times {10}^{-2}$	1.4742	$9.4464\times {10}^{-2}$	1.4744	$4.9373\times {10}^{-1}$	1.4744
0.4	5	$1.4255\times {10}^{-1}$	−	$5.4933\times {10}^{-1}$	−	$2.8688\times {10}^{+0}$	−
	8	$6.8306\times {10}^{-2}$	1.5654	$2.6296\times {10}^{-1}$	1.5675	$1.3743\times {10}^{+0}$	1.5659
	10	$4.9113\times {10}^{-2}$	1.4783	$1.8905\times {10}^{-1}$	1.4788	$9.8804\times {10}^{-1}$	1.4787
	16	$2.4567\times {10}^{-2}$	1.4739	$9.4561\times {10}^{-2}$	1.4740	$4.9416\times {10}^{-1}$	1.4742

and find that the convergence rates in the temporal direction for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ and $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norms) are close to $2-{\alpha }_1$. Based on the above discussion, we know that the MFE scheme (5)–(6) with the temporal graded mesh for solving the multi-term TFRD equations with the initial layer is also feasible and effective.

7. Conclusions

This work presents a Raviart–Thomas MFE method for solving the multi-term TFRD equations with variable coefficients by using the well-known $L1$ formula. The existence, uniqueness, and unconditional stability of the discrete solution are discussed, and the optimal a priori error estimates for u (in the discrete $L^{\infty }\left(L^2(\Omega )\right)$ norm) and $\mathit{\lambda }$ (in the discrete $L^{\infty }\left({\left(L^2(\Omega )\right)}^2\right)$ norm) and the suboptimal a priori error estimate for $\mathit{\lambda }$ (in the discrete $L^{\infty }\left(\mathit{H}(\mathrm{div},\Omega )\right)$ norm) are obtained in this work. In addition, some numerical results are given to demonstrate the effectiveness of the proposed MFE method. In future research, we will try to give theoretical analysis for the MFE method with the temporal graded mesh to solve some FPDEs with the initial layer at $t=0$ and apply the MFE method to solve more FPDEs in scientific and engineering fields.