The unified theory of shifted convolution quadrature for fractional calculus

08/03/2019 ∙ by Yang Liu, et al. ∙ 0

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.



There are no comments yet.


page 13

page 14

page 15

page 16

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


and the Riemann-Liouville differential operator is defined by


where . When , we set , the identity operator.

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


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 [5]. Some other methods such as the shifted Grnwald 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 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 ,


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

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


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


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


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


Another key property of is the homogeneity:


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

Definition 2.2

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

Definition 2.3

A SCQ is convergent of order to if


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, [1]). The consistency of the is a special case of the when , i.e., the is consistent of order for if (see Definition 2.2, [1]).

Remark 2.4

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


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 [1].

Lemma 2.5

(The counterpart of Lemma 3.1, [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 ,


Let , the expression tends to . 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 for a consistent is of the form (see (3.6) in [1])


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, [1]) The is consistent of order if and only if the coefficients in (2.11) satisfy


where are the coefficients of


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


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


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


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 [1], 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 [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 , the following theorem shows that with the starting part , approximates uniformly for bounded .

Theorem 2.9

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

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


the SCQ satisfies


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


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, [1]) 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])


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


with the coefficients defined by .

The fractional BN- method (see [22])

Theorem 2.11

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


where the coefficients satisfy


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


Considering the consistency of the , i.e.,


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


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


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 ,


where the weights satisfy


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


see also (17), (18) in [4]. 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 ,


where are the coefficients of


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


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


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 [16] with replaced by ). Recently, Chen et al. [21] 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


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 shows some superconvergence when (see also [5], 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


is convergent of order provided that


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


whose coefficients can be formulated as


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


One can easily check that the with defined in (2.39) is the WSGL operator (see [6], 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


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


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


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


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. [8] and Ding et al. [11]. Define the generating function by


which can be reformulated as


The generating function (2.46) has been proposed by Gunarathna et al. [18] for the fractional derivative, i.e., for . For the case takes integers and other higher-order formulas, see [32]. 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 [8]) with the generating function as,


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


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


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


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




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


With the relation


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


respectively, where are defined by (2.30), and


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 is defined by,


For simplicity, we define .

Lemma 3.2

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


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 ,


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


where is a real constant. By replacing of with , 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 of a is the set of all complex for which the numerical solutions of (3.5) satisfy


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

Theorem 3.4

The stability regions of the numerical schemes (3.5) are



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