## The unified theory of shifted convolution quadrature for fractional calculus ∗

Yang Liu 1 ∗ , Baoli Yin , Hong Li , Zhimin Zhang 1 1 2 3 ,

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

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

3 Department 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 α 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 θ to cover as many numerical schemes that approximate the operator I α with an integer convergence rate as possible. The constraint on the parameter θ 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 [9,10,22-31,33]. In this paper, we are particularly concerned about the RiemannLiouville calculus operator I α which is defined by

$$I ^ { \alpha } f ( x ) = \frac { 1 } { \Gamma ( \alpha ) } \int _ { 0 } ^ { x } ( x - s ) ^ { \alpha - 1 } f ( s ) \mathrm d s, \quad \text{for } \Re \alpha > 0, & & ( 1. 1 )$$

and the Riemann-Liouville differential operator I -α ( /Rfractur α &gt; 0) is defined by

$$I ^ { - \alpha } f ( x ) = \frac { d ^ { n } } { d x ^ { n } } I ^ { n - \alpha } f ( x ) = \frac { 1 } { \Gamma ( n - \alpha ) } \frac { d ^ { n } } { d x ^ { n } } \int _ { 0 } ^ { x } ( x - s ) ^ { n - \alpha - 1 } f ( s ) d s,$$

where n = /ceilingleft/Rfractur α /ceilingright . When α = 0, we set I 0 = I , the identity operator.

In 1986, Lubich [1] developed the convolution quadrature ( CQ ) theory to approximate the RiemannLiouville calculus with arbitrary α ∈ C at the node x = nh ,

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

∗ Corresponding author.

Email addresses: mathliuyang@imu.edu.cn,

Manuscript submitted to Journal 2019

where h is the mesh size of a uniform grid, ω j denotes the convolution weight and w n,j is the starting weight. For brevity, define x k = kh . For /Rfractur α ≤ 0, the theory requires that f ( j ) (0) = 0( j = 0 , · · · , -/ceilingleft/Rfractur α /ceilingright ). When α takes integers, I α coincides with the traditional integral operator ( α &gt; 0) and differential operator ( α &lt; 0), and the identity operator ( α = 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θ method, see [8, 37]) which approximates the first derivative at the node x n -1 2 and its conterpart for fractional derivative I -α ( α ∈ (0 , 1)) approximated at the node x n -α 2 , with a second-order convergence rate (see the ω ξ ( ) 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¨nwald u 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 , , · · · , 6 are excluded from the CQ theory as well, since all above methods approximate I α at a shifted node.

In this paper, we generalize the CQ theory to cover the methods mentioned above by introducing a shifted parameter θ 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 α ( α ∈ C ) at x = ( n -θ h ) ,

$$( 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 ( j h ).$$

The contributions of the paper are as follows

- · Develop systematic approaches to approximate the fractional calculus I α f ( x )( α ∈ C ) at node x n -θ without assumptions on the regularity of f ( x ). Unlike most papers concerning the approximation at a shifted node with the parameter θ , that pay little attention on the choice of θ , we explore the criterion for this choice from the aspect of generating functions. See Sec. 2.
- · Examine impacts of the parameter θ on the absolute stability regions for different numerical schemes. With a proper shifted parameter θ , we can get a A-stable method, which is superior to others for some problems. See Sec. 3.
- · Construct some new approximation methods (for I α f ( x n -θ )) based on known approximation methods (for I α 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.
- · Generalize the correction technique of the CQ theory to the SCQ theory by introducing a parameter θ , 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.
- · 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 θ with which the fractional Gr¨nwall inequality o 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 θ . 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¨nwall inequality. o 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

$$( 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 ( j h ), \quad x = ( n - \theta ) h,$$

as the shifted convolution quadrature ( SCQ ). Denote Ω α h f ( x ) = h α ∑ n j =0 ω f x j ( + ( θ -j ) h ) as the convolution part and S α h f = h α ∑ s j =0 w n,j f ( jh ) as the starting part of (2.1), respectively. Define the convolution quadrature error by

$$E _ { h } ^ { \alpha } = \Omega _ { h } ^ { \alpha } - I ^ { \alpha }.$$

For the sequence of convolutions weights { ω j } ∞ j =0 we associate with a generating power series ω ξ ( ) = ∑ ∞ j =0 ω ξ j j , and viceversa.

We note that if f is continuous and g is locally integrable, and Ω α h is extended for x ≥ 0 such that

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

then Ω α h commutes with convolution Ω α h ( f ∗ g ) = (Ω α h f ) ∗ g . Hence, the convolution error E α h = Ω α h -I α satisfies (for /Rfractur α ≤ 0 , we require f ( j ) (0) = 0 , j = 0 , · · · , -/ceilingleft/Rfractur α /ceilingright )

$$E _ { h } ^ { \alpha } ( f * g ) = ( E _ { h } ^ { \alpha } f ) * g.$$

Another key property of E α h is the homogeneity:

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

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 α ∈ C ,

Definition 2.1 A SCQ is stable for I α if the convolution weights is

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

Definition 2.2 A SCQ is consistent of order p for I α if the generating function of { ω j } ∞ j =0 satisfies

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

Definition 2.3 A SCQ is convergent of order p to I α if

/negationslash

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

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 θ = 0, i.e., the CQ is consistent of order p for I α if h α ω e ( -h ) = 1 + O h ( p ) (see Definition 2.2, [1]).

Remark 2.4 Combining the homogeneity of E α h and the definition of convergent (2.8), we can get

/negationslash

$$( E _ { h } ^ { \alpha } t ^ { \beta - 1 } ) ( x _ { n } ) = O ( x _ { n } ^ { \alpha - 1 } h ^ { \beta } ) + O ( x _ { n } ^ { \alpha + \beta - 1 - p } h ^ { p } ), \quad \text{for all} \ \beta \in \mathbb { C }, \beta \neq 0, - 1, - 2, \cdots, \quad \ ( 2. 9 )$$

which means for those small β that α + β -1 -p &lt; 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 t α k h -1 )(1) = O h ( k )+ O h ( p ) for k = 1 2 3 , , , · · · , 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 ],

$$e _ { h } ( x ) \coloneqq ( E _ { h } ^ { \alpha } e ^ { t - x } ) ( x ) & = h ^ { \alpha } \sum _ { 0 \leq j h \leq x + \theta h } \omega _ { j } e ^ { ( \theta - j ) h } - ( I ^ { \alpha } e ^ { t - x } ) ( x ) \\ & = h ^ { \alpha } e ^ { \theta h } \sum _ { 0 \leq j h \leq x + \theta h } \omega _ { j } e ^ { - j h } - ( I ^ { \alpha } e ^ { t - x } ) ( x ). \\ \lambda _ { \lambda } - \lambda _ { \lambda } \quad \lambda _ { \lambda } \quad \lambda _ { \lambda }.$$

Let x → ∞ , the expression h e α θh ∑ 0 ≤ jh ≤ x + θh ω e j -jh tends to h e α θh ω 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 /pi1 ( ξ ) for a consistent CQ is of the form (see (3.6) in [1])

$$\varpi ( \xi ) = ( 1 - \xi ) ^ { - \alpha } \left [ c _ { 0 } + c _ { 1 } ( 1 - \xi ) + c _ { 2 } ( 1 - \xi ) ^ { 2 } + \dots + c _ { N - 1 } ( 1 - \xi ) ^ { N - 1 } + ( 1 - \xi ) ^ { N } \tilde { r } ( \xi ) \right ]. \quad ( 2. 1 1 )$$

where ˜( r ξ ) is holomorphic at 1, and constants c j = γ ′ j where γ ′ j are defined by (2.21). We argue that the generating function ω ξ ( ) for SCQ can be expressed similarly by (2.11) with a different definition of the coefficients c i that depend on θ :

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

$$c _ { i } = \gamma _ { i }, \ \ \text{ for } i = 0, 1, \cdots, p - 1,$$

where γ i are the coefficients of

$$\sum _ { i = 0 } ^ { \infty } \gamma _ { i } ( 1 - \xi ) ^ { i } = \xi ^ { \theta } \left ( \frac { - \ln \xi } { 1 - \xi } \right ) ^ { - \alpha }.$$

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

Remark 2.7 We can construct ω ξ ( ) by polynomial functions. Assume ω ξ ( ) has the form

$$\omega ( \xi ) = \left [ \frac { p _ { 1 } ( \xi ) } { p _ { 2 } ( \xi ) } \right ] ^ { \alpha } \frac { p _ { 3 } ( \xi ) } { p _ { 4 } ( \xi ) },$$

where p i ( ξ ) are polynomial functions, and p 2 (1) = 0 . Denote by D the unit disc | ξ | &lt; 1 , and by ¯ D the closed unit disc in the complex plane. Then, for a stable SCQ , which means ω n = O n ( α -1 ) , it holds that ω ξ ( ) is analytic in D . Hence, ω ξ ( ) can be written as

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

/negationslash where ξ j are distinct with | ξ j | = 1 , ν ξ ( ) is analytic on ¯ D , and ν ξ ( j ) = 0 , α j = 0 , -1 , -2 , · · · . It can be shown that /Rfractur α j ≤ /Rfractur α is equivalent to (2.6), see [1]. We limit the choice of the parameter θ by the Conditionω ,

/negationslash

$$C o n d i t i o n \prec \omega \colon \omega ( \xi ) \text{ is analytic without zeros in } D, \text{ and } \Re \alpha _ { j } \leq \Re \alpha.$$

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 ω ξ ( ) satisfying the conditionω 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 , I α h f ( x ) approximates I α f ( x ) uniformly for bounded x .

Theorem 2.9 (See Theorem 2.4 in [1]) Suppose ω ξ ( ) is the generating function of a convergent SCQ . Then we have: For any β = 0 , -1 , -2 , · · · ,

/negationslash

- (i) there exist starting weights w n,j = O n ( α -1 )( n ≥ 0 , j = 0 , · · · , s ) such that for any function

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

the SCQ satisfies

$$I _ { h } ^ { \alpha } f ( x ) - I ^ { \alpha } f ( x ) = O ( h ^ { p } )$$

uniformly for x ∈ [ a, b ] with 0 &lt; a &lt; b &lt; ∞ .

(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

$$( I _ { h } ^ { \alpha } t ^ { q + \beta + 1 } ) ( t _ { n - \theta } ) = ( I ^ { \alpha } t ^ { q + \beta + 1 } ) ( t _ { n - \theta } )$$

for all integer q ≥ 0 such that /Rfractur ( q + β -1) ≤ p -1 . Note that the estimate (2.18) requires that x is bounded away from 0 . The starting weights w n,j in (ii) require additionally those /Rfractur ( q + α + β -1) &lt; 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 /pi1 p ( ξ ) ( p is the convergence order) for the CQ :

- · (Theorem 2.6, [1]) For an implicit linear multistep method which is stable and consistent of order p with the characteristic polynomials ρ and σ , if the zeros of σ ξ ( ) have absolute value less than 1,

then /pi1 p ( ξ ) = [ σ (1 /ξ /ρ ) (1 /ξ ) ] α . Some special cases include (see [1,22])

$$\sum _ { \ } w p ( \Sigma ) = \lfloor 0 \, ( 1 / \Sigma ) / \rho ( 1 / \Sigma ) \rfloor \, \cdots \, \text{ some special cases micro } \{ \sec \, \lfloor 1, \, \sqrt { 4 } \rfloor } \\ \quad \text{The fractional trapezoidal rule: } \varpi _ { 2 } ( \xi ) & = \frac { ( 1 + \xi ) ^ { \alpha } } { 2 ^ { \alpha } ( 1 - \xi ) ^ { \alpha } }, \quad \Re \alpha \geq 0. \\ \quad \text{The fractional BDF-p: } \varpi _ { p } ( \xi ) & = \left [ \sum _ { i = 1 } ^ { p } \frac { 1 } { i } ( 1 - \xi ) ^ { i } \right ] ^ { - \alpha }, \quad 1 \leq p \leq 6. \\ \quad \text{The fractional BT-\theta method: } \varpi _ { 2 } ( \xi ) & = \left [ \frac { 1 - \theta + \theta \xi } { ( 3 / 2 - \theta ) - ( 2 - 2 \theta ) \xi + ( 1 / 2 - \theta ) \xi ^ { 2 } } \right ] ^ { \alpha }, \\ & \quad \text{for } \theta \in ( - \infty, 1 / 2 ), \alpha < 0, \text{ or } \theta \in ( - \infty, 1 / 2 ), \alpha \geq 0. \\ - \, \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$}$} \,, \text{$\text{$\text{$\text{$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$}$} \,, \text{$\text{$\text{$$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$
$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$}\end{array} } }$$

- · The generalized Newton-Gregory formula (see [1])

$$\varpi _ { p } ( \xi ) = ( 1 - \xi ) ^ { - \alpha } \left [ \gamma _ { 0 } ^ { \prime } + \gamma _ { 1 } ^ { \prime } ( 1 - \xi ) + \dots + \gamma _ { p - 1 } ^ { \prime } ( 1 - \xi ) ^ { p - 1 } \right ], \quad & ( 2. 2 1 ) \\. \quad & - \cong \dots \dots \dots - \alpha$$

with the coefficients γ ′ i defined by ∑ ∞ i =0 γ ′ i (1 -ξ ) i = ( ln ξ ξ -1 ) -α . · The fractional BNθ method (see [22])

$$\varpi _ { 2 } ( \xi ) = \frac { 1 - \alpha \theta + \alpha \theta \xi } { \left [ ( 3 / 2 - \theta ) - ( 2 - 2 \theta ) \xi + ( 1 / 2 - \theta ) \xi ^ { 2 } \right ] ^ { a } }, \ \theta \in ( - \infty, 1 ), \ \alpha \in \mathbb { R }, \ \alpha \theta \leq \frac { 1 } { 2 }. \quad \text{(2.22)} \\ \intertext { T h o r o m $ 2 $ 1 $ 1 $ S u n n a e $ \equiv ( \varepsilon ) $ i s t h e n e r e n t i n a $ f u n t i o n $ o $ C O $ w i t h $ c o m e n e n e $ o r d e r $ n $ D e f i n e }$$

Theorem 2.11 Suppose /pi1 p ( ξ ) is the generating function of a CQ with convergence order p . Define

$$\omega ( \xi ) = \varpi _ { p } ( \xi ) \left [ \gamma _ { 0 } ^ { \prime \prime } + \gamma _ { 1 } ^ { \prime \prime } ( 1 - \xi ) + \dots + \gamma _ { p - 1 } ^ { \prime \prime } ( 1 - \xi ) ^ { p - 1 } \right ], & & ( 2. 2 3 ) \\..$$

where the coefficients γ ′′ i satisfy

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

Then the SCQ with convolution weights generated by ω ξ ( ) is convergent of order p .

Proof. For the generating function /pi1 p ( ξ ) = ∑ ∞ j =0 /pi1 j ξ j , we know that /pi1 n = O n ( α -1 ) by the stability of the CQ . Hence, we have ω n = O n ( α -1 ). By replacing ξ with e -h in (2.24), we get

$$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 } ).$$

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

$$h ^ { \alpha } \varpi _ { p } ( e ^ { - h } ) = 1 + O ( h ^ { p } ),$$

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

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

Then by Theorem 2.8 the proof is completed.

A special case for (2.23) that is of vital importance is when α = 0. For any /pi1 p ( ξ ) defined by (2.20)-(2.22), taking α = 0, we obtain the generating function by (2.23) for the approximation of I 0 f ( x ) = f ( x ) at node x n -θ . The coefficients γ ′′ i in (2.24) can be formulated as

$$\gamma _ { i } ^ { \prime \prime } = ( - 1 ) ^ { i } \frac { \Gamma ( \theta + 1 ) } { \Gamma ( i + 1 ) \Gamma ( \theta - i + 1 ) }.$$

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

/negationslash

Corollary 2.12 Let f ( x ) = x β -1 with β = 0 , -1 , -2 , · · · . Then we have the following approximation formula for f ( x n -θ ) of order p ,

$$f ( x _ { n - \theta } ) - \sum _ { j = 0 } ^ { p - 1 } \theta _ { j } f ( x _ { n - j } ) = O ( x _ { n - \theta } ^ { \beta - 1 - p } h ^ { p } ) + O ( x _ { n - \theta } ^ { - 1 } h ^ { \beta } ),$$

where the weights θ j satisfy

$$\sum _ { i = 0 } ^ { p - 1 } \theta _ { i } \xi ^ { i } = \sum _ { i = 0 } ^ { p - 1 } \gamma _ { i } ^ { \prime \prime } ( 1 - \xi ) ^ { i },$$

with coefficients γ ′′ i defined in (2.28).

Remark 2.13 The approximation formula (2.29) holds for f ( x ) = x β -1 g x ( ) with g x ( ) sufficiently differentiable as well. For p = 2 3 , , we can obtain two popular formulas

$$f ( x _ { n - \theta } ) & \approx ( 1 - \theta ) f ( x _ { n } ) + \theta f ( x _ { n - 1 } ), \\ f ( x _ { n - \theta } ) & \approx \frac { 1 } { 2 } ( 1 - \theta ) ( 2 - \theta ) f ( x _ { n } ) + \theta ( 2 - \theta ) f ( x _ { n - 1 } ) + \frac { 1 } { 2 } \theta ( \theta - 1 ) f ( x _ { n - 2 } ),$$

see also (17), (18) in [4]. Considering the conditionω , we limit θ for (2.31) to satisfying θ ≤ 1 2 and θ ≤ 1 -√ 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 α ) of order p ,

$$\omega ( \xi ) = ( 1 - \xi ) ^ { - \alpha } \left [ \gamma _ { 0 } + \gamma _ { 1 } ( 1 - \xi ) + \dots + \gamma _ { p - 1 } ( 1 - \xi ) ^ { p - 1 } \right ],$$

where γ i are the coefficients of

$$\sum _ { i = 0 } ^ { \infty } \gamma _ { i } ( 1 - \xi ) ^ { i } = \xi ^ { \theta } \left ( \frac { - \ln \xi } { 1 - \xi } \right ) ^ { - \alpha }.$$

One can find out that the shift-generalized Newton-Gregory formula reduces to (2.21) when θ = 0. We conclude this section by further exploring some numerical methods that approximate f ( x ) at node x n -θ 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 α ∈ R .

Example 1. (Generalized shifted Gr¨nwald formula) u

Expanding (2.33), we can easily derive that

$$\gamma _ { 0 } = 1, \ \gamma _ { 1 } = - \frac { \alpha } { 2 } - \theta, \ \gamma _ { 2 } = \frac { 1 } { 8 } ( \alpha + 2 \theta ) ^ { 2 } - \frac { 1 } { 2 4 } ( 5 \alpha + 1 2 \theta ).$$

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

$$( 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 ( j h ), \quad \text{with } \omega ( \xi ) = ( 1 - \xi ) ^ { - \alpha }. \quad \text{(2.35)}$$

Now by assuming that f is sufficiently smooth (hence, the starting part can be omitted), and that θ takes nonnegative integers, the relation (2.35) reduces to the shifted Gr¨nwald formula (see [16] with u α replaced by -α ). Recently, Chen et al. [21] has applied the Gr¨nwald formula ( u θ = 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 ω ξ ( ) for the shift-generalized NewtonGregory formula of order p = 2 3 can be formulated as ,

$$\omega ( \xi ) = ( 1 - \xi ) ^ { - \alpha } \left [ 1 - \left ( \frac { \alpha } { 2 } + \theta \right ) ( 1 - \xi ) \right ] \quad \text{of order $2$, with $\theta \leq \frac{1-\alpha}{2}$,} \\ \omega ( \xi ) = ( 1 - \xi ) ^ { - \alpha } \left \{ 1 - \left ( \frac { \alpha } { 2 } + \theta \right ) ( 1 - \xi ) + \left [ \frac { 1 } { 8 } ( \alpha + 2 \theta ) ^ { 2 } - \frac { 1 } { 2 4 } ( 5 \alpha + 1 2 \theta ) \right ] ( 1 - \xi ) ^ { 2 } \right \} \quad \text{of order $3$.}$$

Both of the formula were also proposed by Dimitrov [4] by using the theory developed in [17]. An obvious conclusion is that ω ξ ( ) defined by (2.36) which is of order 2 shows some superconvergence when θ = -α 2 (see also [5], with α replaced by -α ), since generally a CQ with the generating function (1 -ξ ) -α approximates I α 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 ω ξ ( ) defined by (2.36) which is of order 3, that from the discussion above, a SCQ with ω ξ ( ) defined by

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

is convergent of order 3 provided that

$$\theta = \frac { 1 - \alpha } { 2 } - \frac { 1 } { 2 } \sqrt { 1 - \frac { \alpha } { 3 } }, \quad \text{for } \alpha \leq 3.$$

One can check that this choice of θ satisfies the conditionω .

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

Assume two integers p &gt; q , and let θ = , consider the following generating function p

$$\omega ( \xi ) = ( 1 - \xi ) ^ { - \alpha } \left [ \frac { - \alpha - 2 q } { 2 ( p - q ) } + \frac { 2 p + \alpha } { 2 ( p - q ) } \xi ^ { p - q } \right ],$$

whose coefficients can be formulated as

$$\omega _ { k } = \frac { - \alpha - 2 q } { 2 ( p - q ) } \varrho _ { k } + \frac { 2 p + \alpha } { 2 ( p - q ) } \varrho _ { k - ( p - q ) },$$

where /rho1 k are the coefficients of (1 -ξ ) -α and define /rho1 k = 0 for k &lt; 0. Considering the conditionω , we assume ( p, q ) satisfies

$$p + q + \alpha \leq 0.$$

One can easily check that the SCQ with ω ξ ( ) defined in (2.39) is the WSGL operator (see [6], with α replaced by -α ), which is convergent of order 2 by Theorem 2.8. Actually, by (2.40), we have ω n = O n ( α -1 ) since /rho1 n = O n ( α -1 ) (see the proof of Lemma 4.1 with α replaced by -α ). And, by the

Taylor expansion formulas, we can get

$$\text{for expansion formulas, we can get} \\ h ^ { \alpha } e ^ { p h } \omega ( e ^ { - h } ) & = h ^ { \alpha } e ^ { p h } ( 1 - e ^ { - h } ) ^ { - \alpha } \left [ \frac { - \alpha - 2 q } { 2 ( p - q ) } + \frac { 2 p + \alpha } { 2 ( p - q ) } e ^ { h ( q - p ) } \right ] \\ & = h ^ { \alpha } ( 1 - e ^ { - h } ) ^ { - \alpha } \left [ \frac { - \alpha - 2 q } { 2 ( p - q ) } e ^ { h p } + \frac { 2 p + \alpha } { 2 ( p - q ) } e ^ { h q } \right ] \\ & = \left [ 1 - \frac { h } { 2 } + O ( h ^ { 2 } ) \right ] ^ { - \alpha } \left [ 1 - \frac { \alpha } { 2 } h + O ( h ^ { 2 } ) \right ] \\ & = \left [ 1 + \frac { \alpha } { 2 } h + O ( h ^ { 2 } ) \right ] \left [ 1 - \frac { \alpha } { 2 } h + O ( h ^ { 2 } ) \right ] \\ & = 1 + O ( h ^ { 2 } ), \\ \text{h implies the WSGL operator is stable and consistent.  Actually, the second-order WSGL method}$$

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-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

$$\omega ( \xi ) = ( 1 - \xi ) ^ { - \alpha } ( \lambda _ { 1 } + \lambda _ { 2 } \xi ^ { m } ), \quad \lambda _ { 1 } \geq | \lambda _ { 2 } |,$$

where m is a positive integer, can we propose a new formula that is convergent (to I α ) of order 3? Actually, with the following coefficients

$$\lambda _ { 1 } = 1 - \lambda _ { 2 }, \quad \lambda _ { 2 } = \frac { \alpha + 2 \theta } { 2 m }, \quad m = 1 + \frac { \alpha + 2 \theta } { 2 } - \frac { 5 \alpha + 1 2 \theta } { 6 ( \alpha + 2 \theta ) },$$

one can easily check that the SCQ with ω ξ ( ) defined by (2.44) is stable and consistent of order 3. Furthermore, since α +2 θ and m cannot be zero by (2.44), it seems that there is no superconvergence point for the (1 -ξ ) -α to devise a formula of convergence order 3. For the application of (2.43), by considering the condition λ 1 ≥ | λ 2 | (which is derived by the conditionω ), we can take

$$\theta = \frac { m - \alpha } { 2 } - \frac { 1 } { 2 } \sqrt { m ^ { 2 } - \frac { \alpha } { 3 } }, \quad \text{for } \alpha \leq 3 m ^ { 2 }.$$

We remark here that (2.43) generalizes (2.37) from some aspects that, if we take m = 1, under the condition (2 45) (which reduces to (2.38)), then the generating function . ω ξ ( ) defined by (2.43) with (2.44) is exactly the same as the ω ξ ( ) defined by (2.37).

## Example 5. (Generalized BDF2θ method)

In this example we consider the generalized BDF2θ method that generalizes the work of Liu et al. [8] and Ding et al. [11]. Define the generating function ω ξ ( ) by

$$\omega ( \xi ) = \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 }, \text{ \ with } \alpha ( \alpha + \theta ) \geq 0,$$

which can be reformulated as

$$\omega ( \xi ) = \left ( \frac { 3 } { 2 } + \frac { \theta } { \alpha } \right ) ^ { - \alpha } ( 1 - \xi ) ^ { - \alpha } \left ( 1 - \frac { \alpha + 2 \theta } { 3 \alpha + 2 \theta } \xi \right ) ^ { - \alpha }.$$

The generating function (2.46) has been proposed by Gunarathna et al. [18] for the fractional derivative, i.e., for α &lt; 0. For the case θ takes integers and other higher-order formulas, see [32]. By careful derivation, we can check that the SCQ with ω ξ ( ) defined in (2.46) is convergent of order 2. If we take

α = -1, (2.46) is reduced to the BDFθ method (see [8]) with the generating function as,

$$\omega ( \xi ) = \frac { 3 - 2 \theta } { 2 } - ( 2 - 2 \theta ) \xi + \frac { 1 - 2 \theta } { 2 } \xi ^ { 2 }.$$

If we take θ = -1 2 , we get a SCQ that approximates f at the node x n + 1 2 (see [11] with α replaced by -α ),

$$\omega ( \xi ) = \left ( \frac { 3 \alpha - 1 } { 2 \alpha } - \frac { 2 \alpha - 1 } { \alpha } \xi + \frac { \alpha - 1 } { 2 \alpha } \xi ^ { 2 } \right ) ^ { - \alpha }.$$

There are two points deserve discussion:

a) The generating function (2.46) also implies that some superconvergence properties at node x n + α 2 by taking θ = -α 2 .

b) By taking θ = -α , we can obtain a shorter or simpler generating function

$$\omega ( \xi ) = 2 ^ { \alpha } ( 1 - \xi ^ { 2 } ) ^ { - \alpha },$$

and we call the corresponding SCQ the fractional central difference method, for the reason that if we take α = -1, then the generating function ω ξ ( ) = 1 2 (1 -ξ 2 ) implies the classical central difference scheme

$$f ^ { \prime } ( x _ { n - 1 } ) \approx \frac { f ( x _ { n } ) - f ( x _ { n - 2 } ) } { 2 }.$$

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ω can not guarantee a SCQ is A( π 2 )-stable (see Definition 3.3). As is well known, for some problems numerical schemes of A( π 2 )-stable are superior to the others. To be specific, we analyse the following two models

$$C D ^ { \alpha } _ { 0, x } y ( x ) = \lambda y ( x ), \ \ y ( 0 ) = y _ { 0 }, \ \ \alpha \in ( 0, 1 ),$$

and

$$I ^ { \alpha } y ( x ) = \frac { 1 } { \lambda } y ( x ), \quad \alpha \in ( 0, 1 ),$$

respectively. The operator C D α 0 ,x denotes the Caputo fractional derivative operator of order α , which is defined by

$$c D _ { 0, x } ^ { \alpha } y ( x ) = \frac { 1 } { \Gamma ( 1 - \alpha ) } \int _ { 0 } ^ { x } \frac { y ^ { \prime } ( s ) } { ( x - s ) ^ { \alpha } } d s, \quad \alpha \in ( 0, 1 ).$$

With the relation

$$_ { C } D _ { 0, x } ^ { \alpha } y ( x ) = I ^ { - \alpha } ( y ( x ) - y _ { 0 } ), \quad \alpha \in ( 0, 1 ),$$

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

Define z := λh α , the numerical schemes for equation (3.1) and (3.2) are as follows

$$\text{Scheme-I: } \sum _ { k = 0 } ^ { n } \omega _ { n - k } y ^ { k } & = z \sum _ { j = 0 } ^ { p - 1 } \theta _ { j } y ^ { n - j } + f ( x _ { n - \theta } ), \\ \text{Scheme-II: } z \sum _ { k = 0 } ^ { n } \omega _ { n - k } y ^ { k } & = \sum _ { j = 0 } ^ { p - 1 } \theta _ { j } y ^ { n - j },$$

respectively, where θ j are defined by (2.30), and

$$f ( x _ { n - \theta } ) \coloneqq h ^ { \alpha } ( I ^ { - \alpha } y _ { 0 } ) ( x _ { n - \theta } ) = \frac { y _ { 0 } } { ( n - \theta ) ^ { \alpha } \Gamma ( 1 - \alpha ) }.$$

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 α,β ( x ) is defined by,

$$E _ { \alpha, \beta } ( x ) \coloneqq \sum _ { j = 0 } ^ { \infty } \frac { x ^ { j } } { \Gamma ( j \alpha + \beta ) }, \, \text{for those $x$ that make the series converge.}$$

For simplicity, we define E α ( x ) := E α, 1 ( x ).

Lemma 3.2 For | arg λ -π | &lt; (1 -α 2 ) π , the solutions of (3.1) and (3.2) satisfy

$$y ( x ) \to 0 \ \ a s \ \ x \to \infty.$$

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 α ( x ),

$$y ( x ) = y _ { 0 } E _ { \alpha } ( \lambda x ^ { \alpha } ).$$

Note that for α ∈ (0 , 2), we have the asymptotic property (see Theorem 1.6, p.35, [12]) for E α ( x ) that with α 2 π &lt; µ &lt; min { π, απ } ,

$$| E _ { \alpha } ( x ) | \leq \frac { C } { 1 + | x | }, \ \mu \leq | \arg x | \leq \pi, \ \ | x | \geq 0,$$

where C is a real constant. By replacing x of E α ( x ) with λx α , 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 = λh α for which the numerical solutions y n of (3.5) satisfy

$$y ^ { n } \to 0 \ \ a s \ \ n \to \infty.$$

Further, we call a numerical method A δ ( ) -stable if S contains the sector | arg z -π | &lt; δ .

Theorem 3.4 The stability regions of the numerical schemes (3.5) are

$$\text{For Scheme} - I \colon \quad \mathbb { C } \, \left \{ z \, \colon z = \omega ( \xi ) \Big / \sum _ { j = 0 } ^ { p - 1 } \theta _ { j } \xi ^ { j }, | \xi | \leq 1 \right \}, \\ \text{For Scheme} - I I \colon \quad \mathbb { C } \, \left \{ z \, \colon z = \sum _ { j = 0 } ^ { p - 1 } \theta _ { j } \xi ^ { j } \Big / \omega ( \xi ), | \xi | \leq 1 \right \},$$

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 θ ≤ 1 2 , the method is A ( π 2 )-stable. If θ satisfies 1 2 &lt; θ &lt; α +1 2 , the method is conditional stable which means for arbitrary λ &lt; 0, the step size h must be small enough ( λh α &gt; -x 0 ) such that for fixed h , the solution y n tends to 0 as n tends to infinity. If θ = α +1 2 , then for any x 0 &gt; 0, ( -x , 0 0) ∩ S = ∅ , which means for any step size h &gt; 0, y n oscillates or blows up as n tends to infinity. For the scheme-II, by conditionω , we require θ ≤ 1 -α 2 , in which the method is A ( π 2 )-stable. If θ exceeds this value, we will get a conditional stable method, provided θ &lt; 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 ( π 2 )-stable.

In Table 3, we show some stability regions for the SCQ with ω ξ ( ) defined in (2.46). For the Scheme-I under the conditionω , we can conclude that with θ ≤ min { α, 1 2 } , the method is A ( π 2 )-stable. For the Scheme-II, the method is A ( π 2 )-stable if -α ≤ θ ≤ 1 2 . An interesting phenomenon is that when θ &gt; 1 2 , there exists a x 0 &gt; 0, such that for any h &gt; x 0 , the solution y n → 0( n →∞ ).

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 α &gt; 1 2 , since for any x 0 &gt; 0, ( -x , 0 0) ∩ S = ∅ . See also Fig. 1 for the numerical solutions of the Scheme-I with λ = -15, h = 0 2. . Nonetheless, the method is quite suitable for the Scheme-II which is A ( π 2 )-stable for any α ∈ (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¨nwall inequality [13]. o In the spacial direction, the finite element method is adopted to formulate the fully discrete scheme. The equation is

$$C D _ { 0, t } ^ { \alpha } u ( t ) & = \Delta u + g ( x, t ) \quad \text{ for } ( x, t ) \in \Omega \times ( 0, T ] \\ u ( x, 0 ) & = u _ { 0 } ( x ) \quad \text{ for } x \in \Omega \quad \text{ (4.1)} \\ u ( x, t ) & = 0 \quad \text{ for } ( x, t ) \in \partial \Omega \times ( 0, T ),$$

where T &gt; 0, C D α 0 ,t ( α ∈ (0 , 1)) is the Caputo fractional derivative operator defined by (3.3), and Ω is a bounded interval. The temporal interval is separated uniformly with 0 = t 0 &lt; t 1 &lt; · · · &lt; t N = T and define τ = T/N as the step size of the time mesh. For Ω, we denote T h as its shape-regular and quasi-uniform triangulation with the mesh size h . Introduce the subspace X h of H 1 0 (Ω) defined by

$$X _ { h } = \{ \chi \in H _ { 0 } ^ { 1 } ( \Omega ) \, \colon \chi | _ { e } \in \mathbb { P } _ { r } ( x ), e \in \mathcal { T } _ { h } \},$$

where P r ( x ) is defined as the set of polynomials (of x ) with the degree at most r r ( ∈ Z + ).

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

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01136v3-figure-af0df08d-3a23-45a7-9220-ee5ffd08ad41.jpg)

## 4.1 Numerical scheme

We approximate the equation (4.1) at t n -θ by the shift-generalized Newton-Gregory formula of order 2 (which is defined by (2.36) with α replaced by -α ), with the generating function

$$\omega ( \xi ) = ( 1 - \xi ) ^ { \alpha } \left [ 1 + \left ( \frac { \alpha } { 2 } - \theta \right ) ( 1 - \xi ) \right ] \text{ \ for \ max } \{ 0, \frac { \alpha } { 2 } - \frac { 1 - \alpha } { 1 + \alpha } \} \leq \theta \leq \frac { \alpha } { 2 }. \text{ \quad \ \ } ( 4. 3 )$$

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

$$\tau ^ { - \alpha } \sum _ { k = 0 } ^ { n } \omega _ { n - k } v ^ { k } = \Delta v ^ { n - \theta } + g ^ { n - \theta } + \Delta u _ { 0 } + S _ { h } ^ { 0 } \Delta v - S _ { h } ^ { - \alpha } v,$$

where ∆ v n -θ := (1 -θ )∆ v n + ∆ θ v n -1 and g n -θ := g t ( n -θ ). With the space X h , then the fully discrete scheme of (4.1) is to find V n : [0 , T ] → X h such that

$$\tau ^ { - \alpha } \sum _ { k = 0 } ^ { n } \omega _ { n - k } ( V ^ { k }, \chi _ { h } ) + ( \nabla V ^ { n - \theta }, \nabla \chi _ { h } ) = ( g ^ { n - \theta } + \Delta u _ { 0 }, \chi _ { h } ) - S _ { h } ^ { 0 } ( \nabla V, \nabla \chi _ { h } ) - S _ { h } ^ { - \alpha } ( V, \chi _ { h } ) \quad ( 4. 5 )$$

holds for any χ h ∈ X h .

Table 2: The WSGL operator with q = -2 , θ = p

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01136v3-figure-71617ad9-dcef-416a-8e20-8ff36a47b519.jpg)

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

$$\tau ^ { - \alpha } \sum _ { k = 0 } ^ { n } \omega _ { n - k } ( V ^ { k }, \chi _ { h } ) + ( \nabla V ^ { n - \theta }, \nabla \chi _ { h } ) = ( g ^ { n - \theta } + \Delta u _ { 0 }, \chi _ { h } ).$$

## 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¨nwall inequality [13]. In the following analysis, we omit the starting part from our scheme o 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 o ω ξ ( ) defined by (4.3).

For clarity, we denote by ‖ · ‖ the norm of the L 2 (Ω) space, and by ‖ · ‖ s the norm of the Sobolev space H s (Ω).

Lemma 4.1 The coefficients ω k defined by the generating function ω ξ ( ) in (4.3) satisfy

$$\text{emim} \alpha \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{cases} \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text#$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text}$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text $\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text.$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text'$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text`$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text<$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text1$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text"$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$$$

Proof. Denote by κ k the coefficients of (1 -ξ ) α . Then the coefficients κ k satisfy (see [1, 14]),

- (a) κ n = O n ( - -α 1 ),
- (b) κ 0 = 1 , κ n &lt; , κ 0 | n +1 | &lt; κ | n | , n = 1 2 , , · · · ,

$$( c ) \sum _ { \substack { k = 0 \\ \text{$k=0$} } } ^ { \nu \nu _ { 0 } - \ast, \nu _ { n } \succ \nu, | \nu _ { n + 1 | } \succ | \nu _ { n | }, \nu - \ast, \varepsilon, \quad, } \\ ( c ) \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 ω k by κ k as

$$\omega _ { k } = \left ( 1 + \frac { \alpha } { 2 } - \theta \right ) \kappa _ { k } + \left ( \theta - \frac { \alpha } { 2 } \right ) \kappa _ { k - 1 }, \ k = 1, 2, \cdots,$$

and, ω 0 = 1 + α 2 -θ . If k = 1, we have

$$\omega _ { 1 } = \theta ( \alpha + 1 ) - \frac { \alpha } { 2 } ( \alpha + 3 ) \leq - \alpha < 0.$$

Table 3: The SCQ with ω ξ ( ) in (2.46)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01136v3-figure-f84f1684-2fc0-4594-9108-67cda917603f.jpg)

When k ≥ 2, by (4.7) and (b) the condition ω k ≤ 0 is equivalent to the following inequality

$$\frac { \alpha } { 2 } - \theta \leq \frac { \kappa _ { k } } { \kappa _ { k - 1 } - \kappa _ { k } }.$$

With the estimates

$$\frac { \kappa _ { k - 1 } } { \kappa _ { k } } = \frac { k } { k - \alpha - 1 } = \frac { 1 } { 1 - \frac { \alpha + 1 } { k } } \leq \frac { 1 } { 1 - \frac { \alpha + 1 } { 2 } },$$

we can easily check that (4.10) holds for θ ∈ [ α 2 -1 -α 1+ α , α 2 ] . By the definition of ω ξ ( ), we know that ∑ ∞ k =0 ω k = 0, hence by (i), we can get ∑ n k =0 ω k &gt; 0. Now combining (a), (c) with (4.7), we have

$$\sum _ { k = 0 } ^ { n - 1 } \omega _ { k } = \sum _ { k = 0 } ^ { n - 1 } \kappa _ { k } + \left ( \frac { \alpha } { 2 } - \theta \right ) \kappa _ { n - 1 } = \frac { n ^ { - \alpha } } { \Gamma ( 1 - \alpha ) } + O ( n ^ { - 1 - \alpha } ).$$

The proof of the lemma is completed.

Lemma 4.2 Let A k := 1 τ α ∑ k j =0 ω j . Then, A k ≥ A k +1 &gt; 0( k ≥ 0) . Furthermore, there exists a positive π A such that

$$A _ { n - k } \geq \frac { 1 } { \pi _ { A } \tau } \int _ { t _ { k - 1 } } ^ { t _ { k } } \frac { ( t _ { n } - s ) ^ { - \alpha } } { \Gamma ( 1 - \alpha ) } d s, \quad \text{for any } n \geq k \geq 1.$$

Table 4: The fractional central difference method (2.50)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01136v3-figure-50347f2d-4d11-4def-a698-19eeb860277a.jpg)

Proof. By careful derivation one can see that to prove the existence of π A , it is sufficient to demonstrate

$$\pi _ { A } \geq \frac { 1 } { \Gamma ( 2 - \alpha ) } \max _ { n \geq 0 } \Theta _ { n } > 0,$$

where Θ n = ( n +1) 1 -α -n 1 -α ∑ n j =0 ω j . By Lemma 4.1 (ii), Θ n is an increasing sequence, and the limit lim n →∞ Θ 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 ξ ( ) = τ α ω ξ ( ) -1 . We remark that this definition implies

$$\sum _ { j = 0 } ^ { n } P _ { n - j } A _ { j } \equiv 1 \quad \text{for any } n \geq 0.$$

Actually, the sequence { A k } defined in Lemma 4.2 is the coefficients of the function A ξ ( ) = ω ξ ( ) τ α (1 -ξ ) . On the other hand, ∑ n j =0 P n -j A j is the n th coefficient of P ξ A ξ ( ) ( ) = 1 1 -ξ = 1 + ξ + ξ 2 + · · · , 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¨nwall inequality, see [13]. o

Lemma 4.3 Let { φ n } ∞ n =1 and { η l } ∞ l =0 be given nonnegative sequences. Let 0 ≤ θ &lt; 1 . Assume further

Figure 1: h = 0 2, . λ = -15.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01136v3-figure-e6130720-5440-4210-b514-80e165459a87.jpg)

the series ∑ ∞ l =0 η l is bounded with Λ , i.e., ∑ ∞ l =0 η l &lt; Λ , and that the time step size satisfies

$$\tau & \leq \frac { 1 } { \alpha \sqrt { 2 \pi _ { A } \Gamma ( 2 - \alpha ) \Lambda } }. \\ \intertext { c o u n o d } \tilde { L } k \lambda N \quad \text{cosh t h at }$$

Then for any nonnegative sequence { v k } N k =0 such that

$$\sum _ { k = 1 } ^ { n } A _ { n - k } \left [ ( v ^ { k } ) ^ { 2 } - ( v ^ { k - 1 } ) ^ { 2 } \right ] \leq \sum _ { k = 1 } ^ { n } \eta _ { n - k } ( v ^ { k - \theta } ) ^ { 2 } + v ^ { n - \theta } \phi ^ { n } \quad \text{for $1\leq n\leq N$},$$

it holds that

$$v ^ { n } \leq 2 E _ { \alpha } ( 2 \pi _ { A } \Lambda t _ { n } ^ { \alpha } ) \left ( v ^ { 0 } + \sum _ { j = 1 } ^ { n } P _ { n - j } \phi ^ { j } \right ) \ \text{ for } 1 \leq n \leq N.$$

Remark 4.4 A careful examination shows that ∑ n j =1 P n -j is bounded by some constant C which is independent of n . Hence, (4.17) can be simplified to

$$v ^ { n } \leq C ( v ^ { 0 } + \max _ { 1 \leq j \leq n } \phi ^ { j } ) \quad \text{for $1\leq n\leq N$.}$$

We reformulate the numerical scheme (4.6) with the coefficients A k defined in Lemma 4.2 as the following

$$\sum _ { k = 1 } ^ { n } A _ { n - k } ( V ^ { k } - V ^ { k - 1 }, \chi _ { h } ) + ( \nabla V ^ { n - \theta }, \nabla \chi _ { h } ) = & ( g ^ { n - \theta } + \Delta u _ { 0 }, \chi _ { h } ).$$

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 0 h and V n is the numerical solution of (4.6), u 0 h is a proper approximation of u 0 . Then with τ satisfying (4.15), we have the stability estimates as

$$\| U ^ { n } \| \leq \| u _ { 0 h } \| + C \| \Delta u _ { 0 } \| + C \max _ { 1 \leq j \leq n } \| g ^ { n } \|,$$

where the constant C is independent of τ and h .

Theorem 4.6 Let V n be the numerical solution of (4.6). Define U n := V n + u 0 h , where u 0 h is a proper approximation of u 0 . Then U n is a numerical solution of (4.1). Suppose u ∈ C ([0 , T ]; H r +1 (Ω) ∩ H 1 0 (Ω)) ∩ C 3 ([0 , T ]; L 2 (Ω)) , then with τ satisfying (4.15), we have the error estimates as

$$\| U ^ { n } - u ^ { n } \| & \leq C ( \tau ^ { 2 } + h ^ { r + 1 } ),$$

where the constant C is independent of τ 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 τ, h ( ) := max 0 ≤ ≤ n N ‖ u n -U n ‖ . The convergence rate are derived by the formula

$$\text{temporal order } = \log _ { 2 } \frac { E ( 2 \tau, h ) } { E ( \tau, h ) }, \ \text{spatial order } = \log _ { 2 } \frac { E ( \tau, 2 h ) } { E ( \tau, h ) }. \quad \text{(4.22)}$$

Let Ω = (0 1), , T = 1 and divide the interval Ω as 0 = x 0 &lt; x 1 &lt; · · · &lt; 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

$$u ( x, t ) = ( 1 + t ^ { \alpha } + t ^ { 2 \alpha } + t ^ { 3 } ) \sin ( 2 \pi x ),$$

and the corresponding source term g is

$$g ( x, t ) = \sin ( 2 \pi x ) \left [ \Gamma ( 1 + \alpha ) + 4 \pi ^ { 2 } ( 1 + t ^ { \alpha } + t ^ { 2 \alpha } + t ^ { 3 } ) + \frac { 2 \alpha t ^ { \alpha } \Gamma ( 2 \alpha ) } { \Gamma ( \alpha + 1 ) } + \frac { 6 t ^ { 3 - \alpha } } { \Gamma ( 4 - \alpha ) } \right ]. \quad \quad ( 4. 2 4 )$$

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 ( τ, h ) the errors of the scheme with starting part, and by E o ( τ, h ) the errors without the starting part.

In Table 5, with fixed fine space mesh size h = 1 5000 , we take different α and for each α we choose some θ which satisfy max { 0 , α 2 -1 -α 1+ α } ≤ θ ≤ α 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 α is small.

In Table 6, we fix the time mesh size τ = 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 α 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

Table 5: The temporal convergence rate with h = 1 5000

| α   | θ    | τ    |   E c ( τ, h ) | rate   |   E o ( τ, h ) | rate    |
|-----|------|------|----------------|--------|----------------|---------|
| 0.1 |      | 1/10 |    2.35436e-05 | -      |    0.00287659  | -       |
| 0.1 | 0    | 1/20 |    5.47471e-06 | 2.1045 |    0.00279338  | 0.0423  |
| 0.1 |      | 1/40 |    1.13371e-06 | 2.2717 |    0.00271399  | 0.0416  |
| 0.1 |      | 1/80 |    2.55465e-07 | 2.1499 |    0.0026392   | 0.0403  |
| 0.1 |      | 1/10 |    0.000959982 | -      |    0.0422249   | -       |
| 0.1 | 0.05 | 1/20 |    0.000243373 | 1.9798 |    0.0381734   | 0.1455  |
| 0.1 |      | 1/40 |    6.1552e-05  | 1.9833 |    0.0344973   | 0.1461  |
| 0.1 |      | 1/80 |    1.57177e-05 | 1.9694 |    0.0311936   | 0.1452  |
|     |      | 1/10 |    0.00172824  | -      |    0.00564791  | -       |
|     | 0.1  | 1/20 |    0.000438829 | 1.9776 |    0.00256403  | 1.1393  |
|     |      | 1/40 |    0.00011075  | 1.9864 |    0.000453468 | 2.4993  |
| 0.5 |      | 1/80 |    2.80539e-05 | 1.9810 |    0.00089205  | -0.9761 |
|     |      | 1/10 |    0.00322333  | -      |    0.0200719   | -       |
|     | 0.2  | 1/20 |    0.00081512  | 1.9835 |    0.0128495   | 0.6435  |
|     |      | 1/40 |    0.00020512  | 1.9905 |    0.00775159  | 0.7291  |
|     |      | 1/80 |    5.16829e-05 | 1.9887 |    0.00428721  | 0.8545  |
|     |      | 1/10 |    0.00470189  | -      |    0.00608467  | -       |
|     | 0.4  | 1/20 |    0.00121112  | 1.9569 |    0.00152565  | 1.9958  |
|     |      | 1/40 |    0.00030719  | 1.9791 |    0.000699968 | 1.1241  |
| 0.9 |      | 1/80 |    7.75617e-05 | 1.9857 |    0.000601301 | 0.2192  |
|     |      | 1/10 |    0.00490787  | -      |    0.00631974  | -       |
|     | 0.45 | 1/20 |    0.00125984  | 1.9619 |    0.00158234  | 1.9978  |
|     |      | 1/40 |    0.000319039 | 1.9814 |    0.000396145 | 1.9980  |
|     |      | 1/80 |    8.04831e-05 | 1.9870 |    0.000145886 | 1.4412  |

/negationslash 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 (2 14), . 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 ( ξ ) ≡ 1 in (2.14), and it is not hard to construct formulas with p 4 ( ξ ) = 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 θ on the stable regions with changing θ . We emphasize that the shift parameter θ plays an essential role in the stability of a numerical scheme, and careful examination for the choice of θ 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 o 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

Table 6: The spacial convergence rate with τ = 1 1000

| α   | θ    | h         | E c ( τ, h )            | rate     | E ( τ, h )                                        | rate          |
|-----|------|-----------|-------------------------|----------|---------------------------------------------------|---------------|
| 0.3 |      | 1/10 1/20 | 9.90558E-02             | -        | o 9.90559E-02 2.49007E-02 6.23387E-03 4.29525E-03 | -             |
|     | 0    |           | 2.49005E-02             | 1.9921   |                                                   | 1.9921        |
|     |      | 1/40      | 6.23369E-03             | 1.9980   |                                                   | 1.9980        |
|     |      | 1/80 1/10 | 1.55895E-03             | 1.9995   |                                                   | 0.5374        |
|     |      |           | 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       |
|     | 0.3  | 1/10      | 9.85495E-02             | -        | 9.85515E-02                                       | -             |
|     |      | 1/20      | 2.47681E-02             | 1.9924   | 2.47686E-02                                       | 1.9924        |
| 0.8 |      | 1/40      | 6.20048E-03 1.55092E-03 | 1.9980   | 6.20066E-03                                       | 1.9980        |
|     |      | 1/80 1/10 | 9.85493E-02             | 1.9993 - | 1.55102E-03                                       | 1.9992        |
|     | 0.4  |           | 2.47681E-02             | 1.9924   | 9.85515E-02                                       | -             |
|     |      | 1/20 1/40 | 6.20052E-03             | 1.9980   | 2.47687E-02 6.20072E-03                           | 1.9924 1.9980 |
|     |      | 1/80      | 1.55098E-03             | 1.9992   | 1.55108E-03                                       | 1.9992        |

Mongolia Autonomous Region NJYT-17-A07. The work of the fourth author was supported in part by grants NSFC 11871092 and NSAF U1530401.

## References

- [1] C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal., 1986, 17(3): 704-719.
- [2] C. Lubich, A stability analysis of convolution quadraturea for Abel-Volterra integral equations, IMA J. Numer. Anal., 1986, 6(1): 87-101.
- [3] G. Dahlquist, Convergence and stability in the numerical integration of ordinary differential equations, Math. Scand., 1956, 33-53.
- [4] Y. Dimitrov, Numerical approximations for fractional differential equations, Journal of Fractional Calculus and Applications, 2015, 5(3S)(22): 1-45.
- [5] G.H. Gao, H.W. Sun, Z.Z. Sun, Stability and convergence of finite difference schemes for a class of time-fractional sub-diffusion equations based on certain superconvergence, J. Comput. Phys., 2015, 280: 510-528.
- [6] W.Y. Tian, H. Zhou, W.H. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Math. Comput., 2015, 84: 1703-1727.
- [7] F.H. Zeng, Z.Q. Zhang, G.E. Karniadakis, Second-order numerical methods for multi-term fractional differential equations: smooth and non-smooth solutions, Computer Methods in Applied Mechanics and Engineering, 2017, 327: 478-502.
- [8] Y. Liu, Y.W. Du, H. Li, F.W. Liu, Y.J. Wang, Some second-order θ schemes combined with finite element method for nonlinear fractional Cable equation, Numer. Algor., 2019, 80(2): 533-555. https://doi.org/10.1007/s11075-018-0496-0
- [9] J.C. Li, Y.Q. Huang, Y.P. Lin, Developing finite element methods for Maxwell's equations in a Cole-Cole dispersive medium, SIAM J. Sci. Comput., 2011, 33(6): 3153-3174.
- [10] D. Baffet, J.S. Hesthaven, High-order accurate local schemes for fractional differential equations, J. Sci. Comput., 2017, 70(1): 355-385.
- [11] H.F. Ding, C.P. Li, Q. Yi, A new second-order midpoint approximation formula for Riemann-Liouville derivative: algorithm and its application, IMA J. Appl. Math., 2017, 82(5): 909-944.
- [12] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Elsevier, 1998.
- [13] H.L. Liao, W. McLean, J.W. Zhang, A discrete Gr¨ onwall inequality with applications to numerical schemes for subdiffusion problems, SIAM J. Numer. Anal., 2019, 57(1): 218-237.
- [14] F.H. Zeng, C.P. Li, F.W. Liu, I. Turner, The use of finite difference /element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput., 2013, 35(6): A2976-A3000.

| [15]   | Z.B. Wang, S. Vong, Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation, J. Comput. Phys., 2014, 277: 1-15.                     |
|--------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [16]   | M. M. Meerschaert, C. Tadjeran, Finite difference approximations for fractional advection-dispersion flow equations, J. Comput. Appl. Math., 2004, 172(1): 65-77.                                              |
| [17]   | C. Tadjeran, M. M. Meerschaert, H. P. Scheffer, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys., 2006, 213: 205-213.                                   |
| [18]   | W.A. Gunarathna, H.M. Nasir, W.B. Daundasekera, An explicit form for higher order approximations of fractional derivatives, Appl. Numer. Math., 2019.                                                          |
| [19]   | Y.W. Du, Y. Liu, H. Li, Z.C. Fang, S. He, Local discontinuous Galerkin method for a nonlinear time- fractional fourth-order partial differential equation, J. Comput. Phys., 2017, 344: 108-126.               |
| [20]   | Y. Liu, Y.W. Du, H. Li, J.F. Wang, A two-grid finite element approximation for a nonlinear time- fractional Cable equation, Nonlinear Dyn., 2016, 85: 2535-2548.                                               |
| [21]   | H. Chen, F. Holland, M. Stynes, An analysis of the Gr¨nwald-Letnikov u scheme for initial-value problems with weakly singular solutions, Appl. Numer. Math., 2019, 139: 52-61.                                 |
| [22]   | B.L. Yin, Y. Liu, H. Li, Z.M. Zhang, Two families of novel second-order fractional numerical formulas and their applications to fractional differential equations, arXiv preprint arXiv:1906.01242, 2019.      |
| [23]   | B.L. Yin, Y. Liu, H. Li, Z.M. Zhang, Finite element methods based on two families of novel second-order numerical formulas for the fractional Cable model, 2019.                                               |
| [24]   | B.T. Jin, B.Y. Li, Z. Zhou, Correction of high-order BDF convolution quadrature for fractional evolution equations, SIAM J. Sci. Comput., 2017, 39(6): A3129-A3152.                                            |
| [25]   | N. Ford, Y.B. Yan, An approach to construct higher order time discretisation schemes for time fractional partial differential equations with nonsmooth data, Fract. Calc. Appl. Anal., 2017, 20(5): 1076-1105. |
| [26]   | X.J. Yang, General Fractional Derivatives: Theory, Methods and Applications, Chapman and Hall/CRC, 2019.                                                                                                       |
| [27]   | M. Dehghan, M. Abbaszadeh, W. Deng, Fourth-order numerical method for the space-time tempered fractional diffusion-wave equation, Appl. Math. Lett., 2017, 73: 120-127.                                        |
| [28]   | L.B. Feng, F.W. Liu, I. Turner, Q.Q. Yang, P.H. Zhuang, Unstructured mesh finite difference/finite element method for the 2D time-space Riesz fractional diffusion equation on irregular convex domains,       |
|        | Appl. Math. Model., 2018, 59: 441-463. D. Baleanu, K. Diethelm, E. Scalas and J.J. Trujillo, Fractional Calculus: Models and Numerical                                                                         |
| [29]   | Meth- ods, 3 of Series on Complexity, Nonlinearity and Chaos,World Scientific Publishing, New York, NY.USA, 2012.                                                                                              |
| [30]   | M. Zhao, A.J. Cheng, H. Wang, A preconditioned fast Hermite finite element method for space-fractional diffusion equations, Discrete Contin. Dyn. Syst. Ser. B, 2017, 22(9): 3529-3545.                        |
| [31]   | C.P. Li, F.H. Zeng, Numerical methods for fractional calculus, Chapman and Hall/CRC, 2015.                                                                                                                     |
| [32]   | C.P. Li, M. Cai, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations: revisited, Numer. Func. Anal. Opt., 2017, 38(7): 861-890.                                       |
| [33]   | X.J. Yang, F. Gao, H.M. Srivastava, A new computational approach for solving nonlinear local fractional PDEs, J. Comput. Appl. Math., 2018, 339: 285-296.                                                      |
| [34]   | Y. Liu, M. Zhang, H. Li, J.C. Li, High-order local discontinuous Galerkin method combined with WSGD- approximation for a fractional subdiffusion equation, Comput. Math. Appl., 2017, 73(6), 1298-1314.        |
| [35]   | C.C. Ji, Z.Z. Sun, A high-order compact finite difference scheme for the fractional sub-diffusion equation, J. Sci. Comput., 2015, 64(3): 959-985.                                                             |
| [36]   | M.H. Chen, W.H. Deng, Fourth order difference approximations for space Riemann-Liouville derivatives based on weighted and shifted Lubich difference operators, Commun. Comput. Phys., 2014, 16(2): 516- 540.  |
| [37]   | B.L. Yin, Y. Liu, H. Li, S. He, Fast algorithm based on TT-M FE system for space fractional Allen-Cahn equations with smooth and non-smooth solutions, J. Comput. Phys., 2019, 379: 351-372.                   |
Yang Liu Baoli Yin Hong Li Zhimin Zhang
The unified theory of shifted convolution quadrature for fractional calculus
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.