Error Estimates of a Continuous Galerkin Time Stepping Method for Subdiffusion Problem

A continuous Galerkin time stepping method is introduced and analyzed for subdiffusion problem in an abstract setting. The approximate solution will be sought as a continuous piecewise linear function in time t and the test space is based on the discontinuous piecewise constant functions. We prove that the proposed time stepping method has the convergence order O(τ1+α),α∈(0,1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\tau ^{1+ \alpha }), \, \alpha \in (0, 1)$$\end{document} for general sectorial elliptic operators for nonsmooth data by using the Laplace transform method, where τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} is the time step size. This convergence order is higher than the convergence orders of the popular convolution quadrature methods (e.g., Lubich’s convolution methods) and L-type methods (e.g., L1 method), which have only O(τ)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\tau )$$\end{document} convergence for the nonsmooth data. Numerical examples are given to verify the robustness of the time discretization schemes with respect to data regularity.

Note that z α ∈ Σ θ with θ = αθ < π for all z ∈ Σ θ which implies that, by (3), with some π/2 < θ < π, Let S h ⊂ H 1 0 (Ω) denote the piecewise linear finite element space defined on the triangulation of Ω. The finite element method of (1), (2) is to find u h (t) ∈ S h such that for where A h : S h → S h denotes the discrete analogue of the elliptic operator A defined by, see Thomée [37,Chapter 12], where A(·, ·) denotes the bilinear form associated with A. Here P h : L 2 (Ω) → S h denotes the L 2 projection operator defined by Note that the minimal eigenvalue of A h is bounded below by that of A since S h is a subspace of L 2 (Ω), we have, with some π/2 < θ < π, see Lubich et al. [25, page 6], Let v h (t) = u h (t) − u 0h . Then (5), (6) can be written equivalently as In this work, we shall construct and analyze a continuous Galerkin time stepping method for solving (8), (9). It is more convenient for analyzing the time stepping method for (8), (9) than for (5), (6), see Lubich et al. [25]. Many application problems are modeled by using (1), (2), see, e.g., [1,9,10,32], etc. It is not possible to find the analytic solution of (1), (2). Therefore we need to design and analyze the numerical methods for solving (1), (2). Note that the Caputo fractional derivative is a nonlocal operator which makes the numerical analysis of the subdiffusion equation (1), (2) more difficult than the diffusion equation with the integer order derivative. There are two popular ways to approximate the Caputo fractional derivative in literature. One way is to use the convolution quadrature formula to approximate the Caputo fractional derivative, see, e.g., [25,26,46]. Another way is to use the L-type schemes to approximate the Caputo fractional derivative, see, e. g., [16,[20][21][22][23]34,45], etc. Under the assumptions that the solutions of (1), (2) are sufficiently smooth, the numerical methods constructed based on both convolution quadrature and L-type schemes have the optimal convergence orders, see [23,34], etc. However, Stynes et al. [36] and Stynes [35] showed that in general the solutions of (1), (2) have the limited regularities and the solutions behave as O(t α ) even for smooth initial data, which implies that the solutions of (1), (2) are not in C 1 [0, T ], see also [33]. Hence the higher order numerical methods constructed by convolution quadrature and L-type schemes for solving (1), (2) have only first order O(τ ) convergence for both smooth and nonsmooth data, see, e.g., [11,41], etc. To obtain the optimal convergence orders of the higher order numerical methods for solving (1), (2), one may use the corrected schemes to correct the weights of the starting steps of the numerical methods, see, e.g., [12,42], etc., or use the graded meshes to capture the singularities of the solutions of the subdiffusion problems, see, e.g., [16,36,45], etc. Discontinuous Galerkin methods are also well studied for solving fractional subdiffusion equations, see, e.g., [29][30][31] and the references therein. There are other numerical methods for solving fractional partial differential equations, see, e.g., [2][3][4][5][6]8,[13][14][15]24,27,28,39,43,44,46], etc.
Recently, Li et al. [17] analyzed the L1 scheme for solving the superdiffusion problem with the fractional order α ∈ (1, 2) based on the Petrov-Galerkin method. This scheme is first proposed and analyzed by Sun and Wu [34] in the framework of finite difference method under the assumptions that the solution of the problem is sufficiently smooth. Li et al. [17] proved that, without any regularity assumptions for the solution of the problem and without using the corrections of the weights and the graded meshes, the L1 scheme has the convergence order O(τ 3−α ), α ∈ (1, 2) for both smooth and nonsmooth data when the elliptic operator A is assumed to be self-adjoint, positive semidefinite and densely defined in a suitable Hilbert space, see also [18,19], etc.
The purpose of this paper is to consider the continuous Galerkin method for solving the subdiffusion problem with α ∈ (0, 1) by using the similar argument as in Li et al. [17] for the superdiffusion problem with α ∈ (1, 2). We prove that, without any regularity assumptions for the solution of the problem and without using the corrections of the weights and the graded meshes, the proposed time stepping method has the convergence order O(τ 1+α ), α ∈ (0, 1) for general sectorial elliptic operators satisfying the resolvent estimate (3).
The main contributions of this paper are the following: 1. A continuous Galerkin time stepping method for solving subdiffusion problem is introduced for general sectorial elliptic operators satisfying the resolvent estimate (3). 2. The convergence order O(τ 1+α ), α ∈ (0, 1) of the proposed time stepping method for solving homogeneous subdiffusion problem is proved by using the Laplace transform method. 3. The convergence order O(τ 1+α ), α ∈ (0, 1) of the proposed time stepping method for solving inhomogeneous subdiffusion problem is also proved by using the Laplace transform method.
The paper is organized as follows. In Sect. 2, we introduce the continuous Galerkin time stepping method for solving subdiffusion problem. In Sect. 3, the error estimates of the proposed time stepping methods are proved by using discrete Laplace trnsform method for the subdiffusion problems. In Sect. 4, some numerical examples for both fractional ordinary differential equations and subdiffusion problems are given to verify the theoretical results. An Appendix with three lemmas is given in Sect. 5.
By C we denote a positive constant independent of the functions and parameters concerned, but not necessarily the same at different occurrences. We assume that α ∈ (0, 1) in this paper and, for simplicity, we will not explicitly write this assumption in many places.

A Continuous Galerkin Time Stepping Method for (8), (9)
In this section, we shall introduce a continuous Galerkin time stepping method for solving (8), (9).
Let 0 = t 0 < t 1 < · · · < t N = T be a partition of [0, T ] and τ the time step size. We define the following trial space, with k = 0, 1, 2, . . . , N − 1 and q = 2, and the test space respectively. We remark that the trial space W 1 is a continuous function space and the test space is a discontinuous piecewise constant function space with respect to time t. The trial space W 1 is, due to the continuous restriction, of one order higher than the test space W 0 .

Error Estimates of the Time Stepping Method (11), (12)
In this section, we will show the error estimates of the abstract time stepping method (11), (12) by using the Laplace transform method proposed originally by Lubich et al. [25] and developed by Jin et al. [11,12], Yan et al. [42] and Wang et al. [38], etc.

The Homogeneous Case with f = 0 and u 0 = 0
In this subsection, we will consider the error estimates of the time stepping method (11), (12) with f = 0 and u 0 = 0. We thus consider the following homogeneous problem, The abstract time stepping method (11), (12) for solving (14) is now reduced to, with V 0 h = 0, We then have the following theorem: . . , be the solutions of (14) and (15), (16), Step 1: Find the exact solution of (14).
For I 1 , we have, by the resolvent estimate (7), with some constant c > 0, where we use the variable change z = re iθ in the second inequality above.
For I 2 , we have, by Lemma 2, Together these estimates complete the proof of Theorem 1. (28), we follow the approach in Lubich et al. [25, (3.9)] and Jin et al. [11] and obtain the formula for V n h by using the Laplace transform method. By Lemma 1, we may show that V n h is well defined since z α τ ∈ Σ θ for some θ ∈ (π/2, π). An alternative approach for obtaining the formula of V n h is given in Xie et al. [17], where the authors do not apply the Taylor expansion of the analytic function around the origin, and instead directly calculate the inverse Laplace transform V n h = 1 2πi a+iπ a−iπ e nz V h (z) dz for some a > 0, see [17,Lemma 3.6]. In particular, they obtain the stability estimate of V n h , i.e., [17,

Remark 3
In the above estimates for I 1 and I 2 , the constants c and C depend on the angle θ of the integral path Γ . Moreover, the constant C will tend to ∞ as θ → π/2. In our proof of Lemma 1, we require that θ is sufficiently close to π/2, so that the constant in Theorem 1 could be very large. The similar remark is also valid for the error estimates in Theorem 2 in the inhomogeneous case in the next section.

The
Inhomogeneous Case with f = 0 and u 0 = 0 In this subsection, we will consider the error estimates of the time stepping method (11), (12) with f = 0 and u 0 = 0. We thus consider the following inhomogeneous problem, The abstract time stepping method (11), (12) for solving (29) is now reduced to, with V 0 h = 0, We then have the following theorem:
Further we observe that lim z→∞ z α ψ(z) = 0. Hence we get which implies that, by (41), Therefore, by (40), Note that lim a→∞ e (n− j)a a −α = 0 for j ≥ n with any fixed n = 1, 2, . . . , which implies that (39) holds. We next consider the error estimates v h (t n ) − V n h . Subtracting (34) from (32) By Lemma 3, we have, with any small > 0, Together these estimates complete the proof of Theorem 1.

Numerical Simulations
In this section, we will consider some numerical examples for solving both fractional ordinary differential equations and subdiffusion problems by using the time stepping method (11), (12).

Example 1
Our first example is a homogeneous problem. Choose g(t) = 0 in (46) and the initial value y 0 = 1 in (47). In this case the problem has the following exact solution, with α ∈ (0, 1), where E α,1 (z) denotes the Mittag-Leffler function. We choose T = 2 and λ = 1. The exact solution can be calculated by using the MAT-LAB function mlf.m. We obtain the approximate solutions with the different step sizes τ = 1/20, 1/40, 1/80, 1/160. In Table 1, we observe that the experimentally determined convergence order of the time stepping method (48), (49) is about O(τ 2 ) which is better than theoretical convergence order O(τ 1+α ) for small α ∈ (0, 1). In Table 1, we also compare the errors and convergence orders of the time stepping method (48), (49) with the popular L1 scheme [11] and the modified L1 scheme [42]. It is well known that the L1 scheme has only O(τ ) convergence due to the singularity of the solution of the fractional differential equation [11] . After correcting the starting step, the modified L1 scheme has the optimal convergence order O(τ 2−α ) [42]. We observe that the time stepping method (48), (49) indeed captures the singularities of the problem more accurately than the L1 and modified L1 schemes. We remark that the modified L1 scheme has the convergence order O(τ 2−α ), α ∈ (0, 1), however in order to observe this convergence order for small α, we need to consider the error at sufficiently large T . For example, we indeed observe the convergence order O(τ 2−α ) of the modified L1 scheme for α = 0.3 when we choose T ≥ 10. In Table 1, the convergence order of the modified L1 scheme is not close to O(τ 2−α ) for α = 0.3 since T = 2 is not big enough to observe the required convergence order. Since the purpose of this paper is to show the convergence orders of the time stepping method (48), (49), we will not study further the numerical behaviors of the modified L1 scheme.

Example 3
The final example for the fractional ordinary differential equation is an inhomogeneous problem with nonzero initial value. We choose λ = 1 and assume that the exact solution of (46) is y(t) = t β + 1, β > 0. The initial value y 0 = 1 and We choose T = 1 and obtain the approximate solutions with the different step sizes τ = 1/10, 1/20, 1/40, 1/80, 1/160. In Table 3, we show the experimentally determined orders of convergences with β = 1.1 for the time stepping method (48), (49). We also observe that the convergence orders are about O(τ 2 ) which are better than the theoretical convergence order O(τ 1+α ) for small α ∈ (0, 1).

Example 4
In this example, we shall consider a homogeneous subdiffusion problem. We choose f (t, x) = 0 and the initial value u 0 (x) = x(1 − x) in (50), (52). In this case, the exact solution is . In our numerical simulation, we let T = 2. Choose the space step size h = 2 −10 and the different time step sizes τ = 1/10, 1/20, 1/40, 1/80, 1/160, we get the different approximate solutions. The exact solution is calculated by using the MATLAB function mlf.m.
We observe that, in Table 4, the experimentally determined convergence orders are better than the theoretical convergence orders O(τ 1+α ). For α > 0.5, the table shows a second order convergence rate. We also compare the errors and convergence orders of the method (48), (49) with the popular L1 scheme [11] and the modified L1 scheme [42]. We see that the time stepping method (48), (49) captures the singularities of the problem more accurately than the L1 and modified L1 schemes.

Example 6
Consider an inhomogeneous subdiffusion problem with nonzero initial value. Assume that the exact solution of (50)-(52) is u(t, x) = (t β + 1)x(1 − x), β > 0 and The initial value u 0 (x) = x(1 − x). We observe that the experimentally determined convergence orders depend on the smoothness of the solution with respect to the time variable t. In our numerical simulation, we choose β ∈ (1, 2) which implies that u(·, x) ∈ C 1 [0, T ], but u(·, x) / ∈ C 2 [0, T ] for any fixed x. We choose T = 2 and the space step size h = 2 −10 and the different time step sizes τ = 1/10, 1/20, 1/40, 1/80, 1/160 to get the approximate solutions. In Table 6, we show the experimentally determined orders of convergence with β = 1.1, and we see that the convergence orders are consistent with our theoretical results.

Example 7
Consider an inhomogeneous subdiffusion problem in two-dimensional case. With 1], assume that the exact solution of (50)-(52) is We use the same parameters and the time and space step sizes as in Example 6 except in this example we need space partitions in both x 1 and x 2 directions. In Table 7, we show the experimentally determined convergence orders with β = 1.1. We see that the numerical results are consistent with the theoretical results for this example.
Proof We first note that, after a direct calculation, with ζ = e −z , where Li α−2 is the polylogarithm function, see Jin et al. [11]. We will divide the proof of this lemma into the following 3 steps.
Together these estimates complete the proof of Lemma 3.