## 1 Introduction

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

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)]() |

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

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]

which implies that the equivalent form of can be rewritten as

() |

with the Caputo fractional derivative

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,

() |

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

() |

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

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

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

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

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

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

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

(3) |

### 2.1 Discretization schemes

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

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

(7) |

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

(8) |

(9) |

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

By the inverse Laplace transform, there exists [10]

(10) |

with

(11) |

and , .

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

By the inverse Laplace transform, we obtain

(12) |

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

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

Similarly, one has

with . It leads to

(13) |

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

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

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

and

(15) |

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

(16) |

with . The proof is completed.

## 3 Convergence analysis: General source function

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

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

where is sufficiently close to .

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

and Lemma 3.1, it yields and

(17) |

Since

Thus we have

The proof is completed.

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

(18) |

where is sufficiently close to . Let

with

According to Lemma 3.1 and 3.1, we have

and

By the triangle inequality, the desired result is obtained.

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

Let

with

The resolvent estimate (3) and Lemma 3.1 imply directly

(19) |

From (19) and Lemma 3.1, we obtain

Using Lemma 3.1, (19) and the identity

(20) |

we estimate as following

By the triangle inequality, the desired result is obtained.

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

Using identical and

we get

with