# Composite Quadrature Methods for Weakly Singular Convolution Integrals

The well-known Caputo fractional derivative and the corresponding Caputo fractional integral occur naturally in many equations that model physical phenomena under inhomogeneous media. The relationship between the two fractional terms can be readily obtained by applying the Laplace transform to a given equation. We seek to numerically approximate Caputo fractional integrals using a Taylor series expansion for convolution integrals. This naturally extends into being able to approximate convolution integrals for a wider class of convolution integral kernels K(t-s). One of the main advantages under this approach is the ability to numerically approximate weakly singular kernels, which fail to converge under traditional quadrature methods. We provide stability and convergence analysis for these composite quadratures, which offer optimal convergence for approximating functions in C^γ[0,T], where α≤γ≤ 5 and 0<α < 1. For the order γ = 1,2,3,4,5 scheme, the resulting approximation is O(τ^γ) accurate, where τ is the size of the partition of the time domain. By instead utilizing a fractional Taylor series expansion, we are able to obtain for γ∈ (0,5)-{1,2,3,4} order scheme, which yields an approximation of O(τ^γ) with a constant dependent on the kernel function which improves the order of convergence. This allows for a far wider class of functions to be approximated, and by strengthening the regularity assumption, we are able to obtain more accurate results. General convolution integrals exhibit these results up to γ = 2 without the assumption of K being decreasing. Finally, some numerical examples are presented, which validate our findings.

Comments

There are no comments yet.

## Authors

• 1 publication
• 1 publication
07/02/2020

### Approximate solution of the integral equations involving kernel with additional singularity

The paper is devoted to the approximate solutions of the Fredholm integr...
12/25/2020

### Kernel-Independent Sum-of-Exponentials with Application to Convolution Quadrature

We propose an accurate algorithm for a novel sum-of-exponentials (SOE) a...
09/28/2019

### Complete monotonicity-preserving numerical methods for time fractional ODEs

The time fractional ODEs are equivalent to the convolutional Volterra in...
07/08/2021

### Fast accurate approximation of convolutions with weakly singular kernel and its applications

In this article, we present an O(N log N) rapidly convergent algorithm f...
04/10/2020

### Numerical methods for stochastic Volterra integral equations with weakly singular kernels

In this paper, we first establish the existence, uniqueness and Hölder c...
04/25/2020

### Estimates on translations and Taylor expansions in fractional Sobolev spaces

In this paper we study how the (normalised) Gagliardo semi-norms [u]_W^s...
##### 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

We begin by defining the Caputo fractional time-derivative [9,10] of a given function as

 C0Dαtf(t)=1Γ(1−α)∫t0df(s)ds(t−s)−αds,0<α<1

which is a fractional derivative of order . In [2], the use of the Laplace transform was applied to a Caputo fractional derivative term to obtain the fractional integral term

 C0Iαtf(t)=1Γ(α)∫t0(t−s)α−1f(s)ds, (1)

which was studied numerically and convergent schemes were developed for this integral inspired by the works presented in [9-13], [15-16]. The integral (1) can be expressed as the convolution , where . Our work will examine and derive numerical schemes to discretize integrals of the form (1), which has numerous engineering and physics applications, see [11], [12] and [16]. One of the major advantages of applying the Laplace Transform to a fractional derivative term, as seen in [2], is the ability to preserve the same assumption of regularity as required for fractional derivative discretizations and recovering an extra order of accuracy, where denotes the time step size. Further, we now have the ability to relax the regularity assumption from requiring the objective function under the usual L1-method to instead be any , where and . This is achieved by a usual Taylor series expansion to obtain convergence results for whole number values of , and by utilizing a fractional Taylor series expansion to approximate functions with a fractional order of regularity, see [15]. By requiring more regularity, we are able to obtain a higher order of convergence, as seen in Theorems 3.6 and 3.7. This method naturally generalizes to any convolution type-quadrature where the kernel function is positive, decreasing, and satisfies , as seen in Theorems 3.4 and 3.5.

The remainder of the paper is organized as follows. Section 2 will provide discretizations for fractional integrals of the above form, and a general scheme is established for fractional integrals of other forms based on the integral kernel. We obtain general schemes of orders up to order of accuracy. Section 3 establishes all of the necessary consistency, stability, and convergence results for each of these schemes, in addition to a discussion of the implementation of the schemes. We prove optimal order of convergence of our stable schemes, where the order is at least 3. The instability of schemes of order 6 and above are presented as well. The main convergence results are featured in Theorems 3.4 through 3.6. Section 4 presents two fractional integral equations as numerical examples that validate our findings. Future works will consider the application of these methods to fractional order diffusion processes based on the validity of the schemes and the order of accuracy they can provide.

## 2. Discretized numerical schemes

In order to discretize the Caputo fractional integral (1), we must consider the Taylor expansion for a given function at the arbitrary point , the usual time interval over which the integral is considered. That is, the Taylor expansion centered at the point for on , given , is

 f(s)=f(tk)+(s−tk)f′(tk)+(s−tk)22!f′′(tk)+(s−tk)33!f′′′(tk)+... (2)

Define and let and such that . From the above, similar Taylor expansions centered at any given can be constructed for each of the previous points . That is,

 f(tk) =f(tk) (3) f(tk−1) =f(tk)−τkf′(tk)+τ2k2!f′′(tk)−τ3k3!f′′′(tk)+O(τ4k) (4) f(tk−2) =f(tk)−2τkf′(tk)+(2τk)22!f′′(tk)−(2τk)33!f′′′(tk)+O(τ4k) (5) f(tk−3) =f(tk)−3τkf′(tk)+(3τk)22!f′′(tk)−(3τk)33!f′′′(tk)+O(τ4k) (6) ... f(t0) =f(tk)−kτkf′(tk)+(kτk)22!f′′(tk)−(kτk)33!f′′′(tk)+O(τ4k) (7)

Thus, the -th order approximation of for any at the point requires we solve the linear system of equations:

 j−1∑i=0ckif(tk−i) =f(s) (8) =j−1∑i=0(s−tk)ii!f(i)(tk)+O((s−tk)j). (9)

Each of the terms are then replaced by their Taylor expansions about the point and then solved, after neglecting the higher order terms of and . For example, a second order approximation of is provided in [2], Theorem 3.2, which can be recovered by solving the equation

 ck0f(tk)+ck1(f(tk)−τkf′(tk)) =f(tk)+(s−tk)f′(tk).

This equation can be rewritten as a system of equations

 ck0+ck1 =1 −ck1τk =(s−tk),

Which yields the solution , . Therefore, we may numerically approximate the integral as seen in [2] using and as solved for above:

 ∫tn0(tn−s)α−1Γ(α)f(s)ds =n∑k=1∫tktk−1(tn−s)α−1Γ(α)f(s)ds ≈n∑k=1∫tktk−1(tn−s)α−1Γ(α)(ck0f(tk)+ck1f(tk−1))ds,

which recovers the equation that was studied in greater detail in [2]. We remark that under this construction, we satisfy the condition . This directly implies that the coefficients and presented above are positive. We now provide the values of the coefficients for each scheme up to 4th order accuracy. Higher order schemes can be derived using the generalized system of equations (8) and (9). We remark that in general, for each i.

First order accurate:

 ck0 =1, f(s) =f(tk)+O(τk).

Second order accurate:

 ck0 =1−tk−sτk,ck1=tk−sτk, f(s) =(1−tk−sτk)f(tk)+(tk−sτk)f(tk−1)+O(τ2k).

Third order accurate:

 ck0 =(τk+s−tk)(2τk+s−tk)2τ2k,ck1=(tk−s)(2τk+s−tk)τ2k, ck2 =(s−tk)(τk+s−tk)2τ2k f(s) =((τk+s−tk)(2τk+s−tk)2τ2k)f(tk)+((tk−s)(2τk+s−tk)τ2k)f(tk−1) +((s−tk)(τk+s−tk)2τ2k)f(tk−2)+O(τ3k).

Fourth order accurate:

 ck0 =(τk+s−tk)(2τk+s−tk)(3τk+s−tk)6τ3k, ck1 =(tk−s)(2τk+s−tk)(3τk+s−tk)2τ3k, ck2 =(s−tk)(τk+s−tk)(3τk+s−tk)2τ3k, ck3 =(tk−s)(τk+s−tk)(2τk+s−tk)6τ3k f(s) =((τk+s−tk)(2τk+s−tk)(3τk+s−tk)6τ3k)f(tk) +((tk−s)(2τk+s−tk)(3τk+s−tk)2τ3k)f(tk−1) +((s−tk)(τk+s−tk)(3τk+s−tk)2τ3k)f(tk−2) +((tk−s)(τk+s−tk)(2τk+s−tk)6τ3k)f(tk−3)+O(τ4k).

We present the generalized equation to solve for a n-th order accurate approximation to . The above equations may be rewritten by the matrix equation

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣111...10−τk−2τk...−(n−1)τk0(−τk)2(−2τk)2...(−(n−1)τk)2...0(−τk)n−1(−2τk)n−1...(−(n−1)τk)n−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ck0ck1ck2...ckn−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1(s−tk)(s−tk)2...(s−tk)n−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (10)

The first matrix is, in fact, the transpose of the usual Vandermonde matrix [14] where each entry, other than the first row of ones, is a multiple of the time partition . That is, we set , , , … so that

 VTτk→ckn=→ykn, (11)

where is the Vandermonde matrix with each expressed in terms of as above, and

 →ckn=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ck0ck1ck2...ckn−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,→ykn=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣1(s−tk)(s−tk)2...(s−tk)n−1⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

The following lemma asserts the existence of any general n-th order numerical scheme. Since we have

 det(VTτk)=det(Vτk) =∏1≤i

provided that which directly implies that the matrix is invertible under this condition. The following lemma is immediate from the above.

###### Lemma 2.1.

Let . Then, (11) has a unique solution for each .

We now compute the unique solution based on the previous lemma. From [14], we can establish the generalized inverse of the Vandermonde matrix.

###### Theorem 2.2.

Let . Then, (11) has a unique solution for each , with the solution

 [cki]=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∑1≤j≤n(s−tk)j−1(−1)n−i⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝∑1≤m1,m2,...,mn−i,m1,...,mn−1≠jxm1...xmn−i∏1≤i
###### Proof.

From Lemma 2.1, we may invert the matrix to obtain the solution

 →ckn=(VTτk)−1→ykn,

where then from [14], each entry of , is calculated by

 vij=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(−1)n−i⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝∑1≤m1,m2,...,mn−im1,...,mn−1≠jxm1...xmn−i∏1≤i

so we may solve component-wise to find each entry of , , from

 (VTτk)−1→ykn=∑1≤j≤nvijykj,

where , . Thus,

 [cki]=(VTτk)−1→ykn=∑jvijykj (13) (14)

###### Remark 2.3.

By utilizing the fractional Taylor series expansion instead for on , as discussed in [15], we may obtain similar results to those outlined in Theorem 2.2. This can further relax the regularity assumption to .

Using the fractional Taylor series expansion, we can assert that we have an order scheme defined by the following:
order accurate:

 ck0 =1, f(s) =f(tk)+O(ταk).

We now examine the consistency, stability, and convergence of these schemes based on the generalized scheme

 f(s)=n−1∑i=0ckif(tk−i)+O((s−tk)n). (15)

## 3. Numerical Consistency, Stability, and Convergence

### 3.1. Numerical Consistency and Stability

We motivate our discussion of stability by examining the results presented in [1]. The quadrature studied in [1] is of the form

 ∫T0ϕ(s)ds=τN∑j=0wjϕ(jτ)+O(ϵ), (16)

which provides the error estimate given

a sequence of constants, , such that

 ∫T0ϕ(s)ds−τN∑j=0wjϕ(jτ)=R∑l=ρ+1hl(rlcl){ϕ(l−1)(T)−ϕ(l−1)(0)}+O(hR).

We will directly compare these results to the ones established in the previous section to prove stability and assert convergence. Our goal is to decompose the integrand into a convolution integral , where we may relax the continuity assumptions on the kernel function . This goal is motivated in part from the results obtained in [2], where we seek a generalization for the integral kernel. We begin by recalling some basic definitions for quadrature methods. From (1.15) of [1], a quadrature method is said to be consistent if it satisfies

 N∑j=0wj=N

for the global quadrature of the integral (16). We will relate (16) and a generalization of the results provided in [2].

###### Lemma 3.1.

Let for any prescribed . Let denote the order of the desired approximation to the function , let such that and . Then, for an order scheme, as described in Theorem 2.2, we have

 ∫tn0ϕ(s)ds=n∑k=1γ−1∑j=0wkjf(tk−j)+O(ϵ). (17)
###### Proof.

By utilizing the taylor expansion for about the point , we may readily obtain a similar quadrature rule by using Theorem 2.2 and the definition of each defined in (12). By further remarking that for each , then we may write . Thus,

 ∫tn0 ϕ(s)ds=∫tn0f(s)K(tn−s)ds (18) =n∑k=1∫tktk−1f(s)K(tn−s)ds (19) =n∑k=1∫tktk−1(γ−1∑j=0ckj(s)f(tk−j)+O(ϵ))K(tn−s)ds (20) =n∑k=1γ−1∑j=0f(tk−j)∫tktk−1ckj(s)K(tn−s)ds+O(ϵ). (21)

By letting

 wkj=∫tktk−1ckj(s)K(tn−s)ds, (22)

we arrive at the conclusion. ∎

The following remark is a natural extension of the first lemma, which allows for direct comparison to prove stability using the Theorem 3.7 in [1].

###### Remark 3.2.

By expanding the series

 n∑k=1γ−1∑j=0wkjf(tk−j)+O(ϵ),

and by collecting all of the repeating terms for each , we may condense the double summation into a single summation term

 n∑k=1γ−1∑j=0wkjf(tk−j)=n∑k=0(wk0+wk+11+...+wk+γ−1γ−1)f(tk), (23)

where we note that to satisfy the previous lemma. Further, by defining for fixed

 ~wγk=wk0+wk+11+...+wk+γ−1γ−1, (24)

We arrive at a form identical to the generalized quadrature rule posed in [1]

 n∑k=1γ−1∑j=0wkjf(tk−j)+O(ϵ)=n∑k=0~wγkf(tk)+O(ϵ). (25)
###### Theorem 3.3.

The approximation scheme (25) is consistent for any , where is the order of approximation.

###### Proof.

From the consistency requirement in [1], we must show that the scheme (25) satisfies any time step

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∫tn0ϕ(s)ds=τn∑j=0wn−jϕ(jτ)+O(ϵ)n∑j=0wj=n (26)

for any fixed . That is, we have from Remark 3.2

 n∑k=0~wγk =n∑k=0wk0+wk+11+...+wk+γ−1γ−1 (27) =n∑k=1γ−1∑j=0wkj (28) =n∑k=1γ−1∑j=0∫tktk−1ckj(s)K(tn−s)ds (29) =n∑k=1∫tktk−1(γ−1∑j=0ckj(s))K(tn−s)ds. (30)

From (11), the first equation in the Vandermonde matrix requires , hence,

 n∑k=0~wγk=n∑k=1∫tktk−1K(tn−s)ds =∫tn0K(tn−s)ds. (31)

On the other hand, by relabelling the coefficients of (26) and by noting that ,

 ∫tn0f(s)K(tn−s)ds =τn∑j=0wn−jf(jτ)K(tn−jτ)+O(ϵ) (32) =τn∑k=0wn−kf(tk)K(tn−tk)+O(ϵ). (33)

By equating (25) and (33), we have

 τn∑k=0wn−kK(tn−tk)=n∑k=0~wγk =∫tn0K(tn−s)ds. (34)

Since each is arbitrary under this construction, we select to satisfy . Thus, we have for the scheme (25)

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∫tn0f(s)K(tn−s)ds=τn∑k=0wn−kf(tk)K(tn−tk)+O(ϵ)n∑k=0wk=n, (35)

hence the scheme (25) is consistent. For simplicity and for implementation, we take for each k to trivially satisfy these conditions since . ∎

We must further satisfy stability requirements in order to prove the convergence of these schemes for any order . From [1], we have the following theorem asserting stability under arbitrary quadrature rules:

###### Theorem.

(3.7 of [1]) The stability polynomial

 Σ(μ;λτ)= (1−λτw0K(0))μN−λτw1K(τ)μN−1−... −λτwNK(nτ) (36)

is Schur, if . Assuming each and satisfy , the recurrence for

 y(t)=f(t)+λ∫tntn−TK(tn−s)y(s)ds

when for is stable whenever , given .

We remark that under these results, we must simply satisfy the requirement that each in (25) to satisfy a similar stability criterion for this generalized quadrature. This leads immediately to two stability results:

###### Theorem 3.4.

Let be a positive kernel on and let for all k. Then, the approximation scheme (25) is stable for and , where is the order of approximation.

###### Proof.

The case where is immediate since , hence . For , then

 ~wk2 =wk0+wk+11 (37) =∫tktk−1ck0(s)K(s)ds+∫tk+1tkck+11(s)K(s)ds (38) ≥mins∈[t1,T]|K(s)|(∫tktk−1ck0(s)ds+∫tk+1tkck+11(s)ds) (39) =mins∈[t1,T]|K(s)|(∫tk+1tkck+10(s)ds+∫tk+1tkck+11(s)ds) (40) =mins∈[t1,T]|K(s)|(∫tktk−1ck0(s)+ck+11(s)ds) (41) =mins∈[t1,T]|K(s)|τk+1≥0. (42)

Using similar analysis we are able to come to the same conclusion for and , given . Therefore, when , the scheme (25) is stable. ∎

We require additional assumptions on the integral kernel to ensure that the scheme is stable in the case where the order of approximation to (25) is any order .

###### Theorem 3.5.

Let be a positive, nonincreasing kernel on and let for all k. The approximation scheme (25) is stable for any order of accuracy.

###### Proof.

We begin by showing that for each . That is, we use the relationship established in Remark 3.2. We will present the argument for the cases where and deduce the pattern from there. We remark that under the construction found in Theorem 2.2 that for then , provided . Therefore, when , we have

 ~wk3 =wk0+wk+11+wk+22 =∫tktk−1ck0(s)K(s)ds+∫tk+1tkck+11(s)K(s)ds+∫tk+2tk+1ck+22(s)K(s)ds ≥K(tk+1)(∫tktk−1ck0(s)ds+∫tk+1tkck+11(s)ds+∫tk+2tk+1ck+22(s)ds)