Higher Order Time Stepping Methods for Subdiffusion Problems Based on Weighted and Shifted Grünwald–Letnikov Formulae with Nonsmooth Data

Abstract

Two higher order time stepping methods for solving subdiffusion problems are studied in this paper. The Caputo time fractional derivatives are approximated by using the weighted and shifted Grünwald–Letnikov formulae introduced in Tian et al. (Math Comput 84:2703–2727, 2015). After correcting a few starting steps, the proposed time stepping methods have the optimal convergence orders $O(k^2)$ and $O(k^3)$, respectively for any fixed time t for both smooth and nonsmooth data. The error estimates are proved by directly bounding the approximation errors of the kernel functions. Moreover, we also present briefly the applicabilities of our time stepping schemes to various other fractional evolution equations. Finally, some numerical examples are given to show that the numerical results are consistent with the proven theoretical results.

1 Introduction

Two higher order time stepping methods based on the weighted and shifted Grünwald-Letnikov formulae in Tian et al. [34] are introduced and analyzed for the following subdiffusion problem, with $0< \alpha <1$,

$$\begin{aligned} \, _{0}^{C} D_{t}^{\alpha } u(t) + A u(t) = f(t), \quad \text{ for } \; 0 < t \le T \quad \text{ with } \; u(0) =u_{0}, \end{aligned}$$

(1)

where $A=-\varDelta$ and $\varDelta$ denotes the Laplacian defined on a regular domain $\varOmega \subset {\mathbb {R}}^{d},d=1,2,3$ with smooth boundary $\partial \varOmega$ and $D(A)=H_{0}^{1}(\varOmega ) \cap H^{2}(\varOmega )$. Here the initial value $u_{0} \in L_{2}(\varOmega )$ and the smoothness of the source term f is described in Theorems 1, 2 in Sect. 2 and 3, respectively. The time fractional derivative $\, _{0}^{C} D_{t}^{\alpha } u(t), \, 0< \alpha <1$ is defined in the sense of Caputo, see, e.g., Diethelm [10],

$$\begin{aligned} \, _{0}^{C} D_{t}^{\alpha } u(t) = \frac{1}{\varGamma (1- \alpha )} \int _{0}^{t} (t- s)^{-\alpha } u_{t}(s) \, ds. \end{aligned}$$

More generally the operator A in (1) could be any linear, selfadjoint, positive definite operator with compact inverse, defined in $D(A) \subset L_{2}(\varOmega )$, and satisfies the following resolvent estimates, with $\pi /2< \theta _{0} < \pi$, see, e.g., Lubich et al. [23] and Thomée [33],

$$\begin{aligned} \Vert (zI + A)^{-1} \Vert \le C | z|^{-1} \quad \text{ for } \; z \in \varSigma _{\theta _{0}} = \{ z \ne 0: \, | \text{ arg } \text{ z }| < \theta _{0} \}. \end{aligned}$$

(2)

It is easy to see that for any $z \in \varSigma _{\theta }$ with $\theta \in (\pi /2, \theta _{0})$, we have $z^{\alpha } \in \varSigma _{\theta _{0}}$ since, with $0< \alpha <1$,

$$\begin{aligned} |\arg ( z^{\alpha })| = |\alpha \arg (z) |< \alpha \theta< \theta < \theta _{0}, \end{aligned}$$

which implies that, by (2), see, e.g., Jin et al. [14, (2.3)],

$$\begin{aligned} \Vert (z^{\alpha } I + A)^{-1} \Vert \le C |z |^{-\alpha }, \quad \forall \, z \in \varSigma _{\theta } = \{ z \ne 0: \, | \text{ arg } \text{ z }| < \theta \}. \end{aligned}$$

(3)

In Sects. 2 and 3, with some suitable approximation $z_{k}$ of z, we shall choose $\theta \in (\pi /2, \theta _{0})$ sufficiently close to $\pi /2$ such that $z_{k}^{\alpha } \in \varSigma _{\theta _{0}}$ which guarantees that $(z_{k}^{\alpha } I + A)^{-1}$ exists.

Recently, Meerschaert et al. [27] and Tian et al. [34] introduced the weighted and shifted Gr$\ddot{u}$nwald-Letnikov difference operators to approximate the Riemann-Liouville fractional derivative and applied such difference operators to solve space fractional partial differential equations under the assumptions that the solution is sufficiently smooth and satisfies the homogeneous boundary conditions. To our knowledge, we have not seen any works in literature to apply such weighted and shifted Grünwald-Letnikov difference operator to construct higher order time discretization schemes for solving the subdiffusion equation (1). One of the possible reasons for lacking such works may be that the solution of (1) is not sufficiently smooth and it has the singularity near $t=0$. For example, in the homogeneous case of (1) with $f=0$, one has the following stability estimate, [29], with $\Vert \cdot \Vert$ the norm in $L_{2}(\varOmega )$,

$$\begin{aligned} \Vert \, _{0}^{C}D^{\alpha }_{t} u\Vert \le Ct^{-\alpha }\Vert u_{0}\Vert , \end{aligned}$$

(4)

which shows that the $\alpha$-th order Caputo derivative of the solution of (1) becomes unbounded as $t\rightarrow 0$. Hence, the $C^{2}$-regularity assumption, generally, does not hold for the exact solution of (1). Numerical experiments indicate that the convergence orders of some numerical methods for solving (1) actually do not hold uniformly in t even for the smooth data $u_{0},$ see, e.g., Jin et al. [15], Stynes et al. [31] and Stynes [30]. Therefore, an attempt has been made in this paper to consider the higher order time discretization schemes for (1), based on the weighted and shifted Grünwald-Letnikov schemes developed in Tian et al. [34] and prove the optimal convergence orders of the proposed schemes with both smooth and nonsmooth data.

There are several approaches to improve the convergence orders of the numerical methods for solving (1), when the solution is not sufficiently smooth. One approach is to correct some weights of the numerical methods in order to capture the singularity of the solution. This idea was first introduced by Lubich et al. [23] for second order time stepping scheme applied to an evolution equation with positive memory term. After correcting some weights of the numerical methods, Lubich et al. [23] proved optimal convergence of the corrected numerical method for both smooth and nonsmooth data. Jin et al. [16] derived some higher order numerical methods in time for the problem (1), where the fractional derivatives are approximated by using the convolution quadrature generated by using the backward difference formulae. By correcting some starting steps of the numerical methods, Jin et al. [16] established that the corrected numerical methods have optimal order of convergence for any fixed time t for both smooth and nonsmooth data, see also [14, 17]. Subsequently, Yan et al. [36] corrected the starting steps of the L1 scheme for solving (1) and proved that the modified L1 scheme has the optimal convergence order $O (k^{2-\alpha }), 0< \alpha <1$. More recently, Xing and Yan [35] analyzed a numerical method for solving (1), where the Caputo fractional derivative is expressed by using the Hadamard finite-part integral which is then approximated by using the quadratic interpolation polynomials. After correcting some starting steps and some weights of the high-order numerical methods, Xing and Yan [35] derived the optimal convergence order $O(k^{3- \alpha }), 0< \alpha <1$ of the corrected numerical methods for both smooth and nonsmooth data. For the recent development of the corrections of numerical methods for (1), we refer the readers to the survey paper [12], see also [37]. For other numerical methods for solving time fractional diffusion equation, we refer to [3, 5,6,7,8, 11, 15, 18,19,20,21,22, 24, 25, 28, 32, 38,39,40], etc.

The aim of this paper is to prove that the proposed numerical methods have the optimal convergence orders $O(k^{2})$ and $O(k^{3})$, respectively, by correcting a few starting steps of the numerical methods for both smooth and nonsmooth datas. Compared to other higher order time stepping methods in the literature for solving time diffusion problem (1), the proposed methods have two advantages: (i) The weights of our numerical methods are much simpler than those obtained by approximating the fractional derivative with the quadratic interpolating polynomials, see, e.g., in Xing and Yan [35] and further, these weights have a special structure as mentioned in Tian et al. [34], which may be useful for constructing some fast algorithms and also for proving their stability and error analyses; (ii) The weights of the proposed numerical schemes are related not only to the order of the fractional derivative, but also to the shifted numbers, which imply that our methods are more related to the equation itself, see, e.g., Tian et al. [34].

The main contributions of this paper are as follows.

• 1. Based on the weighted and shifted Gr$\ddot{u}$nwald -Letnikov schemes proposed in Tian et al. [34], two new corrected higher order time discretization methods are introduced and the convergence orders are shown to be of $O(k^2)$ and $O(k^3),$ respectively for both smooth and nonsmooth data.

• 2. The error estimates of the corrected numerical methods are proved in both homogeneous and inhomogeneous cases.

• 3. With the help of Laplace transform techniques, it is shown that the error estimates are even suitable for more general elliptic operator A, which satisfies the resolvent estimate (2).

The paper is organized as follows. In Sect. 2, we consider the error estimates of the time discretization scheme for (1) with the convergence order $O(k^{2})$ for both smooth and nonsmooth data. In Sect. 3, we derive the error estimates of the time discretization scheme with the convergence order $O(k^{3})$ again for both smooth and nonsmooth data. Finally in Sect. 4, numerical examples are presented to show that numerical results are consistent with the theoretical results.

By C,  we denote a positive constant independent of discretization parameter k, but not necessarily the same at different occurrences.

2 Second Order Time Stepping Scheme

In this section, we analyze a second order time discretization scheme for approximating the solution of the problem (1). After correcting some starting steps of the scheme, the optimal order of convergence is derived for the problem with both smooth and nonsmooth data.

Based on Tian et al. [34], we shall introduce a scheme to approximate the Riemann-Liouville fractional derivative $_{0} ^{R}D^{\alpha }_{t}\phi (t)$. Let $0=t_{0}<t_{1}<\cdots <t_{N}=T$ be a time partition of [0, T] and k be the step size. We define the following numerical scheme to approximate $_{0} ^{R}D^{\alpha }_{t} \phi (t)$ at $t=t_{n},n\ge 1$

$$\begin{aligned} \, _{0}^{R}D_{t}^{\alpha } \phi (t_{n}) \approx \, _{L} D^{\alpha }_{k,p,q} \phi (t_{n}) :=\frac{\alpha -2q}{2(p-q)} B_{k,p}^{\alpha } \phi (t_{n})+\frac{2p-\alpha }{2(p-q)} B_{k,q}^{\alpha } \phi (t_{n}), \end{aligned}$$

(5)

where

$$\begin{aligned} B_{k,p}^{\alpha } \phi (t_{n}) =k^{-\alpha }\sum _{j=0}^{n+p}g_{j}^{(\alpha )}\phi (t_{n-j+p}). \end{aligned}$$

(6)

Here $g_{j}^{(\alpha )},j=0,1,2,\ldots ,$ are generated by the generating function $\delta _{1}(\zeta )=(1-\zeta )$, that is,

$$\begin{aligned} (\delta _{1} (\zeta ))^{\alpha } = (1- \zeta )^{\alpha } =\sum _{j=0}^{\infty }g_{j}^{(\alpha )}\zeta ^{n} \; \text{ with } \; g_{j}^{(\alpha )}=(-1)^{j} \left( \begin{array}{c} \alpha \\ j\\ \end{array} \right) . \end{aligned}$$

(7)

When $p=0, \, q=-1$ or $q=0, \, p=-1$, the equation (5) leads to

$$\begin{aligned} \, _{L} D^{\alpha }_{k,p,q} \phi (t_{n}) = k^{-\alpha }\sum _{j=0}^{n} w^{(\alpha )}_{n-j} \phi (t_{j}), \end{aligned}$$

(8)

where

$$\begin{aligned}&w^{(\alpha )}_{0}=\frac{\alpha +2}{2}g_{0}^{(\alpha )}, \quad w^{(\alpha )}_{j}=\frac{\alpha +2}{2}g_{j}^{(\alpha )}- \frac{\alpha }{2}g_{j-1}^{(\alpha )}, \quad j=1,2,\ldots ,n. \end{aligned}$$

(9)

In Table 1, we show the differences numerically between the weights generated by (9) and by BDF2, i.e., backward difference formula with convergence order 2, in Jin et al. [16] with $n=5$. Here $w_{j}^{(\alpha )}$ and $b_{j}^{(\alpha )}, 0< \alpha <1, \, j=0, 1, 2, \ldots , n$ denote the weights generated by (9) and by BDF2, respectively with respect to the different $0<\alpha <1$.

Table 1 Comparison of the weights generated by (9) and BDF2 with $n=5$

The corresponding generating function $\delta (\zeta )$ of the weights $w^{(\alpha )}_{j}, j=0, 1, 2, \ldots$ in (9) satisfies

$$\begin{aligned} (\delta (\zeta ))^{\alpha }&= \sum _{j=0}^{\infty } w^{(\alpha )}_{j}\zeta ^{j}= \frac{\alpha +2}{2}g_{0}^{(\alpha )}+ \left( \frac{\alpha +2}{2}g_{1}^{(\alpha )}- \frac{\alpha }{2}g_{0}^{(\alpha )} \right) \zeta \nonumber \\&\quad + \left( \frac{\alpha +2}{2}g_{2}^{(\alpha )}- \frac{\alpha }{2}g_{1}^{(\alpha )} \right) \zeta ^{2} + \left( \frac{\alpha +2}{2}g_{3}^{(\alpha )}- \frac{\alpha }{2}g_{2}^{(\alpha )} \right) \zeta ^{3} +\cdots \nonumber \\&=\frac{\alpha +2}{2} \big (g_{0}^{(\alpha )}+g_{1}^{(\alpha )}\zeta +\cdots +g_{n}^{(\alpha )}\zeta ^{n} \big ) -\frac{\alpha }{2}\zeta \big (g_{0}^{(\alpha )}+g_{1}^{(\alpha )}\zeta +\cdots +g_{n}^{(\alpha )}\zeta ^{n} \big )\nonumber \\&=\left( \frac{\alpha +2}{2}-\frac{\alpha }{2}\zeta \right) (1-\zeta )^{\alpha }. \end{aligned}$$

(10)

As in Lubich et al. [23], we denote the discrete Laplace transform of the sequence $w_{j}^{(\alpha )}, j=0, 1, 2, \ldots$ by ${\tilde{w}} (\zeta ) = \sum _{j=0}^{\infty } w_{j}^{(\alpha )} \zeta ^{j}$. We then have, by (10),

$$\begin{aligned} {\tilde{w}} (\zeta ) = \sum _{j=0}^{\infty } w_{j}^{(\alpha )} \zeta ^{j}= \delta (\zeta )^{\alpha } = \left( \frac{\alpha +2}{2}-\frac{\alpha }{2}\zeta \right) (1-\zeta )^{\alpha }. \end{aligned}$$

(11)

The following lemma gives the series expansion of ${\tilde{w}} (\zeta )$ in (11) in terms of $(1- \zeta )$, see, the similar argument used in [16, (16)].

Lemma 1

Let ${\tilde{w}} (\zeta )$ be defined by (11). Then, the following expansion holds:

$$\begin{aligned} {\tilde{w}}(\zeta )^{1/\alpha } = \left( 1+ \frac{1}{2} (1- \zeta ) + \frac{1-\alpha }{8} (1- \zeta )^2 + \cdots \right) (1- \zeta ) \quad \text{ as } \; \zeta \rightarrow 1. \end{aligned}$$

Proof

From (10), we observe that

$$\begin{aligned} {\tilde{w}}(\zeta )^{1/\alpha }&= \left( \frac{\alpha +2}{2}-\frac{\alpha }{2}\zeta \right) ^{1/\alpha } (1-\zeta ) \nonumber \\&= \left( 1+ \frac{1}{2} (1- \zeta ) + \frac{1-\alpha }{8} (1- \zeta )^2 + \cdots \right) (1- \zeta ) \quad \text{ as } \; \zeta \rightarrow 1, \end{aligned}$$

(12)

where we have used the following binomial expansion, with $\beta \in {\mathbb {R}}$,

$$\begin{aligned} (1+z)^{\beta } = 1+ \beta z + \frac{\beta (\beta +1)}{2} z^2 + \cdots , \quad \text{ as } \; z \rightarrow 0. \end{aligned}$$

(13)

This completes the proof of the Lemma 1. $\square$

We now introduce a fully discrete scheme for solving (1). Let ${\mathcal {T}}_{h}$ denote a triangulation of $\varOmega$ with h the maximal length of the sides on ${\mathcal {T}}_{h}$. Let $S_{h}\subset H^{1}_0(\varOmega )$ denote the piecewise continuous linear finite element space.

For any fixed $t \in (0, T]$, the finite element method of (1) is to find $u_{h} (t) \in S_{h}$ such that

$$\begin{aligned} \, ^{C}_{0}D^{\alpha }_{t}u_{h}(t)+A_{h}u_{h}(t)=f_{h}(t), \quad \text{ for }\, 0<t\le T \quad \text{ with }\,u_{h}(0)=u_{0h}, \end{aligned}$$

(14)

where $A_{h}:S_{h}\rightarrow S_{h}$ denotes the discrete analogue of A defined with some suitable bilinear form $a(\cdot ,\cdot )$ defined on $H^1_0(\varOmega ) \times H^1_0(\varOmega )$ associated with the operator A, by

$$\begin{aligned} (A_{h}u_{h},\chi )= a ( u_h,\chi ), \quad \forall \chi \in S_{h}, \end{aligned}$$

and $f_{h}=P_{h} f$, where $P_{h}:L_{2}(\varOmega )\rightarrow S_{h}$ denotes the $L_{2}$ project operator given by

$$\begin{aligned} (P_{h}v,\chi )=(v,\chi ), \quad \forall \chi \in S_{h}. \end{aligned}$$

Here $u_{0h} \in S_{h}$ denotes some approximation of $u_{0}\in L_{2}(\varOmega )$. When $u_{0}$ is nonsmooth, we choose $u_{0h}=P_{h}u_{0}$ and when $u_{0}$ is smooth, that is $u_{0} \in D(A)$, we may choose $u_{0h}=\mathbf{{R}}_{h}u_{0}$, where $\mathbf{{R}}_{h}:H_{0}^{1}(\varOmega )\rightarrow S_{h}$ denotes the Ritz projection or elliptic projection defined by

$$\begin{aligned} a ( \mathbf{{R}}_{h}v, \chi )= a ( v, \chi ), \quad \forall \chi \in S_{h}. \end{aligned}$$

Let $V_{h}(t)=u_{h}(t)-u_{0h}.$ Then, the equation (14) is equivalent to

$$\begin{aligned} _{0} ^{C}D^{\alpha }_{t} V_{h}(t)+A_{h} V_{h}(t) =f_{h}(t) -A_{h}u_{0h} \; \text{ with } \; V_{h}(0)=0. \end{aligned}$$

(15)

Since $\, _{0} ^{C}D^{\alpha }_{t} V_{h}(t)=\, _{0} ^{R}D^{\alpha }_{t}( V_{h}(t)-V_{h}(0)),0<\alpha <1$, now (15) is rewritten as

$$\begin{aligned} _{0} ^{R}D^{\alpha }_{t} V_{h}(t)+A_{h} V_{h}(t)=f_{h}(t) -A_{h}u_{0h} \; \text{ with } \; V_{h}(0)=0. \end{aligned}$$

(16)

An application of Taylor’s expansion as in Jin et al. [17] yields

$$\begin{aligned} f_{h}(t) = f_{h}(0) + R_{h}(t), \quad R_{h}(t) = t f^{\prime }_{h}(0) + \big (t*f^{\prime \prime }_{h} \big )(t), \end{aligned}$$

where $g*h$ denotes the convolution of g and h.

By the Laplace transform method, we obtain, with ${\hat{g}} (z)$ denoting the Laplace transform of g(t),

$$\begin{aligned} z^{\alpha } {\hat{V}} (z) + A_{h} {\hat{V}}_{h}(z)&= (f_{h}(0)- A_{h} u_{0h} ) z^{-1} + {\hat{R}}_{h}(z). \end{aligned}$$

A use of the inverse Laplace transform shows at $t=t_{n}$

$$\begin{aligned} V_{h}(t_{n})&= \frac{1}{2 \pi i} \int _{\varGamma } e^{zt} (z^{\alpha } + A_{h})^{-1} z^{-1}(f_{h}(0)- A_{h}u_{0h}) \,dz\nonumber \\&\quad +\frac{1}{2 \pi i} \int _{\varGamma } e^{zt} (z^{\alpha } + A_{h})^{-1}{\hat{R}}_{h}(z) \, dz, \end{aligned}$$

(17)

where $\varGamma$ is defined by, see, e.g., Lubich et al. [23], with some $\theta \in (\pi /2, \theta _{0})$,

$$\begin{aligned} \varGamma = \varGamma _{\theta } := \{ z \in {\mathbb {C}}: \, | \arg z| = \theta \}. \end{aligned}$$

(18)

Now we shall consider the time discretization scheme of (1). To improve the accuracy near $t=0$, we follow the approach in Lubich et al. [23] to correct the values of the first step in the time discretization.

Let $V^{n}\approx V_{h}(t_{n})$ be the approximation of $V_{h} (t_{n})$. We define the following time discretization scheme for approximating (1)

$$\begin{aligned} k^{-\alpha }\sum _{j=0}^{n}w^{(\alpha )}_{n-j}V^{j}+A_{h}V^{n}=a_{n} (f_{h}(0)-A_{h}u_{0h}) +R_{h}(t_{n}) \; \text{ with } \; V^{0}=0, \end{aligned}$$

(19)

where, with $c_{0}=1/2$,

$$\begin{aligned} a_{n}= {\left\{ \begin{array}{ll} 1+ c_{0}, \quad &{}n=1, \\ 1, \quad &{} n \ge 2. \end{array}\right. } \end{aligned}$$

Applying the discrete Laplace transform in both sides of (19), we obtain

$$\begin{aligned}&\sum _{n=1}^{\infty } \Big ( k^{-\alpha } \sum _{j=1}^{n} w^{(\alpha )}_{n-j} V^{j} \Big ) \zeta ^{n} + \sum _{n=1}^{\infty } (A_{h}V^{n}) \zeta ^{n} \nonumber \\&\quad =\left( \frac{\zeta }{1-\zeta } + c_{0} \zeta \right) (f_{h}(0)-A_{h}u_{0h} ) + \sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n}. \end{aligned}$$

(20)

It is easy to see with ${\tilde{V}} (\zeta ) = \sum _{j=0}^{\infty } V^{j} \zeta ^{j}$

$$\begin{aligned} \sum _{n=1}^{\infty }\left( \sum _{j=1}^{n} w^{(\alpha )}_{n-j}V^{j} \right) \zeta ^{n}=\left( \sum _{j=0}^{\infty }w^{(\alpha )}_{j}\zeta ^{j} \right) (V^{1} \zeta +V^{2}\zeta ^{2}+\cdots )={\tilde{w}}(\zeta ){\widetilde{V}}(\zeta ), \end{aligned}$$

which implies by (20) that

$$\begin{aligned}&{\widetilde{V}}(\zeta )=(k^{-\alpha }{\tilde{w}}(\zeta )+A_{h})^{-1}\left( \left( \frac{\zeta }{1-\zeta } + c_{0} \zeta \right) (f_{h}(0)-A_{h}u_{0h} ) +\sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n}\right) . \end{aligned}$$

Further, with ${\tilde{w}} (\zeta )$ given in (11), we set

$$\begin{aligned} z_{k} =k^{-1} {\tilde{w}}(\zeta )^{\frac{1}{\alpha }}, \end{aligned}$$

(21)

and

$$\begin{aligned} \mu (\zeta )=k z_{k}\Big ( \frac{\zeta }{1-\zeta }+c_{0} \zeta \Big ) = {\tilde{w}}(\zeta )^{\frac{1}{\alpha }}\Big ( \frac{\zeta }{1-\zeta }+c_{0} \zeta \Big ). \end{aligned}$$

(22)

By the inverse discrete Laplace transform, it follows for $n \ge 1$ and using the variable change $\zeta = e^{-zk}$ that

$$\begin{aligned} V^{n}&=\frac{1}{2 \pi i} \int _{|\zeta |=\rho }\zeta ^{-n-1}{\widetilde{V}}(\zeta )d\zeta =\frac{k}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} {\widetilde{V}}(e^{-zk}) \, dz\nonumber \\&= \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} z_{k}^{-1} \mu (\zeta ) \big ( f_{h}(0) -A_{h}u_{0h}\big ) \, dz\nonumber \\&\quad +\frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} k \big ( \sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n} \big ) \, dz, \end{aligned}$$

(23)

where, see Lubich et al. [23], with $\varGamma$ defined in (18),

$$\begin{aligned} \varGamma _{k} = \{ z \in \varGamma : \, | \mathfrak {I}z | \le \pi /k \}. \end{aligned}$$

(24)

Below, we state our main result in this section whose proof will be provided subsequently.

Theorem 1

Let $V_{h}(t_{n})$ and $V^{n}$ be defined in (17) and (23), respectively. Assume that $u_{0} \in L_{2}(\varOmega )$ and $f\in C^{1}([0,T]; L_{2}(\varOmega ))$ and $\int _{0}^{t_{n}} (t_{n}-s)^{ \alpha -1} \Vert f^{\prime \prime } (s) \Vert \, ds<\infty$ for $t_{n} \in (0,T]$. Then, there is a positive constant C, independent of k,  such that for $0< \alpha <1$

$$\begin{aligned} \Vert V_{h}(t_{n}) - V^{n} \Vert&\le C k^{2} \Big ( t_{n}^{ -2} \Vert u_{0}\Vert + t_{n}^{ \alpha -2} \Vert f(0) \Vert \\&\quad + t_{n}^{ \alpha -1} \Vert f^{\prime } (0) \Vert + \int _{0}^{t_{n}} (t_{n}-s)^{ \alpha -1} \Vert f^{\prime \prime } (s) \Vert \, ds \Big ). \end{aligned}$$

To prove Theorem 1, the following two lemmas will be useful.

Lemma 2

Let $z_{k}$ and $\mu (\zeta )$ with $\zeta = e^{-zk}$ be defined by(21) and (22), respectively. Assume that

$$\begin{aligned} K_{1}(z)=z^{-1}(z^{\alpha }+A_h)^{-1}A_h, \quad K_{2}(z)=z^{-1}(z^{\alpha }+A_h)^{-1}. \end{aligned}$$

(25)

Then with $\varGamma _{k}$ defined by (24), the following estimates hold:

$$\begin{aligned}&| \mu (e^{- z k}) -1| \le C |zk|^{2}, \quad z \in \varGamma _{k}, \end{aligned}$$

(26)

$$\begin{aligned}&C|z| \le | z_{k} | \le C|z|, \quad z \in \varGamma _{k}, \end{aligned}$$

(27)

$$\begin{aligned}&\big \Vert \mu (\zeta ) K_{1}(z_{k}) - K_{1}(z) \big \Vert \le C k^{2} |z|, \quad z \in \varGamma _{k}, \end{aligned}$$

(28)

$$\begin{aligned}&\big \Vert \mu (\zeta ) K_{2}(z_{k}) - K_{2}(z) \big \Vert \le C k^{2} |z|^{1-\alpha }, \quad z \in \varGamma _{k}. \end{aligned}$$

(29)

Proof

We first show (26). Now, zk is uniformly bounded for any $z \in \varGamma _{k}$ since, with $\theta \in (\pi /2, \theta _{0})$,

$$\begin{aligned} |zk| =|z|k=\frac{|\mathfrak {I}z|}{\sin \theta }k\le \frac{\frac{\pi }{k}}{\sin \theta }k=\frac{\pi }{\sin \theta }=\text{ const. }, \quad \text{ for } \; z \in \varGamma _{k}. \end{aligned}$$

Further we note that, by (22) and Lemma 1 with $c_{0}=1/2$,

$$\begin{aligned} \mu (\zeta ) -1&= {\tilde{w}} (z)^{1/\alpha } \Big (\frac{\zeta }{1-\zeta }+c_{0} \zeta \Big ) -1 \\&= \Big (1+ \frac{1}{2} (1- \zeta ) + \frac{1-\alpha }{8} (1- \zeta )^2 + \cdots \Big ) (1- \zeta ) \Big ( \frac{\zeta }{1-\zeta } + c_{0} \zeta \Big ) -1 \\&= \Big (1+ \frac{1}{2} (1- \zeta ) + \frac{1-\alpha }{8} (1- \zeta )^2 + \cdots \Big ) \big ( \zeta + c_{0} \zeta (1- \zeta ) \big ) -1 \\&= \Big (1+ \frac{1}{2} (1- \zeta ) + \frac{1-\alpha }{8} (1- \zeta )^2 + \cdots \Big ) \Big ( 1+ (c_{0}-1) (1- \zeta ) - c_{0} (1- \zeta )^2 \Big ) -1 \\&= O \big ( (1-\zeta )^2 \big ) \quad \text{ as } \; \zeta \rightarrow 1, \end{aligned}$$

and this implies that

$$\begin{aligned} \mu (e^{-z k}) - 1 = O \big ( (1- e^{-zk})^2 \big ) = O \big ( (zk)^2 \big ) \quad \text{ as } \; zk \rightarrow 0. \end{aligned}$$

Hence, there exists $\delta _{0}>0$ with $0< \delta _{0} \le \frac{\pi }{\sin \theta }$ such that (26) holds for $0 \le |zk| \le \delta _{0}, \, z \in \varGamma _{k}$.

For large zk with $\delta _{0} \le |zk| \le \frac{\pi }{\sin \theta }$, we now note by (22) that

$$\begin{aligned} \mu (\zeta )-1 = {\tilde{w}}(\zeta )^{\frac{1}{\alpha }} \Big ( \frac{\zeta }{1-\zeta } + c_{0} \zeta \Big ) = \Big ( \frac{\alpha +2}{2} - \frac{\alpha }{2} \Big )^{\frac{1}{\alpha }} (1- \zeta ) \Big ( \frac{\zeta }{1-\zeta } + c_{0} \zeta \Big ). \end{aligned}$$

This is continuous at any $\zeta \ne 1$, which implies that $\mu (e^{-zk}) -1$ is continuous at any $z \ne 0$. Since every continuous function is bounded on the closed and bounded domain, therefore, $\mu (e^{-zk})-1$ is bounded on $\delta _{0} \le |zk| \le \frac{\pi }{\sin \theta }, z \in \varGamma _{k}$. Then,

$$\begin{aligned} |\mu (e^{-zk})-1| \le C = { C \delta _{0}^{-2} \delta _{0}^{2} \le C \delta _{0}^{-2} |zk|^{2} \le C |zk|^{2}. } \end{aligned}$$

(30)

Hence, (26) also holds for $\delta _{0} \le |zk| \le \frac{\pi }{\sin \theta }, \, z \in \varGamma _{k}$ which completes the estimate (26).

In order to prove the estimate (27), it suffices to show $|\frac{z}{z_{k}}|$ is bounded for any $z \in \varGamma _{k}$. Now, a use of (21) yields

$$\begin{aligned} \frac{|z|}{|z_{k}|} = \frac{|z k|}{|{\tilde{w}} (e^{- z k})^{\frac{1}{\alpha }}|}, \; z \in \varGamma _{k}. \end{aligned}$$

To show $|\frac{z}{z_{k}}|$ is bounded for any $z \in \varGamma _{k}$, we consider two cases: one for the small zk and the other for the large zk.

For the small zk, observe that

$$\begin{aligned} \lim _{x \rightarrow 0} \frac{x}{{\tilde{w}} (e^{-x})^{\frac{1}{\alpha }}}&=\lim _{x \rightarrow 0} \frac{x}{(\sum _{j=0}^{\infty }w_{j}^{(\alpha )}(e^{- x})^{j})^{\frac{1}{\alpha }}} =\lim _{x \rightarrow 0} \frac{x}{(x^{\alpha }+d_{1}x^{2+\alpha }+d_{2}x^{3+\alpha }+\cdots )^{\frac{1}{\alpha }}} \\&=\lim _{x \rightarrow 0} \frac{1}{(1+d_{1}x^{2}+\cdots )^{\frac{1}{\alpha }}} =1, \end{aligned}$$

which implies that $\frac{|z|}{|z_{k}|}$ is bounded for $0\le |zk| \le \delta _{0},z\in \varGamma _{k}$ with some suitable $\delta _{0}>0, \; 0< \delta _{0} \le \frac{\pi }{\sin \theta }$.

For the large zk, we note that $\frac{|z|}{|z_{k}|}$ is continuous at any $w=zk\ne 0,z\in \varGamma _{k}$ which implies the boundedness of $\frac{|z|}{|z_{k}|}$ for large |zk| with $\delta _{0}\le |zk|\le \frac{\pi }{\sin \theta },z\in \varGamma _{k}$. Thus, we complete the estimate (27). Similarly, it is easy to show that $|\frac{z_{k}}{z}|$ is also bounded for any $z \in \varGamma _{k}$.

For (28), we first observe with $\zeta =e^{-zk}$ that

$$\begin{aligned} z_{k}-z&= \frac{ \left( \sum _{j=0}^{\infty }w_{j}^{(\alpha )}(e^{- zk})^{j} \right) ^{\frac{1}{\alpha }}-zk}{k} = \frac{(zk)(1+d_{1}(zk)^{2}+\cdots )^{\frac{1}{\alpha }}-zk}{k} \\&= \frac{(zk)(1+\frac{d_{1}}{\alpha }(zk)^{2}+\cdots )-zk}{k} = k^{-1} O \big ( (zk)^{3} \big ), \quad \text{ as } \; zk \rightarrow 0. \end{aligned}$$

This implies that there exists $0<\delta _{0} \le \frac{\pi }{\sin \theta }$ such that

$$\begin{aligned} |z_{k}-z|\le Ck^{2}|z^{3}|, \quad \text{ for } \; 0\le |zk|\le \delta _{0}, z \in \varGamma _{k}. \end{aligned}$$

For large |zk| with $\delta _{0}\le |zk|\le \frac{\pi }{\sin \theta },z\in \varGamma _{k}$, an application of (27) shows

$$\begin{aligned} |z_{k}-z|\le |z_{k}|+|z|\le C |z|\le C( k^{2}|z|^{3})\frac{1}{|zk|^{2}}\le C ( k^{2}|z|^{3})\frac{1}{\delta ^{2}_{0}}\le C( k^{2}|z|^{3}). \end{aligned}$$

Hence,

$$\begin{aligned} |z_{k}-z|\le C ( k^{2}|z|^{3}), \quad \text{ for } \; z \in \varGamma _{k}. \end{aligned}$$

(31)

Following the idea of the proof in Lubich et al. [23, (4.6)] and noting that $\Vert K_{1}^{\prime } (z) \Vert \le C |z|^{-2}$ in [23, (3.12)], we obtain by mean value theorem and using (31),

$$\begin{aligned} \Vert K _{1}(z_{k}) - K_{1}(z) \Vert&\le { C |z|^{-2} k^{2} |z|^{3}} \le C k^{2} |z|, \quad \text{ for } \; z \in \varGamma _{k}. \end{aligned}$$

(32)

As in the proof of [36, Lemma 3.12], and noting that $| K_{1}(z_{k}) | \le C |z_{k}|^{-1} \le C |z|^{-1}, \, z \in \varGamma _{k}$, we now arrive by (32) and (26) at

$$\begin{aligned} \big \Vert \mu (\zeta ) K_{1}(z_{k}) - K_{1}(z) \big \Vert&\le \big \Vert \big ( \mu (\zeta ) -1 \big ) K_{1}(z_{k}) \big \Vert + \big \Vert K_{1}(z_{k}) - K_{1}(z) \big \Vert \\&\le C | z k |^{2} |z|^{ -1} +C k^{2} |z| \le C k^{2} | z|, \; z \in \varGamma _{k}, \end{aligned}$$

which completes the proof of (28).

Finally in order to estimate (29), a use of the mean value theorem, (31) with $K_{2}(z)=z^{-1}(z^{\alpha }+A_{h})^{-1}$ and $\Vert K_{2}^{\prime } (z) \Vert \le C |z|^{-2-\alpha }$ yields

$$\begin{aligned} \Vert K_{2} (z_{k}) - K_{2}(z) \Vert&\le { C |z|^{-2-\alpha } k^{2} |z|^{3}} \le C k^{2} |z|^{1-\alpha }, \quad z \in \varGamma _{k}. \end{aligned}$$

(33)

Further, noting that $| K_{2}(z_{k}) | \le C |z_{k}|^{-1-\alpha } \le C |z|^{-1-\alpha }, \, z \in \varGamma _{k}$, we obtain, by (33) and (26),

$$\begin{aligned} \big \Vert \mu (\zeta ) K_{2}(z_{k}) - K_{2}(z) \big \Vert&\le \big \Vert \big ( \mu (\zeta ) -1 \big ) K_{2}(z_{k}) \big \Vert + \big \Vert K_{2}(z_{k}) - K_{2}(z) \big \Vert \\&\le | z k |^{2} C |z|^{ -1-\alpha } +C k^{2} |z|^{1-\alpha } \le C k^{2} | z|^{1-\alpha }, \; z \in \varGamma _{k}, \end{aligned}$$

which shows (29).

Altogether, it concludes the proof of the Lemma 2. $\square$

In the following lemma, with $z_{k}$ defined in (21), we claim that $z_{k}^{\alpha } \in \varSigma _{\theta _{0}}$ for some $\theta _{0}\in (\frac{\pi }{2},\pi )$.

Lemma 3

Let $\theta >\pi /2$ be sufficiently close to $\pi /2$. Let $z_{k}$ be defined by (21). Then we have, for some $\theta _{0} \in (\pi /2, \pi )$, with $\varGamma _{k}$ defined by (24),

$$\begin{aligned} z_{k}^{\alpha } \in \varSigma _{\theta _{0}}, \quad \text{ for } \; z \in \varGamma _{k}. \end{aligned}$$

(34)

Proof

By the definition of $z_{k}$ in (21) and the expression of ${\tilde{w}} (\zeta )$ in Lemma 1, we obtain

$$\begin{aligned} z_{k}^{\alpha }=k^{-\alpha } \Big ( \frac{\alpha +2}{2} - \frac{\alpha }{2} \zeta \Big )(1-\zeta )^{\alpha }=k^{-\alpha }\Big ( \frac{\alpha }{2}(1-\zeta )^{\alpha +1}+(1-\zeta )^{\alpha }\Big ). \end{aligned}$$

(35)

It suffices to show that both $\frac{\alpha }{2}(1-\zeta )^{\alpha +1}$ and $(1-\zeta )^{\alpha }$ lie in $\varSigma _{\theta _{0}}$ for all $z\in \varGamma _{k}$. Recall that $z\in \varGamma _{k}$ satisfies $\mathfrak {I}(zk)\in (0,\pi ]$ and $\arg (z)=\theta$ with $\theta >\frac{\pi }{2}$. Note that $z_{k}^{\alpha }$ depends on z continuously [14, proof of Lemma 3.6]. It suffices to consider the case for z with $\arg (z)=\frac{\pi }{2}$ and $\mathfrak {I}(zk)\in (0,\pi ]$. In other words, suppose that we can prove $\frac{\alpha }{2}(1-\zeta )^{\alpha +1}\in \varSigma _{\theta _{0}}$ for all z with $\arg (z)=\frac{\pi }{2}$ and $\mathfrak {I}(zk)\in (0,\pi ]$, then since $z_{k}$ depends on z continuously, there exist $\theta _{0}\in (0,\pi )$ such that $\frac{\alpha }{2}(1-\zeta )^{\alpha +1}\in \varSigma _{\theta _{0}}$. Note that $\arg (z)=\frac{\pi }{2}$ with $\mathfrak {I}(zk)\in (0,\pi ]$ implies that $\zeta =e^{-zk}=e^{-i\varphi },\varphi \in (0,\pi ]$. We next show $\frac{\alpha }{2}(1-\zeta )^{\alpha +1}\in \varSigma _{\theta _{0}}$ and $(1-\zeta )^{\alpha }\in \varSigma _{\theta _{0}}$ for $\zeta =e^{-i\varphi },\varphi \in (0,\pi ]$. Write

$$\begin{aligned} 1-\zeta =1-e^{-i\varphi }=1-\cos \varphi +i\sin \varphi , \end{aligned}$$

then, $0<\mathfrak {R}(1-\zeta )=1-\cos \varphi <1$ and $0<\mathfrak {I}(1-\zeta )=\sin \varphi <1$ for $\varphi \in (0,\pi ]$. Thus, we arrive at $\arg \frac{\alpha }{2}(1-\zeta )^{\alpha +1}\in (0,\pi )$ with $0<\alpha <1$ and $\arg (1-\zeta )^{\alpha }\in (0,\frac{\pi }{2})$ with $0<\alpha <1$. Then, there exists $\theta _{0}\in (\frac{\pi }{2},\pi )$ such that $z_{k}^{\alpha }\in \varSigma _{\theta _{0}}$ and this completes the proof. $\square$

Remark 1

In Jin et al. [14, Lemma 3.7], the authors have proved that for all $-\pi \le \theta < \pi$, there exists $\theta _{0} \in (\pi /2, \pi )$ such that $z_{k}^{\alpha } \in \varSigma _{\theta _{0}}$ for all $z \in \varSigma _{\theta }$. Actually in our analysis, we only need to show $z_{k}^{\alpha } \in \varSigma _{\theta _{0}}$ for all $z \in \varGamma _{k}$ for some $\theta >\pi /2$ sufficiently close to $\pi /2$.

Lemma 4

Let $z_{k}$ be defined as in (21), then, there is a positive constant C independent of k such that

$$\begin{aligned} \Big \Vert (z^{\alpha } +A_h)^{-1} z^{-2} - (z_{k}^{\alpha } +A_h)^{-1} \Big (k \sum _{n=1}^{\infty } t_{n}\zeta ^{n} \Big ) \Big \Vert \le C k^{2} | z |^{-\alpha }. \end{aligned}$$

Proof

Apply Lemmas 2, 3 to arrive at

$$\begin{aligned}&\Big \Vert (z^{\alpha } +A_h)^{-1} z^{-2} - (z_{k}^{\alpha } +A_h)^{-1} \Big (k \sum _{n=1}^{\infty } t_{n}\zeta ^{n} \Big ) \Big \Vert \\&\quad \le \Vert (z^{\alpha } +A_h)^{-1} z^{-2}- (z_{k}^{\alpha } +A_h)^{-1} z_{k}^{-2}\Vert + \Big \Vert (z_{k}^{\alpha } +A_h)^{-1} z_{k}^{-2} \Big ( 1- z_{k}^{2} k \sum _{n=1}^{\infty }t_{n} \zeta ^{n} \Big ) \Big \Vert . \\&\quad \le C \Big ( \big \Vert (z^{\alpha } +A_h)^{-2} z^{\alpha -3} \big \Vert + \big \Vert (z^{\alpha } +A_h)^{-1}z^{-3} \big \Vert \Big ) \Vert z_{k}-z \Vert \\&\qquad + \big \Vert (z_{k}^{\alpha } +A_h)^{-1}\big \Vert \, |z_{k}|^{-2} \, \Big |1- {\tilde{w}} (\zeta )^{\frac{2}{\alpha }} \frac{\zeta }{(1-\zeta )^{2}}\Big | \\&\quad \le Ck^{2}|z|^{-\alpha }. \end{aligned}$$

This completes the rest of the proof. $\square$

Now we turn to the proof of Theorem 1.

Proof

(Proof of Theorem 1) Subtracting (17) from (23), we arrive at

$$\begin{aligned} V_{h}(t_{n}) - V^{n} =I_{1} + I_{2}, \end{aligned}$$

where, with $K_{2}(z)$ defined by (25),

$$\begin{aligned} I_{1}&= \frac{1}{2 \pi i} \int _{ \varGamma / \varGamma _{k}} e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} z^{-1} ( f_{h}(0) -A_{h} u_{0h} ) \, dz \\&\quad + \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} \Big (K_{2}(z)- \mu (e^{-zk})K_{2}(z_{k}) \Big ) ( f_{h}(0) -A_{h} u_{0h} ) \, dz \\&= I_{11}+I_{12}, \\ I_{2}&= \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} {\hat{R}}_{h}(z) \, dz \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big (k \sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n} \Big ) \, dz = I_{21}+I_{22}. \end{aligned}$$

For $I_{1}$, apply the bound $\Vert (z^{\alpha }+A_h)^{-1}A_h\Vert \le C$, (3), (28) and (29) to obtain

$$\begin{aligned} \Vert I_{1} \Vert \le C k^{2} t_{n}^{-2} \Vert u_{0h} \Vert + Ck^{2} t_{n}^{\alpha -2} \Vert f_{h}(0) \Vert . \end{aligned}$$

For $I_{2}$, we note that

$$\begin{aligned} I_{21}&= \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}}(z^{\alpha } +A_{h})^{-1} \hat{R^{1}_{h}} (z) \, dz \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty } R^{1}_{h} (t_{n}) \zeta ^{n} \Big ) \, dz, \\ I_{22}&= \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} \hat{R^{2}_{h}} (z) \, dz \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty } R^{2}_{h} (t_{n}) \zeta ^{n} \Big ) \, dz, \end{aligned}$$

where

$$\begin{aligned} R_{h}(t) = t f^{ \prime }_{h} (0) + \big ( t*f^{ \prime \prime } _{h} \big )(t) =: R^{1}_{h}(t) + R^{2}_{h}(t). \end{aligned}$$

For $I_{21}$, we easily bound it as

$$\begin{aligned} \Vert I_{21}\Vert&=\Big \Vert \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} z^{-2} \, dz f^{\prime }_{h}(0) \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \left( k \sum _{n=1}^{\infty } R^{1}_{h} (t_{n}) \zeta ^{n} \right) \, dz \Big \Vert \\&=\Big \Vert \frac{1}{2 \pi i} \int _{ \varGamma /\varGamma _{k} } e^{z t_{n}}(z^{\alpha } +A_{h})^{-1} z^{-2} \, dz f^{\prime }_{h} (0) \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} \Big ( (z^{\alpha } +A_{h})^{-1} z^{-2}-(z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty } t_n \zeta ^{n} \Big ) \Big ) \, dz f^{\prime }_{h} (0) \Big \Vert . \end{aligned}$$

An application of Lemma 4 yields

$$\begin{aligned} \Vert I_{21} \Vert \le C k^{2} t_{n}^{ \alpha -1} \Vert f^{\prime }_{h} (0) \Vert . \end{aligned}$$

(36)

For $I_{22}$, following the arguments as in Jin et al. [14, 17], we arrive at

$$\begin{aligned} \Vert I_{22}\Vert \le C k^{2} \int _{0}^{t_{n}} (t_{n}-s)^{ \alpha -1 } \Vert f^{\prime \prime }_{h} (s) \Vert \, ds. \end{aligned}$$

Together these estimates complete the proof of Theorem 1. $\square$

3 Third Order Time Discretization Method

In this section, we introduce a third order time discretization scheme for solving (1) based on the weighted and shifted Grünwald-Letnikov difference operator introduced in Tian et al. [34].

Let us define the following weighted and shifted Grünwald-Letnikov difference operator $\, _{L} D^{\alpha }_{k,p,q,r}$ to approximate the Riemann-Liouville fractional derivative operator $\, _{0}^{R}D_{t}^{\alpha }$ by

$$\begin{aligned} \, _{0}^{R}D_{t}^{\alpha } \phi (t) \approx \, _{L} D^{\alpha }_{k,p,q,r} \phi (t) :=\lambda _{1}B_{k,p}^{\alpha }u(t)+\lambda _{2}B_{k,q}^{\alpha } \phi (t)+\lambda _{3}B_{k,r}^{\alpha }\phi (t), \end{aligned}$$

(37)

where p, q, r are integers and mutually nonequal, and $B_{k,p}^{\alpha } \phi (t_{n})$ are defined by (6) and

$$\begin{aligned}&\lambda _{1}=\frac{12qr-(6q+6r+1)\alpha +3\alpha ^{2}}{12(qr-pq-pr+p^{2})}, \\&\lambda _{2}=\frac{12pr-(6p+6r+1)\alpha +3\alpha ^{2}}{12(pr-pq-qr+q^{2})}, \\&\lambda _{3}=\frac{12pq-(6p+6q+1)\alpha +3\alpha ^{2}}{12(pq-pr-qr+r^{2})}. \end{aligned}$$

When $p=0, \, q=-1, \, r=-2$, we obtain

$$\begin{aligned}&\lambda _{1}=\frac{24+17\alpha +3\alpha ^{2}}{24}, \quad \lambda _{2}=\frac{-22\alpha -6\alpha ^{2}}{24}, \quad \lambda _{3}=\frac{5\alpha +3\alpha ^{2}}{24}, \end{aligned}$$

(38)

and we then arrive for $n\ge 2$ at

$$\begin{aligned} _{0} D^{\alpha }_{t} \phi (t_n) \approx \, _{L} D^{\alpha }_{k,p,q,r} \phi (t_{n}):= \lambda _{1}B_{k,0}^{\alpha } \phi (t_{n})+\lambda _{2}B_{k,-1}^{\alpha } \phi (t_{n})+\lambda _{3}B_{k,-2}^{\alpha } \phi (t_{n}). \end{aligned}$$

(39)

Thus, we obtain

$$\begin{aligned} \, _{L} D^{\alpha }_{k,p,q,r} \phi (t_{n})&= \lambda _{1}k^{-\alpha }\sum _{j=0}^{n}g_{j}^{(\alpha )} \phi (t_{n-j}) + \lambda _{2}k^{-\alpha }\sum _{j=0}^{n-1}g_{j}^{(\alpha )} \phi (t_{n-j-1}) \nonumber \\&\quad +\lambda _{3}k^{-\alpha }\sum _{j=0}^{n-2}g_{j}^{(\alpha )} \phi (t_{n-j-2}) \nonumber \\&= k^{-\alpha } \lambda _{1} g_{0}^{(\alpha )} \phi (t_{n})+ k^{-\alpha } (\lambda _{1} g_{1}^{(\alpha )}+\lambda _{2} g_{0}^{(\alpha )}) \phi (t_{n-1})\nonumber \\&\quad +k^{-\alpha }(\lambda _{1}g_{2}^{(\alpha )}+\lambda _{2}g_{1}^{(\alpha )}+\lambda _{3}g_{0}^{(\alpha )}) \phi (t_{n-2})\nonumber \\&\quad +k^{-\alpha }(\lambda _{1}g_{3}^{(\alpha )}+\lambda _{2}g_{2}^{(\alpha )}+\lambda _{3}g_{1}^{(\alpha )}) \phi (t_{n-3})+\cdots \nonumber \\&= k^{-\alpha }\sum _{j=0}^{n}w^{(\alpha )}_{n-j} \phi (t_{j}), \end{aligned}$$

(40)

where

$$\begin{aligned} w_{j}^{(\alpha )}= {\left\{ \begin{array}{ll} \lambda _{1}g_{0}^{(\alpha )}, &{}\quad j=0, \\ \lambda _{1}g_{1}^{(\alpha )}+\lambda _{2}g_{0}^{(\alpha )}, &{}\quad j=1, \\ \lambda _{1}g_{j}^{(\alpha )}+\lambda _{2}g_{j-1}^{(\alpha )}+\lambda _{3}g_{j-2}^{(\alpha )}, &{}\quad j=2,3,\ldots ,n. \end{array}\right. } \end{aligned}$$

The discrete Laplace transform of $w_{j}^{(\alpha )}, j=0, 1, 2, \ldots$ is given by

$$\begin{aligned} {\tilde{w}} (\zeta )&=\sum _{j=0}^{\infty } w^{(\alpha )}_{j}\zeta ^{j} =\lambda _{1}g_{0}^{(\alpha )}+\left( \lambda _{1}g_{1}^{(\alpha )}+\lambda _{2}g_{0}^{(\alpha )}\right) \zeta \nonumber \\&\quad +\left( \lambda _{1}g_{2}^{(\alpha )}+\lambda _{2}g_{1}^{(\alpha )} +\lambda _{3}g_{0}^{(\alpha )}\right) \zeta ^{2} +\left( \lambda _{1}g^{(\alpha )}_{3}+\lambda _{2}g^{(\alpha )}_{2} +\lambda _{3}g_{0}^{(\alpha )}\right) \zeta ^{3}+\cdots \nonumber \\&=(\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2})(1-\zeta )^{\alpha }. \end{aligned}$$

(41)

Below, following the idea in Lemma 1, we consider the series expansion of the function ${\tilde{w}}(\zeta )$ in (41).

Lemma 5

Let ${\tilde{w}} (\zeta )$ be defined by (41). Then,

$$\begin{aligned} {\tilde{w}}(\zeta )^{1/\alpha } = \Big (1- \frac{\lambda _{2}+2 \lambda _{3}}{\alpha } (1- \zeta ) + \cdots \Big ) (1- \zeta ) \quad \text{ as } \; \zeta \rightarrow 1. \end{aligned}$$

(42)

Proof

From (41), the binomial expansion (13) and $\lambda _{1} + \lambda _{2} + \lambda _{3}=1$, it now follows that

$$\begin{aligned} {\tilde{w}}(\zeta )^{1/\alpha }&= \big ( \lambda _{1} + \lambda _{2} \zeta + \lambda _{3} \zeta ^2 \big )^{1/\alpha } (1-\zeta ) \nonumber \\&= \Big ( 1- (\lambda _{2} + 2 \lambda _{3} ) (1- \zeta ) + \lambda _{3} (1- \zeta )^2 \Big )^{1/\alpha } (1- \zeta ) \nonumber \\&=\Big \{ 1+ \frac{1}{\alpha } \big ( - (\lambda _{2} + 2 \lambda _{3}) (1- \zeta ) + \lambda _{3} (1- \zeta )^2 \big ) \nonumber \\&\quad + \frac{\frac{1}{\alpha } \left( \frac{1}{\alpha }-1\right) }{2} \Big ( - (\lambda _{2}+ 2 \lambda _{3}) (1- \zeta ) + \lambda _{3} (1- \zeta )^2 \Big )^2 + \cdots \Big \} (1- \zeta ) \nonumber \\&= \Big (1- \frac{\lambda _{2}+2 \lambda _{3}}{\alpha } (1- \zeta ) + \cdots \Big ) (1- \zeta ) \quad \text{ as } \; \zeta \rightarrow 1. \end{aligned}$$

(43)

This concludes the rest of the proof. $\square$

Next, we turn to the solution of (16) when $f_{h}(t)$ is written as

$$\begin{aligned} f_{h}(t)=f_{h}(0)+f'_{h}(0)t+R_{h}(t), \quad R_{h}(t)=\frac{t^{2}}{2!}f''_{h}(0)+(\frac{t^{2}}{2!}*f'''_{h} )(t), \end{aligned}$$

(44)

where $g*h$ denotes the convolution of g and h.

An application of the Laplace transform to (16) with respect to the time variable t yields

$$\begin{aligned} \hat{V_{h}}(t)=(z^{\alpha }+A_{h})^{-1} \big ( (f_{h}(0)-A_{h} u_{0h}) z^{-1} +f'_{h}(0)z^{-2}+{\hat{R}}_{h} (z) \big ). \end{aligned}$$

(45)

By the inverse Laplace transform, the solution of (16) takes the following form at $t=t_{n}$,

$$\begin{aligned} V_{h}(t_{n})=&\frac{1}{2\pi i} \int _{\varGamma }e^{zt}(z^{\alpha }+A_{h})^{-1} z^{-1} (f_{h}(0)-A_{h} u_{0h}) \,dz \nonumber \\&\quad +\frac{1}{2\pi i} \int _{\varGamma }e^{zt}\Big ((z^{\alpha }+A_{h})^{-1}z^{-2}f'_{h}(0)+(z^{\alpha }+A_{h})^{-1} {\hat{R}}_{h}(z)\Big )\,dz. \end{aligned}$$

(46)

Let $V^{n}\approx V_{h}(t_{n}),n=0,1,2,\ldots ,N$ denote the approximate solution of the following time discretization scheme for solving (16), with $V^{0}=0$,

$$\begin{aligned} k^{-\alpha } \sum _{j=0}^{n}w_{n-j}^{(\alpha )}V^{j}+A_{h}V^{n} =&f_{h}(0) -A_{h}u_{0h} +f^{\prime }_{h}(0)t_{n}+R_{h}(t_{n})\nonumber \\&\quad +a_{1} (f_{h}(0)-A_{h}u_{0h}) +b_{1}kf^{\prime }_{h}(0), \quad n=1, \end{aligned}$$

(47)

$$\begin{aligned} k^{-\alpha } \sum _{j=0}^{n} w_{n-j}^{(\alpha )}V^{j}+A_{h}V^{n} =&f_{h}(0)-A_{h}u_{0h} +f^{\prime }_{h}(0)t_{n}+R_{h}(t_{n}) \nonumber \\&\quad +a_{2} (f_{h}(0)-A_{h}u_{0h}) +b_{2}kf^{\prime }_{h}(0), \quad n=2, \end{aligned}$$

(48)

$$\begin{aligned} k^{-\alpha } \sum _{j=0}^{n}w_{n-j}^{(\alpha )}V^{j}+A_{h}V^{n}&= f_{h}(0)-A_{h}u_{0h} +f^{\prime }_{h}(0)t_{n}+R_{h}(t_{n}), \quad n \ge 3, \end{aligned}$$

(49)

where $w_{j}^{(\alpha )},j=0,1,2,\ldots$ are defined by (40) and the coefficients $a_{1}, a_{2}, b_{1}, b_{2}$ satisfy

$$\begin{aligned} a_{1}=\frac{11}{12}, \quad a_{2}=-\frac{5}{12}, \end{aligned}$$

(50)

and

$$\begin{aligned} b_{1}+b_{2}=\frac{1}{12}. \end{aligned}$$

(51)

Now we come to the following main theorem in this section

Theorem 2

Let $V_{h}(t_{n})$ and $V^{n}$ be the solutions of (16) and (47)–(49), respectively. Assume that $u_{0} \in L_{2} (\varOmega )$ and $f\in C^{2}([0,T]; L_{2}(\varOmega ))$ and $\int _{0}^{t_{n}} (t_{n}-s)^{ \alpha -1} \Vert f^{\prime \prime \prime } (s) \Vert \, ds<\infty$ for $t_{n} \in (0,T]$. Let $u_{0h}=P_{h}u_{0}$, then there exists a positive constant C independent of k such that

$$\begin{aligned} \Vert V_{h}(t_{n}) - V^{n} \Vert&\le C k^{3} \Big ( t_{n}^{-3} \Vert u_{0}\Vert + t_{n}^{ \alpha -3} \Vert f(0) \Vert + t_{n}^{ \alpha -2} \Vert f^{\prime } (0) \Vert + t_{n}^{ \alpha -1} \Vert f^{\prime \prime } (0) \Vert \\&\quad + \int _{0}^{t_{n}} (t_{n}-s)^{ \alpha -1} \Vert f^{\prime \prime \prime } (s) \Vert \, ds \Big ). \end{aligned}$$

To prove Theorem 2, we need the following lemmas.

Lemma 6

With ${\tilde{w}} (\zeta )$ given by (41), let $z_{k}$ and $\mu (\zeta )$ with $\zeta = e^{-zk}$ be defined by

$$\begin{aligned} z_{k}=k^{-1} {\tilde{w}}(\zeta )^{\frac{1}{\alpha }}, \end{aligned}$$

(52)

and

$$\begin{aligned} \mu (\zeta )=k z_{k}\Big (\frac{\zeta }{1-\zeta }+\sum _{j=1}^{2}a_{j} \zeta ^{j}\Big ) = {\tilde{w}}(\zeta )^{\frac{1}{\alpha }} \Big (\frac{\zeta }{1-\zeta }+\sum _{j=1}^{2}a_{j} \zeta ^{j}\Big ). \end{aligned}$$

(53)

Further, let $K_{1}(z)$ and $K_{2} (z)$ be given by (25). Then, with $\varGamma _{k}$ as in (24), the following estimates hold:

$$\begin{aligned}&| \mu (e^{- z k}) -1| \le C |zk|^{3}, \quad z \in \varGamma _{k}, \end{aligned}$$

(54)

$$\begin{aligned}&C|z| \le | z_{k} | \le C|z|, \quad z \in \varGamma _{k}, \end{aligned}$$

(55)

$$\begin{aligned}&\big \Vert \mu (\zeta ) K_{1}(z_{k}) - K_{1}(z) \big \Vert \le C k^{3} |z|^2, \quad z \in \varGamma _{k}, \end{aligned}$$

(56)

$$\begin{aligned}&\big \Vert \mu (\zeta ) K_{2}(z_{k}) - K_{2}(z) \big \Vert \le C k^{3} |z|^{2-\alpha }, \quad z \in \varGamma _{k}. \end{aligned}$$

(57)

Proof

The proof is similar to the proof of Lemma 2. For each inequality, again two cases such as the small zk and the large zk are considered.

For (54), a use of Lemma 5 yields, with $a_{1}= 11/12, \, a_{2}=-5/12$ by (50),

$$\begin{aligned} \mu (\zeta ) -1&= {\tilde{w}} (z)^{1/\alpha } \Big (\frac{\zeta }{1-\zeta }+ a_{1} \zeta + a_{2} \zeta ^2 \Big ) -1 \\&= \Big (1- \frac{\lambda _{2} + 2 \lambda _{3}}{\alpha } (1- \zeta ) + \cdots \Big ) (1- \zeta ) \Big ( \frac{\zeta }{1-\zeta } + a_{1} \zeta + a_{2} \zeta ^2 \Big ) -1 \\&= \Big (1- \frac{\lambda _{2} + 2 \lambda _{3}}{\alpha } (1- \zeta ) + \cdots \Big ) \Big ( \zeta + a_{1} \zeta (1- \zeta ) + a_{2} \zeta ^2 (1- \zeta ) \Big ) -1 \\&= \Big (1- \frac{\lambda _{2} + 2 \lambda _{3}}{\alpha } (1- \zeta ) + \cdots \Big ) \\&\quad \cdot \Big ( 1+ (a_{1}+a_{2}-1) (1-\zeta ) -(a_{1}+2 a_{2}) (1-\zeta )^2 + a_{2} (1- \zeta )^3 \Big ) -1 \\&= O \big ( (1-\zeta )^3 \big ) \quad \text{ as } \; \zeta \rightarrow 1, \end{aligned}$$

and then this implies that

$$\begin{aligned} \mu (e^{-z k}) - 1 = O \big ( (1- e^{-zk})^3 \big ) = O \big ( (zk)^3 \big ) \quad \text{ as } \; zk \rightarrow 0. \end{aligned}$$

Hence (54) holds for small zk with $0 \le |zk| \le \delta _{0}$ with some positive $\delta _{0} >0, \, 0 \le \delta _{0} < \frac{\pi }{\sin (\theta )}, \, \theta \in (\pi /2, \theta _{0})$. Note that $\mu (\zeta )$ is continuous at any point except $\zeta \ne 1$, which implies that, following the argument as in (30), $\mu (e^{-zk})-1$ is bounded for large zk with $\delta _{0} \le |zk| \le \frac{\pi }{\sin (\theta )}$. Hence (54) follows.

We next estimate (55). With ${\tilde{w}} (\zeta )$ defined by (41), we arrive at

$$\begin{aligned} \frac{|z|}{|z_{k}|} = \frac{|z k|}{| {\tilde{w}}(e^{- z k})^{\frac{1}{\alpha }}|}. \end{aligned}$$

Hence,

$$\begin{aligned} \lim _{x \rightarrow 0} \frac{x}{{\tilde{w}}(e^{- x})^{\frac{1}{\alpha }}}&=\lim _{x \rightarrow 0} \frac{x}{(\sum _{j=0}^{\infty }w^{(\alpha )}_{j}(e^{- x})^{j})^{\frac{1}{\alpha }}} \\&=\lim _{x \rightarrow 0} \frac{x}{(x^{\alpha }+d_{1}x^{3+\alpha }+d_{2}x^{4+\alpha }+\cdots )^{\frac{1}{\alpha }}} \\&=\lim _{x \rightarrow 0} \frac{1}{(1+d_{1}x^{3}+\cdots )^{\frac{1}{\alpha }}} =1, \end{aligned}$$

which implies that $\frac{|z|}{|z_{k}|}$ is bounded for small |zk| with $0\le |zk| \le \delta _{0}, z\in \varGamma _{k}$ for some suitable $0< \delta _{0} \le \frac{\pi }{\sin \theta }$. Note also that $\frac{|z|}{|z_{k}|}$ is continuous at any $w=zk\ne 0, z\in \varGamma _{k}$, which implies that $\frac{|z|}{|z_{k}|}$ is also bounded for large |zk| with $\delta _{0}\le |zk|\le \frac{\pi }{\sin \theta },z\in \varGamma _{k}$. Hence, we proved the boundedness of $\frac{|z|}{|z_{k}|}$ for any $z \in \varGamma _{k}$. Similarly, we may show that $\frac{|z_{k}|}{|z|}$ is also bounded for any $z \in \varGamma _{k}$. Thus, we derive the estimate (55).

To estimate (56), we observe that

$$\begin{aligned} z_{k}-z&= \frac{\big ( \sum _{j=0}^{\infty }w_{j}(e^{- zk})^{j} \big )^{\frac{1}{\alpha }}-zk}{k} = \frac{(zk)(1+d_{1}(zk)^{3}+\cdots )^{\frac{1}{\alpha }}-zk}{k} \\&= \frac{(zk)(1+\frac{d_{1}}{\alpha }(zk)^{3}+\cdots )-zk}{k} = k^{-1} O \big ( (zk)^{4} \big ), \quad \text{ as } \; zk \rightarrow 0, \end{aligned}$$

which implies that there exists $0<\delta _{0} \le \frac{\pi }{\sin \theta }$ such that

$$\begin{aligned} |z_{k}-z|\le Ck^{3}|z^{4}|\quad \text{ for } \; 0\le |zk|\le \delta _{0},z \in \varGamma _{k}. \end{aligned}$$

For large |zk|, with $\delta _{0}\le |zk|\le \frac{\pi }{\sin \theta }, \, z\in \varGamma _{k}$, we have, by (55),

$$\begin{aligned} |z_{k}-z|\le |z_{k}|+|z|\le C |z|\le C( k^{3}|z|^{4})\frac{1}{|zk|^{3}}\le C ( k^{3}|z|^{4})\frac{1}{\delta ^{3}_{0}}\le C( k^{3}|z|^{4}). \end{aligned}$$

Thus, we obtain

$$\begin{aligned} |z_{k}-z|\le C ( k^{3}|z|^{4}), \; z\in \varGamma _{k}. \end{aligned}$$

(58)

Then an application of the mean-value theorem with (58) and $\Vert K_{1}^{\prime } (z) \Vert \le C |z|^{-2}$ shows

$$\begin{aligned} \Vert K_{1}(z_{k}) - K_{1}(z) \Vert&\le { C |z|^{-2} k^{3} |z|^{4}} \le C k^{3} |z|^{2}, \; z\in \varGamma _{k}. \end{aligned}$$

(59)

Following the same line of proof of (28), we arrive from (54), (59), and $| K_{1}(z_{k}) | \le C |z_{k}|^{-1} \le C |z|^{-1}, z \in \varGamma _{k}$ at

$$\begin{aligned} \big \Vert \mu (\zeta ) K_{1}(z_{k}) - K_{1}(z) \big \Vert&\le \big \Vert \big ( \mu (\zeta ) -1 \big ) K_{1}(z_{k}) \big \Vert + \big \Vert K_{1}(z_{k}) - K_{1}(z) \big \Vert \\&\le | z k |^{3} C |z|^{ -1} +C k^{3} |z|^{2} \le C k^{3} | z|^{2}, \; z\in \varGamma _{k}, \end{aligned}$$

which shows (56).

Finally in order to show (57), we note that $K_{2}(z)=z^{-1}(z^{\alpha }+A_{h})^{-1}$ and $\Vert K_{2}^{\prime } (z) \Vert \le C |z|^{-2-\alpha }$. Then, by mean-value theorem and (58), we obtain

$$\begin{aligned} \Vert K_{2} (z_{k}) - K_{2}(z) \Vert =\Vert K'_{2}(z)\Vert |z_{k}-z| \le { C |z|^{-2-\alpha } k^{3} |z|^{4}} \le C k^{3} |z|^{2-\alpha }, \; z\in \varGamma _{k}. \end{aligned}$$

Following the same arguments as the proof of (29), a use of $|K_{2}(z_{k}) | \le C |z_{k}|^{-1-\alpha } \le C |z|^{-1-\alpha }$ yields

$$\begin{aligned} \big \Vert \mu (\zeta ) K_{2}(z_{k}) - K_{2}(z) \big \Vert&\le \big \Vert \big ( \mu (\zeta ) -1 \big ) K_{2}(z_{k}) \big \Vert + \big \Vert K_{2}(z_{k}) - K_{2}(z) \big \Vert \\&\le | z k |^{3} C |z|^{ -1-\alpha } +C k^{3} |z|^{2-\alpha } \le C k^{3} | z|^{2-\alpha }, \; z\in \varGamma _{k}. \end{aligned}$$

Hence, we prove (57).

Together with these estimates, we complete the proof of Lemma 6. $\square$

In the following lemma, with $z_{k}$ defined by (52), we will show that $z_{k}^{\alpha } \in \varSigma _{\theta _{0}}$ for some $\theta _{0} \in (\pi /2, \pi )$.

Lemma 7

Let $\theta >\pi /2$ be sufficiently close to $\pi /2$. Let $z_{k}$ be defined by (52). Then, there exists $\theta _{0} \in (\pi /2, \pi )$ such that

$$\begin{aligned} z_{k}^{\alpha } \in \varSigma _{\theta _{0}} \quad \text{ for } \text{ all } \; z \in \varGamma _{k}, \end{aligned}$$

(60)

where $\varGamma _{k}$ is defined by (24).

Proof

By the definition of $z_{k}$ in (52) and the expression of ${\tilde{w}}(\zeta )$ in Lemma 5, we rewrite

$$\begin{aligned} z_{k}^{\alpha }=k^{-\alpha }(\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2})(1-\zeta )^{\alpha }, \end{aligned}$$

(61)

where $\lambda _{1}=\frac{24+17\alpha +3\alpha ^{2}}{24}$, $\lambda _{2}=\frac{-22\alpha -6\alpha ^{2}}{24}$, $\lambda _{3}=\frac{5\alpha +3\alpha ^{2}}{24}$. It suffices to show both $\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2}$ and $(1-\zeta )^{\alpha }$ lie in $\varSigma _{\theta _{0}}$ for all $z\in \varGamma _{k}$. Recall that $z\in \varGamma _{k}$ satisfies $\mathfrak {I}(zk)\in (0,\pi ]$ and $\arg (z)=\theta$ with $\theta >\frac{\pi }{2}.$ Note that $z_{k}^{\alpha }$ depends on z continuously [14, proof of Lemma 3.6], and $\arg (z)=\frac{\pi }{2}$ with $\mathfrak {I}(zk)\in (0,\pi ]$ implies that $\zeta =e^{-zk}=e^{-i\varphi },\varphi \in (0,\pi ]$. Hence, we only need to show $(\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2})\in \varSigma _{\theta _{0}}$ for $\zeta =e^{-i\varphi },\varphi \in (0,\pi ].$ Observe that

$$\begin{aligned} \lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2}&=\lambda _{1}+\lambda _{2}e^{-i\varphi }+\lambda _{3}e^{-i2\varphi } \\&=(\lambda _{1}+\lambda _{2}\cos \varphi +\lambda _{3}\cos 2\varphi ) -i(\lambda _{2}\sin \varphi +\lambda _{3}\sin 2\varphi ), \end{aligned}$$

and therefore, we arrive at

$$\begin{aligned}&\mathfrak {R}(\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2})\\&\quad =\lambda _{1}+\lambda _{2}\cos \varphi +2\lambda _{3}\cos ^{2}\varphi -\lambda _{3} \\&\quad =\frac{5\alpha +3\alpha ^{2}}{12}\cos ^{2}\varphi -\frac{22\alpha +6\alpha ^{2}}{24}\cos \varphi +1+\frac{\alpha }{2}\\&\quad =\frac{\alpha ^{2}}{4} \left( \cos \varphi -\frac{1}{2}\right) ^{2} +\frac{5\alpha }{12}(\cos \varphi -1)^{2}+\frac{\alpha }{12} (1-\cos \varphi )+\left( 1-\frac{\alpha ^{2}}{16}\right) . \end{aligned}$$

Then, $\mathfrak {R}(\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2})>0$ for $\varphi \in (0,\pi )$ and $0<\alpha <1$.

Similarly,

$$\begin{aligned} \mathfrak {I}(\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2})&=\frac{11\alpha +3\alpha ^{2}}{12}\sin \varphi -\frac{5\alpha +3\alpha ^{2}}{12}\sin \varphi \cos \varphi \nonumber \\&=\frac{\alpha ^{2}}{4}\sin \varphi (1-\cos \varphi )+\frac{5\alpha }{12}\sin \varphi (1-\cos \varphi )+\frac{\alpha }{2}\sin \varphi >0. \end{aligned}$$

(62)

Thus, $\arg (\lambda _{1}+\lambda _{2}\zeta +\lambda _{3}\zeta ^{2}) \in (0,\frac{\pi }{2})$ for $\varphi \in (0,\pi )$ and $\arg (1-\zeta )^{\alpha } \in (\frac{\pi }{2},\pi )$ with $0<\alpha <1$ in Lemma 3. Then, we can infer the existence of $\theta _{0}\in (0,\frac{\pi }{2})$ such that $z_{k}^{\alpha }\in \varSigma _{\theta _{0}}$ and this concludes the rest of the proof. $\square$

Lemma 8

Let $z_{k}$ be defined as in (52). Let $b_{1},b_{2}$ be defined as in (51). Then, the following estimate with $\zeta = e^{-zk}$, holds:

$$\begin{aligned} \Big \Vert z_{k}^{-2} - \Big ( \sum _{n=1}^{\infty } n \zeta ^{n} +\sum _{j=1}^{2} b_{j} \zeta ^{j} \Big ) k^{2} \Big \Vert \le C k^{3} | z |. \end{aligned}$$

Proof

Let $\zeta =e^{-zk}$ and $x=zk$, then by Lemma 5, we have

$$\begin{aligned}&z_{k}^{2} \Big ( \sum _{n=1}^{\infty } n \zeta ^{n} +\sum _{j=1}^{2} b_{j} \zeta ^{j} \Big ) k^{2}-1 \\&\quad = \big ({\tilde{w}}(e^{-x})\big )^{\frac{2}{\alpha }} \Big ( \sum _{n=1}^{\infty } n (e^{-x})^{n} +\sum _{j=1}^{2} b_{j} (e^{-x})^{j} \Big ) -1 \\&\quad =\Big (\frac{e^{-x}}{(1-e^{-x})^{2}}+\sum _{j=1}^{2}b_{j}(e^{-x})^{j}\Big )\big ({\tilde{w}}(e^{-x})\big )^{\frac{2}{\alpha }}-1\\&\quad = \Big ( \frac{1}{x^{2}}(1- 1/12 x^{2}+0x^{3}+\cdots ) +b_{1}-b_{1}x+ 1/2 b_{1}x^{2}-\frac{1}{3!}b_{1}x^{3}\\&\qquad +b_{2}-2b_{2}x+2b_{2}x^{2}-\frac{8}{3!}b_{2}x^{3} \Big ) (x^{\alpha }+d_{1}x^{\alpha +3}+d_{2}x^{\alpha +4}+\cdots )^{\frac{2}{\alpha }}-1\\&\quad =(b_{1}+b_{2}- 1/12 )x^{2}+d_{1}x^{3}+d_{2}x^{4+\alpha }+\cdots , \end{aligned}$$

for some suitable positive constants $d_{1}, d_{2}$.

Combining this with (51), we obtain

$$\begin{aligned} \Big \Vert (z_{k})^{-2} - \Big ( \sum _{n=1}^{\infty } n \zeta ^{n} +\sum _{j=1}^{2} b_{j} \zeta ^{j} \Big ) k^{2} \Big \Vert \le C k^{3} | z |. \end{aligned}$$

This completes the proof of Lemma 8. $\square$

Remark 2

Observe that following the arguments as in Jin et al. [16], we can choose suitable coefficients $b_{1}=\frac{1}{12},b_{2}=0.$

Lemma 9

Let $z_{k}$ be defined as in (52). Let $b_{1},b_{2}$ be defined as in (51). Then, there holds

$$\begin{aligned} \Big \Vert (z^{\alpha } +A)^{-1} z^{-2} - (z_{k}^{\alpha } +A)^{-1} \Big ( \sum _{n=1}^{\infty } n \zeta ^{n} + \sum _{j=1}^{2} b_{j} \zeta ^{j} \Big ) k^{2} \Big \Vert \le C k^{3} | z |^{1- \alpha }. \end{aligned}$$

Proof

The proof is similar to the proof of Jin et al. [16, Lemma C.1.] and hence, we omit the proof here. $\square$

Lemma 10

Let $z_{k}$ be defined as in (52), then, there exists a positive constant C independent of k such that

$$\begin{aligned} \Big \Vert (z^{\alpha } +A)^{-1} z^{-3} - (z_{k}^{\alpha } +A)^{-1} \Big (k \sum _{n=1}^{\infty } \frac{t_{n}^2}{2!} \zeta ^{n} \Big ) \Big \Vert \le C k^{3} | z |^{-\alpha }. \end{aligned}$$

Proof

From Lemma 6, it follows that

$$\begin{aligned}&\Big \Vert (z^{\alpha } +A)^{-1} z^{-3} - (z_{k}^{\alpha } +A)^{-1} \Big (k \sum _{n=1}^{\infty } \frac{t_{n}^2}{2!} \zeta ^{n} \Big ) \Big \Vert \\&\quad \le \Vert (z^{\alpha } +A)^{-1} z^{-3}- (z_{k}^{\alpha } +A)^{-1} z_{k}^{-3}\Vert + \Big \Vert (z_{k}^{\alpha } +A)^{-1} \Big ( z_{k}^{-3}- k \sum _{n=1}^{\infty }\frac{t_{n}^2}{2!} \zeta ^{n} \Big ) \Big \Vert . \\&\quad \le \Big ( C \Vert (z^{\alpha } +A)^{-2}z^{\alpha -4} \Vert +C \Vert (z^{\alpha } +A)^{-1}z^{-4} \Vert \Big ) | z_{k}-z | \\&\qquad + \Vert (z_{k}^{\alpha } +A)^{-1} \Vert \Big |z_{k}^{-3}- k \sum _{n=1}^{\infty } \frac{t_{n}^2}{2!} \zeta ^{n}\Big | \le Ck^{3}|z|^{-\alpha }. \end{aligned}$$

This concludes the proof. $\square$

We are now ready to prove the main theorem of this section.

Proof

(Proof of Theorem2) We now calculate the approximate solution $V^{n}$ defined in (47)–(49). Taking the discrete Laplace transform in (47)–(49), we arrive at

$$\begin{aligned}&\sum _{n=1}^{\infty } \Big ( k^{-\alpha } \sum _{j=1}^{n} w^{(\alpha )}_{n-j} V^{j} \Big ) \zeta ^{n} + \sum _{n=1}^{\infty } (A_{h} V^{n}) \zeta ^{n} \\&\quad = (f_{h}(0) -A_{h} u_{0h}) \Big ( a_{1} \zeta + a_{2} \zeta ^{2}+ \frac{\zeta }{1-\zeta } \Big ) \\&\qquad + \sum _{n=1}^{\infty }t_{n} f^{\prime }_{h}(0)\zeta ^{n}+b_{1}k \zeta f^{\prime }_{h}(0) +b_{2}k\zeta ^{2} f^{\prime }_{h}(0) + \sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n}. \end{aligned}$$

Note that

$$\begin{aligned} \sum _{n=1}^{\infty }\left( \sum _{j=1}^{n} w^{(\alpha )}_{n-j}V^{j} \right) \zeta ^{n}=\left( \sum _{j=0}^{\infty } w^{(\alpha )}_{j}\zeta ^{j} \right) (V^{1} \zeta +V^{2}\zeta ^{2}+\cdots )={\tilde{w}}(\zeta ){\widetilde{V}}(\zeta ), \end{aligned}$$

and hence

$$\begin{aligned} {\widetilde{V}}(\zeta )&= (k^{-\alpha }{\tilde{w}}(\zeta )+A_{h})^{-1} \Big ( \big (f_{h}(0) -A_{h} u_{0h} \big ) \big (a_{1} \zeta + a_{2} \zeta ^{2}+ \frac{\zeta }{1-\zeta }\big ) \\&\quad +\left( \sum _{n=1}^{\infty }n\zeta ^{n} +b_{1}\zeta +b_{2}\zeta ^{2} \right) k f^{\prime }_{h}(0)+\sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n} \Big ). \end{aligned}$$

A use of the inverse discrete Laplace transform yields, with $\mu (\zeta )$ defined by (53), and $\varGamma _{k}$ as in (24).

$$\begin{aligned} V^{n}&=\frac{1}{2 \pi i} \int _{|\zeta |=\rho }\zeta ^{-n-1}{\widetilde{V}}(\zeta )d\zeta =\frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}}e^{zk}{\widetilde{V}}(e^{-zk})e^{-zk}k dz\nonumber \\&= \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} z_{k}^{-1} \mu (\zeta ) \big (f_{h}(0) -A_{h} u_{0h}\big ) \, dz\nonumber \\&\quad +\frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \left( \sum _{n=1}^{\infty }n\zeta ^{n} +b_{1} \zeta +b_{2}\zeta ^{2} \right) k^{2} f^{\prime }_{h}(0) \, dz \nonumber \\&\quad +\frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} k \big ( \sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n} \big ) \, dz. \end{aligned}$$

(63)

Now, subtracting (46) from (63), we arrive at

$$\begin{aligned} V_{h}(t_{n}) - V^{n} =I_{1} + I_{2} + I_{3}, \end{aligned}$$

where, with $K_{2}(z)$ defined by (25),

$$\begin{aligned} I_{1}&=\frac{1}{2 \pi i} \int _{ \varGamma / \varGamma _{k}} e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} z^{-1} (f_{h}(0)-A_{h} u_{0h} ) \, dz \\&\quad + \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} \Big (K_{2}(z)-\mu (e^{-zk})K_{2}(z_{k}) \Big ) (f_{h}(0)-A_{h} u_{0h} ) \, dz =I_{11}+I_{12}, \end{aligned}$$

and

$$\begin{aligned}&I_{2}=\frac{1}{2 \pi i} \int _{ \varGamma / \varGamma _{k}} e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} z^{-2} f^{\prime }_{h} (0) \, dz \\&\qquad + \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} \Big \{ (z^{\alpha } +A_{h})^{-1}z^{-2} f^{\prime }_{h}(0) \\&\qquad -(z_{k}^{\alpha } +A_{h})^{-1}\Big (\sum _{n=1}^{\infty }n\zeta ^{n} + \sum _{j=1}^{2}b_{j}\zeta ^{j} \Big )k^{2} f^{\prime }_{h}(0) \Big \} \, dz =I_{21}+I_{22}, \end{aligned}$$

and

$$\begin{aligned} I_{3}&=\frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} {\hat{R}}_{h}(z) \, dz \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big (k \sum _{n=1}^{\infty } R_{h}(t_{n}) \zeta ^{n} \Big ) \, dz =I_{31}+I_{32}. \end{aligned}$$

For $I_{1}$, apply the bound $\Vert (z^{\alpha }+A_h)^{-1}A_h\Vert \le C$, (3), (56) and (57) to obtain

$$\begin{aligned} \Vert I_{1} \Vert \le C k^{3} t_{n}^{-3} \Vert u_{0h} \Vert + Ck^{3} t_{n}^{\alpha -3} \Vert f_{h}(0) \Vert . \end{aligned}$$

For $I_{21}$, by (3), it follows that

$$\begin{aligned} \Vert I_{21} \Vert\le & {} C\Vert \int _{ \varGamma /\varGamma _{k}} e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} z^{-2} f^{\prime }_{h} (0) \, dz \Vert \\\le & {} C k^{3} t_{n}^{\alpha -2} \Vert f^{\prime } _{h}(0) \Vert . \end{aligned}$$

For $I_{22}$, a use of the Lemma 9 shows

$$\begin{aligned} \Vert I_{22}\Vert \le&C \Big \Vert \int _{\varGamma _{k}} e^{zt_{n}} (z^{\alpha } +A_{h})^{-1}z^{-2} f^{\prime }_{h}(0)\;dz \\&- \int _{\varGamma _{k}} e^{zt_{n}} (z_{k}^{\alpha } +A_{h})^{-1}\big [ \big (\sum _{n=1}^{\infty }n\zeta ^{n}+\sum _{j=1}^{2}b_{j}\zeta ^{j}\big )k^{2}f^{\prime }_{h}(0)\big ]\;dz\Big \Vert \\ \le&C k^{3} \int _{\varGamma _{k}} e^{-c t_{n} |z|}|z|^{1-\alpha }\, |dz| \Vert f^{\prime }_{h}(0) \Vert \le C k^{3} t_{n}^{\alpha -2} \Vert f^{\prime }_{h} (0) \Vert . \end{aligned}$$

Thus, we obtain

$$\begin{aligned} \Vert I_{2} \Vert \le C k^{3} t_{n}^{\alpha -2} \Vert f^{\prime }_{h} (0) \Vert . \end{aligned}$$

For $I_{3}$, we observe that

$$\begin{aligned} I_{31}&= \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}}(z^{\alpha } +A_{h})^{-1} \hat{R^{1}_{h}} (z) \, dz \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty } R^{1}_{h} (t_{n}) \zeta ^{n} \Big ) \, dz, \end{aligned}$$

and

$$\begin{aligned} I_{32}&= \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} \hat{R^{2}_{h}} (z) \, dz \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty } R^{2}_{h} (t_{n}) \zeta ^{n} \Big ) \, dz, \end{aligned}$$

where

$$\begin{aligned} R_{h}(t) = \frac{t^{2}}{2!} f^{\prime \prime }_{h} (0) + \big (\frac{t^{2}}{2!}*f^{\prime \prime \prime } _{h} \big )(t) =: R^{1}_{h}(t) + R^{2}_{h}(t). \end{aligned}$$

For $I_{31}$, we easily bound it as

$$\begin{aligned} \Vert I_{31}\Vert&=\Big \Vert \frac{1}{2 \pi i} \int _{ \varGamma } e^{z t_{n}} (z^{\alpha } +A_{h})^{-1} z^{-3} \, dz f^{\prime \prime }_{h}(0) \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} (z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty } R^{1}_{h} (t_{n}) \zeta ^{n} \Big ) \, dz \Big \Vert \\&=\Big \Vert \frac{1}{2 \pi i} \int _{ \varGamma /\varGamma _{k} } e^{z t_{n}}(z^{\alpha } +A_{h})^{-1} z^{-3} \, dz f^{\prime \prime }_{h} (0) \\&\quad - \frac{1}{2 \pi i} \int _{\varGamma _{k}} e^{z t_{n}} \Big ( (z^{\alpha } +A_{h})^{-1} z^{-3}-(z_{k}^{\alpha } +A_{h})^{-1} \Big ( k \sum _{n=1}^{\infty }\frac{t_{n}^{2}}{2!} \zeta ^{n} \Big ) \Big ) \, dz f^{\prime \prime }_{h} (0) \Big \Vert . \end{aligned}$$

An application of Lemma 10 yields

$$\begin{aligned} \Vert I_{31} \Vert \le C k^{3} t_{n}^{ \alpha -1} \Vert f^{\prime \prime }_{h} (0) \Vert . \end{aligned}$$

(64)

For $I_{32}$, following the arguments as in Jin et al. [14, 17], we now arrive at

$$\begin{aligned} \Vert I_{32}\Vert \le C k^{3} \int _{0}^{t_{n}} (t_{n}-s)^{ \alpha -1 } \Vert f^{\prime \prime \prime }_{h} (s) \Vert \, ds. \end{aligned}$$

Together these estimates complete the proof of Theorem 2. $\square$

In order to prove the error analysis in the completely discrete schemes, we now recall the error estimates of the semidiscrete scheme as is developed in [13] for the problem with $f=0$,

$$\begin{aligned} \Vert u(t_n)-u_h(t_n)\Vert \le C\;t^{-\alpha }\, h^2 \Vert u_0\Vert . \end{aligned}$$

(65)

Denote $U^{n}= V^{n} + P_{h} u_{0}$ with $V^{n}$ defined by (19) or (47)–(49), respectively. Using (65), we then have the following fully discrete error estimate

$$\begin{aligned} \Vert u(t_n)-U^n\Vert \le C \Big (t_n^{-\alpha } h^2 + t_n^{-m} k^m\Big )\; \Vert u_0\Vert , \end{aligned}$$

where $m=2$ or 3, respectively, for second order and third order time stepping schemes.

Since our main objective in this article is to derive higher order time stepping schemes, therefore, we may generalize our present results to include mass lumping scheme as discussed in [13].

4 Some Generalizations

This section is devoted to some generalization of the present methods to various other problems, which includes evolution equations with positive memory.

For instance, our generalizations include the following type of problems:

• 1. Evolution equations with positive memory called time diffusion-wave equation, as in [23],

  $$\begin{aligned} u'(x,t)+ \, _{0}^{R}D^{-\alpha }_{t} A u(x,t)=0, \quad \alpha \in (0,1), \end{aligned}$$

  (66)

  where $\, _{0}^{R}D^{-\alpha }_{t}$ denotes the Riemann-Liouville fractional integral operator.

• 2. The parabolic integro-differential equation with singular kernel, see, [26]

  $$\begin{aligned} u'(x,t)+\big (I + \, _{0}^{R}D^{-\alpha }_{t} \big ) A u(x,t)=0, \quad \alpha \in (0,1). \end{aligned}$$

  (67)

• 3. The Rayleigh-Stokes problem described by the time-fractional differential equation as in [4]

  $$\begin{aligned} u'(x,t)+ \big (I+\gamma \, _{0}^{R}D^{-\alpha }_{t} \big ) A u(x,t)=0,\quad \alpha \in (0,1), \end{aligned}$$

  (68)

  where $\gamma$ is a positive constant. In order to unify problems (66)–(68), we define ${\mathcal {J}}^\alpha$ denoting a time integral/differential operator and consider the unified problem by

  $$\begin{aligned} u'(x,t)+{\mathcal {J}}^\alpha A u(x,t)=0. \end{aligned}$$

  (69)

Now, a use of Laplace transforms in (69) yields

$$\begin{aligned} z{\hat{u}}+h(z)A{\hat{u}}=v, \end{aligned}$$

with some function h(z) depending on $\alpha$. Hence, with $\beta (z)=h(z)^{-1}$ we formally write the representation of solution as ${\hat{u}}=\beta (z) (z \beta (z) I+ A)^{-1} v =: {\hat{E}}_h(z)v.$ Here, for the problem (66), note that $\beta (z)=z^{\alpha },$ for the problem (67), $\beta (z)= z^\alpha /(1+z^{\alpha }),$ and for (68), $\beta (z)=1/(1+\gamma z^\alpha )$.

Assume that one can properly choose $\theta$ in $(\pi /2,\theta _0)$ with $\theta _0 \in (\pi /2,\pi )$ such that $z_k\beta (z_k)\in \varSigma _{\theta _0}$, where $z_k$ is an approximation of $z\in \varSigma _{\theta }.$ This is indeed possible in all given examples. With this, the resolvent estimate yields

$$\begin{aligned} \Vert (z\beta (z)I+ A )^{-1}\Vert \le \frac{M_{\theta _{0}}}{|z\beta (z)|}\quad \forall z\in \varSigma _{\theta }, \end{aligned}$$

(70)

where $M_{\theta _0}=1/\sin (\pi -\theta _0)$.

With corresponding ${A}_h$ as semidiscrete approximation of A ( as in Sect. 2), the semidiscrete method becomes: find ${\bar{u}}_h(t)\in V_h$ such that

$$\begin{aligned} \bar{u}_h'+ {\mathcal {J}}^\alpha A_h\bar{u}_h= 0\quad \quad t\in (0,T], \quad \bar{u}_h(0)=v_h. \end{aligned}$$

(71)

Following the arguments of [9, 23], one easily deduces property (70), when A is replaced by $A_h$.

Therefore, our two time stepping methods can be appropriately applied and the desired estimates can be easily derived.

Other generalizations also include the fractional cable equations with mass lumping by Al-Maskari and Karaa [2] and references therein. Further, it may include Fokker-Plank spatial discretization as our analysis does not depend on selfadjointness of the operator A and for complete discrete scheme, we may use the proposed time stepping schemes.

Finally, the present schemes can be applied to the Eq. (1), when A is a fractional Laplacian like $A= (-\varDelta )^{s},\;s\in (0,1)$, that is, with $0< \alpha , \beta <1$ and some positive constant $\delta >0$,

$$\begin{aligned} u' + \, _{0}^{R}D^{\alpha }_{t} \big (-\varDelta \big )^{s}\; u + \,\delta \, _{0}^{R}D^{\beta }_{t} u= & {} f \quad \text{ in } \; \varOmega \times (0,T], \end{aligned}$$

(72)

$$\begin{aligned} u= & {} 0 \quad \text{ in } \; \varOmega ^c \times (0,T], \end{aligned}$$

(73)

$$\begin{aligned} u(0)= & {} u_0 \quad \text{ in } \; \varOmega . \end{aligned}$$

(74)

An appropriate modification of the arguments in [1], we obtain the semidiscrete error estimate as:

$$\begin{aligned} \Vert u(t)-u_h(t)\Vert \le C t^{-\alpha }\;h^{s+\min \{s,1/2-\epsilon \}}\; \Vert u_0\Vert , \end{aligned}$$

with $\epsilon > 0$ small. Therefore, when it is combined with our proposed time stepping schemes, the final error estimate reads as: with $U^{n} \approx u_{h} (t_{n})$,

$$\begin{aligned} \Vert u(t_n)-U^n\Vert \le C \Big (t_n^{-\alpha }\;h^{s+\min \{s,1/2-\epsilon \}} + t_n^{-m} k^m\Big )\; \Vert u_0\Vert ,\quad m=2, 3. \end{aligned}$$

5 Numerical Simulations

In this section, we present five numerical examples to show that the numerical results are consistent with the theoretical results obtained in this paper. The first three examples are solved by using the numerical method (19) for both homogeneous and inhomogeneous problems in one- and two-dimensional cases. The last two examples are computed by using the numerical method (47)–(49) for both homogeneous and inhomogeneous problems in one-dimensional case.

Example 1

Consider, with $0<\alpha <1$,

$$\begin{aligned}&\, _0^C D_{t}^{\alpha } u(x,t) - \frac{\partial ^2 u (x, t)}{\partial x^2} = 0, \quad 0< x<1, \quad 0< t \le T, \\&u(0, t) = u(1, t) =0, \\&u(x, 0) = u_{0} (x), \end{aligned}$$

where (a) $u_{0}(x) = x(1-x)$ (smooth data) and (b) $u_{0}(x)= \chi _{[0, 1/2]}$ (nonsmooth data).

Let $N_{h}$ be a positive integer. Let $0=x_{0}< x_{1}< x_{2}< \cdots < x_{N_{h}} =1$ be the space partition and h the space step size. We shall use the piecewise linear finite element method to consider the space discretization.

Let $0<t_{0}< t_{1}< \cdots < t_{N} =T$ be the time partition and k the time step size. To observe the convergence order of the numerical method, we first need to calculate the reference solution $u_{ref}(t)$ at some fixed time T with very small step sizes $h_{ref} = 2^{-6}$ and $k_{ref} = 2^{-10}$.

By Theorem 1, we see that the numerical method (19) has the second order convergence. To see this convergence order, we shall calculate the approximate solution of u(T) at $T=1$ with the space step size $h=2^{-6}$ and the different time step sizes $k =\kappa * k_{ref}$ with $\kappa = [2^2, 2^3, 2^4,2^5, 2^6]$. In Table 2, we observe the convergence orders $O(k^{2})$ of the corrected scheme (19), where the rows with (a) denote the errors and the experimentally determined convergence orders in the smooth data case and the rows with (b) denote the errors and the orders in the nonsmooth data case. For each $\alpha$, we choose the average convergence order of the computed orders obtained by using the different time step sizes.

Table 2 Time convergence orders for the corrected scheme (19) in Example 1 at $T=1$

For the numerical method (19) with $c_{0}=0$, that is, for the uncorrected scheme, we observe that, in Table 3, the experimentally determined convergence order is only O(k) with both smooth and nonsmooth data.

Table 3 Time convergence orders for the uncorrected scheme (19) with $c_{0}=0$ in Example 1 at $T=1$

The second example is an inhomogeneous problem with zero initial value and the source term f which is smooth in time.

Example 2

Consider, with $0<\alpha <1$,

$$\begin{aligned}&\, _0^C D_{t}^{\alpha } u(x,t) - \frac{\partial ^2 u (x, t)}{\partial x^2} = f(x,t), \quad 0< x<1, \quad 0< t \le T, \\&u(0, t) = u(1, t) =0, \\&u(x, 0) = 0, \end{aligned}$$

where $f(x, t)= (\cos (t)+\sin (t)) ( 1 + \chi _{(0, 1/2)}(x))$.

We use the same parameters as in the numerical simulations in Example 1. In Table 4, we observe that the experimentally determined convergence order of the corrected scheme (19) indeed is $O(k^{2})$ for all $0< \alpha <1$ (Table 5).

Table 4 Time convergence orders for the corrected scheme (19) in Example 2 at $T=1$

Table 5 Time convergence orders for the uncorrected scheme (19) with $c_{0}=0$ in Example 2 at $T=1$

The third example is a two-dimensional example and we shall consider an inhomogeneous problem with nonsmooth initial data and the source term f which is smooth in time.

Example 3

Consider

$$\begin{aligned}&_0^C D_{t}^{\alpha } u (x,y, t) - \varDelta u (x, y, t) = f(x,y, t), \quad t \in (0, T], \; (x, y) \in \varOmega , \end{aligned}$$

(75)

$$\begin{aligned}&u(x, y, 0) = u_{0}(x, y), \quad (x, y) \in \varOmega , \end{aligned}$$

(76)

$$\begin{aligned}&u(x,y, t)=0, \quad t \in (0, T], \; (x, y) \in \partial \varOmega , \end{aligned}$$

(77)

where $\varOmega = (0, 1) \times (0, 1)$, $u_{0}(x, y)= \chi _{[0, 1/2]}(x) \chi _{[0, 1/2]}(y)$ and $f(x, y, t)= (\cos (t)+\sin (t)) ( 1 + \chi _{(0, 1/2)}(x)) ( 1 + \chi _{(0, 1/2)}(y))$. Note that f is smooth with respect to the time variable t.

Let $N_{h}$ be a positive integer. Let $0=x_{0}< x_{1}< x_{2}< \cdots < x_{N_{h}} =1$ and $0=y_{0}< y_{1}< y_{2}< \cdots < y_{N_{h}} =1$ be the partition of $\varOmega$. We divide $\varOmega$ into some triangles with the same sizes and let h be the maximal length of the sides of the triangle. We shall use the piecewise linear finite element method to consider the space discretization on the triangulation of $\varOmega$.

Let $0<t_{0}< t_{1}< \cdots < t_{N} =T$ be the time partition and k the time step size. We shall use the very small space step size $h_{ref} = 2^{-6}$ and the time step size $k_{ref} = 2^{-10}$ to calculate the reference solution at time T.

We shall choose $T=1$ in our simulation. We calculate the approximate solutions with the space step size $h=2^{-6}$ and the time step sizes $k =\kappa * k_{ref}$ with $\kappa = [2^2, 2^3, 2^4,2^5, 2^6]$. In Table 6, the experimentally determined convergence orders $O(k^{2})$ are observed as we expected.

Table 6 Time convergence orders for the corrected scheme (19) in Example 3 at $T=1$

In Table 7, we observe the experimentally determined convergence orders O(k) for the uncorrected scheme (19) with $c_{0}=0$ as we expected.

Table 7 Time convergence orders for the uncorrected scheme (19) with $c_{0}=0$ in Example 3 at $T=1$

In the next two examples, we shall consider the experimentally determined convergence orders of the numerical method (47)–(49).

Example 4

In this example, we shall use the numerical method (47)–(49) to solve Example 1. We use the same parameters as in the numerical simulation in Example 1. In Table 8, we observe that the corrected scheme (47)–(49) has the convergence orders $O(k^{3})$ with both smooth and nonsmooth data as expected.

Table 8 Time convergence orders for the corrected scheme (47)–(49) in Example 4 at $T=1$

For the uncorrected scheme (47)–(49), that is, $a_{1}=a_{2}=b_{1}=b_{2}=0$ in (47)–(49), we observe that, in Table 9, the experimentally determined convergence order is only O(k) with both smooth and nonsmooth data.

Table 9 Time convergence orders for the uncorrected scheme (47)–(49) with $a_{1}=a_{2}=b_{1}=b_{2}=0$ in Example  4 at $T=1$

Example 5

In this example, we shall use the numerical method (47)–(49) to solve Example 2. We use the same parameters as in the numerical simulation in Example 1. In Table 10, we also observe the convergence orders $O(k^{3})$ of the corrected scheme (47)–(49) in the inhomogeneous case.

In Table 11, the experimentally determined convergence order of the uncorrected scheme (47)–(49) with $a_{1}=a_{2}=b_{1}=b_{2}=0$ has only convergence order O(k).

Table 10 Time convergence orders for the corrected scheme (47)–(49) in Example 5 at $T=1$

Table 11 Time convergence orders for the uncorrected scheme (47)–(49) with $a_{1}=a_{2}=b_{1}=b_{2}=0$ in Example 5 at $T=1$

Remark 3

In Tables 8 and 10, we observe that the experimentally determined convergence orders are slightly better than the theoretical orders. The possible reason may be that the proposed numerical methods involve both fractional orders and the shifted numbers. These combinations which are more related to the equation may be instrumental in helping us to provide possibly more accurate computational results. But, we do not have a theory yet to establish it. Therefore, in future, we shall continue to investigate this interesting observation.

6 Declarations

Data sharing not applicable to this article as no datasets were generated or analysed during the current study.