Efficient Temporal Third/Fourth-Order Finite Element Method for a Time-Fractional Mobile/Immobile Transport Equation with Smooth and Nonsmooth Data

Abstract

In recent years, the numerical theory of fractional models has received more and more attention from researchers, due to the broad and important applications in materials and mechanics, anomalous diffusion processes and other physical phenomena. In this paper, we propose two efficient finite element schemes based on convolution quadrature for solving the time-fractional mobile/immobile transport equation with the smooth and nonsmooth data. In order to deal with the weak singularity of solution near $t=0$, we choose suitable corrections for the derived schemes to restore the third/fourth-order accuracy in time. Error estimates of the two fully discrete schemes are presented with respect to data regularity. Numerical examples are given to illustrate the effectiveness of the schemes.

1. Introduction

In this paper, we are interested in the robust high-order numerical methods for the following time-fractional mobile/immobile transport equation [1]:

$$\begin{array}{ccc}\hfill ({\kappa }_1{\partial }_t+{\kappa }_2{}_{RL}{D}_{0,t}^{\alpha })(u(x,t)-v)& =& \mu \Delta u(x,t)+f(x,t),\mathrm{in}\Omega \times (0,T],\hfill \end{array}$$

(1)

which is subject to the boundary condition $u=0$ on $\partial \Omega \times (0,T]$ and the initial condition $u(x,0)=v\left(x\right)$ in $\overline{\Omega }$. Here, $\Omega \subset {\mathbb{R}}^2$ is a bounded convex polygonal domain and the parameters ${\kappa }_1,{\kappa }_2$, and $\mu$ are positive constants. The source term f and the initial condition v are both given functions. The notation ${}_{RL}{D}_{0,t}^{\alpha }$ with $\alpha \in (0,1)$ denotes the Riemann–Liouville derivative defined by

$${}_{RL}{D}_{0,t}^{\alpha }u(·,t)=\frac{1}{\Gamma (1-\alpha )}\frac{\partial }{\partial t}{\int }_0^t{(t-s)}^{-\alpha }u(·,s)\mathrm{d}s.$$

In recent years, Equation (1) has received extensive attention, and it has significant applications in diverse fields, especially in the modeling of groundwater solutes [2,3,4]. There are some numerical studies for Equation (1) or its variants; see [1,5,6,7,8,9] and the references therein.

In our recent paper [1], we studied the solution regularity of Equation (1) and proposed two efficient finite element schemes by employing convolution quadrature based on the backward Euler and the second-order backward difference methods. It is proved that the proposed schemes are robust with respect to data regularity, including the nonsmooth initial data, i.e., $v\in L^2(\Omega )$. One of the distinct features of the two schemes is that they do not impose any compatibility conditions on the source term f, so this avoids some unrealistic assumption for the solution regularity, i.e., the assumption of the sufficiently smooth solution. When dealing with the nonsmooth solution problem, we observe numerically that the convolution quadrature generated by the high-order backward difference formula (denoted by BDF) is only of first-order accuracy when no corrections are added. This motivates us to investigate how to choose some suitable corrections in order to restore the desired high-order accuracy.

Since the convolution quadrature has a nice stable discrete structure, for which it is easy to perform the analysis of the resulting numerical scheme [10,11], it seems that it is suitable for developing high-order numerical schemes of the fractional model with the Riemann–Liouville derivative. In order to restore the high-order accuracy, one usually needs to add some corrections to the original scheme. There are a few papers on this direction [12,13,14,15,16,17,18], just to name a few. One may refer to two review papers [19,20] for more details. In [14,16,17], the convolution quadratures generated by the backward Euler and the second-order backward difference methods were applied to obtain the robust scheme for the time-fractional Cattaneo equation, the two-term time-fractional diffusion-wave equation, and the modified anomalous subdiffusion equation, respectively. In [12], the authors constructed and analyzed an efficient fractional Crank–Nicolson finite element scheme with two corrections at the starting time levels. They proved that the proposed scheme is of second-order accuracy in time for both smooth and nonsmooth data. Recently, Wang et al. improved the results by adding only one single correction at the starting time level and presented the optimal error estimates of their resulting scheme [18]. In [13], Jin et al. developed the BDFk (with the k denoting the order of BDF) convolution quadrature with suitable correction terms to deal with the temporal discretization of the fractional evolution equation. Later on, Shi and Chen extended the BDFk convolution quadrature to numerically solve the fractional Feynman–Kac equation with Lévy flight [15]. As far as we know, there is no convolution quadrature generated by high-order BDF to solve Equation (1) with smooth and nonsmooth problem data. The goal of this paper is to fill this gap.

The contributions of this paper are listed below. First, based on the convolution quadrature generated by BDF3/BDF4, we develop the robust temporal third/fourth-order finite element scheme for problem (1) by carefully choosing the suitable corrections at the starting two/three steps. Second, the error estimates of the resulting schemes are expressed in term of data regularity, and are proved rigorously, cf. Theorems 2–5. Third, the designed numerical tests show that the two schemes indeed yield the temporal third/fourth-order theoretical accuracy, so this further verifies numerically the correctness of the error estimates and illustrates the significant improvement in the temporal accuracy of our schemes compared to that of [1], cf. Tables 4–6 in Section 6.

The rest of this paper is organized as follows. In Section 2, we propose the temporal third-order finite element scheme based on the convolution quadrature generated by BDF3. In Section 3, the method to obtain the suitable correction terms is discussed in detail. In Section 4, we prove that the BDF3 is of third-order accuracy, with respect to smooth and nonsmooth data. The generalization of the BDF3 is given in Section 5, and numerical examples are provided in Section 6 to verify the effectiveness of the two proposed numerical schemes. Finally, we present the conclusion in the last section, i.e., Section 7. The symbol c in this paper denotes a positive constant that is independent of the temporal and spatial step sizes.

2. The Temporal Third-Order Numerical Scheme

In this section, we propose a temporal third-order fully discrete scheme for Equation (1) based on convolution quadrature in time and the standard Galerkin finite element method in space.

2.1. Spatial Semidiscrete Scheme

We first introduce the semidiscrete scheme by using the finite element method in space and recall the corresponding error estimates [1].

Denote the continuous piecewise linear finite element space $V_h$ as follows:

$$V_h=\left\{v_h\in H_0^1(\Omega ):v_h{|}_K\mathrm{is}\mathrm{the}\mathrm{linear}\mathrm{function}\mathrm{and}\forall K\in {\mathcal{T}}_h\right\},$$

where ${\mathcal{T}}_h$ is the partition of the domain $\Omega$ with $h={max}_{K\in {\mathcal{T}}_h}{h}_K$, and $h_K$ is the diameter of the triangles K [21]. The semidiscrete scheme for (1) is given by the following: Find $u_h\left(t\right)=u_h(·,t)\in V_h$ such that the following holds:

$$(({\kappa }_1{\partial }_t+{\kappa }_2{}_{RL}{D}_{0,t}^{\alpha })(u_h-v_h),\chi )+\mu (\nabla u_h,\nabla \chi )=(f,\chi ),\forall \chi \in V_h,$$

with the initial value condition $u_h\left(0\right)=v_h\in V_h$. Here, $v_h$ is the proper approximation to the function v. In this paper, we choose $v_h=P_{h}v$ for $v\in L^2(\Omega )$ and $v_h=R_{h}v$ when $v\in D(\Delta )=H_0^1(\Omega )\cap H^2(\Omega )$ [21]. Here, $P_h$ and $R_h$ are the $L^2(\Omega )$ orthogonal projection operator and Ritz projection operator, respectively. By using the discrete Laplacian $-{\Delta }_h:V_h\to V_h$, defined by

$$-\left({\Delta }_h\phi ,\chi \right)=(\nabla \phi ,\nabla \chi ),\forall \phi ,\chi \in V_h,$$

one further has

$$({\kappa }_1{\partial }_t+{\kappa }_2{}_{RL}{D}_{0,t}^{\alpha })(u_h\left(t\right)-v_h)-\mu {\Delta }_h{u}_h\left(t\right)=f_h\left(t\right),t>0,$$

(2)

where $f_h\left(t\right)=P_{h}f\left(t\right)$.

The error estimates for the semidiscrete scheme (2) were investigated in [1]; we collect the corresponding results as follows.

Theorem 1.

Denote the space $\mathbf{S}:=L^{\infty }(0,T;L^2(\Omega ))$. Then, for the semidiscrete scheme (2), we have

(a) 

If $v\in L^2(\Omega )$, then $\parallel u\left(t\right)-u_h\left(t\right)\parallel \le c{h}^2(t^{-\alpha }\parallel v\parallel +|ln(t^{\alpha }/h^2){|\parallel f\parallel }_{\mathbf{S}}).$

(b) 

If $v\in D(\Delta )$, then $\parallel u\left(t\right)-u_h\left(t\right)\parallel \le c{h}^2(\parallel \Delta v\parallel +|ln(t^{\alpha }/h^2){|\parallel f\parallel }_{\mathbf{S}}).$

2.2. Fully Discrete Schemes

We divide the time interval $[0,T]$ into uniform grids ${\{t_n=n\tau \}}_{n=0}^N$ with $\tau =T/N$. Here, N is a given positive integer. Based on the convolution quadrature generated by BDFk, the numerical approximation of ${}_{RL}{D}_{0,t}^{\alpha }{u}_h\left(t\right)$ at $t=t_n$ is described by

$${\partial }_{\tau }^{\alpha }{u}_h^n=\frac{1}{{\tau }^{\alpha }}\sum _{j=0}^n{w}_{n-j}{u}_h^j,$$

in which $u_h^n=u_h\left(t_n\right)$ and the weights $w_j$ are provided by

$$\frac{1}{{\tau }^{\alpha }}\sum _{j=0}^{\infty }{w}_j{\xi }^j={\delta }_{\tau }{\left(\xi \right)}^{\alpha },$$

(3)

with ${\delta }_{\tau }\left(\xi \right)=\frac{1}{\tau }{\sum }_{j=1}^k\frac{1}{j}{(1-\xi )}^j$.

Applying the convolution quadrature generated by BDF3/BDF4 for solving the semidiscrete scheme (2), one may obtain the following fully discrete finite element scheme: Find the numerical approximation $U_h^n,n\ge 1$ of $u\left(t_n\right)$, such that the following holds:

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })(U_h^n-v_h)-\mu {\Delta }_h{U}_h^n=F_h^n,$$

(4)

where $n=1,2,\cdots ,N,U_h^0=v_h$, and $F_h^n=f_h\left(t_n\right)$. It is known that the scheme (4) generally exhibits only temporal first-order accuracy, due to the low regularity of the solution near $t=0$. Such a fact is also verified by the numerical examples in Section 6.

In what follows of this section, we consider the scheme (4) generated by BDF3 and use the idea presented in [13] to derive a modified scheme which can restore the temporal third-order accuracy. To this end, we define ${\tilde{f}}_h\left(t\right)=f_h\left(t\right)-f_h\left(0\right)-t{f}_h^{\prime }\left(0\right)$ and denote $w_h\left(t\right)=u_h\left(t\right)-v$. The semidiscrete scheme (2) can be rewritten as follows:

$$({\kappa }_1{\partial }_t+{\kappa }_2{}_{RL}{D}_{0,t}^{\alpha })w_h\left(t\right)-\mu {\Delta }_h{w}_h\left(t\right)=\mu {\Delta }_h{v}_h+f_h\left(0\right)+t{f}_h^{\prime }\left(0\right)+{\tilde{f}}_h\left(t\right).$$

(5)

In view of the construction of BDF2 in [1], we may naturally consider the following two corrections for (4) at the starting two steps:

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })W_h^1-\mu {\Delta }_h{W}_h^1={\lambda }_{11}(\mu {\Delta }_h{v}_h+f_h\left(0\right))+{\lambda }_{12}\tau f_h^{\prime }\left(0\right)+{\tilde{f}}_h^1,$$

(6)

and

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })W_h^2-\mu {\Delta }_h{W}_h^2={\lambda }_{21}(\mu {\Delta }_h{v}_h+f_h\left(0\right))+2{\lambda }_{22}\tau f_h^{\prime }\left(0\right)+{\tilde{f}}_h^2,$$

(7)

where $W_h^n=U_h^n-v,n\ge 0$. The four parameters ${\lambda }_{ij}(i,j=1,2)$ are unknown and need to be determined.

3. Determination of the Two Corrections

In this part, we determine the four parameters ${\lambda }_{ij}(i,j=1,2)$ in the two corrections (6) and (7) based on the comparison of the semidiscrete solution in (5) and the fully discrete solution in (4). Let the notation $\widehat{(·)}$ be the Laplace transform of $(·)$, i.e., $\widehat{g}\left(z\right):={\int }_0^{\infty }{e}^{-zt}g\left(t\right)\mathrm{d}t$. Denote the kernel $K\left(z\right)$ as follows:

$$K\left(z\right)=z^{-1}{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}.$$

(8)

Define the generating function $\overline{g}\left(\xi \right)={\sum }_{n=0}^{\infty }{g}^n{\xi }^n$ with a given sequence ${\left(g^n\right)}_{n=0}^{\infty }$. We have the following results for the integral representations.

Lemma 1.

The semidiscrete solution in (5) and the fully discrete solution in (4) (with $W_h^n$ replacing $U_h^n$) have the following integral representations:

$$w_h\left(t\right)=\frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }}{e}^{zt}(K\left(z\right)(\mu {\Delta }_h{v}_h+f_h\left(0\right))+z^{-1}K\left(z\right)f_h^{\prime }\left(0\right)+zK\left(z\right){\widehat{\tilde{f}}}_h\left(z\right))dz,$$

(9)

and

$$\begin{array}{ccc}\hfill W_h^n& =& \frac{1}{2\pi \mathrm{i}}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)(\tau {\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right)(\mu {\Delta }_h{v}_h+f_h\left(0\right))\hfill \\ & & +{\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_2\left(e^{-z\tau }\right){\tau }^2{f}_h^{\prime }\left(0\right)+\tau {\delta }_{\tau }\left(e^{-z\tau }\right){\overline{\tilde{F}}}_h\left(e^{-z\tau }\right))\mathrm{d}z,\hfill \end{array}$$

(10)

where

$${\varpi }_1\left(z\right)=\frac{z}{1-z}+({\lambda }_{11}-1)z+({\lambda }_{12}-1)z^2,$$

and

$${\varpi }_2\left(z\right)=\frac{z}{{(1-z)}^2}+({\lambda }_{12}-1)z+(2{\lambda }_{22}-1)z^2.$$

Here, the contours ${\Gamma }_{\theta ,\delta }:=\{r{e}^{\pm i\theta }:r\ge \delta \}\cup \{\delta e^{i\varphi }:\left|\varphi \right|\le \theta \}$ with $\delta >0$ and $\theta \in (\pi /2,\pi )$ and ${\Gamma }_{\theta ,\delta }^{\tau }:=\{z\in {\Gamma }_{\theta ,\delta }:|Im\left(z\right)|\le \pi /\tau \}$.

Proof. 

By using the Laplace transform technology, we can easily obtain the integral representation (9) from (5). One may refer to [1] for further details.

It remains to derive the fully discrete solution in (4) based on discrete Laplace transform, i.e., generating function.

Let ${\tilde{F}}_h^n=F_h^n-f_h^0-t_n{f}_h^{\prime }\left(0\right)$. The scheme (4) has the following equivalent form:

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })W_h^n-\mu {\Delta }_h{W}_h^n=\mu {\Delta }_h{v}_h+f_h\left(0\right)+t_n{f}_h^{\prime }\left(0\right)+{\tilde{F}}_h^n,n\ge 3,$$

(11)

with the two initial corrections (6) and (7). Multiplying ${\xi }^n$ on both sides of (11), and summing up for $n=1,2,\cdots$, we obtain the following:

$$\begin{array}{ccc}\hfill & & \sum _{n=1}^{\infty }({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })W_h^n{\xi }^n-\mu \sum _{n=1}^{\infty }{\Delta }_h{W}_h^n{\xi }^n\hfill \\ & =& (\sum _{n=3}^{\infty }{\xi }^n+{\lambda }_{11}\xi +{\lambda }_{12}{\xi }^2)(\mu {\Delta }_h{v}_h+f_h\left(0\right))\hfill \\ & & +(\sum _{n=3}^{\infty }n{\xi }^n+{\lambda }_{12}\xi +2{\lambda }_{22}{\xi }^2)\tau f_h^{\prime }\left(0\right)+{\overline{\tilde{F}}}_h\left(\xi \right).\hfill \end{array}$$

(12)

Noting $W_h^0=0$, one has

$$\sum _{n=1}^{\infty }({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })W_h^n{\xi }^n=({\kappa }_1{\delta }_{\tau }\left(\xi \right)+{\kappa }_2{\delta }_{\tau }{\left(\xi \right)}^{\alpha }){\overline{W}}_h\left(\xi \right).$$

From (12), we have

$$\begin{array}{ccc}\hfill & & ({\kappa }_1{\delta }_{\tau }\left(\xi \right)+{\kappa }_2{\delta }_{\tau }{\left(\xi \right)}^{\alpha }-\mu {\Delta }_h){\overline{W}}_h\left(\xi \right)\hfill \\ & =& \left(\frac{\xi }{1-\xi }+({\lambda }_{11}-1)\xi +({\lambda }_{12}-1){\xi }^2\right)(\mu {\Delta }_h{v}_h+f_h\left(0\right))\hfill \\ & & +\left(\frac{\xi }{{(1-\xi )}^2}+({\lambda }_{12}-1)\xi +(2{\lambda }_{22}-1){\xi }^2\right)\tau f_h^{\prime }\left(0\right)+{\overline{\tilde{F}}}_h\left(\xi \right).\hfill \end{array}$$

(13)

Thus,

$$\begin{array}{ccc}\hfill {\overline{W}}_h\left(\xi \right)& =& {({\kappa }_1{\delta }_{\tau }\left(\xi \right)+{\kappa }_2{\delta }_{\tau }{\left(\xi \right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\hfill \\ & & ·((\frac{\xi }{1-\xi }+({\lambda }_{11}-1)\xi +({\lambda }_{12}-1){\xi }^2)(\mu {\Delta }_h{v}_h+f_h\left(0\right))\hfill \\ & & +(\frac{\xi }{{(1-\xi )}^2}+({\lambda }_{12}-1)\xi +(2{\lambda }_{22}-1){\xi }^2)\tau f_h^{\prime }\left(0\right)+{\overline{\tilde{F}}}_h\left(\xi \right)),\hfill \end{array}$$

from which we have

$$W_h^n=\frac{1}{2\pi \mathrm{i}}{\int }_{\left|\xi \right|=r}{\xi }^{-n-1}{\overline{W}}_h\left(\xi \right)\mathrm{d}\xi =\frac{\tau }{2\pi \mathrm{i}}{\int }_{{\Gamma }_r^{\tau }}{e}^{z{t}_n}{\overline{W}}_h\left(e^{-z\tau }\right)\mathrm{d}z,$$

(14)

with $r\in (0,1)$ and the notation ${\Gamma }_r^{\tau }:=\{z=-ln\left(r\right)/\tau +iy:y\in \mathbb{R}\mathrm{and}|y|\le \pi /\tau \}$. Here, we have used the change of variables $\xi =e^{-z\tau }$ in the second equality of (14). Denote the complex region as $\Sigma$ enclosed by ${\Gamma }_r^{\tau },{\Gamma }_{\theta ,\delta }^{\tau }$, and ${\Gamma }_{\pm }^{\tau }$, where ${\Gamma }_{\pm }^{\tau }:=\mathbb{R}\pm i\pi /\tau$. Note that the integrand $e^{z{t}_n}{\overline{W}}_h\left(e^{-z\tau }\right)$ in (14) is analytic with respect to z in $\Sigma$, so applying the Cauchy’s theorem, we have

$$\frac{\tau }{2\pi \mathrm{i}}\left({\int }_{{\Gamma }_r^{\tau }}-{\int }_{{\Gamma }_{+}^{\tau }}-{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}+{\int }_{{\Gamma }_{-}^{\tau }}\right)e^{z{t}_n}{\overline{W}}_h\left(e^{-z\tau }\right)\mathrm{d}z=0.$$

This leads to the desired result (10) since there holds

$${\int }_{{\Gamma }_{+}^{\tau }}{e}^{z{t}_n}{\overline{W}}_h\left(e^{-z\tau }\right)\mathrm{d}z={\int }_{{\Gamma }_{-}^{\tau }}{e}^{z{t}_n}{\overline{W}}_h\left(e^{-z\tau }\right)\mathrm{d}z.$$

All this completes the proof of the lemma.  □

From the representations of (9) and (10), it is intuitive to select ${\lambda }_{ij}(i,j=1,2)$ in (10) such that

$$\parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right)\tau -K\left(z\right)\parallel \le c{\tau }^3,$$

(15)

and

$$\parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_2\left(e^{-z\tau }\right){\tau }^2-z^{-1}K\left(z\right)\parallel \le c{\tau }^3,$$

(16)

for any given $z\in {\Gamma }_{\theta ,\delta }^{\tau }$.

For the left hand sides of (15) and (16), the triangle inequality leads to

$$\begin{array}{ccc}\hfill & & \parallel ({\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right)\tau -1)K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)+K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-K\left(z\right)\parallel \hfill \\ & \le & |{\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right)\tau -1|\parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)\parallel +\parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-K\left(z\right)\parallel ,\hfill \end{array}$$

(17)

and

$$\begin{array}{ccc}\hfill & & \parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}({\varpi }_2\left(e^{-z\tau }\right){\tau }^2-{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2})\hfill \\ & & +{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-1}K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-z^{-1}K\left(z\right)\parallel \hfill \\ & \le & \parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel |{\varpi }_2\left(e^{-z\tau }\right){\tau }^2-{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2}|\hfill \\ & & +\parallel {\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-1}K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-z^{-1}K\left(z\right)\parallel .\hfill \end{array}$$

(18)

We shall need the following lemma.

Lemma 2.

Let $K\left(z\right)$ and ${\delta }_{\tau }\left(\xi \right)$ be defined as in (8) and (3) with $k=3$, respectively. Then, for any $z\in {\Gamma }_{\theta ,\delta }^{\tau }$, we have

$$\parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-K\left(z\right)\parallel \le c{\tau }^3{\left|z\right|}^{2-\alpha }.$$

Proof. 

Since $K\left(z\right)=z^{-1}{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}$, we have

$$\begin{array}{ccc}\hfill & & \parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-K\left(z\right)\parallel \hfill \\ & \le & \parallel {\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-1}{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}-z^{-1}{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & \le & |{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-1}-z^{-1}|\parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & & +{\left|z\right|}^{-1}\parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}-{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel .\hfill \end{array}$$

In view of Lemma B.1. in [13], we derive that

$$\begin{array}{ccc}\hfill |{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-1}-z^{-1}|& \le & c{\tau }^3{\left|z\right|}^2,\hfill \\ \hfill \parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel & \le & cmin\left\{\right|{\delta }_{\tau }\left(e^{-z\tau }\right){|}^{-1},\left|{\delta }_{\tau }\left(e^{-z\tau }\right){|}^{-\alpha }\right\}\hfill \end{array}$$

(19)

$$\begin{array}{ccc}& \le & cmin\left\{\right|z{|}^{-1},\left|z{|}^{-\alpha }\right\},\hfill \end{array}$$

(20)

and

$$\begin{array}{ccc}& & \parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}-{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & \le & |{\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-{\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)-{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }|\hfill \\ & & ·\parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \parallel {({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & \le & c{\tau }^3{\left(\right|z|}^4+{\left|z\right|}^{3+\alpha }{)·min\{\left|z\right|}^{-1}{,\left|z\right|}^{-\alpha }\}.\hfill \end{array}$$

Thus, we further deduce that

$$\begin{array}{ccc}& & \parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-K\left(z\right)\parallel \hfill \\ & \le & c{\tau }^3{\left|z\right|}^2{min\left\{\right|z|}^{-1}{,\left|z\right|}^{-\alpha }\}+c{\tau }^3{\left(\right|z|}^3+{\left|z\right|}^{2+\alpha }{)·min\{\left|z\right|}^{-1}{,\left|z\right|}^{-\alpha }\}\hfill \\ & \le & c{\tau }^3{\left(\right|z|}^2+{\left|z\right|}^3+{\left|z\right|}^{2+\alpha }{)min\{\left|z\right|}^{-1}{,\left|z\right|}^{-\alpha }\}\le c{\tau }^3{\left|z\right|}^{2-\alpha },\hfill \end{array}$$

for $z\in {\Gamma }_{\theta ,\delta }^{\tau }$, where the last inequality holds since $\left|z\right|<1$.  □

By Lemma 2, we only need to choose

$$|{\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right){\tau -1|\le c|z|}^3{\tau }^3,$$

(21)

and

$$|{\varpi }_2\left(e^{-z\tau }\right){\tau }^2-{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2}|\le c|z|{\tau }^3.$$

(22)

It is derived in ([13], cf. (16) and (17)) that the above error estimates can be valid by choosing ${\lambda }_{11}=23/12,{\lambda }_{12}=13/12,{\lambda }_{21}=7/12$, and ${\lambda }_{22}=1/2$. So, the suitable corrections in (6) and (7) can be written as follows:

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^1-\mu {\Delta }_h{U}_h^1=11/12(\mu {\Delta }_h{v}_h+f_h\left(0\right))+1/12\tau f_h^{\prime }\left(0\right)+F_h^1,$$

(23)

and

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^2-\mu {\Delta }_h{U}_h^2=-5/12(\mu {\Delta }_h{v}_h+f_h\left(0\right))+F_h^2.$$

(24)

The fully discrete scheme

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^n-\mu {\Delta }_h{U}_h^n=F_h^n,n\ge 3,$$

(25)

with the two corrections (23) and (24) as the starting two steps is referred to as the BDF3, and can restore third-order accuracy in time, which will be proved in the next section.

4. Error Estimates of the BDF3

In this part, we present the error estimates for the modified scheme (25). We first consider the error estimate for a homogenous problem with $v\in D(\Delta )$ and $f=0$. In the proof of error estimates, we may explicitly or implicitly employ the identity ${\Delta }_h{R}_h=P_h\Delta$ and the $L^2$-stability of the orthogonal projection $P_h$ when needed.

Theorem 2.

Let u and $U_h^n$ be the solutions of (1) and (25), respectively. Assuming $v\in D(\Delta )$ and $f=0$. Then, for $v_h:=R_{h}v$, we have

$$\parallel u\left(t_n\right)-U_h^n\parallel \le c(h^2+{\tau }^3{t}_n^{\alpha -3})\parallel \Delta v\parallel .$$

Proof. 

It suffices to consider the bound of $\parallel u_h\left(t_n\right)-U_h^n\parallel$ in view of Theorem 1. From (9) and (10) in Lemma 1, we can obtain

$$\begin{array}{ccc}\hfill \parallel w_h\left(t_n\right)-W_h^n\parallel & \le & \parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}(K\left(z\right)-K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)\tau {\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right))\left(\mu {\Delta }_h{v}_h\right)\mathrm{d}z\parallel \hfill \\ & & +\parallel \frac{1}{2\pi \mathrm{i}}{\int }_{{\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}K\left(z\right)\left(\mu {\Delta }_h{v}_h\right)\mathrm{d}z\parallel =:I_1+I_2.\hfill \end{array}$$

(26)

Since $\parallel K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\parallel \le cmin\left\{\right|z|}^{-1}{,\left|z\right|}^{-\alpha }\}$, in view of (17) with Lemma 2 and letting $\delta \le 1/t_n$, we have

$$\begin{array}{ccc}\hfill I_1& \le & c{\tau }^3\parallel {\Delta }_h{v}_h\parallel \left({\int }_{\delta }^{\pi /(\tau sin\theta )}{e}^{r{t}_{n}cos\theta }{r}^{2-\alpha }\mathrm{d}r+{\int }_{-\theta }^{\theta }{e}^{\delta t_n|cos\eta |}{\delta }^{3-\alpha }\mathrm{d}\eta \right)\hfill \\ & \le & c{\tau }^3{t}_n^{\alpha -3}\parallel \Delta v\parallel .\hfill \end{array}$$

For $I_2$, we get

$$I_2\le c\parallel {\Delta }_h{v}_h\parallel {\int }_{\pi /(\tau sin\theta )}^{\infty }{e}^{t_{n}rcos\theta }{r}^{-\alpha -1}\mathrm{d}r.$$

Note that $r\ge \pi /(\tau sin\theta )$, one has $r^{-1}\le r^2{\tau }^3{(sin\theta /\pi )}^3$. So, we have

$$\begin{array}{ccc}\hfill I_2& \le & c{\tau }^3\parallel {\Delta }_h{v}_h\parallel {\int }_{\pi /(\tau sin\theta )}^{\infty }{e}^{t_{n}rcos\theta }{r}^{2-\alpha }\mathrm{d}r\hfill \\ & \le & c{t}_n^{\alpha -3}{\tau }^3\parallel {\Delta }_h{v}_h\parallel {\int }_0^{\infty }{e}^{scos\theta }{s}^{2-\alpha }\mathrm{d}s\hfill \\ & \le & c{t}_n^{\alpha -3}{\tau }^3\parallel \Delta v\parallel .\hfill \end{array}$$

Putting the error estimates of the two parts, $I_1$ and $I_2$, into (26), we complete the proof of the theorem. □

Next, we state the error estimate for the homogenous problem with $v\in L^2(\Omega )$ and $f=0$.

Theorem 3.

Let u and $U_h^n$ be the solutions of (1) and (25), respectively. Assuming $v\in L^2(\Omega )$ and $f=0$. Then, for $v_h:=P_{h}v$ we have

$$\parallel u\left(t_n\right)-U_h^n\parallel \le c(h^2{t}_n^{-\alpha }+{\tau }^3{t}_n^{-3})\parallel v\parallel .$$

Proof. 

Let $\tilde{K}\left(z\right)=\mu {({\kappa }_{1}z+{\kappa }_1{z}^{\alpha }-\mu {\Delta }_h)}^{-1}{\Delta }_h$. In view of (26), we have

$$\begin{array}{ccc}& & \parallel w_h\left(t_n\right)-W_h^n\parallel \hfill \\ & \le & \parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}\tilde{\tilde{K}}\left(z\right)v_h\mathrm{d}z\parallel \hfill \\ & & +\parallel \frac{1}{2\pi \mathrm{i}}{\int }_{{\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}{z}^{-1}\tilde{K}\left(z\right)v_h\mathrm{d}z\parallel =:I_1+I_2,\hfill \end{array}$$

where $\tilde{\tilde{K}}\left(z\right)=z^{-1}\tilde{K}\left(z\right)-\tau {\varpi }_1\left(e^{-z\tau }\right)\tilde{K}\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)$. By the triangle inequality, the term $\tilde{\tilde{K}}\left(z\right)$ has the following error estimate:

$$\parallel \tilde{\tilde{K}}\left(z\right)\parallel \le |\tau {\varpi }_1\left(e^{-z\tau }\right)-z^{-1}|\parallel \tilde{K}\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\parallel +\left|z\right|}^{-1}\parallel \tilde{K}\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)-\tilde{K}\left(z\right)\parallel =:I_{11}+I_{12}.$$

For the first term $I_{11}$, using the result in (19) and the criterion (21), we have

$$\begin{array}{ccc}& & |\tau {\varpi }_1\left(e^{-z\tau }\right)-z^{-1}|\hfill \\ & \le & |\tau {\varpi }_1\left(e^{-z\tau }\right){\delta }_{\tau }\left(e^{-z\tau }\right)({\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-1}-z^{-1})|\hfill \\ & & +|(\tau {\varpi }_1\left(e^{-z\tau }\right){\delta }_{\tau }\left(e^{-z\tau }\right)-1){\left|\right|z|}^{-1}\le c{\tau }^3{\left|z\right|}^2.\hfill \end{array}$$

So,

$$\begin{array}{ccc}\hfill I_{11}& \le & c{\tau }^3{\left|z\right|}^2\parallel \tilde{K}\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)\parallel \hfill \\ & \le & c{\tau }^3{\left|z\right|}^2\parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha })-I\parallel \hfill \\ & \le & c{\tau }^3{\left|z\right|}^2.\hfill \end{array}$$

For the second term $I_{12}$, we have

$$\begin{array}{ccc}\hfill \left|z\right|I_{12}& =& \parallel ({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }){({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\hfill \\ & & -({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }){({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & \le & \parallel \left[{\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-({\kappa }_{1}z+{\kappa }_2{z}^{\alpha })\right]{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & & +\parallel ({\kappa }_{1}z+{\kappa }_2{z}^{\alpha })\left[{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}-{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\right]\parallel \hfill \\ & =:& I_{121}+I_{122}.\hfill \end{array}$$

Applying the error estimates in Lemma B.1. of [13] again, we derive that

$$I_{121}\le c{\tau }^3{\left(\right|z|}^4+{\left|z\right|}^{3+\alpha }{)min\{\left|z\right|}^{-1}{,\left|z\right|}^{-\alpha }\}\le c{\tau }^3{\left|z\right|}^3.$$

Since

$$\begin{array}{ccc}\hfill I_{122}& =& \parallel ({\kappa }_{1}z+{\kappa }_2{z}^{\alpha })\left[{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}-{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\right]\parallel \hfill \\ & \le & \parallel ({\kappa }_{1}z+{\kappa }_2{z}^{\alpha })\left[{\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-({\kappa }_{1}z+{\kappa }_2{z}^{\alpha })\right]\hfill \\ & & ·{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & =& |{\kappa }_{1}z+{\kappa }_2{z}^{\alpha }|\parallel {({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel ·I_{121}\le c{\tau }^3{\left|z\right|}^3,\hfill \end{array}$$

we obtain the error estimate for $I_{12}$: $I_{12}\le c{\tau }^3{\left|z\right|}^2.$

Consequently, by choosing $\delta \le 1/t_n$, we obtain

$$\begin{array}{ccc}\hfill I_1& \le & c{\tau }^3\parallel v_h\parallel \left({\int }_{\delta }^{\pi /(\tau sin\theta )}{e}^{r{t}_{n}cos\theta }{r}^2\mathrm{d}r+{\int }_{-\theta }^{\theta }{e}^{\delta t_n|cos\eta |}{\delta }^3\mathrm{d}\eta \right)\hfill \\ & \le & c{t}_n^{-3}{\tau }^3\parallel v\parallel .\hfill \end{array}$$

For the second term $I_2$, similar to the derivation of $I_2$ in Theorem 2, we can easily derive the following inequality:

$$\begin{array}{ccc}\hfill I_2& \le & c\parallel v_h\parallel {\int }_{\pi /(\tau sin\theta )}^{\infty }{e}^{r{t}_{n}cos\theta }{r}^{-1}\mathrm{d}r\hfill \\ & \le & c{\tau }^3\parallel v_h\parallel {\int }_0^{\infty }{e}^{scos\theta }{s}^2\mathrm{d}s\le c{t}_n^{-3}{\tau }^3\parallel v\parallel ,\hfill \end{array}$$

where the first inequality is valid since $\parallel \tilde{K}\left(z\right)\parallel =\parallel ({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }){({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}-I\parallel \le c$. This, together with the error estimate of $I_1$ and Theorem 1, ends the proof. □

Now, we consider the inhomogenous problem, i.e., $f\ne 0$ and $v=0$.

Theorem 4.

Let u and $U_h^n$ be the solutions of (1) and (25) with $v=0$, respectively. Then, we have

$$\begin{array}{ccc}\hfill \parallel u\left(t_n\right)-U_h^n\parallel & \le & c(h^2|ln(t_n^{\alpha }/h^2){|\parallel f\parallel }_{\mathbf{S}}+{\tau }^3{t}_n^{\alpha -3}\parallel f\left(0\right)\parallel \hfill \\ & & +{\tau }^3{t}_n^{\alpha -2}\parallel f^{\prime }\left(0\right)\parallel +{\tau }^3{t}_n^{\alpha -1}\parallel f^{″}\left(0\right)\parallel +{\tau }^3{\int }_0^{t_n}{(t_n-s)}^{\alpha -1}\parallel f^{‴}\left(s\right)\parallel \mathrm{d}s).\hfill \end{array}$$

Proof. 

It follows from (9) and (10) in Lemma 1 that

$$\parallel w_h\left(t_n\right)-W_h^n\parallel \le I_1+I_2+I_3,$$

where

$$\begin{array}{ccc}\hfill I_1& =& \parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}\left(K\left(z\right)-K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)\tau {\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_1\left(e^{-z\tau }\right)\right)f_h\left(0\right)\mathrm{d}z\parallel \hfill \\ & & +\parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}K\left(z\right)f_h\left(0\right)\mathrm{d}z\parallel ,\hfill \\ \hfill I_2& =& \parallel \frac{1}{2\pi \mathrm{i}}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}\left(z^{-1}K\left(z\right)-K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_2\left(e^{-z\tau }\right){\tau }^2\right)f_h^{\prime }\left(0\right)\mathrm{d}z\parallel \hfill \\ & & +\parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }\setminus {\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}K\left(z\right)z^{-1}{f}_h^{\prime }\left(0\right)\mathrm{d}z\parallel =:I_{21}+I_{22},\hfill \\ \hfill I_3& =& \parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }}{e}^{z{t}_n}zK\left(z\right){\widehat{\tilde{f}}}_h\left(z\right)\mathrm{d}z-\frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right)\tau {\delta }_{\tau }\left(e^{-z\tau }\right){\overline{\tilde{F}}}_h\left(e^{-z\tau }\right)\mathrm{d}z\parallel .\hfill \end{array}$$

From the derivation in Theorem 2, we can easily obtain the following error estimate of $I_1$:

$$I_1\le c{\tau }^3{t}_n^{\alpha -3}\parallel f\left(0\right)\parallel .$$

Next, we consider the two terms $I_2$ and $I_3$, respectively. Through the triangle inequality, we can obtain

$$\begin{array}{ccc}& & \parallel z^{-1}K\left(z\right)-K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\delta }_{\tau }\left(e^{-z\tau }\right){\varpi }_2\left(e^{-z\tau }\right){\tau }^2\parallel \hfill \\ & \le & \parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\left[{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2}-{\varpi }_2\left(e^{-z\tau }\right){\tau }^2\right]\parallel \hfill \\ & & +\parallel z^{-2}{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}-{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2}{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel \hfill \\ & =:& I_{21}+I_{22}.\hfill \end{array}$$

Since $\parallel {({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}{\parallel \le cmin\left\{\right|z|}^{-1}{,\left|z\right|}^{-\alpha }\}$, using the criterion (16), we have $I_{21}\le c{\tau }^3{\left|z\right|}^{1-\alpha }$. For $I_{22}$, the following inequality holds:

$$\begin{array}{ccc}\hfill I_{22}& \le & \parallel {\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2}\left[{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}-{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\right]\parallel \hfill \\ & & +\parallel (z^{-2}-{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{-2}){({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\parallel .\hfill \end{array}$$

Analogous to the derivation of the error estimates of $I_{121}$ and $I_{122}$ for $I_{12}$ in Theorem 3, we can obtain $I_{22}\le c{\tau }^3{\left|z\right|}^{1-\alpha }$. So, we choose $\delta =1/t_n$ and get

$$I_2\le c{\tau }^3\parallel f_h^{\prime }\left(0\right)\parallel \left[{\int }_{\delta }^{\pi /(\tau sin\theta )}{e}^{r{t}_{n}cos\theta }{r}^{1-\alpha }\mathrm{d}r+{\int }_{-\theta }^{\theta }{e}^{\delta t_n|cos\eta |}{\delta }^{2-\alpha }\mathrm{d}\eta \right]\le c{\tau }^3{t}_n^{\alpha -2}\parallel f^{\prime }\left(0\right)\parallel .$$

For the last term $I_3$, by using the Taylor expansion of $f_h\left(t\right)$ at $t=0$, one has $f_h\left(t\right)=f_h\left(0\right)+t{f}_h^{\prime }\left(0\right)+\frac{t^2}{2}{f}_h^{″}\left(0\right)+\frac{t^2}{2}\ast f_h^{‴}$. So ${\tilde{f}}_h\left(t\right)=\frac{t^2}{2}{f}_h^{″}\left(0\right)+\frac{t^2}{2}\ast f_h^{‴}=:f_1\left(t\right)+f_2\left(t\right)$. It suffices to analyze the error estimates for source terms of the form $f_1\left(t\right)$ and $f_2\left(t\right)$, respectively. Assume that ${\tilde{f}}_h\left(t\right)\equiv f_1\left(t\right)$. Then the term $I_3$ has the following form:

$$\begin{array}{ccc}\hfill I_3& =& \parallel \frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }}{e}^{z{t}_n}{z}^{-2}K\left(z\right)f_h^{″}\left(0\right)\mathrm{d}z\hfill \\ & & -\frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}K\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\tau }^3{\varpi }_{21}\left({\delta }_{\tau }\left(e^{-z\tau }\right)\right){\delta }_{\tau }\left(e^{-z\tau }\right)f_h^{″}\left(0\right)\mathrm{d}z\parallel ,\hfill \end{array}$$

where ${\varpi }_{21}\left(\xi \right)=\xi (1+\xi )/{(1-\xi )}^3$. By repeating the similar argument in the above error estimate for $I_2$, one can derive that $I_3\le c{\tau }^3{t}_n^{\alpha -1}\parallel f^{″}\left(0\right)\parallel$. When ${\tilde{f}}_h\left(t\right)\equiv f_2\left(t\right)$, the corresponding semidiscrete Galerkin solution in (9) can be written as follows:

$$u_h\left(t_n\right)=(E\ast {\tilde{f}}_h)\left(t_n\right)=(E\ast (\frac{t^2}{2}\ast f_h^{‴}))\left(t_n\right)=((E\ast \frac{t^2}{2})\ast f_h^{‴})\left(t_n\right),$$

where the operator $E\left(t\right)$ is given by

$$E\left(t\right)=\frac{1}{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }}{e}^{zt}{({\kappa }_{1}z+{\kappa }_2{z}^{\alpha }-\mu {\Delta }_h)}^{-1}\mathrm{d}z.$$

Furthermore, by (13), we have

$${\overline{U}}_h\left(\xi \right)={({\kappa }_1{\delta }_{\tau }\left(\xi \right)+{\kappa }_2{\delta }_{\tau }{\left(\xi \right)}^{\alpha }-\mu {\Delta }_h)}^{-1}{\overline{f}}_2\left(\xi \right)=:\overline{E}\left({\delta }_{\tau }\left(\xi \right)\right){\overline{f}}_2\left(\xi \right),$$

from which we can derive the representation of the fully discrete solution $U_h^n$ below:

$$U_h^n=\sum _{k=0}^n{E}_{\tau }^{n-k}{f}_2\left(t_k\right),$$

(27)

with $\overline{E}\left({\delta }_{\tau }\left(\xi \right)\right)={\sum }_{n=0}^{\infty }{E}_{\tau }^n{\xi }^n$. In view of (14), we get

$$E_{\tau }^n=\frac{\tau }{2\pi i}{\int }_{{\Gamma }_{\theta ,\delta }^{\tau }}{e}^{z{t}_n}{({\kappa }_1{\delta }_{\tau }\left(e^{-z\tau }\right)+{\kappa }_2{\delta }_{\tau }{\left(e^{-z\tau }\right)}^{\alpha }-\mu {\Delta }_h)}^{-1}\mathrm{d}z.$$

Denote $E_{\tau }\left(t\right)={\sum }_{n=0}^{\infty }{E}_{\tau }^n{\delta }_{t_n}\left(t\right)$ with the Dirac delta function ${\delta }_{t_n}\left(t\right)$ at $t_n$, we can rewrite (27) as

$$U_h^n=(E_{\tau }\ast f_2)\left(t_n\right)=((E_{\tau }\ast \frac{t^2}{2})\ast f_h^{‴})\left(t_n\right).$$

Since

$$\overline{E_{\tau }\ast \frac{t^2}{2}}\left(\xi \right)=\sum _{n=0}^{\infty }\sum _{j=0}^n{E}_{\tau }^{n-j}\frac{t_j^2}{2}{\xi }^n=\overline{E}\left({\delta }_{\tau }\left(\xi \right)\right){\omega }_{21}\left(\xi \right){\tau }^2,$$

one can derive that

$$\parallel ((E_{\tau }-E)\ast \frac{t^2}{2})\left(t_n\right)\parallel \le c{\tau }^3{t}_n^{\alpha -1}.$$

We show that the above inequality is valid for $\forall t\in (t_{n-1},t_n)$. Indeed, by using the Taylor series expansion of $E\left(t\right)$ at $t=t_n$, we have

$$(E\ast \frac{t^2}{2})\left(t\right)=(E\ast \frac{t^2}{2})\left(t_n\right)+(t-t_n)(E\ast t)\left(t_n\right)+\frac{{(t-t_n)}^2}{2}(E\ast 1)\left(t_n\right)+\frac{1}{2}{\int }_{t_n}^t{(t-s)}^{2}E\left(s\right)\mathrm{d}s.$$

This Taylor expansion also holds for $(E_{\tau }\ast \frac{t^2}{2})\left(t\right)$. One can easily obtain

$$\parallel ((E_{\tau }-E)\ast t)\left(t_n\right)\parallel \le c{\tau }^2{t}_n^{\alpha -1},$$

and

$$\parallel ((E_{\tau }-E)\ast 1)\left(t_n\right)\parallel \le c\tau t_n^{\alpha -1}.$$

Besides, there holds the following estimates:

$$\parallel {\int }_{t_n}^t(t-s)E\left(s\right)\mathrm{d}s\parallel \le c{\int }_{t_n}^t{(t-s)}^2{s}^{\alpha -1}\mathrm{d}s\parallel \le c{\tau }^3{t}^{\alpha -1},$$

and

$$\parallel {\int }_{t_n}^t(t-s)E_{\tau }\left(s\right)\mathrm{d}s\parallel \le c{(t_n-t)}^2\parallel E_{\tau }^n\parallel \le c{\tau }^3{t}_n^{\alpha -1}\le c{\tau }^3{t}^{\alpha -1},$$

with $t\in (t_{n-1},t_n)$. Thus, combining with Theorem 1, we complete the proof. □

5. Extending to the Temporal Fourth-Order Scheme

In this section, we consider a higher-order numerical scheme based on the idea discussed in derivation of BDF3. For sake of simplicity, we focus on the temporal fourth-order scheme. By performing the preceding argument, we introduce temporal fourth-order finite element scheme (denoted by BDF4) based on convolution quadrature generated by BDF4: Find $U_h^n,n\ge 1$, such that the following holds:

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^n-\mu {\Delta }_h{U}_h^n=F_h^n,n\ge 4,$$

(28)

with the following three corrections:

$$\begin{array}{ccc}\hfill ({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^1-\mu {\Delta }_h{U}_h^1& =& 31/24(\mu {\Delta }_h{v}_h+f_h\left(0\right))+1/6\tau f_h^{\prime }\left(0\right)+F_h^1,\hfill \\ \hfill ({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^2-\mu {\Delta }_h{U}_h^2& =& -7/6(\mu {\Delta }_h{v}_h+f_h\left(0\right))-1/12\tau f_h^{\prime }\left(0\right)+F_h^2,\hfill \end{array}$$

and

$$({\kappa }_1{\partial }_{\tau }+{\kappa }_2{\partial }_{\tau }^{\alpha })U_h^3-\mu {\Delta }_h{U}_h^3=3/8(\mu {\Delta }_h{v}_h+f_h\left(0\right))+F_h^3.$$

The corresponding error estimates for BDF4 is described by the following theorem.

Theorem 5.

Let u and $U_h^n$ be the solutions of (1) and (28), respectively. Then, we have

$$\begin{array}{ccc}\hfill \parallel u\left(t_n\right)-U_h^n\parallel & \le & c(h^2|ln(t_n^{\alpha }/h^2){|\parallel f\parallel }_{\mathbf{S}}+(h^2{t}_n^{\beta -\alpha }+{\tau }^4{t}_n^{\beta -4})\parallel \tilde{v}\parallel +{\tau }^4{t}_n^{\alpha -4}\parallel f\left(0\right)\parallel \hfill \\ & & +{\tau }^4{t}_n^{\alpha -3}\parallel f^{\prime }\left(0\right)\parallel +{\tau }^4{t}_n^{\alpha -2}\parallel f^{″}\left(0\right)\parallel +{\tau }^4{t}_n^{\alpha -1}\parallel f^{‴}\left(0\right)\parallel \hfill \\ & & +{\tau }^4{\int }_0^{t_n}{(t_n-s)}^{\alpha -1}\parallel f^{‴\prime }\left(s\right)\parallel \mathrm{d}s),\hfill \end{array}$$

where $\beta =0$ and $\tilde{v}=v$ if $v\in L^2(\Omega )$ and $\beta =\alpha$ and $\tilde{v}=\Delta v$ if $v\in D(\Delta )$.

Proof. 

Similar to the idea discussed in the error estimates of Theorems 2–4, one can easily deduce the desired error estimate results for BDF4, and thus we omit the details here. □

6. Numerical Examples

Now, we perform the numerical examples to test the error estimates of BDF3 (25) and BDF4 (28). Let $\Omega ={(0,1)}^2$ and $T=0.1$. We consider the following three types of problem data:

(a)

$v(x,y)=xy(1-x)(1-y)$ and $f=0$;

(b)

$v(x,y)={\chi }_{(0,3/4]\times (0,1)}(x,y)$ and $f=0$;

(c)

$v=0$ and $f=(1+t^{4.1}){\chi }_{(0,3/4]\times (0,1)}(x,y)$.

Here, ${\chi }_S$ is denoted as the characteristic function of the set S.

We remark that the three types of data considered above cover all cases of error theories, and therefore this is sufficient for our numerical tests. Moreover, we focus on the temporal convergence rate since the spatial one is well understood. Since the analytical solution of equation with the above data is hard to obtain, we use reference solution instead and the reference solution is computed by BDF4 with fixing $h=1/10$ and $\tau =T/4096$. The $L^2(\Omega )$ norm errors are measured by $\parallel u\left(t_n\right)-U_h^n\parallel$. The numerical results are demonstrated in

Table 1. The temporal $L^2$ norm errors for case (a) with $h=1/10$, by original schemes.

Scheme	N	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
		${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate
	16	$2.52\times {10}^{-4}$	-	2.67$\times {10}^{-4}$	-	3.52$\times {10}^{-4}$	-
BDF3 (25)	32	1.26$\times {10}^{-4}$	1.00	1.33$\times {10}^{-4}$	1.00	1.76$\times {10}^{-4}$	1.00
	64	6.27$\times {10}^{-5}$	1.01	6.64$\times {10}^{-5}$	1.00	8.78$\times {10}^{-5}$	1.00
	128	3.13$\times {10}^{-5}$	1.00	3.32$\times {10}^{-5}$	1.00	4.38$\times {10}^{-5}$	1.00
	16	2.55$\times {10}^{-4}$	-	2.70$\times {10}^{-4}$	-	3.55$\times {10}^{-4}$	-
BDF4 (28)	32	1.26$\times {10}^{-4}$	1.02	1.34$\times {10}^{-4}$	1.01	1.76$\times {10}^{-4}$	1.01
	64	6.27$\times {10}^{-5}$	1.01	6.65$\times {10}^{-5}$	1.01	8.78$\times {10}^{-5}$	1.00
	128	3.13$\times {10}^{-5}$	1.00	3.32$\times {10}^{-5}$	1.00	4.38$\times {10}^{-5}$	1.00

,

Table 2. The temporal $L^2$ norm errors for case (b) with $h=1/10$, by original schemes.

Scheme	N	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
		${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate
	16	5.25$\times {10}^{-3}$	-	5.57$\times {10}^{-3}$	-	7.39$\times {10}^{-3}$	-
BDF3 (25)	32	2.62$\times {10}^{-3}$	1.00	2.78$\times {10}^{-3}$	1.00	3.69$\times {10}^{-3}$	1.00
	64	1.30$\times {10}^{-3}$	1.01	1.39$\times {10}^{-3}$	1.00	1.84$\times {10}^{-3}$	1.00
	128	6.51$\times {10}^{-4}$	1.00	6.91$\times {10}^{-4}$	1.00	9.19$\times {10}^{-4}$	1.00
	16	5.31$\times {10}^{-3}$	-	5.62$\times {10}^{-3}$	-	7.44$\times {10}^{-3}$	-
BDF4 (28)	32	2.63$\times {10}^{-3}$	1.02	2.78$\times {10}^{-3}$	1.01	3.70$\times {10}^{-3}$	1.01
	64	1.31$\times {10}^{-3}$	1.01	1.39$\times {10}^{-3}$	1.01	1.84$\times {10}^{-3}$	1.00
	128	6.51$\times {10}^{-4}$	1.00	6.91$\times {10}^{-4}$	1.00	9.19$\times {10}^{-4}$	1.00

,

Table 3. The temporal $L^2$ norm errors for case (c) with $h=1/10$, by original schemes.

Scheme	N	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
		${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate
	16	2.59$\times {10}^{-4}$	-	2.75$\times {10}^{-4}$	-	3.63$\times {10}^{-4}$	-
BDF3 (25)	32	1.29$\times {10}^{-4}$	1.00	1.37$\times {10}^{-4}$	1.00	1.81$\times {10}^{-4}$	1.00
	64	6.45$\times {10}^{-5}$	1.01	6.83$\times {10}^{-5}$	1.00	9.04$\times {10}^{-5}$	1.00
	128	3.22$\times {10}^{-5}$	1.00	3.41$\times {10}^{-5}$	1.00	4.51$\times {10}^{-5}$	1.00
	16	2.63$\times {10}^{-4}$	-	2.77$\times {10}^{-4}$	-	3.65$\times {10}^{-4}$	-
BDF4 (28)	32	1.30$\times {10}^{-4}$	1.02	1.37$\times {10}^{-4}$	1.01	1.81$\times {10}^{-4}$	1.01
	64	6.45$\times {10}^{-5}$	1.01	6.84$\times {10}^{-5}$	1.01	9.04$\times {10}^{-5}$	1.00
	128	3.22$\times {10}^{-5}$	1.00	3.41$\times {10}^{-5}$	1.00	4.51$\times {10}^{-5}$	1.00

,

Table 4. The temporal $L^2$ norm errors for case (a) with $h=1/10$, by modified schemes.

Scheme	N	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
		${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate
	16	$3.03\times {10}^{-5}$	-	2.34$\times {10}^{-5}$	-	8.34$\times {10}^{-6}$	-
BDF2 [1]	32	7.56$\times {10}^{-6}$	2.00	5.81$\times {10}^{-6}$	2.01	2.17$\times {10}^{-6}$	1.94
	64	1.89$\times {10}^{-6}$	2.00	1.45$\times {10}^{-6}$	2.00	5.52$\times {10}^{-7}$	1.98
	128	4.71$\times {10}^{-7}$	2.00	3.62$\times {10}^{-7}$	2.00	1.39$\times {10}^{-7}$	1.99
	16	1.21$\times {10}^{-6}$	-	1.10$\times {10}^{-6}$	-	3.15$\times {10}^{-7}$	-
BDF3 (25)	32	1.67$\times {10}^{-7}$	2.86	1.45$\times {10}^{-7}$	2.93	3.18$\times {10}^{-8}$	3.31
	64	2.17$\times {10}^{-8}$	2.94	1.85$\times {10}^{-8}$	2.97	3.55$\times {10}^{-9}$	3.16
	128	2.75$\times {10}^{-9}$	2.98	2.33$\times {10}^{-9}$	2.99	4.20$\times {10}^{-10}$	3.08
	16	2.32$\times {10}^{-7}$	-	1.34$\times {10}^{-7}$	-	1.46$\times {10}^{-7}$	-
BDF4 (28)	32	3.59$\times {10}^{-9}$	6.02	2.00$\times {10}^{-9}$	6.07	4.66$\times {10}^{-9}$	4.97
	64	1.51$\times {10}^{-10}$	4.57	1.34$\times {10}^{-10}$	3.90	2.61$\times {10}^{-10}$	4.16
	128	7.50$\times {10}^{-12}$	4.33	8.78$\times {10}^{-12}$	3.93	1.55$\times {10}^{-11}$	4.08

,

Table 5. The temporal $L^2$ norm errors for case (b) with $h=1/10$, by modified schemes.

Scheme	N	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
		${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate
	16	6.60$\times {10}^{-4}$	-	5.23$\times {10}^{-4}$	-	3.32$\times {10}^{-4}$	-
BDF2 [1]	32	1.63$\times {10}^{-4}$	2.02	1.29$\times {10}^{-4}$	2.02	8.08$\times {10}^{-5}$	2.04
	64	4.06$\times {10}^{-5}$	2.01	3.20$\times {10}^{-5}$	2.01	2.00$\times {10}^{-5}$	2.02
	128	1.01$\times {10}^{-5}$	2.00	7.96$\times {10}^{-6}$	2.01	4.96$\times {10}^{-6}$	2.01
	16	5.07$\times {10}^{-5}$	-	4.43$\times {10}^{-5}$	-	4.43$\times {10}^{-5}$	-
BDF3 (25)	32	5.94$\times {10}^{-6}$	3.10	5.15$\times {10}^{-6}$	3.11	4.90$\times {10}^{-6}$	3.18
	64	7.20$\times {10}^{-7}$	3.04	6.22$\times {10}^{-7}$	3.05	5.76$\times {10}^{-7}$	3.09
	128	8.88$\times {10}^{-8}$	3.02	7.64$\times {10}^{-8}$	3.02	6.99$\times {10}^{-8}$	3.04
	16	6.38$\times {10}^{-5}$	-	5.24$\times {10}^{-5}$	-	4.62$\times {10}^{-5}$	-
BDF4 (28)	32	5.59$\times {10}^{-7}$	6.84	4.62$\times {10}^{-7}$	6.82	5.59$\times {10}^{-7}$	6.37
	64	3.15$\times {10}^{-8}$	4.15	2.61$\times {10}^{-8}$	4.15	3.11$\times {10}^{-8}$	4.17
	128	1.87$\times {10}^{-9}$	4.07	1.55$\times {10}^{-9}$	4.07	1.83$\times {10}^{-9}$	4.08

and

Table 6. The temporal $L^2$ norm errors for case (c) with $h=1/10$, by modified schemes.

Scheme	N	$\mathit{\alpha }=0.1$		$\mathit{\alpha }=0.5$		$\mathit{\alpha }=0.9$
		${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate	${\mathit{L}}^2$ Error	Rate
	16	3.14$\times {10}^{-5}$	-	2.43$\times {10}^{-5}$	-	9.82$\times {10}^{-6}$	-
BDF2 [1]	32	7.82$\times {10}^{-6}$	2.01	6.04$\times {10}^{-6}$	2.01	2.51$\times {10}^{-6}$	1.97
	64	1.95$\times {10}^{-6}$	2.00	1.50$\times {10}^{-6}$	2.01	6.34$\times {10}^{-7}$	1.99
	128	4.86$\times {10}^{-7}$	2.00	3.75$\times {10}^{-7}$	2.00	1.59$\times {10}^{-7}$	1.99
	16	1.50$\times {10}^{-6}$	-	1.33$\times {10}^{-6}$	-	6.06$\times {10}^{-7}$	-
BDF3 (25)	32	1.94$\times {10}^{-7}$	2.95	1.68$\times {10}^{-7}$	2.99	6.91$\times {10}^{-8}$	3.13
	64	2.47$\times {10}^{-8}$	2.98	2.10$\times {10}^{-8}$	3.00	8.26$\times {10}^{-9}$	3.06
	128	3.11$\times {10}^{-9}$	2.99	2.63$\times {10}^{-9}$	3.00	1.01$\times {10}^{-9}$	3.03
	16	3.84$\times {10}^{-7}$	-	2.79$\times {10}^{-7}$	-	1.73$\times {10}^{-7}$	-
BDF4 (28)	32	1.08$\times {10}^{-8}$	5.15	8.20$\times {10}^{-9}$	5.09	6.16$\times {10}^{-9}$	4.81
	64	6.04$\times {10}^{-10}$	4.16	4.77$\times {10}^{-10}$	4.10	3.52$\times {10}^{-10}$	4.13
	128	3.57$\times {10}^{-11}$	4.08	2.88$\times {10}^{-11}$	4.05	2.11$\times {10}^{-11}$	4.06

. From these tables, we can see that when no correction terms are added to BDF3 and BDF4, the temporal convergence accuracy is only first-order accuracy for both (see Table 1, Table 2 and Table 3), while after adding suitable correction terms (cf. Table 4, Table 5 and Table 6), we observe the theoretical temporal third-/fourth-order accuracy, which numerically verifies the importance of adding correction terms and the correctness of our constructed correction terms in BDF3/BDF4.

In addition, we also numerically compare the BDF3 and BDF4 with the BDF2 proposed in [1]; see Table 4, Table 5 and Table 6. The numerical results further indicate the correctness of our error estimates presented in Theorems 2–5 and reveal that our schemes greatly improve the temporal convergence accuracy of the BDF2.

7. Conclusions

In this paper, we developed an efficient temporal third/fourth-order finite element scheme for model (1) with smooth and nonsmooth data. The schemes were obtained by using the Galerkin finite element method in space and the convolution quadrature generated by the BDF3 and BDF4 in time, so this is why we refer to the resulting schemes as BDF3 and BDF4, respectively. In order to restore the desired accuracy for both smooth and nonsmooth data, we carefully chose the suitable corrections for the BDF3 and BDF4 at the starting two/three steps. Error estimates of the two modified schemes were established, with respect to the data regularity. Finally, the numerical examples confirmed the robustness and accuracy of the proposed schemes. With the BDF3/BDF4 modified at the starting two/three steps, we greatly improved the convergence results presented in [1].