    # 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^α 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.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 which is defined by

 Iαf(x)=1Γ(α)∫x0(x−s)α−1f(s)ds,for Rα>0, (1.1)

and the Riemann-Liouville differential operator is defined by

 I−αf(x)=dndxnIn−αf(x)=1Γ(n−α)dndxn∫x0(x−s)n−α−1f(s)ds, (1.2)

where . When , we set , the identity operator.

In 1986, Lubich  developed the convolution quadrature () theory to approximate the Riemann-Liouville calculus with arbitrary at the node ,

 (Iαf)(x)≈hαn∑j=0ωjf(x−jh)+hαs∑j=0wn,jf(jh), (1.3)

where is the mesh size of a uniform grid, denotes the convolution weight and is the starting weight. For brevity, define . For , the theory requires that . When takes integers, coincides with the traditional integral operator and differential operator , and the identity operator . From this aspect, we can conclude that the 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 and its conterpart for fractional derivative approximated at the node , 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 . Some other methods such as the shifted Grnwald formula , the WSGL operator , the method developed by Ding et al.  and recently proposed higher-order approximation formulas  that generalize the fractional BDFp for are excluded from the theory as well, since all above methods approximate at a shifted node.

In this paper, we generalize the theory to cover the methods mentioned above by introducing a shifted parameter to develop the shifted convolution quadrature () theory and further design some new methods with higher-order convergence rate. Specifically, we approximate the operator at ,

 (Iαf)(x)≈hαn∑j=0ωjf(x+(θ−j)h)+hαs∑j=0wn,jf(jh). (1.4)

The contributions of the paper are as follows

Develop systematic approaches to approximate the fractional calculus at node without assumptions on the regularity of . 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 ) based on known approximation methods (for ), 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 theory to the 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 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 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 and develop some definitions, lemmas and theorems for the . Some existing or newly proposed schemes are discussed by several examples. In Sec. 3, we analyse the stability regions for some 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 theory, with the purpose of the easy application of the discrete fractional Grnwall 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, ) which extends the classical theorem of Dahlquist  on linear multistep methods to fractional ones. For convenience in the subsequent analysis we introduce some notations and definitions. Define

 (Iαhf)(x)=hαn∑j=0ωjf(x+(θ−j)h)+hαs∑j=0wn,jf(jh),x=(n−θ)h, (2.1)

as the shifted convolution quadrature (). Denote as the convolution part and as the starting part of (2.1), respectively. Define the convolution quadrature error by

 Eαh=Ωαh−Iα. (2.2)

For the sequence of convolutions weights we associate with a generating power series , and viceversa.

We note that if is continuous and is locally integrable, and is extended for such that

 Ωαhf(x)=hα∑0≤jh≤x+θhωjf(x+(θ−j)h), (2.3)

then commutes with convolution . Hence, the convolution error satisfies (for , we require )

 Eαh(f∗g)=(Eαhf)∗g. (2.4)

Another key property of is the homogeneity:

 (Eαhtβ−1)(x)=xα+β−1(Eαh/xtβ−1)(1). (2.5)

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 ,

###### Definition 2.1

A SCQ is stable for if the convolution weights is

 ωn=O(nα−1). (2.6)
###### Definition 2.2

A SCQ is consistent of order for if the generating function of satisfies

 hαeθhω(e−h)=1+O(hp). (2.7)
###### Definition 2.3

A SCQ is convergent of order to if

 (Eαhtβ−1)(1)=O(hβ)+O(hp),for all β∈C,β≠0,−1,−2,⋯. (2.8)

We take note of the fact that the definition of stability and convergence of the are the same as the corresponding definitions of the (see Definition 2.1 and 2.3, ). The consistency of the is a special case of the when , i.e., the is consistent of order for if (see Definition 2.2, ).

###### Remark 2.4

Combining the homogeneity of and the definition of convergent (2.8), we can get

 (Eαhtβ−1)(xn)=O(xα−1nhβ)+O(xα+β−1−pnhp),for all β∈C,β≠0,−1,−2,⋯, (2.9)

which means for those small that , the convergence rate will be much lower than . 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 , generalize the arguments for the in .

###### Lemma 2.5

(The counterpart of Lemma 3.1, ) If for then the is consistent of order .

Proof. First we examine the convolution error for the function with respect to on the interval ,

 eh(x):=(Eαhet−x)(x)=hα∑0≤jh≤x+θhωje(θ−j)h−(Iαet−x)(x)=hαeθh∑0≤jh≤x+θhωje−jh−(Iαet−x)(x). (2.10)

Let , the expression tends to . The rest argument of the proof is exactly the same as that of the Lemma 3.1 in , which is omitted here. The proof of the lemma is completed.

As is pointed out in , the structure of the generating function for a consistent is of the form (see (3.6) in )

 ϖ(ξ)=(1−ξ)−α[c0+c1(1−ξ)+c2(1−ξ)2+⋯+cN−1(1−ξ)N−1+(1−ξ)N~r(ξ)]. (2.11)

where is holomorphic at , and constants where are defined by (2.21). We argue that the generating function for can be expressed similarly by (2.11) with a different definition of the coefficients that depend on :

###### Lemma 2.6

(The counterpart of Lemma 3.2, ) The is consistent of order if and only if the coefficients in (2.11) satisfy

 ci=γi,for i=0,1,⋯,p−1, (2.12)

where are the coefficients of

 ∞∑i=0γi(1−ξ)i=ξθ(−lnξ1−ξ)−α. (2.13)

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

###### Remark 2.7

We can construct by polynomial functions. Assume has the form

 ω(ξ)=[p1(ξ)p2(ξ)]αp3(ξ)p4(ξ), (2.14)

where are polynomial functions, and . Denote by the unit disc , and by the closed unit disc in the complex plane. Then, for a stable , which means , it holds that is analytic in . Hence, can be written as

 ω(ξ)=ν(ξ)ℓ∏j=0(ξ−ξj)−αj,ξ0=1,α0=α, (2.15)

where are distinct with , is analytic on , and , . It can be shown that is equivalent to (2.6), see . We limit the choice of the parameter by the Condition-,

 Condition-ω: ω(ξ) is analytic without zeros in D, and Rαj≤Rα. (2.16)

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 :

###### Theorem 2.8

A with convolution weights defined by a generating function satisfying the condition- is convergent of order if and only if it is stable and consistent of order .

Proof. The theorem is a generalization of the Theorem 2.4 in , whose proof consists of several lemmas. Since the definition of stability and convergence of the are the same as those of the , we merely generalize the lemmas in  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 , the following theorem shows that with the starting part , approximates uniformly for bounded .

###### Theorem 2.9

(See Theorem 2.4 in ) Suppose is the generating function of a convergent . Then we have: For any ,

(i) there exist starting weights such that for any function

 f(x)=xβ−1g(x) with g % sufficiently differentiable, (2.17)

the SCQ satisfies

 Iαhf(x)−Iαf(x)=O(hp) (2.18)

uniformly for with .

(ii) there exist bounded starting weights such that for any function (2.17), the estimate (2.18) holds uniformly for bounded .

###### Remark 2.10

The starting weights in (i) are derived by letting

 (Iαhtq+β+1)(tn−θ)=(Iαtq+β+1)(tn−θ) (2.19)

for all integer such that . Note that the estimate (2.18) requires that is bounded away from . The starting weights in (ii) require additionally those such that (2.18) holds for bounded .

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

(Theorem 2.6, ) For an implicit linear multistep method which is stable and consistent of order with the characteristic polynomials and , if the zeros of have absolute value less than , then . Some special cases include (see [1, 22])

 (2.20)

The generalized Newton-Gregory formula (see )

 ϖp(ξ)=(1−ξ)−α[γ′0+γ′1(1−ξ)+⋯+γ′p−1(1−ξ)p−1], (2.21)

with the coefficients defined by .

The fractional BN- method (see )

 ϖ2(ξ)=1−αθ+αθξ[(3/2−θ)−(2−2θ)ξ+(1/2−θ)ξ2]α,θ∈(−∞,1], α∈R, αθ≤12. (2.22)
###### Theorem 2.11

Suppose is the generating function of a with convergence order . Define

 ω(ξ)=ϖp(ξ)[γ′′0+γ′′1(1−ξ)+⋯+γ′′p−1(1−ξ)p−1], (2.23)

where the coefficients satisfy

 ∞∑i=0γ′′i(1−ξ)i=ξθ. (2.24)

Then the SCQ with convolution weights generated by is convergent of order .

Proof. For the generating function , we know that by the stability of the . Hence, we have . By replacing with in (2.24), we get

 eθhp−1∑i=0γ′′i(1−e−h)i=1−eθh∞∑i=pγ′′i(1−e−h)i=1+O(hp). (2.25)

Considering the consistency of the , i.e.,

 hαϖp(e−h)=1+O(hp), (2.26)

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

 hαeθhω(e−h)=1+O(hp). (2.27)

Then by Theorem 2.8 the proof is completed.

A special case for (2.23) that is of vital importance is when . For any defined by (2.20)-(2.22), taking , we obtain the generating function by (2.23) for the approximation of at node . The coefficients in (2.24) can be formulated as

 γ′′i=(−1)iΓ(θ+1)Γ(i+1)Γ(θ−i+1). (2.28)

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

###### Corollary 2.12

Let with . Then we have the following approximation formula for of order ,

 f(xn−θ)−p−1∑j=0θjf(xn−j)=O(xβ−1−pn−θhp)+O(x−1n−θhβ), (2.29)

where the weights satisfy

 p−1∑i=0θiξi=p−1∑i=0γ′′i(1−ξ)i, (2.30)

with coefficients defined in (2.28).

###### Remark 2.13

The approximation formula (2.29) holds for with sufficiently differentiable as well. For , we can obtain two popular formulas

 f(xn−θ)≈(1−θ)f(xn)+θf(xn−1),f(xn−θ)≈12(1−θ)(2−θ)f(xn)+θ(2−θ)f(xn−1)+12θ(θ−1)f(xn−2), (2.31)

see also (17), (18) in . Considering the condition-, we limit for (2.31) to satisfying and , respectively. As one can see from (2.29), to get a convergence order of for bounded if is not so regular, the staring part is needed according to the Theorem 2.9.

Another family of generating functions for the 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 ) of order ,

 ω(ξ)=(1−ξ)−α[γ0+γ1(1−ξ)+⋯+γp−1(1−ξ)p−1], (2.32)

where are the coefficients of

 ∞∑i=0γi(1−ξ)i=ξθ(−lnξ1−ξ)−α. (2.33)

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

Example 1. (Generalized shifted Grnwald formula)

Expanding (2.33), we can easily derive that

 γ0=1,γ1=−α2−θ,γ2=18(α+2θ)2−124(5α+12θ). (2.34)

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

 (Iαf)(x)≈hαn∑j=0ωjf(x+(θ−j)h)+hαs∑j=0wn,jf(jh),with ω(ξ)=(1−ξ)−α. (2.35)

Now by assuming that is sufficiently smooth (hence, the starting part can be omitted), and that takes nonnegative integers, the relation (2.35) reduces to the shifted Grnwald formula (see  with replaced by ). Recently, Chen et al.  has applied the Grnwald formula () 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 Newton-Gregory formula of order can be formulated as

 ω(ξ)=(1−ξ)−α[1−(α2+θ)(1−ξ)]of order 2, with θ≤1−α2,ω(ξ)=(1−ξ)−α{1−(α2+θ)(1−ξ)+[18(α+2θ)2−124(5α+12θ)](1−ξ)2}of order % 3. (2.36)

Both of the formula were also proposed by Dimitrov  by using the theory developed in . An obvious conclusion is that defined by (2.36) which is of order shows some superconvergence when (see also , with replaced by ), since generally a with the generating function approximates 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 , that from the discussion above, a with defined by

 ω(ξ)=(1−ξ)−α[1−(α2+θ)(1−ξ)] (2.37)

is convergent of order provided that

 θ=1−α2−12√1−α3,for α≤3. (2.38)

One can check that this choice of satisfies the condition-.

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

Assume two integers , and let , consider the following generating function

 ω(ξ)=(1−ξ)−α[−α−2q2(p−q)+2p+α2(p−q)ξp−q], (2.39)

whose coefficients can be formulated as

 ωk=−α−2q2(p−q)ϱk+2p+α2(p−q)ϱk−(p−q), (2.40)

where are the coefficients of and define for . Considering the condition-, we assume satisfies

 p+q+α≤0. (2.41)

One can easily check that the with defined in (2.39) is the WSGL operator (see , with replaced by ), which is convergent of order by Theorem 2.8. Actually, by (2.40), we have since (see the proof of Lemma 4.1 with replaced by ). And, by the Taylor expansion formulas, we can get

 hαephω(e−h)=hαeph(1−e−h)−α[−α−2q2(p−q)+2p+α2(p−q)eh(q−p)]=hα(1−e−h)−α[−α−2q2(p−q)ehp+2p+α2(p−q)ehq]=[1−h2+O(h2)]−α[1−α2h+O(h2)]=[1+α2h+O(h2)][1−α2h+O(h2)]=1+O(h2), (2.42)

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

 ω(ξ)=(1−ξ)−α(λ1+λ2ξm),λ1≥|λ2|, (2.43)

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

 λ1=1−λ2,λ2=α+2θ2m,m=1+α+2θ2−5α+12θ6(α+2θ), (2.44)

one can easily check that the with defined by (2.44) is stable and consistent of order . Furthermore, since and cannot be zero by (2.44), it seems that there is no superconvergence point for the to devise a formula of convergence order . For the application of (2.43), by considering the condition (which is derived by the condition-), we can take

 θ=m−α2−12√m2−α3,for α≤3m2. (2.45)

We remark here that (2.43) generalizes (2.37) from some aspects that, if we take , under the condition (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.  and Ding et al. . Define the generating function by

 (2.46)

which can be reformulated as

 ω(ξ)=(32+θα)−α(1−ξ)−α(1−α+2θ3α+2θξ)−α. (2.47)

The generating function (2.46) has been proposed by Gunarathna et al.  for the fractional derivative, i.e., for . For the case takes integers and other higher-order formulas, see . By careful derivation, we can check that the with defined in (2.46) is convergent of order . If we take , (2.46) is reduced to the BDF- method (see ) with the generating function as,

 ω(ξ)=3−2θ2−(2−2θ)ξ+1−2θ2ξ2. (2.48)

If we take , we get a that approximates at the node (see  with replaced by ),

 ω(ξ)=(3α−12α−2α−1αξ+α−12αξ2)−α. (2.49)

There are two points deserve discussion:
a) The generating function (2.46) also implies that some superconvergence properties at node by taking .
b) By taking , we can obtain a shorter or simpler generating function

 ω(ξ)=2α(1−ξ2)−α, (2.50)

and we call the corresponding the fractional central difference method, for the reason that if we take , then the generating function implies the classical central difference scheme

 f′(xn−1)≈f(xn)−f(xn−2)2. (2.51)

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

 CDα0,xy(x)=λy(x),y(0)=y0,α∈(0,1), (3.1)

and

 Iαy(x)=1λy(x),α∈(0,1), (3.2)

respectively. The operator denotes the Caputo fractional derivative operator of order , which is defined by

 CDα0,xy(x)=1Γ(1−α)∫x0y′(s)(x−s)αds,α∈(0,1). (3.3)

With the relation

 CDα0,xy(x)=I−α(y(x)−y0),α∈(0,1), (3.4)

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

Define , the numerical schemes for equation (3.1) and (3.2) are as follows

 Scheme-I: n∑k=0ωn−kyk=zp−1∑j=0θjyn−j+f(xn−θ),Scheme-II: zn∑k=0ωn−kyk=p−1∑j=0θjyn−j, (3.5)

respectively, where are defined by (2.30), and

 f(xn−θ):=hα(I−αy0)(xn−θ)=y0(n−θ)αΓ(1−α). (3.6)

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 ) The two-parameter Mittag-Leffler function is defined by,

 Eα,β(x):=∞∑j=0xjΓ(jα+β), for those x that make the series converge. (3.7)

For simplicity, we define .

###### Lemma 3.2

For , the solutions of (3.1) and (3.2) satisfy

 y(x)→0asx→∞. (3.8)

Proof. For the equation (3.2), which is a special case of (1.2) in , 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 ,

 y(x)=y0Eα(λxα). (3.9)

Note that for , we have the asymptotic property (see Theorem 1.6, p.35, ) for that with ,

 |Eα(x)|≤C1+|x|,μ≤|argx|≤π,|x|≥0, (3.10)

where is a real constant. By replacing of with , the proof of the lemma is completed.

###### Definition 3.3

The stability region of a is the set of all complex for which the numerical solutions of (3.5) satisfy

 yn→0asn→∞. (3.11)

Further, we call a numerical method -stable if contains the sector .

###### Theorem 3.4

The stability regions of the numerical schemes (3.5) are

 For Scheme-I: C∖{z:z=ω(ξ)/p−1∑j=0θjξj,|ξ|≤1},For Scheme-II: C∖{z:z=p−1∑j=0θjξj/ω(ξ),|ξ|≤1}, (3.12)

respectively.

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

We conclude this section by illustrating the stability regions of some s and make some notations. Generally, we require that there exists a positive such that the interval is contained in the stability region .

Table 1 illustrates the stability regions for the shift-generalized Newton-Gregory formula of order . For the Scheme-I with , the method is -stable. If satisfies , the method is conditional stable which means for arbitrary , the step size must be small enough () such that for fixed , the solution tends to as tends to infinity. If , then for any