The unified theory of shifted convolution quadrature for fractional calculus

TL;DR

This paper introduces the shifted convolution quadrature ($SCQ$) theory, extending existing methods for approximating fractional operators by incorporating a shift parameter, analyzing superconvergence, and validating with numerical tests.

Contribution

The paper develops a generalized $SCQ$ theory with a shifted parameter, broadening the scope of convolution quadrature methods for fractional calculus and analyzing superconvergence phenomena.

Findings

Superconvergence observed in some schemes.

New formulas designed for stability and convergence analysis.

Numerical tests confirm theoretical predictions.

Abstract

The convolution quadrature theory is a systematic approach to analyse the approximation of the Riemann-Liouville fractional operator $I^{\alpha}$ at node $x_{n}$. In this paper, we develop the shifted convolution quadrature ($SCQ$) theory which generalizes the theory of convolution quadrature by introducing a shifted parameter $\theta$ to cover as many numerical schemes that approximate the operator $I^{\alpha}$ with an integer convergence rate as possible. The constraint on the parameter $\theta$ is discussed in detail and the phenomenon of superconvergence for some schemes is examined from a new perspective. For some technique purposes when analysing the stability or convergence estimates of a method applied to PDEs, we design some novel formulas with desired properties under the framework of the $SCQ$. Finally, we conduct some numerical tests with nonsmooth solutions to further confirm our theory.

Tables

Table 1. Table 1: The shift-generalized Newton-Gregory formula of order 2 2 2

Scheme-I	Scheme-II
$\omega \left(\xi \right)={\left(1-\xi \right)}^{\alpha }\left[1+\left(\frac{\alpha }{2}-\theta \right)\left(1-\xi \right)\right]$	$\omega \left(\xi \right)={\left(1-\xi \right)}^{-\alpha }\left[1-\left(\frac{\alpha }{2}+\theta \right)\left(1-\xi \right)\right]$

Table 2. Table 2: The WSGL operator with q = − 2 , θ = p formulae-sequence 𝑞 2 𝜃 𝑝 q=-2,\theta=p

Scheme-I	Scheme-II
$\omega \left(\xi \right)={\left(1-\xi \right)}^{\alpha }\left[\frac{4+\alpha }{2\left(p+2\right)}+\frac{2p-\alpha }{2\left(p+2\right)}{\xi }^{p+2}\right]$	$\omega \left(\xi \right)={\left(1-\xi \right)}^{-\alpha }\left[\frac{4-\alpha }{2\left(p+2\right)}+\frac{2p+\alpha }{2\left(p+2\right)}{\xi }^{p+2}\right]$

Table 3. Table 3: The S ​ C ​ Q 𝑆 𝐶 𝑄 SCQ with ω ​ ( ξ ) 𝜔 𝜉 \omega(\xi) in ( 2.46 )

Scheme-I	Scheme-II
$\omega \left(\xi \right)={\left(\frac{3\alpha -2\theta }{2\alpha }-\frac{2\alpha -2\theta }{\alpha }\xi +\frac{\alpha -2\theta }{2\alpha }{\xi }^2\right)}^{\alpha }$	$\omega \left(\xi \right)={\left(\frac{3\alpha +2\theta }{2\alpha }-\frac{2\alpha +2\theta }{\alpha }\xi +\frac{\alpha +2\theta }{2\alpha }{\xi }^2\right)}^{-\alpha }$

Table 4. Table 4: The fractional central difference method ( 2.50 )

Scheme-I	Scheme-II
$\omega \left(\xi \right)=2^{-\alpha }{\left(1-{\xi }^2\right)}^{\alpha }$, with $\theta =\alpha$	$\omega \left(\xi \right)=2^{\alpha }{\left(1-{\xi }^2\right)}^{-\alpha }$, with $\theta =-\alpha$

Table 5. Table 5: The temporal convergence rate with h = 1 5000 ℎ 1 5000 h=\frac{1}{5000}

$\alpha$	$\theta$	$\tau$	$E_c(\tau ,h)$	rate	$E_o(\tau ,h)$	rate
		1/10	2.35436E-05	—	2.87659E-03	—
	0	1/20	5.47471E-06	2.1045	2.79338E-03	0.0423
		1/40	1.13371E-06	2.2717	2.71399E-03	0.0416
0.1		1/80	2.55465E-07	2.1499	2.63920E-03	0.0403
		1/10	9.59982E-04	—	4.22249E-02	—
	0.05	1/20	2.43373E-04	1.9798	3.81734E-02	0.1455
		1/40	6.15520E-05	1.9833	3.44973E-02	0.1461
		1/80	1.57177E-05	1.9694	3.11936E-02	0.1452
		1/10	1.72824E-03	—	5.64791E-03	—
	0.1	1/20	4.38829E-04	1.9776	2.56403E-03	1.1393
		1/40	1.10750E-04	1.9864	4.53468E-04	2.4993
0.5		1/80	2.80539E-05	1.9810	8.92050E-04	-0.9761
		1/10	3.22333E-03	—	2.00719E-02	—
	0.2	1/20	8.15120E-04	1.9835	1.28495E-02	0.6435
		1/40	2.05120E-04	1.9905	7.75159E-03	0.7291
		1/80	5.16829E-05	1.9887	4.28721E-03	0.8545
		1/10	4.70189E-03	—	6.08467E-03	—
	0.4	1/20	1.21112E-03	1.9569	1.52565E-03	1.9958
		1/40	3.07190E-04	1.9791	6.99968E-04	1.1241
0.9		1/80	7.75617E-05	1.9857	6.01301E-04	0.2192
		1/10	4.90787E-03	—	6.31974E-03	—
	0.45	1/20	1.25984E-03	1.9619	1.58234E-03	1.9978
		1/40	3.19039E-04	1.9814	3.96145E-04	1.9980
		1/80	8.04831E-05	1.9870	1.45886E-04	1.4412

Table 6. Table 6: The spacial convergence rate with τ = 1 1000 𝜏 1 1000 \tau=\frac{1}{1000}

$\alpha$	$\theta$	$h$	$E_c(\tau ,h)$	rate	$E_o(\tau ,h)$	rate
		1/10	9.90558E-02	—	9.90559E-02	—
	0	1/20	2.49005E-02	1.9921	2.49007E-02	1.9921
		1/40	6.23369E-03	1.9980	6.23387E-03	1.9980
0.3		1/80	1.55895E-03	1.9995	4.29525E-03	0.5374
		1/10	9.90560E-02	—	9.90561E-02	—
	0.15	1/20	2.49008E-02	1.9921	2.49009E-02	1.9920
		1/40	6.23393E-03	1.9980	6.23403E-03	1.9980
		1/80	1.55920E-03	1.9993	7.03535E-03	-0.1745
		1/10	9.85495E-02	—	9.85515E-02	—
	0.3	1/20	2.47681E-02	1.9924	2.47686E-02	1.9924
		1/40	6.20048E-03	1.9980	6.20066E-03	1.9980
0.8		1/80	1.55092E-03	1.9993	1.55102E-03	1.9992
		1/10	9.85493E-02	—	9.85515E-02	—
	0.4	1/20	2.47681E-02	1.9924	2.47687E-02	1.9924
		1/40	6.20052E-03	1.9980	6.20072E-03	1.9980
		1/80	1.55098E-03	1.9992	1.55108E-03	1.9992

Equations

$$\begin{array}{r}{I}^{\alpha }f(x)=\frac{1}{\Gamma (\alpha )}{\int }_0^x(x-s{)}^{\alpha -1}f(s)\mathrm{d}s,\text{for }\Re \alpha >0,\end{array}$$

$$\begin{array}{r}{I}^{-\alpha }f(x)=\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}{I}^{n-\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}{\int }_0^x(x-s{)}^{n-\alpha -1}f(s)\mathrm{d}s,\end{array}$$

$$(I^{\alpha }f)(x)\approx h^{\alpha }\sum _{j=0}^n{\omega }_{j}f(x-jh)+h^{\alpha }\sum _{j=0}^s{w}_{n,j}f(jh),$$

$$(I^{\alpha }f)(x)\approx h^{\alpha }\sum _{j=0}^n{\omega }_{j}f(x+(\theta -j)h)+h^{\alpha }\sum _{j=0}^s{w}_{n,j}f(jh).$$

$$(I_h^{\alpha }f)(x)=h^{\alpha }\sum _{j=0}^n{\omega }_{j}f(x+(\theta -j)h)+h^{\alpha }\sum _{j=0}^s{w}_{n,j}f(jh),x=(n-\theta )h,$$

$$E_h^{\alpha }={\Omega }_h^{\alpha }-I^{\alpha }.$$

$${\Omega }_h^{\alpha }f(x)=h^{\alpha }\sum _{0\le jh\le x+\theta h}{\omega }_{j}f(x+(\theta -j)h),$$

$$E_h^{\alpha }(f\ast g)=(E_h^{\alpha }f)\ast g.$$

$$(E_h^{\alpha }{t}^{\beta -1})(x)=x^{\alpha +\beta -1}(E_{h/x}^{\alpha }{t}^{\beta -1})(1).$$

$${\omega }_n=O(n^{\alpha -1}).$$

$$h^{\alpha }{e}^{\theta h}\omega (e^{-h})=1+O(h^p).$$

$$(E_h^{\alpha }{t}^{\beta -1})(1)=O(h^{\beta })+O(h^p),\text{for all }\beta \in \mathbb{C},\beta \ne 0,-1,-2,\cdots .$$

$$(E_h^{\alpha }{t}^{\beta -1})(x_n)=O(x_n^{\alpha -1}{h}^{\beta })+O(x_n^{\alpha +\beta -1-p}{h}^p),\text{for all }\beta \in \mathbb{C},\beta \ne 0,-1,-2,\cdots ,$$

$$\begin{array}{rl}{e}_h(x):=(E_h^{\alpha }{e}^{t-x})(x)& =h^{\alpha }\sum _{0\le jh\le x+\theta h}{\omega }_j{e}^{(\theta -j)h}-(I^{\alpha }{e}^{t-x})(x)\\ & =h^{\alpha }{e}^{\theta h}\sum _{0\le jh\le x+\theta h}{\omega }_j{e}^{-jh}-(I^{\alpha }{e}^{t-x})(x).\end{array}$$

$$\begin{split}\varpi(\xi)=(1-\xi)^{-\alpha}\big{[}c_{0}+c_{1}(1-\xi)+c_{2}(1-\xi)^{2}+\cdots+c_{N-1}(1-\xi)^{N-1}+(1-\xi)^{N}\tilde{r}(\xi)\big{]}.\end{split}$$

$$\begin{array}{r}{c}_i={\gamma }_i,\text{for }i=0,1,\cdots ,p-1,\end{array}$$

$$\begin{split}\sum_{i=0}^{\infty}\gamma_{i}(1-\xi)^{i}=\xi^{\theta}\bigg{(}\frac{-\ln\xi}{1-\xi}\bigg{)}^{-\alpha}.\end{split}$$

$$\omega(\xi)=\bigg{[}\frac{p_{1}(\xi)}{p_{2}(\xi)}\bigg{]}^{\alpha}\frac{p_{3}(\xi)}{p_{4}(\xi)},$$

$$\omega (\xi )=\nu (\xi )\prod _{j=0}^{\ell }(\xi -{\xi }_j{)}^{-{\alpha }_j},{\xi }_0=1,{\alpha }_0=\alpha ,$$

$$\text{Condition-ω: }\omega (\xi )\text{ is analytic without zeros in D, and }\Re {\alpha }_j\le \Re \alpha .$$

$$\begin{array}{r}f(x)=x^{\beta -1}g(x)\text{with }g\text{ sufficiently differentiable},\end{array}$$

$$\begin{array}{r}{I}_h^{\alpha }f(x)-I^{\alpha }f(x)=O(h^p)\end{array}$$

$$\begin{array}{r}(I_h^{\alpha }{t}^{q+\beta +1})(t_{n-\theta })=(I^{\alpha }{t}^{q+\beta +1})(t_{n-\theta })\end{array}$$

$$\begin{split}\text{The fractional trapezoidal rule: }\varpi_{2}(\xi)&=\frac{(1+\xi)^{\alpha}}{2^{\alpha}(1-\xi)^{\alpha}},\quad\Re\alpha\geq 0.\\ \text{The fractional BDF-p: }\varpi_{p}(\xi)&=\bigg{[}\sum_{i=1}^{p}\frac{1}{i}(1-\xi)^{i}\bigg{]}^{-\alpha},\quad 1\leq p\leq 6.\\ \text{The fractional BT-$\theta$ method: }\varpi_{2}(\xi)&=\bigg{[}\frac{1-\theta+\theta\xi}{(3/2-\theta)-(2-2\theta)\xi+(1/2-\theta)\xi^{2}}\bigg{]}^{\alpha},\\ &\quad\text{ for }\theta\in(-\infty,1/2),\alpha<0,\text{ or }\theta\in(-\infty,1/2],\alpha\geq 0.\end{split}$$

$$\begin{split}\varpi_{p}(\xi)=(1-\xi)^{-\alpha}\big{[}\gamma^{\prime}_{0}+\gamma^{\prime}_{1}(1-\xi)+\cdots+\gamma^{\prime}_{p-1}(1-\xi)^{p-1}\big{]},\end{split}$$

$$\begin{split}\varpi_{2}(\xi)=\frac{1-\alpha\theta+\alpha\theta\xi}{\big{[}(3/2-\theta)-(2-2\theta)\xi+(1/2-\theta)\xi^{2}\big{]}^{\alpha}},\quad\theta\in(-\infty,1],~{}\alpha\in\mathbb{R},~{}\alpha\theta\leq\frac{1}{2}.\end{split}$$

$$\begin{split}\omega(\xi)=\varpi_{p}(\xi)\big{[}\gamma^{\prime\prime}_{0}+\gamma^{\prime\prime}_{1}(1-\xi)+\cdots+\gamma^{\prime\prime}_{p-1}(1-\xi)^{p-1}\big{]},\end{split}$$

$$\begin{array}{r}\sum _{i=0}^{\infty }{\gamma }_i^{\prime \prime }(1-\xi {)}^i={\xi }^{\theta }.\end{array}$$

$$\begin{array}{r}{e}^{\theta h}\sum _{i=0}^{p-1}{\gamma }_i^{\prime \prime }(1-e^{-h}{)}^i=1-e^{\theta h}\sum _{i=p}^{\infty }{\gamma }_i^{\prime \prime }(1-e^{-h}{)}^i=1+O(h^p).\end{array}$$

$$\begin{array}{r}{h}^{\alpha }{\varpi }_p(e^{-h})=1+O(h^p),\end{array}$$

Full text

The unified theory of shifted convolution quadrature for fractional calculus

††thanks: Corresponding author.

Email addresses: mathliuyang@imu.edu.cn,

Manuscript submitted to Journal  2019

Yang Liu1∗, Baoli Yin1, Hong Li1, Zhimin Zhang2,3

1School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China;

2Beijing Computational Science Research Center, Beijing 100193, China;

3Department of Mathematics, Wayne State University, Detroit, MI 48202, USA

Abstract: The convolution quadrature theory is a systematic approach to analyse the approximation of the Riemann-Liouville fractional operator $I^{\alpha}$ at node $x_{n}$. In this paper, we develop the shifted convolution quadrature ($SCQ$) theory which generalizes the theory of convolution quadrature by introducing a shifted parameter $\theta$ to cover as many numerical schemes that approximate the operator $I^{\alpha}$ with an integer convergence rate as possible. The constraint on the parameter $\theta$ is discussed in detail and the phenomenon of superconvergence for some schemes is examined from a new perspective. For some technique purposes when analysing the stability or convergence estimates of a method applied to PDEs, we design some novel formulas with desired properties under the framework of the $SCQ$. Finally, we conduct some numerical tests with nonsmooth solutions to further confirm our theory.

Keywords: shifted convolution quadrature, generating functions, Riemann-Liouville fractional calculus operator, stability regions

1 Introduction

The fractional calculus has drawn much attention in recent years for its wide applications and theoretical interests, see [31, 24, 25, 27, 28, 29, 30, 23, 22, 33, 26, 9, 10]. In this paper, we are particularly concerned about the Riemann-Liouville calculus operator $I^{\alpha}$ which is defined by

[TABLE]

and the Riemann-Liouville differential operator $I^{-\alpha}(\Re\alpha>0)$ is defined by

[TABLE]

where $n=\lceil\Re\alpha\rceil$. When $\alpha=0$, we set $I^{0}=I$, the identity operator.

In 1986, Lubich [1] developed the convolution quadrature ($CQ$) theory to approximate the Riemann-Liouville calculus with arbitrary $\alpha\in\mathbb{C}$ at the node $x=nh$,

[TABLE]

where $h$ is the mesh size of a uniform grid, $\omega_{j}$ denotes the convolution weight and $w_{n,j}$ is the starting weight. For brevity, define $x_{k}=kh$. For $\Re\alpha\leq 0$, the theory requires that $f^{(j)}(0)=0(j=0,\cdots,-\lceil\Re\alpha\rceil)$. When $\alpha$ takes integers, $I^{\alpha}$ coincides with the traditional integral operator $(\alpha>0)$ and differential operator $(\alpha<0)$, and the identity operator $(\alpha=0)$. From this aspect, we can conclude that the $CQ$ theory generalized the traditional methods for approximating calculus operators with integer orders to fractional, i.e., arbitrary orders.

Nonetheless, some classical methods are still excluded from the theory, such as the Crank-Nicolson scheme (a special case of the BDF2-$\theta$ method, see [8, 37]) which approximates the first derivative at the node $x_{n-\frac{1}{2}}$ and its conterpart for fractional derivative $I^{-\alpha}(\alpha\in(0,1))$ approximated at the node $x_{n-\frac{\alpha}{2}}$, with a second-order convergence rate (see the $\omega(\xi)$ of order 2 in (2.36)). Such kind of superconvergence schemes are important for the numerical analysis for PDEs since the resulted schemes possess some good characters, see [5]. Some other methods such as the shifted Gr$\ddot{\text{u}}$nwald formula [16], the WSGL operator [6], the method developed by Ding et al. [11] and recently proposed higher-order approximation formulas [18] that generalize the fractional BDFp for $p=1,2,\cdots,6$ are excluded from the $CQ$ theory as well, since all above methods approximate $I^{\alpha}$ at a shifted node.

In this paper, we generalize the $CQ$ theory to cover the methods mentioned above by introducing a shifted parameter $\theta$ to develop the shifted convolution quadrature ($SCQ$) theory and further design some new methods with higher-order convergence rate. Specifically, we approximate the operator $I^{\alpha}(\alpha\in\mathbb{C})$ at $x=(n-\theta)h$,

[TABLE]

The contributions of the paper are as follows

$\bullet$ Develop systematic approaches to approximate the fractional calculus $I^{\alpha}f(x)(\alpha\in\mathbb{C})$ at node $x_{n-\theta}$ without assumptions on the regularity of $f(x)$. Unlike most papers concerning the approximation at a shifted node with the parameter $\theta$, that pay little attention on the choice of $\theta$, we explore the criterion for this choice from the aspect of generating functions. See Sec. 2.

$\bullet$ Examine impacts of the parameter $\theta$ on the absolute stability regions for different numerical schemes. With a proper shifted parameter $\theta$, we can get a A-stable method, which is superior to others for some problems. See Sec. 3.

$\bullet$ Construct some new approximation methods (for $I^{\alpha}f(x_{n-\theta})$) based on known approximation methods (for $I^{\alpha}f(x_{n})$), see the Theorem 2.11 and Example 4. Reveal some facts about the superconvergence for the numerical schemes (known or newly developed) from a new perspective, see the Example 2.

$\bullet$ Generalize the correction technique of the $CQ$ theory to the $SCQ$ theory by introducing a parameter $\theta$, see the Theorem 2.9 and Remark 2.10. This generalization is important since solutions of fractional calculus equations generally show some singularity at initial node. Now with the correction technique, all methods that belong to the framework of the $SCQ$ can be modified to obtain the optimal convergence rate.

$\bullet$ Apply a novel class of second-order shift-generalized Newton-Gregory formula (2.32) to the time fractional diffusion equation with stability analysis and error estimate, see Sec. 4. Based on the analysis in foregoing sections we now pick on a set of $\theta$ with which the fractional Grönwall inequality can be employed in the subsequent numerical analysis.

We organize the rest of the paper as follows: In Sec. 2, we generalize the theory of the $CQ$ and develop some definitions, lemmas and theorems for the $SCQ$. Some existing or newly proposed schemes are discussed by several examples. In Sec. 3, we analyse the stability regions for some $SCQ$s aforementioned, discuss the impact on the regions for different $\theta$. In Sec. 4, we devise a novel numerical scheme for the time-fractional diffusion equation by the $SCQ$ theory, with the purpose of the easy application of the discrete fractional Gr$\ddot{\text{o}}$nwall inequality. Some lemmas are proved and the stability estimates as well as the optimal convergence order are derived. In the end of the section, we conduct some numerical tests to further confirm our theoretical analysis. Finally, in Sec. 5, we make some conclusions and discuss some approaches we may take in the future work.

2 Stability, consistency and convergence of the $SCQ$

In this section, we mainly generalize the equivalence theorem developed by Lubich (Theorem 2.5, [1]) which extends the classical theorem of Dahlquist [3] on linear multistep methods to fractional ones. For convenience in the subsequent analysis we introduce some notations and definitions. Define

[TABLE]

as the shifted convolution quadrature ($SCQ$). Denote $\Omega_{h}^{\alpha}f(x)=h^{\alpha}\sum_{j=0}^{n}\omega_{j}f(x+(\theta-j)h)$ as the convolution part and $S_{h}^{\alpha}f=h^{\alpha}\sum_{j=0}^{s}w_{n,j}f(jh)$ as the starting part of (2.1), respectively. Define the convolution quadrature error by

[TABLE]

For the sequence of convolutions weights $\{\omega_{j}\}_{j=0}^{\infty}$ we associate with a generating power series $\omega(\xi)=\sum_{j=0}^{\infty}\omega_{j}\xi^{j}$, and viceversa.

We note that if $f$ is continuous and $g$ is locally integrable, and $\Omega_{h}^{\alpha}$ is extended for $x\geq 0$ such that

[TABLE]

then $\Omega_{h}^{\alpha}$ commutes with convolution $\Omega_{h}^{\alpha}(f\ast g)=(\Omega_{h}^{\alpha}f)\ast g$. Hence, the convolution error $E_{h}^{\alpha}=\Omega_{h}^{\alpha}-I^{\alpha}$ satisfies (for $\Re\alpha\leq 0$ , we require $f^{(j)}(0)=0,j=0,\cdots,-\lceil\Re\alpha\rceil$)

[TABLE]

Another key property of $E_{h}^{\alpha}$ is the homogeneity:

[TABLE]

We remark that (2.4) and (2.5) are crucial for Theorem 2.8. The proof of (2.3) and (2.4) are omitted here since their correctness can be directly checked. Next, we introduce three definitions that are closely connected in the subsequent Theorem 2.8:

For arbitrary $\alpha\in\mathbb{C}$,

Definition 2.1

A SCQ is stable for $I^{\alpha}$ if the convolution weights is

[TABLE]

Definition 2.2

A SCQ is consistent of order $p$ for $I^{\alpha}$ if the generating function of $\{\omega_{j}\}_{j=0}^{\infty}$ satisfies

[TABLE]

Definition 2.3

A SCQ is convergent of order $p$ to $I^{\alpha}$ if

[TABLE]

We take note of the fact that the definition of stability and convergence of the $SCQ$ are the same as the corresponding definitions of the $CQ$ (see Definition 2.1 and 2.3, [1]). The consistency of the $CQ$ is a special case of the $SCQ$ when $\theta=0$, i.e., the $CQ$ is consistent of order $p$ for $I^{\alpha}$ if $h^{\alpha}\omega(e^{-h})=1+O(h^{p})$ (see Definition 2.2, [1]).

Remark 2.4

Combining the homogeneity of $E_{h}^{\alpha}$ and the definition of convergent (2.8), we can get

[TABLE]

which means for those small $\beta$ that $\alpha+\beta-1-p<0$, the convergence rate will be much lower than $p$. Actually for fractional calculus equations, the solution is generally of weak regular at initial node. We shall cope with such problem in the Theorem 2.9.

The following two lemmas which reveal some facts about the consistency of the $SCQ$, generalize the arguments for the $CQ$ in [1].

Lemma 2.5

(The counterpart of Lemma 3.1, [1]) If $(E_{h}^{\alpha}t^{k-1})(1)=O(h^{k})+O(h^{p})$ for $k=1,2,3,\cdots,$ then the $SCQ$ is consistent of order $p$.

**Proof. ** First we examine the convolution error for the function $e^{t-x}$ with respect to $t$ on the interval $[0,x]$,

[TABLE]

Let $x\to\infty$, the expression $h^{\alpha}e^{\theta h}\sum_{0\leq jh\leq x+\theta h}\omega_{j}e^{-jh}$ tends to $h^{\alpha}e^{\theta h}\omega(e^{-h})$. The rest argument of the proof is exactly the same as that of the Lemma 3.1 in [1], which is omitted here. The proof of the lemma is completed.

As is pointed out in [1], the structure of the generating function $\varpi(\xi)$ for a consistent $CQ$ is of the form (see (3.6) in [1])

[TABLE]

where $\tilde{r}(\xi)$ is holomorphic at $1$, and constants $c_{j}=\gamma^{\prime}_{j}$ where $\gamma^{\prime}_{j}$ are defined by (2.21). We argue that the generating function $\omega(\xi)$ for $SCQ$ can be expressed similarly by (2.11) with a different definition of the coefficients $c_{i}$ that depend on $\theta$:

Lemma 2.6

(The counterpart of Lemma 3.2, [1]) The $SCQ$ is consistent of order $p$ if and only if the coefficients $c_{i}$ in (2.11) satisfy

[TABLE]

where $\gamma_{i}$ are the coefficients of

[TABLE]

Proof. The proof is almost the same as Lemma 3.2 in [1], and here is omitted.

Remark 2.7

We can construct $\omega(\xi)$ by polynomial functions. Assume $\omega(\xi)$ has the form

[TABLE]

where $p_{i}(\xi)$ are polynomial functions, and $p_{2}(1)=0$. Denote by $D$ the unit disc $|\xi|<1$, and by $\bar{D}$ the closed unit disc in the complex plane. Then, for a stable $SCQ$, which means $\omega_{n}=O(n^{\alpha-1})$, it holds that $\omega(\xi)$ is analytic in $D$. Hence, $\omega(\xi)$ can be written as

[TABLE]

where $\xi_{j}$ are distinct with $|\xi_{j}|=1$, $\nu(\xi)$ is analytic on $\bar{D}$, and $\nu(\xi_{j})\neq 0$, $\alpha_{j}\neq 0,-1,-2,\cdots$. It can be shown that $\Re\alpha_{j}\leq\Re\alpha$ is equivalent to (2.6), see [1]. We limit the choice of the parameter $\theta$ by the Condition-$\omega$,

[TABLE]

With the above analysis we can establish the main theorem in this paper that connects the definitions of being stable, consistent and convergent for the $SCQ$:

Theorem 2.8

A $SCQ$ with convolution weights defined by a generating function $\omega(\xi)$ satisfying the condition-$\omega$ is convergent of order $p$ if and only if it is stable and consistent of order $p$.

Proof. The theorem is a generalization of the Theorem 2.4 in [1], whose proof consists of several lemmas. Since the definition of stability and convergence of the $SCQ$ are the same as those of the $CQ$, we merely generalize the lemmas in [1] concerning the consistency, i.e., the Lemma 3.1 and Lemma 3.2 therein. Now with Lemma 2.5 and Lemma 2.6, we complete the proof of the theorem.

For a convergent $SCQ$, the following theorem shows that with the starting part $S_{h}^{\alpha}$, $I_{h}^{\alpha}f(x)$ approximates $I^{\alpha}f(x)$ uniformly for bounded $x$.

Theorem 2.9

(See Theorem 2.4 in [1]) Suppose $\omega(\xi)$ is the generating function of a convergent $SCQ$. Then we have: For any $\beta\neq 0,-1,-2,\cdots$,

(i) there exist starting weights $w_{n,j}=O(n^{\alpha-1})(n\geq 0,j=0,\cdots,s)$ such that for any function

[TABLE]

the SCQ satisfies

[TABLE]

uniformly for $x\in[a,b]$ with $0<a<b<\infty$.

(ii) there exist bounded starting weights $w_{n,j}$ such that for any function (2.17), the estimate (2.18) holds uniformly for bounded $x$.

Remark 2.10

The starting weights $w_{n,j}$ in (i) are derived by letting

[TABLE]

for all integer $q\geq 0$ such that $\Re(q+\beta-1)\leq p-1$. Note that the estimate (2.18) requires that $x$ is bounded away from [math]. The starting weights $w_{n,j}$ in (ii) require additionally those $\Re(q+\alpha+\beta-1)<p$ such that (2.18) holds for bounded $x$.

In the rest of the section we shall collect or devise as many generating functions for the $SCQ$ as possible. Indeed any generating function for the $CQ$ can be transformed into the one for the $SCQ$, as is proved in Theorem 2.11. So let us first recall some classical or newly developed generating functions $\varpi_{p}(\xi)$ ($p$ is the convergence order) for the $CQ$:

$\bullet$ (Theorem 2.6, [1]) For an implicit linear multistep method which is stable and consistent of order $p$ with the characteristic polynomials $\rho$ and $\sigma$, if the zeros of $\sigma(\xi)$ have absolute value less than $1$, then \varpi_{p}(\xi)=\big{[}\sigma(1/\xi)/\rho(1/\xi)\big{]}^{\alpha}. Some special cases include (see [1, 22])

[TABLE]

$\bullet$ The generalized Newton-Gregory formula (see [1])

[TABLE]

with the coefficients $\gamma^{\prime}_{i}$ defined by \sum_{i=0}^{\infty}\gamma^{\prime}_{i}(1-\xi)^{i}=\big{(}\frac{\ln\xi}{\xi-1}\big{)}^{-\alpha}.

$\bullet$ The fractional BN-$\theta$ method (see [22])

[TABLE]

Theorem 2.11

Suppose $\varpi_{p}(\xi)$ is the generating function of a $CQ$ with convergence order $p$. Define

[TABLE]

where the coefficients $\gamma^{\prime\prime}_{i}$ satisfy

[TABLE]

Then the SCQ with convolution weights generated by $\omega(\xi)$ is convergent of order $p$.

Proof. For the generating function $\varpi_{p}(\xi)=\sum_{j=0}^{\infty}\varpi_{j}\xi^{j}$, we know that $\varpi_{n}=O(n^{\alpha-1})$ by the stability of the $CQ$. Hence, we have $\omega_{n}=O(n^{\alpha-1})$. By replacing $\xi$ with $e^{-h}$ in (2.24), we get

[TABLE]

Considering the consistency of the $CQ$, i.e.,

[TABLE]

by combining (2.23), (2.25) and (2.26), we have

[TABLE]

Then by Theorem 2.8 the proof is completed.

A special case for (2.23) that is of vital importance is when $\alpha=0$. For any $\varpi_{p}(\xi)$ defined by (2.20)-(2.22), taking $\alpha=0$, we obtain the generating function by (2.23) for the approximation of $I^{0}f(x)=f(x)$ at node $x_{n-\theta}$. The coefficients $\gamma^{\prime\prime}_{i}$ in (2.24) can be formulated as

[TABLE]

Combining Theorem 2.11, (2.5) and (2.8), we have the following result

Corollary 2.12

Let $f(x)=x^{\beta-1}$ with $\beta\neq 0,-1,-2,\cdots$. Then we have the following approximation formula for $f(x_{n-\theta})$ of order $p$,

[TABLE]

where the weights $\theta_{j}$ satisfy

[TABLE]

with coefficients $\gamma^{\prime\prime}_{i}$ defined in (2.28).

Remark 2.13

The approximation formula (2.29) holds for $f(x)=x^{\beta-1}g(x)$ with $g(x)$ sufficiently differentiable as well. For $p=2,3$, we can obtain two popular formulas

[TABLE]

see also (17), (18) in [4]. Considering the condition-$\omega$, we limit $\theta$ for (2.31) to satisfying $\theta\leq\frac{1}{2}$ and $\theta\leq 1-\frac{\sqrt{2}}{2}$, respectively. As one can see from (2.29), to get a convergence order of $p$ for bounded $x$ if $f(x)$ is not so regular, the staring part is needed according to the Theorem 2.9.

Another family of generating functions for the $SCQ$ is the shift-generalized Newton-Gregory formula which is a natural result from the analysis of Lemma 2.6:

Corollary 2.14

(The shift-generalized Newton-Gregory formula) The SCQ with weights generated by the following generating function is convergent (to $I^{\alpha}$) of order $p$,

[TABLE]

where $\gamma_{i}$ are the coefficients of

[TABLE]

One can find out that the shift-generalized Newton-Gregory formula reduces to (2.21) when $\theta=0$. We conclude this section by further exploring some numerical methods that approximate $f(x)$ at node $x_{n-\theta}$ in the following examples. All methods mentioned in the examples belong to the SCQ, hence correction technique can be applied to the methods if the solution is not regular enough. For simplicity, we assume $\alpha\in\mathbb{R}$.

**Example 1. **(Generalized shifted Gr$\ddot{\text{u}}$nwald formula)

Expanding (2.33), we can easily derive that

[TABLE]

For the case with convergence order $p=1$, we have the approximation formula as

[TABLE]

Now by assuming that $f$ is sufficiently smooth (hence, the starting part can be omitted), and that $\theta$ takes nonnegative integers, the relation (2.35) reduces to the shifted Gr$\ddot{\text{u}}$nwald formula (see [16] with $\alpha$ replaced by $-\alpha$). Recently, Chen et al. [21] has applied the Gr$\ddot{\text{u}}$nwald formula ($\theta=0$) to the time fractional PDEs and derived a sharp convergence rate which is in line with (2.9).

**Example 2. **(Discussion on the superconvergence)

With the coefficients in (2.34), the generating function $\omega(\xi)$ for the shift-generalized Newton-Gregory formula of order $p=2,3$ can be formulated as

[TABLE]

Both of the formula were also proposed by Dimitrov [4] by using the theory developed in [17]. An obvious conclusion is that $\omega(\xi)$ defined by (2.36) which is of order $2$ shows some superconvergence when $\theta=-\frac{\alpha}{2}$ (see also [5], with $\alpha$ replaced by $-\alpha$), since generally a $CQ$ with the generating function $(1-\xi)^{-\alpha}$ approximates $I^{\alpha}f(x_{n})$ with a lower order. From the perspective of generating function (2.32), we can always find some superconvergence points that effective numerical methods can be proposed. To the best knowledge of the authors, there is no literature exploring the superconvergence property of the scheme with $\omega(\xi)$ defined by (2.36) which is of order $3$, that from the discussion above, a $SCQ$ with $\omega(\xi)$ defined by

[TABLE]

is convergent of order $3$ provided that

[TABLE]

One can check that this choice of $\theta$ satisfies the condition-$\omega$.

**Example 3. **(The WSGL method, see [6, 7, 15, 20, 19, 34, 35, 36])

Assume two integers $p>q$, and let $\theta=p$, consider the following generating function

[TABLE]

whose coefficients can be formulated as

[TABLE]

where $\varrho_{k}$ are the coefficients of $(1-\xi)^{-\alpha}$ and define $\varrho_{k}=0$ for $k<0$. Considering the condition-$\omega$, we assume $(p,q)$ satisfies

[TABLE]

One can easily check that the $SCQ$ with $\omega(\xi)$ defined in (2.39) is the WSGL operator (see [6], with $\alpha$ replaced by $-\alpha$), which is convergent of order $2$ by Theorem 2.8. Actually, by (2.40), we have $\omega_{n}=O(n^{\alpha-1})$ since $\varrho_{n}=O(n^{\alpha-1})$ (see the proof of Lemma 4.1 with $\alpha$ replaced by $-\alpha$). And, by the Taylor expansion formulas, we can get

[TABLE]

which implies the WSGL operator is stable and consistent. Actually, the second-order WSGL method is constructed by the first-order fractional BDF with specially designed weights. For some other numerical methods constructed by the fractional BDF but with higher-order convergence rates, see [6, 34, 35, 36].

**Example 4. **(Further discussion on the WSGL method)

A interesting consideration is that with the structure of the generating function (2.39), i.e., with

[TABLE]

where $m$ is a positive integer, can we propose a new formula that is convergent (to $I^{\alpha}$) of order $3$? Actually, with the following coefficients

[TABLE]

one can easily check that the $SCQ$ with $\omega(\xi)$ defined by (2.44) is stable and consistent of order $3$. Furthermore, since $\alpha+2\theta$ and $m$ cannot be zero by (2.44), it seems that there is no superconvergence point for the $(1-\xi)^{-\alpha}$ to devise a formula of convergence order $3$. For the application of (2.43), by considering the condition $\lambda_{1}\geq|\lambda_{2}|$ (which is derived by the condition-$\omega$), we can take

[TABLE]

We remark here that (2.43) generalizes (2.37) from some aspects that, if we take $m=1$, under the condition $(\ref{S.36})$ (which reduces to (2.38)), then the generating function $\omega(\xi)$ defined by (2.43) with (2.44) is exactly the same as the $\omega(\xi)$ defined by (2.37).

**Example 5. **(Generalized BDF2-$\theta$ method)

In this example we consider the generalized BDF2-$\theta$ method that generalizes the work of Liu et al. [8] and Ding et al. [11]. Define the generating function $\omega(\xi)$ by

[TABLE]

which can be reformulated as

[TABLE]

The generating function (2.46) has been proposed by Gunarathna et al. [18] for the fractional derivative, i.e., for $\alpha<0$. For the case $\theta$ takes integers and other higher-order formulas, see [32]. By careful derivation, we can check that the $SCQ$ with $\omega(\xi)$ defined in (2.46) is convergent of order $2$. If we take $\alpha=-1$, (2.46) is reduced to the BDF-$\theta$ method (see [8]) with the generating function as,

[TABLE]

If we take $\theta=-\frac{1}{2}$, we get a $SCQ$ that approximates $f$ at the node $x_{n+\frac{1}{2}}$ (see [11] with $\alpha$ replaced by $-\alpha$),

[TABLE]

There are two points deserve discussion:

a) The generating function (2.46) also implies that some superconvergence properties at node $x_{n+\frac{\alpha}{2}}$ by taking $\theta=-\frac{\alpha}{2}$.

b) By taking $\theta=-\alpha$, we can obtain a shorter or simpler generating function

[TABLE]

and we call the corresponding $SCQ$ the fractional central difference method, for the reason that if we take $\alpha=-1$, then the generating function $\omega(\xi)=\frac{1}{2}(1-\xi^{2})$ implies the classical central difference scheme

[TABLE]

However, the application of the fractional central difference method is limited. See table 4 in section 3 and the notation therein.

3 Stability regions

In this section we pay special attention on the stability regions (see Definition 3.3) of the $SCQ$ when applied to the fractional equations after omitting the starting part. This work is motivated by the fact that the condition-$\omega$ can not guarantee a $SCQ$ is A($\frac{\pi}{2}$)-stable (see Definition 3.3). As is well known, for some problems numerical schemes of A($\frac{\pi}{2}$)-stable are superior to the others. To be specific, we analyse the following two models

[TABLE]

and

[TABLE]

respectively. The operator ${}_{C}D_{0,x}^{\alpha}$ denotes the Caputo fractional derivative operator of order $\alpha$, which is defined by

[TABLE]

With the relation

[TABLE]

we can approximate the Caputo type derivatives by the $SCQ$ developed in the Sec. 2.

Define $z:=\lambda h^{\alpha}$, the numerical schemes for equation (3.1) and (3.2) are as follows

[TABLE]

respectively, where $\theta_{j}$ are defined by (2.30), and

[TABLE]

Before examine the stability region of the numerical schemes (3.5), we state some properties of the analytic solutions of the model equation (3.1) and (3.2).

Definition 3.1

(See [12]) The two-parameter Mittag-Leffler function $E_{\alpha,\beta}(x)$ is defined by,

[TABLE]

For simplicity, we define $E_{\alpha}(x):=E_{\alpha,1}(x)$.

Lemma 3.2

For $|\arg\lambda-\pi|<(1-\frac{\alpha}{2})\pi$, the solutions of (3.1) and (3.2) satisfy

[TABLE]

**Proof. ** For the equation (3.2), which is a special case of (1.2) in [2], see the proof therein. We mainly focus on the differential equation (3.1). By Laplace transform, we can express the analytic solution by the Mittag-Leffler function $E_{\alpha}(x)$,

[TABLE]

Note that for $\alpha\in(0,2)$, we have the asymptotic property (see Theorem 1.6, p.35, [12]) for $E_{\alpha}(x)$ that with $\frac{\alpha}{2}\pi<\mu<\min\{\pi,\alpha\pi\}$,

[TABLE]

where $C$ is a real constant. By replacing $x$ of $E_{\alpha}(x)$ with $\lambda x^{\alpha}$, the proof of the lemma is completed.

The Lemma 3.2 naturally leads to the following definition (see also [2]),

Definition 3.3

The stability region $S$ of a $SCQ$ is the set of all complex $z=\lambda h^{\alpha}$ for which the numerical solutions $y^{n}$ of (3.5) satisfy

[TABLE]

Further, we call a numerical method $A(\delta)$-stable if $S$ contains the sector $|\arg z-\pi|<\delta$.

Theorem 3.4

The stability regions of the numerical schemes (3.5) are

[TABLE]

respectively.

Proof. The technique used in this theorem is the same as the Theorem 2.1 in [2]. We omit the proof here.

We conclude this section by illustrating the stability regions of some $SCQ$s and make some notations. Generally, we require that there exists a positive $x_{0}$ such that the interval $(-x_{0},0)$ is contained in the stability region $S$.

Table 1 illustrates the stability regions for the shift-generalized Newton-Gregory formula of order $2$. For the Scheme-I with $\theta\leq\frac{1}{2}$, the method is $A(\frac{\pi}{2})$-stable. If $\theta$ satisfies $\frac{1}{2}<\theta<\frac{\alpha+1}{2}$, the method is conditional stable which means for arbitrary $\lambda<0$, the step size $h$ must be small enough ($\lambda h^{\alpha}>-x_{0}$) such that for fixed $h$, the solution $y^{n}$ tends to [math] as $n$ tends to infinity. If $\theta=\frac{\alpha+1}{2}$, then for any $x_{0}>0$, $(-x_{0},0)\cap S=\varnothing$, which means for any step size $h>0$, $y^{n}$ oscillates or blows up as $n$ tends to infinity. For the scheme-II, by condition-$\omega$, we require $\theta\leq\frac{1-\alpha}{2}$, in which the method is $A(\frac{\pi}{2})$-stable. If $\theta$ exceeds this value, we will get a conditional stable method, provided $\theta<\frac{1}{2}$.

In Table 2, we briefly depict the stability regions for WSGL operators, with the choice of pairs $(p,q)=(0,-2)$ or $(1,-2)$ for the Scheme-I, and $(p,q)=(0,-2)$ or $(-1,-2)$ for the Scheme-II. One can find that only with $(p,q)=(1,-2)$ the method is conditional stable and for other three pairs the method is $A(\frac{\pi}{2})$-stable.

In Table 3, we show some stability regions for the $SCQ$ with $\omega(\xi)$ defined in (2.46). For the Scheme-I under the condition-$\omega$, we can conclude that with $\theta\leq\min\{\alpha,\frac{1}{2}\}$, the method is $A(\frac{\pi}{2})$-stable. For the Scheme-II, the method is $A(\frac{\pi}{2})$-stable if $-\alpha\leq\theta\leq\frac{1}{2}$. An interesting phenomenon is that when $\theta>\frac{1}{2}$, there exists a $x_{0}>0$, such that for any $h>x_{0}$, the solution $y^{n}\to 0(n\to\infty)$.

In Table 4 we analyse the reason why the application for the fractional central difference method is limited. For the Scheme-I, one can find that the solution $y^{n}$ blows up for $\alpha>\frac{1}{2}$, since for any $x_{0}>0$, $(-x_{0},0)\cap S=\varnothing$. See also Fig. 1 for the numerical solutions of the Scheme-I with $\lambda=-15$, $h=0.2$. Nonetheless, the method is quite suitable for the Scheme-II which is $A(\frac{\pi}{2})$-stable for any $\alpha\in(0,1)$.

4 Applications to PDEs

In this section we apply a class of novel numerical schemes to the time-fractional diffusion equation (4.1). To facilitate the numerical analysis below, we require the scheme satisfies the assumptions of the discrete fractional Gr$\ddot{\text{o}}$nwall inequality [13]. In the spacial direction, the finite element method is adopted to formulate the fully discrete scheme. The equation is

[TABLE]

where $T>0$, ${}_{C}D_{0,t}^{\alpha}(\alpha\in(0,1))$ is the Caputo fractional derivative operator defined by (3.3), and $\Omega$ is a bounded interval. The temporal interval is separated uniformly with $0=t_{0}<t_{1}<\cdots<t_{N}=T$ and define $\tau=T/N$ as the step size of the time mesh. For $\Omega$, we denote $\mathcal{T}_{h}$ as its shape-regular and quasi-uniform triangulation with the mesh size $h$. Introduce the subspace $X_{h}$ of $H_{0}^{1}(\Omega)$ defined by

[TABLE]

where $\mathbb{P}_{r}(x)$ is defined as the set of polynomials (of $x$) with the degree at most r$$(r\in\mathbb{Z}^{+}).

4.1 Numerical scheme

We approximate the equation (4.1) at $t_{n-\theta}$ by the shift-generalized Newton-Gregory formula of order 2 (which is defined by (2.36) with $\alpha$ replaced by $-\alpha$), with the generating function

[TABLE]

Considering the relation (3.4), we define $v(t):=u(t)-u_{0}$. By Theorem 2.9, we can easily formulate the temporal semi-discrete scheme for (4.1) as follows

[TABLE]

where $\Delta v^{n-\theta}:=(1-\theta)\Delta v^{n}+\theta\Delta v^{n-1}$ and $g^{n-\theta}:=g(t_{n-\theta})$.

With the space $X_{h}$, then the fully discrete scheme of (4.1) is to find $V^{n}:[0,T]\to X_{h}$ such that

[TABLE]

holds for any $\chi_{h}\in X_{h}$.

If we omit the starting part in (4.5), we get the scheme

[TABLE]

4.2 Stability and convergence analysis

To analyse the stability and convergence of the scheme (4.1), we employ the key tool of the discrete fractional Gr$\ddot{\text{o}}$nwall inequality [13]. In the following analysis, we omit the starting part from our scheme as results will not be affected under suitable conditions. To meet the assumptions of the discrete fractional Grönwall inequality, we first in Lemmas 4.1-4.2 prove some properties of the coefficients of $\omega(\xi)$ defined by (4.3).

For clarity, we denote by $\|\cdot\|$ the norm of the $L^{2}(\Omega)$ space, and by $\|\cdot\|_{s}$ the norm of the Sobolev space $H^{s}(\Omega)$.

Lemma 4.1

The coefficients $\omega_{k}$ defined by the generating function $\omega(\xi)$ in (4.3) satisfy

(i) $\omega_{0}=1+\frac{\alpha}{2}-\theta>0$, $\omega_{k}\leq 0,~{}\text{for }k=1,2,\cdots$, if \theta\in\big{[}\frac{\alpha}{2}-\frac{1-\alpha}{1+\alpha},\frac{\alpha}{2}\big{]}

(ii) $\sum_{k=0}^{n}\omega_{k}>0$, and $\sum_{k=0}^{n}\omega_{k}=\frac{n^{-\alpha}}{\Gamma(1-\alpha)}+O(n^{-1-\alpha})$.

**Proof. ** Denote by $\kappa_{k}$ the coefficients of $(1-\xi)^{\alpha}$. Then the coefficients $\kappa_{k}$ satisfy (see [1, 14]),

(a) $\kappa_{n}=O(n^{-\alpha-1})$,

(b) $\kappa_{0}=1,\kappa_{n}<0,|\kappa_{n+1}|<|\kappa_{n}|,n=1,2,\cdots$,

(c) $\displaystyle\sum_{k=0}^{n-1}\kappa_{k}=\frac{\Gamma(n-\alpha)}{\Gamma(1-\alpha)\Gamma(n)}=\frac{n^{-\alpha}}{\Gamma(1-\alpha)}+O(n^{-1-\alpha}),n=1,2,\cdots.$

On the other hand, we can express $\omega_{k}$ by $\kappa_{k}$ as

[TABLE]

and, $\omega_{0}=1+\frac{\alpha}{2}-\theta$. If $k=1$, we have

[TABLE]

When $k\geq 2$, by (4.7) and (b) the condition $\omega_{k}\leq 0$ is equivalent to the following inequality

[TABLE]

With the estimates

[TABLE]

we can easily check that (4.10) holds for \theta\in\big{[}\frac{\alpha}{2}-\frac{1-\alpha}{1+\alpha},\frac{\alpha}{2}\big{]}.

By the definition of $\omega(\xi)$, we know that $\sum_{k=0}^{\infty}\omega_{k}=0$, hence by (i), we can get $\sum_{k=0}^{n}\omega_{k}>0$. Now combining (a), (c) with (4.7), we have

[TABLE]

The proof of the lemma is completed.

Lemma 4.2

Let $A_{k}:=\frac{1}{\tau^{\alpha}}\sum_{j=0}^{k}\omega_{j}$. Then, $A_{k}\geq A_{k+1}>0(k\geq 0)$. Furthermore, there exists a positive $\pi_{A}$ such that

[TABLE]

**Proof. ** By careful derivation one can see that to prove the existence of $\pi_{A}$, it is sufficient to demonstrate

[TABLE]

where $\Theta_{n}=\frac{(n+1)^{1-\alpha}-n^{1-\alpha}}{\sum_{j=0}^{n}\omega_{j}}$. By Lemma 4.1 (ii), $\Theta_{n}$ is an increasing sequence, and the limit $\lim_{n\to\infty}\Theta_{n}$ exists. The proof is completed.

The complementary discrete convolution kernels $P_{j}$ is essential to the development of the inequality (4.17), which is defined as the coefficients of the function $P(\xi)=\tau^{\alpha}\omega(\xi)^{-1}$. We remark that this definition implies

[TABLE]

Actually, the sequence $\{A_{k}\}$ defined in Lemma 4.2 is the coefficients of the function $A(\xi)=\frac{\omega(\xi)}{\tau^{\alpha}(1-\xi)}$. On the other hand, $\sum_{j=0}^{n}P_{n-j}A_{j}$ is the $n$th coefficient of $P(\xi)A(\xi)=\frac{1}{1-\xi}=1+\xi+\xi^{2}+\cdots$, which means (4.14).

Now based on Lemmas 4.1-4.2, for coefficients $A_{k}$ defined in Lemma 4.2, we have the following discrete fractional Gr$\ddot{\text{o}}$nwall inequality, see [13].

Lemma 4.3

Let $\{\phi^{n}\}_{n=1}^{\infty}$ and $\{\eta_{l}\}_{l=0}^{\infty}$ be given nonnegative sequences. Let $0\leq\theta<1$. Assume further the series $\sum_{l=0}^{\infty}\eta_{l}$ is bounded with $\Lambda$, i.e., $\sum_{l=0}^{\infty}\eta_{l}<\Lambda$, and that the time step size satisfies

[TABLE]

Then for any nonnegative sequence $\{v^{k}\}_{k=0}^{N}$ such that

[TABLE]

it holds that

[TABLE]

Remark 4.4

A careful examination shows that $\sum_{j=1}^{n}P_{n-j}$ is bounded by some constant $C$ which is independent of $n$. Hence, (4.17) can be simplified to

[TABLE]

We reformulate the numerical scheme (4.6) with the coefficients $A_{k}$ defined in Lemma 4.2 as the following

[TABLE]

Now we are in a position to analyse the stability and error estimates of the scheme by the same approaches taken in [13], hence, we omit the proof of the following theorems.

Theorem 4.5

Define $U^{n}:=V^{n}+u_{0h}$ and $V^{n}$ is the numerical solution of (4.6), $u_{0h}$ is a proper approximation of $u_{0}$. Then with $\tau$ satisfying (4.15), we have the stability estimates as

[TABLE]

where the constant $C$ is independent of $\tau$ and $h$.

Theorem 4.6

Let $V^{n}$ be the numerical solution of (4.6). Define $U^{n}:=V^{n}+u_{0h}$, where $u_{0h}$ is a proper approximation of $u_{0}$. Then $U^{n}$ is a numerical solution of (4.1). Suppose $u\in C([0,T];H^{r+1}(\Omega)\cap H_{0}^{1}(\Omega))\cap C^{3}([0,T];L^{2}(\Omega))$, then with $\tau$ satisfying (4.15), we have the error estimates as

[TABLE]

where the constant $C$ is independent of $\tau$ and $h$.

4.3 Numerical experiments

In this subsection, we implement some numerical experiments to further confirm our convergence estimates (4.21). To this end, define the error $E(\tau,h):=\max_{0\leq n\leq N}\|u^{n}-U^{n}\|$. The convergence rate are derived by the formula

[TABLE]

Let $\Omega=(0,1)$, $T=1$ and divide the interval $\Omega$ as $0=x_{0}<x_{1}<\cdots<x_{M}=1$ with the mesh size $h=1/M$. Define the finite element space $X_{h}$ as the set of piecewise linear polynomials with $r=1$. The exact solution is taken as

[TABLE]

and the corresponding source term $g$ is

[TABLE]

One can see that there is some weak singularity for $u$ at initial value, hence, by Theorem 2.9 we add the starting part to obtain a second-order convergence in time. In the following tables, we denoted by $E_{c}(\tau,h)$ the errors of the scheme with starting part, and by $E_{o}(\tau,h)$ the errors without the starting part.

In Table 5, with fixed fine space mesh size $h=\frac{1}{5000}$, we take different $\alpha$ and for each $\alpha$ we choose some $\theta$ which satisfy \max\big{\{}0,\frac{\alpha}{2}-\frac{1-\alpha}{1+\alpha}\big{\}}\leq\theta\leq\frac{\alpha}{2}. Now one can easily check that with the starting part, we have obtained a second-order convergence in time. For the scheme without the starting part, the convergence rate is much lower, especially when $\alpha$ is small.

In Table 6, we fix the time mesh size $\tau=\frac{1}{1000}$, and similarly, compare the spacial convergence rate for the scheme with or without the starting part. We can see that the convergence rate is optimal when the starting part is added, and for those $\alpha$ close to zero (which means there is a stronger singularity for the solution), the rate becomes much lower if the starting part is omitted.

5 Conclusion

In this paper, the shifted convolution quadrature theory is developed based on the extensible framework established by Lubich. The definition of consistency is generalized and the equivalent theorem is established for the $SCQ$ theory. The fertility of the generalized framework is demonstrated by being able to transform generating functions of the $CQ$ to those of the $SCQ$ (Theorem 2.11), develop the shift-generalized Newton-Gregory formula (Corollary 2.14), include as many existing approximation formula with integer convergence rate as possible (Sec. 2) and easy to design new formulas with desired structure (Example 4). We shall point out that the allowable structure of a generating function is not limited to the case $(\ref{S.2.1})$, as Lubich [1] said ”condition (2.14) can be considerably relaxed, however, the class (2.14) is probably large enough for all practical applications”. Indeed, all formulas mentioned in this paper have further assumed $p_{4}(\xi)\equiv 1$ in (2.14), and it is not hard to construct formulas with $p_{4}(\xi)\neq 1$ identically which are stable and consistent. Another merit of the $SCQ$ is that it inherits the correction techniques from the $CQ$, by which the high-order convergence rate can be obtained numerically despite the weak regularity for solutions at initial node. To further explore the stability properties of the schemes proposed, we analyse the impact of the parameter $\theta$ on the stable regions with changing $\theta$. We emphasize that the shift parameter $\theta$ plays an essential role in the stability of a numerical scheme, and careful examination for the choice of $\theta$ is of vital importance before employing new schemes. For some special designed generating functions, we apply the resulted schemes to PDEs with the technique of the fractional Grönwall inequality to analyse the stability and convergence. The results of the numerical experiments further confirm our theory analysis.

Authors think there are at least two approaches in our future work: a) Propose some simple generating functions that meet with the techniques already developed for the numerical analysis for different PDEs, and b) develop new techniques that are suitable for as many different types of generating functions as possible.

Acknowledgements

The work of the first author was supported in part by the NSFC grant 11661058. The work of the third author was supported in part by the NSFC grant 11761053, the NSF of Inner Mongolia 2017MS0107, and the program for Young Talents of Science and Technology in Universities of Inner Mongolia Autonomous Region NJYT-17-A07. The work of the fourth author was supported in part by grants NSFC 11871092 and NSAF U1530401.