Higher Order Time Stepping Methods for Subdiffusion Problems Based on Weighted and Shifted Grünwald–Letnikov Formulae with Nonsmooth Data

Two higher order time stepping methods for solving subdiffusion problems are studied in this paper. The Caputo time fractional derivatives are approximated by using the weighted and shifted Grünwald–Letnikov formulae introduced in Tian et al. (Math Comput 84:2703–2727, 2015). After correcting a few starting steps, the proposed time stepping methods have the optimal convergence orders O(k2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(k^2)$$\end{document} and O(k3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ O(k^3)$$\end{document}, respectively for any fixed time t for both smooth and nonsmooth data. The error estimates are proved by directly bounding the approximation errors of the kernel functions. Moreover, we also present briefly the applicabilities of our time stepping schemes to various other fractional evolution equations. Finally, some numerical examples are given to show that the numerical results are consistent with the proven theoretical results.

In Sects. 2 and 3, with some suitable approximation z k of z, we shall choose θ ∈ (π/2, θ 0 ) sufficiently close to π/2 such that z α k ∈ Σ θ 0 which guarantees that (z α k I + A) −1 exists. Recently, Meerschaert et al. [27] and Tian et al. [34] introduced the weighted and shifted Grünwald-Letnikov difference operators to approximate the Riemann-Liouville fractional derivative and applied such difference operators to solve space fractional partial differential equations under the assumptions that the solution is sufficiently smooth and satisfies the homogeneous boundary conditions. To our knowledge, we have not seen any works in literature to apply such weighted and shifted Grünwald-Letnikov difference operator to construct higher order time discretization schemes for solving the subdiffusion equation (1). One of the possible reasons for lacking such works may be that the solution of (1) is not sufficiently smooth and it has the singularity near t = 0. For example, in the homogeneous case of (1) with f = 0, one has the following stability estimate, [29], with · the norm in L 2 (Ω), which shows that the α-th order Caputo derivative of the solution of (1) becomes unbounded as t → 0. Hence, the C 2 -regularity assumption, generally, does not hold for the exact solution of (1). Numerical experiments indicate that the convergence orders of some numerical methods for solving (1) actually do not hold uniformly in t even for the smooth data u 0 , see, e.g., Jin et al. [15], Stynes et al. [31] and Stynes [30]. Therefore, an attempt has been made in this paper to consider the higher order time discretization schemes for (1), based on the weighted and shifted Grünwald-Letnikov schemes developed in Tian et al. [34] and prove the optimal convergence orders of the proposed schemes with both smooth and nonsmooth data. There are several approaches to improve the convergence orders of the numerical methods for solving (1), when the solution is not sufficiently smooth. One approach is to correct some weights of the numerical methods in order to capture the singularity of the solution. This idea was first introduced by Lubich et al. [23] for second order time stepping scheme applied to an evolution equation with positive memory term. After correcting some weights of the numerical methods, Lubich et al. [23] proved optimal convergence of the corrected numerical method for both smooth and nonsmooth data. Jin et al. [16] derived some higher order numerical methods in time for the problem (1), where the fractional derivatives are approximated by using the convolution quadrature generated by using the backward difference formulae. By correcting some starting steps of the numerical methods, Jin et al. [16] established that the corrected numerical methods have optimal order of convergence for any fixed time t for both smooth and nonsmooth data, see also [14,17]. Subsequently, Yan et al. [36] corrected the starting steps of the L1 scheme for solving (1) and proved that the modified L1 scheme has the optimal convergence order O(k 2−α ), 0 < α < 1. More recently, Xing and Yan [35] analyzed a numerical method for solving (1), where the Caputo fractional derivative is expressed by using the Hadamard finite-part integral which is then approximated by using the quadratic interpolation polynomials. After correcting some starting steps and some weights of the high-order numerical methods, Xing and Yan [35] derived the optimal convergence order O(k 3−α ), 0 < α < 1 of the corrected numerical methods for both smooth and nonsmooth data. For the recent development of the corrections of numerical methods for (1), we refer the readers to the survey paper [12], see also [37]. For other numerical methods for solving time fractional diffusion equation, we refer to [3,[5][6][7][8]11,15,[18][19][20][21][22]24,25,28,32,[38][39][40], etc.
The aim of this paper is to prove that the proposed numerical methods have the optimal convergence orders O(k 2 ) and O(k 3 ), respectively, by correcting a few starting steps of the numerical methods for both smooth and nonsmooth datas. Compared to other higher order time stepping methods in the literature for solving time diffusion problem (1), the proposed methods have two advantages: (i) The weights of our numerical methods are much simpler than those obtained by approximating the fractional derivative with the quadratic interpolating polynomials, see, e.g., in Xing and Yan [35] and further, these weights have a special structure as mentioned in Tian et al. [34], which may be useful for constructing some fast algorithms and also for proving their stability and error analyses; (ii) The weights of the proposed numerical schemes are related not only to the order of the fractional derivative, but also to the shifted numbers, which imply that our methods are more related to the equation itself, see, e.g., Tian et al. [34].
The main contributions of this paper are as follows.
1. Based on the weighted and shifted Grünwald -Letnikov schemes proposed in Tian et al. [34], two new corrected higher order time discretization methods are introduced and the convergence orders are shown to be of O(k 2 ) and O(k 3 ), respectively for both smooth and nonsmooth data. 2. The error estimates of the corrected numerical methods are proved in both homogeneous and inhomogeneous cases. 3. With the help of Laplace transform techniques, it is shown that the error estimates are even suitable for more general elliptic operator A, which satisfies the resolvent estimate (2).
The paper is organized as follows. In Sect. 2, we consider the error estimates of the time discretization scheme for (1) with the convergence order O(k 2 ) for both smooth and nonsmooth data. In Sect. 3, we derive the error estimates of the time discretization scheme with the convergence order O(k 3 ) again for both smooth and nonsmooth data. Finally in Sect. 4, numerical examples are presented to show that numerical results are consistent with the theoretical results.
By C, we denote a positive constant independent of discretization parameter k, but not necessarily the same at different occurrences.

Second Order Time Stepping Scheme
In this section, we analyze a second order time discretization scheme for approximating the solution of the problem (1). After correcting some starting steps of the scheme, the optimal order of convergence is derived for the problem with both smooth and nonsmooth data.

Lemma 1 Letw(ζ ) be defined by
Proof From (10), we observe that where we have used the following binomial expansion, with β ∈ R, This completes the proof of the Lemma 1.
We now introduce a fully discrete scheme for solving (1). Let T h denote a triangulation of Ω with h the maximal length of the sides on T h . Let S h ⊂ H 1 0 (Ω) denote the piecewise continuous linear finite element space.
For any fixed t ∈ (0, T ], the finite element method of (1) is to find where A h : S h → S h denotes the discrete analogue of A defined with some suitable bilinear form a(·, ·) defined on H 1 0 (Ω) × H 1 0 (Ω) associated with the operator A, by and f h = P h f , where P h : L 2 (Ω) → S h denotes the L 2 project operator given by Here u 0h ∈ S h denotes some approximation of u 0 ∈ L 2 (Ω). When u 0 is nonsmooth, we choose u 0h = P h u 0 and when u 0 is smooth, that is u 0 ∈ D(A), we may choose u 0h = R h u 0 , where R h : H 1 0 (Ω) → S h denotes the Ritz projection or elliptic projection defined by Then, the equation (14) is equivalent to An application of Taylor's expansion as in Jin et al. [17] yields where g * h denotes the convolution of g and h.
By the Laplace transform method, we obtain, withĝ(z) denoting the Laplace transform of g(t), where Γ is defined by, see, e.g., Lubich et al. [23], with some θ ∈ (π/2, θ 0 ), Now we shall consider the time discretization scheme of (1). To improve the accuracy near t = 0, we follow the approach in Lubich et al. [23] to correct the values of the first step in the time discretization.
Let V n ≈ V h (t n ) be the approximation of V h (t n ). We define the following time discretization scheme for approximating (1) where, with c 0 = 1/2, Applying the discrete Laplace transform in both sides of (19), we obtain It is easy to see withṼ which implies by (20) that Further, withw(ζ ) given in (11), we set and By the inverse discrete Laplace transform, it follows for n ≥ 1 and using the variable change ζ = e −zk that where, see Lubich et al. [23], with Γ defined in (18), Below, we state our main result in this section whose proof will be provided subsequently.

Lemma 2
Let z k and μ(ζ ) with ζ = e −zk be defined by (21) and (22), respectively. Assume that Then with Γ k defined by (24), the following estimates hold: Proof We first show (26). Now, zk is uniformly bounded for any z ∈ Γ k since, with θ ∈ (π/2, θ 0 ), Further we note that, by (22) and Lemma 1 with c 0 = 1/2, and this implies that Hence, there exists δ 0 > 0 with 0 < δ 0 ≤ π sin θ such that (26) For large zk with δ 0 ≤ |zk| ≤ π sin θ , we now note by (22) that This is continuous at any ζ = 1, which implies that μ(e −zk ) − 1 is continuous at any z = 0. Since every continuous function is bounded on the closed and bounded domain, therefore, Hence, (26) also holds for δ 0 ≤ |zk| ≤ π sin θ , z ∈ Γ k which completes the estimate (26). In order to prove the estimate (27), it suffices to show | z z k | is bounded for any z ∈ Γ k . Now, a use of (21) yields To show | z z k | is bounded for any z ∈ Γ k , we consider two cases: one for the small zk and the other for the large zk.
For the small zk, observe that For the large zk, we note that |z| |z k | is continuous at any w = zk = 0, z ∈ Γ k which implies the boundedness of |z| |z k | for large |zk| with δ 0 ≤ |zk| ≤ π sin θ , z ∈ Γ k . Thus, we complete the estimate (27). Similarly, it is easy to show that | z k z | is also bounded for any z ∈ Γ k . For (28), we first observe with ζ = e −zk that This implies that there exists 0 < δ 0 ≤ π sin θ such that For large |zk| with δ 0 ≤ |zk| ≤ π sin θ , z ∈ Γ k , an application of (27) shows Hence, Following the idea of the proof in Lubich et al. [23, (4.6)] and noting that K 1 (z) ≤ C|z| −2 in [23, (3.12)], we obtain by mean value theorem and using (31), As in the proof of [36,Lemma 3.12], and noting that we now arrive by (32) and (26) at which completes the proof of (28). Finally in order to estimate (29), a use of the mean value theorem, (33) and (26), which shows (29). Altogether, it concludes the proof of the Lemma 2.

Lemma 4 Let z k be defined as in
Proof Apply Lemmas 2, 3 to arrive at This completes the rest of the proof. Now we turn to the proof of Theorem 1.

An application of Lemma 4 yields
For I 22 , following the arguments as in Jin et al. [14,17], we arrive at Together these estimates complete the proof of Theorem 1.

Third Order Time Discretization Method
In this section, we introduce a third order time discretization scheme for solving (1) based on the weighted and shifted Grünwald-Letnikov difference operator introduced in Tian et al. [34]. Let us define the following weighted and shifted Grünwald-Letnikov difference operator L D α k, p,q,r to approximate the Riemann-Liouville fractional derivative operator R where p, q, r are integers and mutually nonequal, and B α k, p φ(t n ) are defined by (6) and When p = 0, q = −1, r = −2, we obtain and we then arrive for n ≥ 2 at Thus, we obtain where The discrete Laplace transform of w (α) j , j = 0, 1, 2, . . . is given bỹ Below, following the idea in Lemma 1, we consider the series expansion of the functioñ w(ζ ) in (41).
Next, we turn to the solution of (16) when f h (t) is written as where g * h denotes the convolution of g and h.
An application of the Laplace transform to (16) with respect to the time variable t yieldŝ By the inverse Laplace transform, the solution of (16) takes the following form at t = t n , Let V n ≈ V h (t n ), n = 0, 1, 2, . . . , N denote the approximate solution of the following time discretization scheme for solving (16), with V 0 = 0, where w (α) j , j = 0, 1, 2, . . . are defined by (40) and the coefficients a 1 , a 2 , b 1 , b 2 satisfy and Now we come to the following main theorem in this section Theorem 2 Let V h (t n ) and V n be the solutions of (16) and (47)-(49), respectively. Assume that u 0 ∈ L 2 (Ω) and f ∈ C 2 ([0, T ]; L 2 (Ω)) and t n 0 (t n − s) α−1 f (s) ds < ∞ for t n ∈ (0, T ]. Let u 0h = P h u 0 , then there exists a positive constant C independent of k such that To prove Theorem 2, we need the following lemmas. (41), let z k and μ(ζ ) with ζ = e −zk be defined by

Lemma 6 Withw(ζ ) given by
and Further, let K 1 (z) and K 2 (z) be given by (25). Then, with Γ k as in (24), the following estimates hold: Proof The proof is similar to the proof of Lemma 2. For each inequality, again two cases such as the small zk and the large zk are considered. For (54), a use of Lemma 5 yields, with a 1 = 11/12, a 2 = −5/12 by (50), and then this implies that 3 as zk → 0.
Together with these estimates, we complete the proof of Lemma 6.

Lemma 8 Let z k be defined as in
Proof Let ζ = e −zk and x = zk, then by Lemma 5, we have for some suitable positive constants d 1 , d 2 .
Combining this with (51), we obtain This completes the proof of Lemma 8.

Lemma 9
Let z k be defined as in (52). Let b 1 , b 2 be defined as in (51). Then, there holds Proof The proof is similar to the proof of Jin et al. [16,Lemma C.1.] and hence, we omit the proof here. (52), then, there exists a positive constant C independent of k such that

Lemma 10 Let z k be defined as in
Proof From Lemma 6, it follows that This concludes the proof.
We are now ready to prove the main theorem of this section.
Proof (Proof of Theorem 2) We now calculate the approximate solution V n defined in (47)-(49). Taking the discrete Laplace transform in (47)-(49), we arrive at A use of the inverse discrete Laplace transform yields, with μ(ζ ) defined by (53), and Γ k as in (24).
Now, subtracting (46) from (63), we arrive at where, with K 2 (z) defined by (25), and and For I 1 , apply the bound (z α + A h ) −1 A h ≤ C, (3), (56) and (57) to obtain For I 21 , by (3), it follows that For I 22 , a use of the Lemma 9 shows Thus, we obtain For I 3 , we observe that and For I 31 , we easily bound it as An application of Lemma 10 yields For I 32 , following the arguments as in Jin et al. [14,17], we now arrive at Together these estimates complete the proof of Theorem 2.
where M θ 0 = 1/ sin(π − θ 0 ). With corresponding A h as semidiscrete approximation of A ( as in Sect. 2), the semidiscrete method becomes: findū h (t) ∈ V h such that Following the arguments of [9,23], one easily deduces property (70), when A is replaced by A h . Therefore, our two time stepping methods can be appropriately applied and the desired estimates can be easily derived.
Other generalizations also include the fractional cable equations with mass lumping by Al-Maskari and Karaa [2] and references therein. Further, it may include Fokker-Plank spatial discretization as our analysis does not depend on selfadjointness of the operator A and for complete discrete scheme, we may use the proposed time stepping schemes.
Finally, the present schemes can be applied to the Eq. (1), when A is a fractional Laplacian like A = (−Δ) s , s ∈ (0, 1), that is, with 0 < α, β < 1 and some positive constant δ > 0, An appropriate modification of the arguments in [1], we obtain the semidiscrete error estimate as: with > 0 small. Therefore, when it is combined with our proposed time stepping schemes, the final error estimate reads as: with U n ≈ u h (t n ), u(t n ) − U n ≤ C t −α n h s+min{s,1/2− } + t −m n k m u 0 , m = 2, 3.

Numerical Simulations
In this section, we present five numerical examples to show that the numerical results are consistent with the theoretical results obtained in this paper. The first three examples are solved by using the numerical method (19) for both homogeneous and inhomogeneous problems in one-and two-dimensional cases. The last two examples are computed by using the numerical method (47)-(49) for both homogeneous and inhomogeneous problems in one-dimensional case.
where (a) u 0 (x) = x(1 − x) (smooth data) and (b) u 0 (x) = χ [0,1/2] (nonsmooth data).  (19) in Example 1 at T = 1  Table 3 Time convergence orders for the uncorrected scheme (19) with c 0 = 0 in Example 1 at T = 1 Let N h be a positive integer. Let 0 = x 0 < x 1 < x 2 < · · · < x N h = 1 be the space partition and h the space step size. We shall use the piecewise linear finite element method to consider the space discretization.
Let 0 < t 0 < t 1 < · · · < t N = T be the time partition and k the time step size. To observe the convergence order of the numerical method, we first need to calculate the reference solution u re f (t) at some fixed time T with very small step sizes h re f = 2 −6 and k re f = 2 −10 .
By Theorem 1, we see that the numerical method (19) has the second order convergence. To see this convergence order, we shall calculate the approximate solution of u(T ) at T = 1 with the space step size h = 2 −6 and the different time step sizes k = κ * k re f with κ = [2 2 , 2 3 , 2 4 , 2 5 , 2 6 ]. In Table 2, we observe the convergence orders O(k 2 ) of the corrected scheme (19), where the rows with (a) denote the errors and the experimentally determined convergence orders in the smooth data case and the rows with (b) denote the errors and the orders in the nonsmooth data case. For each α, we choose the average convergence order of the computed orders obtained by using the different time step sizes.
For the numerical method (19) with c 0 = 0, that is, for the uncorrected scheme, we observe that, in Table 3, the experimentally determined convergence order is only O(k) with both smooth and nonsmooth data.
The second example is an inhomogeneous problem with zero initial value and the source term f which is smooth in time.  (19) in Example 2 at T = 1  Table 5 Time convergence orders for the uncorrected scheme (19) with c 0 = 0 in Example 2 at T = 1 where f (x, t) = (cos(t) + sin(t))(1 + χ (0,1/2) (x)).
We use the same parameters as in the numerical simulations in Example 1. In Table 4, we observe that the experimentally determined convergence order of the corrected scheme (19) indeed is O(k 2 ) for all 0 < α < 1 ( Table 5).
The third example is a two-dimensional example and we shall consider an inhomogeneous problem with nonsmooth initial data and the source term f which is smooth in time.
Let N h be a positive integer. Let 0 = x 0 < x 1 < x 2 < · · · < x N h = 1 and 0 = y 0 < y 1 < y 2 < · · · < y N h = 1 be the partition of Ω. We divide Ω into some triangles with the same sizes and let h be the maximal length of the sides of the triangle. We shall use the piecewise linear finite element method to consider the space discretization on the triangulation of Ω.
Let 0 < t 0 < t 1 < · · · < t N = T be the time partition and k the time step size. We shall use the very small space step size h re f = 2 −6 and the time step size k re f = 2 −10 to calculate the reference solution at time T . Table 6 Time convergence orders for the corrected scheme (19) in Example 3 at T = 1  Table 7 Time convergence orders for the uncorrected scheme (19) with c 0 = 0 in Example 3 at T = 1 We shall choose T = 1 in our simulation. We calculate the approximate solutions with the space step size h = 2 −6 and the time step sizes k = κ * k re f with κ = [2 2 , 2 3 , 2 4 , 2 5 , 2 6 ]. In Table 6, the experimentally determined convergence orders O(k 2 ) are observed as we expected.
In Table 7, we observe the experimentally determined convergence orders O(k) for the uncorrected scheme (19) with c 0 = 0 as we expected.
In the next two examples, we shall consider the experimentally determined convergence orders of the numerical method (47)-(49).

Example 4
In this example, we shall use the numerical method (47)-(49) to solve Example 1. We use the same parameters as in the numerical simulation in Example 1. In Table 8, we observe that the corrected scheme (47)-(49) has the convergence orders O(k 3 ) with both smooth and nonsmooth data as expected.

Example 5
In this example, we shall use the numerical method (47)-(49) to solve Example 2. We use the same parameters as in the numerical simulation in Example 1. In Table 10, we also observe the convergence orders O(k 3 ) of the corrected scheme (47)-(49) in the inhomogeneous case.
In Table 11, the experimentally determined convergence order of the uncorrected scheme (47)-(49) with a 1 = a 2 = b 1 = b 2 = 0 has only convergence order O(k). Tables 8 and 10, we observe that the experimentally determined convergence orders are slightly better than the theoretical orders. The possible reason may be that the proposed numerical methods involve both fractional orders and the shifted numbers. These combinations which are more related to the equation may be instrumental in helping us to provide possibly more accurate computational results. But, we do not have a theory yet to establish it. Therefore, in future, we shall continue to investigate this interesting observation.