# Modified BDF2 schemes for subdiffusion models with a singular source term

The aim of this paper is to study the time stepping scheme for approximately solving the subdiffusion equation with a weakly singular source term. In this case, many popular time stepping schemes, including the correction of high-order BDF methods, may lose their high-order accuracy. To fill in this gap, in this paper, we develop a novel time stepping scheme, where the source term is regularized by using a k-fold integral-derivative and the equation is discretized by using a modified BDF2 convolution quadrature. We prove that the proposed time stepping scheme is second-order, even if the source term is nonsmooth in time and incompatible with the initial data. Numerical results are presented to support the theoretical results.

• 32 publications
• 5 publications
• 45 publications
03/07/2020

### High-order Time Stepping Schemes for Semilinear Subdiffusion Equations

The aim of this paper is to develop and analyze high-order time stepping...
11/25/2021

### Computing weakly singular and near-singular integrals in high-order boundary elements

We present algorithms for computing weakly singular and near-singular in...
08/28/2020

08/20/2022

### A Modified Trapezoidal Rule for a Class of Weakly Singular Integrals in n Dimensions

In this paper we propose and analyze a general arbitrarily high-order mo...
07/24/2022

### Integral Representations and Quadrature Schemes for the Modified Hilbert Transformation

We present quadrature schemes to calculate matrices, where the so-called...
12/08/2020

### Numerical analysis of a wave equation for lossy media obeying a frequency power law

We study a wave equation with a nonlocal time fractional damping term th...
12/27/2021

### Correction of high-order L_k approximation for subdiffusion

The subdiffusion equations with a Caputo fractional derivative of order ...

## 1 Introduction

For anomalous, non-Brownian diffusion, a mean squared displacement often follows the following power-law

 ⟨x2(t)⟩≃Kαtα.

Prominent examples for subdiffusion include the classical charge carrier transport in amorphous semiconductors, tracer diffusion in subsurface aquifers, porous systems, dynamics of a bead in a polymeric network, or the motion of passive tracers in living biological cells [18, 19]

. Subdiffusion of this type is characterised by a long-tailed waiting time probability density function

, corresponding to the time-fractional diffusion equation with and without an external force field [19, Eq. (88)]

 (♠) ∂tu(x,t)−∂1−αtAu(x,t)=f(x,t),  0<α<1.

Here is a given source function, and the operator denotes Laplacian on a polyhedral domain () with a homogenous Dirichlet boundary condition. The fractional derivative is taken in the Riemann-Liouville sense, that is, with the fractional integration operator

 Jαf(t)=1Γ(α)∫t0(t−τ)α−1f(τ)dτ=1Γ(α)tα−1∗f(t),

and denotes the Laplace convolution: .

Since the Riemann-Liouvile fractional derivative and the Caputo fractional derivative can be written in the form [22, p. 76]

 ∂αtu(x,t)=CDαtu(x,t)+1Γ(1−α)t−αu(x,0),

which implies that the equivalent form of can be rewritten as

 (♡) ∂tu(x,t)−CD1−αtAu(x,t)=f(x,t)+Au(x,0)Γ(α)t−(1−α),  0<α<1

with the Caputo fractional derivative

 CDαtu(t)=1Γ(1−α)∫t0(t−s)−αu′(s)ds,  0

Applying the fractional integration operator to both sides of , we obtain the equivalent form of as, see [17, Eq. (1.6)] or [26, Eq. (2.3)], namely,

 (♣) CDαtu(x,t)−Au(x,t)=1Γ(1−α)t−α∗f(x,t)−JαAu(x,t)|t=0Γ(1−α)t−α, 0<α<1.

As another example, the fractal mobile/immobile models for solute transport associate with power law decay PDF describing random waiting times in the immobile zone, leads to the following models [24, Eq. (15)]

 (♢) ∂tu(x,t)+CDαtu(x,t)−Au(x,t)=−1Γ(1−α)t−αu(x,0), 0<α<1.

Note that the right hand side in aforementioned PDE models ()-() might be nonsmooth in the time variable. In this paper, we consider the subdiffusion model with weakly singular source term:

 (1) CDαtu(x,t)−Au(x,t)=g(x,t):=tμ∘f(x,t)

with the initial condition , and the homogeneous Dirichlet boundary conditions. The symbol can be either the convolution or the product, and is a parameter such that

 μ>−1  if ∘ denotes convolution,  and  μ≥−α  if ∘ denotes product.

The well-posedness could be proved using the separation of variables and Mittag–Leffler functions, see e.g. [23, Eq. (2.11)].

Note that many existing time stepping schemes may lose their high-order accuracy when the source term is nonsmooth in the time variable. As an example, it was reported in [10, Section 4.1] that the convolution quadrature generated by step BDF method (with initial correction) converges with order , provided that the source term behaves like , , see Lemma 3.2 in [31], also see Table 1. The aim of this paper is to fill in this gap.

It is well-known that the smoothness of all the data of (1) (e.g., ) do not imply the smoothness of the solution which has an initial layer at (i.e., unbounded near ) [22, 23, 28]. There are already two predominant discretization techniques in time direction to restore the desired convergence rate for subdiffusion under appropriate regularity source function. The first type is that the nonuniform time meshes/graded meshes are employed to compensate/capture the singularity of the continuous solution near under the appropriate regularity source function and initial data, see [3, 11, 13, 16, 20, 21, 28]. See also spectral method with specially designed basis functions [4, 8, 33]. The second type is that, based on correction of high-order BDF or approximation, the desired high-order convergence rates can be restored even for nonsmooth initial data. For fractional ODEs, one idea is to use starting quadrature weights to correct the fractional integrals [14] (or fractional substantial calculus [1])

 Jαg(t)=1Γ(α)∫t0(t−τ)α−1g(τ)dτ  with  g(t)=tμf(t),  μ>−1,

where the algorithms rely on expanding the solution into power series of . For fractional PDEs, a common practice is to split the source term into

 g(t)=g(0)+k−1∑l=1tll!∂ltg(0)+tk−1(k−1)!∗∂ktg.

Then approximating by may to a modified BDF2 scheme with correction in the first step [5]. The correction of high-order BDF or convolution quadrature are well developed in [10, 27, 32] when the source term sufficiently smooth in the time variable. Performing the integral on both sides for (1), e.g, approximate by , a second-order time-stepping schemes are given in [34], where the singular source function is with a spatially dependent function . How to deal with a more general source term, which might be nonsmooth in the time variable, is still unavailable in the literature.

In this paper, we develop a novel second-order time stepping scheme (IDk-BDF2) for solving the subdiffusion (1) with a weakly singular source term, where the low regularly source term is regularized by using a -fold integral-derivative (IDk) and the equation is discretized by using a modified BDF2 convolution quadrature. We prove that the proposed time stepping scheme is second-order, even if the source term is nonsmooth in time and incompatible with the initial data. Numerical results are presented to support the theoretical results.

The paper is organized as follows. In Section 2, we introduce the development of the IDk-BDF2 scheme for model (1). In Section 3 and 4, based on operational calculus, the detailed convergence analysis of IDk-BDF2 are provided, respectively, for general source function and certain form . Then the desired results with the low regularity source term are obtained in Section 5. To show the effectiveness of the presented schemes, the results of numerical experiments are reported in Section 6. Finally, we conclude the paper with some remarks in the last section.

## 2 IDk-BDF2 Method

In this section, we first provide IDk-BDF2 method for solving subdiffusion (1) if the source term possess the mild regularity. Let with . Then the model (1) can be rewritten as

 (2) ∂αtV(t)−AV(t)=Av+g(t),0

From [15] and [29], we know that the operator satisfies the following resolvent estimate

 ∥∥(z−A)−1∥∥≤cϕ|z|−1∀z∈Σϕ

for all , where is a sector of the complex plane . Hence, with for all . Then, there exists a positive constant such that

 (3) ∥∥(zα−A)−1∥∥≤c|z|−α∀z∈Σθ.

### 2.1 Discretization schemes

Let and . By first fundamental theorem of calculus, we may rewrite (2) as

 (4) ID1 Method:  ∂αtV(t)−AV(t)=∂t(tAv+G(t)),0
 (5)

Let , be a uniform partition of the time interval with the step size , and let denote the approximation of and . The convolution quadrature generated by BDF2 approximates the Riemann-Liouville fractional derivative by

 (6) ∂ατφn:=1ταn∑j=0ωjφn−j

with . Here the weights are the coefficients in the series expansion

 (7) δατ(ξ)=1τα∞∑j=0ωjξjwithδτ(ξ):=1τ(32−2ξ+12ξ2).

Then IDk-BDF2 method for (4) and (5) are, respectively, designed by

 (8) ID1−BDF2 Method:    ∂ατVn−AVn=∂τ(tnAv+Gn).
 (9) ID2−BDF2 Method:    ∂ατVn−AVn=∂2τ(t2n2Av+Gn).
###### Remark 2.1

In the time semidiscrete approximation (8) and (9), we require , i.e., the initial data is reasonably smooth. However one can use the schemes (8) and (9) to prove the error estimates with the nonsmooth data , see Theorems 5.2 and 5.3. Here, we mainly focus on the time semidiscrete approximation (8) and (9), since the spatial discretization is well understood. For example, we choose if and if following [29, 30].

### 2.2 Solution representation for (4) and (5)

Taking the Laplace transform in both sides of (4), it leads to

 ˆV(z)=(zα−A)−1(z−1Av+zˆG(z)).

By the inverse Laplace transform, there exists [10]

 (10) V(t)=12πi∫Γθ,κezt(zα−A)−1(z−1Av+zˆG(z))dz

with

 (11) Γθ,κ={z∈C:|z|=κ,|argz|≤θ}∪{z∈C:z=re±iθ,r≥κ}

and , .

Similarly, applying the Laplace transform in both sides of (5), it yields

 ˆV(z)=(zα−A)−1(z−1Av+z2ˆG(z)).

By the inverse Laplace transform, we obtain

 (12) V(t)=12πi∫Γθ,κezt(zα−A)−1(z−1Av+z2ˆG(z))dz.

### 2.3 Discrete solution representation for (8) and (9)

Given a sequence and take to be its generating power series. Let be given in (7) and , . Then the discrete solution of (8) is represented by

 Vn=12πi∫Γτθ,κeztn(δατ(e−zτ)−A)−1δτ(e−zτ)τ(γ1(e−zτ)τAv+˜G(e−zτ))dz

with . Multiplying the (8) by and summing over with , we obtain

 ∞∑n=1∂ατVnξn−∞∑n=1AVnξn=∞∑n=1∂τ(tnAv+Gn)ξn.

From (6) and (7), we have

 ∞∑n=1∂ατVnξn=∞∑n=11ταn∑j=0ωjVn−jξn=∞∑n=01ταn∑j=0ωjVn−jξn=∞∑j=01τα∞∑n=jωjVn−jξn=∞∑j=01τα∞∑n=0ωjVnξn+j=1τα∞∑j=0ωjξj∞∑n=0Vnξn=δατ(ξ)˜V(ξ).

Similarly, one has

 ∞∑n=1∂τtnAvξn=δτ(ξ)γ1(ξ)τAv,∞∑n=1∂τGnξn=δτ(ξ)˜G(ξ)

 (13) ˜V(ξ)=(δατ(ξ)−A)−1δτ(ξ)(γ1(ξ)τAv+˜G(ξ)).

According to Cauchy’s integral formula, and the change of variables , and Cauchy’s theorem, one has [10]

 (14)

with . The proof is completed.

Let be given in (7) and , . Then the discrete solution of (9) is represented by

 Vn=τ2πi∫Γτθ,κeztn(δατ(e−zτ)−A)−1δ2τ(e−zτ)(γ2(e−zτ)2τ2Av+˜G(e−zτ))dz

with . Multiplying the (9) by and summing over with , we obtain

 ∞∑n=1∂ατVnξn−∞∑n=1AVnξn=∞∑n=1∂2τ(t2n2Av+Gn)ξn.

The similar arguments can be performed as Lemma 2.3, it yields

 ∞∑n=1∂ατVnξn=δατ(ξ)˜V(ξ),∞∑n=1∂2τt2nAvξn=δ2τ(ξ)γ2(ξ)τ2Av,∞∑n=1∂2τGnξn=δ2τ(ξ)˜G(ξ),  γ2(ξ)=ξ+ξ2(1−ξ)3,

and

 (15) ˜V(ξ)=(δατ(ξ)−A)−1δ2τ(ξ)(γ2(ξ)2τ2Av+˜G(ξ)).

Using Cauchy’s integral formula, and the change of variables , and Cauchy’s theorem, one has

 (16) Vn=τ2πi∫Γτθ,κeztn(δατ(e−zτ)−A)−1δ2τ(e−zτ)(γ2(ξ)2τ2Av+˜G(e−zτ))dz

with . The proof is completed.

## 3 Convergence analysis: General source function g(x,t)

In this section, we provide the detailed convergence analysis of ID1-BDF2 in (8) approximation for the subdiffusion (4), and ID2-BDF2 can be similarly augmented.

### 3.1 A few technical lemmas

First, we give some lemmas that will be used. [10] Let be given in (7). Then there exist the positive constants , and with such that

 c1|z|≤|δτ(e−zτ)|≤c2|z|,|δτ(e−zτ)−z|≤cτ2|z|3,|δατ(e−zτ)−zα|≤cτ2|z|2+α, δτ(e−zτ)∈Σπ/2+ε∀z∈Γτθ,κ.

Let be given in (7) and with . Then there exist a positive constants such that

 ∣∣∣γl(e−zτ)l!τl+1−z−l−1∣∣∣≤cτl+1,∀z∈Γτθ,κ,

where is sufficiently close to .

The arguments can be performed in [27] for . For , using

 11−e−zττ−z−1=z−(1−e−zτ)τ−1(1−e−zτ)τ−1z,

and Lemma 3.1, it yields and

 (17) ∣∣(1−e−zτ)τ−1z∣∣≥c|z|2  ∀z∈Γτθ,κ.

Since

Thus we have

 ∣∣∣11−e−zττ−z−1∣∣∣≤cτ.

The proof is completed.

Let be given in (7) and with . Then there exist a positive constants such that

 (18) ∣∣∣δτ(e−zτ)γl(e−zτ)l!τl+1−z−l∣∣∣≤cτl+1|z|+cτ2|z|2−l,∀z∈Γτθ,κ,

where is sufficiently close to . Let

 δτ(e−zτ)γl(e−zτ)l!τl+1−z−l=J1+J2

with

 J1=δτ(e−zτ)γl(e−zτ)l!τl+1−δτ(e−zτ)z−l−1andJ2=δτ(e−zτ)z−l−1−z−l.

According to Lemma 3.1 and 3.1, we have

 |J1|=∣∣ ∣∣δτ(e−zτ)(γl(e−zτ)l!τl+1−z−l−1)∣∣ ∣∣≤c2|z|cτl+1≤cτl+1|z|

and

 |J2|=∣∣(δτ(e−zτ)−z)z−l−1∣∣≤cτ2|z|2−l.

By the triangle inequality, the desired result is obtained.

Let be given by (7) and with . Then there exist a positive constants such that

 ∥∥∥(δατ(e−zτ)−A)−1δτ(e−zτ)γl(e−zτ)l!τl+1−(zα−A)−1z−l∥∥∥≤cτl+1|z|1−α+cτ2|z|2−l−α.

Let

 (δατ(e−zτ)−A)−1δτ(e−zτ)γl(e−zτ)l!τl+1−(zα−A)−1z−l=I+II

with

 I=(δατ(e−zτ)−A)−1[δτ(e−zτ)γl(e−zτ)l!τl+1−z−l],II=[(δατ(e−zτ)−A)−1−(zα−A)−1]z−l.

The resolvent estimate (3) and Lemma 3.1 imply directly

 (19) ∥(δατ(e−zτ)−A)−1∥≤c|z|−α.

From (19) and Lemma 3.1, we obtain

 ∥I∥≤cτl+1|z|1−α+cτ2|z|2−l−α.

Using Lemma 3.1, (19) and the identity

 (20)

we estimate as following

 ∥II∥≤cτ2|z|2+αc|z|−αc|z|−α|z|−l≤cτ2|z|2−l−α.

By the triangle inequality, the desired result is obtained.

Let be given by (7) and . Then there exist a positive constants such that

 ∥∥(δατ(e−zτ)−A)−1δτ(e−zτ)γ1(e−zτ)τ2A−(zα−A)−1z−1A∥∥≤cτ2|z|.

Using identical and

 (δατ(e−zτ)−A)−1δτ(e−zτ)A=−δτ(e−zτ)+(δατ(e−zτ)−A)−1δατ(e−zτ)δτ(e−zτ)A,

we get

 (δατ(e−zτ)−A)−1δτ(e−zτ)γ1(e−zτ)τ2A−(zα−A)−1z−1A=J1+J2+J3+J4

with

 J1=(δατ(e−zτ)−A)−1δατ(e−zτ)(δτ(e−zτ)γ1(e−zτ)τl+1−z−1),J2=(δατ(e−zτ)−A)−1(δατ(e−zτ)−zα)z−1,J3=((δατ(e−zτ)−A)−1−(zα−A)−1)zα−1,J4=z−1−δτ(e−zτ)γ1(e−zτ)τ2.

According to (19) and Lemmas 3.1, 3.1 with , we estimate , and as following

 ∥J1∥≤c|z|−α|z|ατ2|z|≤cτ2|z|,∥J2∥≤c|z|−ατ2|z|2+α|z|−1≤cτ2|z|,∥J4∥≤cτ2|z|.

From Lemma 3.1, (19) and the identity (20), we estimate as following

 ∥J3∥≤cτ2|z|2+α|z|−α|z|−α|z|α−1≤cτ2|z|.

By the triangle inequality, the desired result is obtained.

### 3.2 Error analysis for general source function g(x,t)

From , the Taylor expansion of source function with the remainder term in integral form:

 1∗g(t)=G(t)=G(0)+tG′(0)+t22G′′(0)+t22∗G′′′(t)=J1g(0)+tg(0)+t22g′(0)+t22∗g′′(t).

Then we obtain the following results with . Let and be the solutions of (4) and (8), respectively. Let and with . Then

 (21) ∥V(tn)−Vn∥≤(cτl+1tα−2n+cτ2tα+l−3n)∥∥g(l−1)(0)∥∥.

Using (10) and (14), there exist

 V(tn)=12πi∫Γθ,κeztn(zα−A)−11zlg(l−1)(0)dz,

and

 Vn=12πi∫Γτθ,κeztn(δατ(e−zτ)−A)−1δτ(e−zτ)γl(e−zτ)l!τl+1g(l−1)(0)dz,

where is sufficiently close to , and . Let

 V(tn)−Vn=J1+J2

with

 J1=12πi∫Γτθ,κeztn[(zα−A)−1zl−(δατ(e−zτ)−A)−1δτ(e−zτ)γ