Efficient finite element numerical solution of the variable coefficient fractional subdiffusion equation

Based on the weighted and shifted Grünwald formula, a fully discrete finite element scheme is derived for the variable coefficient time-fractional subdiffusion equation. Firstly, the unconditional stable and convergent of the fully discrete scheme in L1(H1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L^{1}(H^{1})$\end{document}-norm is proved. Secondly, through a new estimate approach, the superclose properties are obtained. The global superconvergence order O(τ2+hm+1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{O}(\tau ^{2}+h^{m+1})$\end{document} is deduced with the help of interpolation postprocessing technique. Finally, some numerical results are provided to verify the theoretical analysis.


Introduction
In the recent few decades, the remarkable applications of fractional calculus in diverse engineering fields have been gradually realized and, meanwhile, the discussion of the related fractional partial differential equations (FPDEs) becomes a hot topic of many scholars. For instance, the reader may refer to the work [1][2][3][4][5][6][7][8][9][10][11]. Like traditional partial differential equations, most commonly the exact solutions of the FPDEs are not available. Even if their solutions can be found, they are usually in the form of series, which are difficult to evaluate. So the numerical investigation of the FPDEs has been a vital topic in recent years.
Fractional derivatives are nonlocal operators and they have the character of history dependence. When solving approximation solutions for time FPDEs on the current time layer, one needs to retain the information about all the previous time layers, which makes the storage expensive. Based on this aspect, developing high-order numerical methods for solving time FPDEs are quite valuable. In the present studies, we propose high-order methods by first improving the spatial accuracy with the compact difference operator for time FPDEs, which needs few grid points to produce a highly accurate solution. Gao and Sun [12] applied the L1 approximation for the time-fractional derivative and developed a compact finite difference scheme for the fractional sub-diffusion equation. The stability and L ∞ convergence are proved by the energy method. Wang and Vong [13,14] estab-lished high-order schemes for the fractional diffusion-wave equation and the modified anomalous fractional subdiffusion equation and fractional Klein-Gordon equation, respectively.
Another way for developing a high-order method for solving time FPDEs is to improve the temporal accuracy. Recently, a weighted and shifted Grünwald difference (WSGD) operator was developed for solving space fractional diffusion equations in [15]. The WSGD operator to approximate the Riemann-Liouville derivative can achieve second-order accuracy. Gao et al. [16] presented a new numerical differentiation formula for the Caputo fractional derivative, called the L1-2 formula. Then the difference schemes for the time-fractional sub-diffusion equations in bounded and unbounded spatial domains were constructed based on the L1-2 formula. However, they did not give the analysis on the stability and convergence of the obtained schemes. Later, Gao et al. [17] proposed a second-order difference scheme based on the certain superconvergence at some particular points of the fractional derivative by the first-order GL formula. The obtained scheme can achieve the global second-order accuracy in time. Recently, Alikhanov [18] proposed a new difference analog of the Caputo fractional derivative, called the L2-1 σ formula. The stability and convergence of these schemes in L 2 -norm were analyzed by the energy method.
The finite element method is one of the effective numerical methods for solving classical PDEs. For the FPDEs, finite element method also can be a useful and effective numerical method. In recent years, some valuable papers were concerned with the finite element method for the FPDEs. Roop and Ervin [19][20][21] investigated the theoretical framework for the Galerkin finite element approximation to some kinds of the FPDEs. Jiang and Ma [22] analyzed the finite element method and showed that the optimal convergent rate O(τ 2-α + N -m ) can be obtained, where m measures the regularity of the solution in space. Jin et al. [23][24][25] constructed some finite element schemes for solving some kinds of the FPDEs. They proved that those schemes were stable and the numerical solution was convergent. In [26] a fully discrete finite element scheme for solving the two-dimensional space-and time-fractional Bloch-Torrey equations is developed and the stability and convergence estimate of the fully discrete scheme are proved. Wang and Yang [27] studied the wellposedness of a variable coefficient conservative fractional elliptic differential equation and its weak formulation. Recently, Ren et al. [28,29] presented two fully discrete schemes for the diffusion equations and diffusion-wave equations, respectively. The unconditional stability and superconvergence error estimates of the obtained schemes are investigated using the integral identities and postprocessing techniques. The optimal time accuracy O(τ 2-α ) (0 < α < 1) and O(τ 3-α ) (1 < α < 2) is obtained, respectively. However, according to current knowledge of the authors, there are few studies on the numerical treatment of FPDEs, high efficient numerical methods of time and space superconvergence analysis are still limited.
Many researchers have observed that for certain classes of problems the rate of convergence of the finite element solution and/or its derivatives at some special points exceeds the best global rate. This phenomenon has been termed "superconvergence" and has been analyzed mathematically because of its practical importance in engineering computations. This motivates us to consider the superclose and global superconvergence properties of the FPDEs. The main goal of this paper is to construct a high-order fully discrete finite element scheme and establish the corresponding superconvergence error estimate. The method follows the idea of the weighted and shifted GL difference operators [14,15]. Based on the equivalence of Riemann-Liouville derivative and Caputo derivative under some regularity assumptions, we derive a second-order accuracy formula to approximate Caputo fractional derivative and the standard Galerkin finite element method approach for the spatial discretization. The unconditional stability and L 1 (H 1 )-norm convergence are proved rigorously.
The outline of this paper is as follows. In Sect. 2, some preliminary numerical formulas and useful lemmas are prepared. In Sect. 3, a fully discrete scheme for the variable coefficient fractional subdiffusion equation is developed and the unconditional stability of the obtained scheme is proved. In Sect. 4, the superclose and superconvergence analysis for the scheme are presented. In Sect. 5, some numerical examples are presented to verify our theoretical results. Some conclusions are given in the last section. Throughout this paper, the notation c denotes a generic constant, which may differ at different occurrences, but it is always independent of the mesh size h, the time step size τ . Let W k,p (Ω) the standard Sobolev space of k-differential functions in L p (Ω), its norm by · k,p , and the norm of H k (Ω) by · k . When k = 0, we let L 2 (Ω) denote the corresponding space defined on Ω with norm · .

Preliminaries
In this section, some useful notations, lemmas and formulas will be prepared for the forthcoming work.
For temporal discretization, we divide the interval [0, T] into N -subintervals with τ = T/N and t k = kτ , k = 0, 1, 2, . . . , N . Let u n = u(x, y, t n ). In order to develop a secondorder approximation of the Riemann-Liouville fractional derivative, we introduce the shifted Grünwald approximation to the Riemann-Liouville fractional derivative given by (cf. [30]) where r is an integer and g (α) . Based on the method introduced by [14,15], we suppose that f ∈ C 2+α (R), then it hold that uniformly holds in t ∈ R as τ → 0. We have By using a simple calculation the left-hand side of (2) is equivalent to the following equation: (3)

Lemma 1 ([14, 17]) Let λ (α) k be defined by (3), then, for any positive integer m and real vector
In order to establish fully discrete finite element scheme, we introduce some notations as follows. We assume that Ω can be partitioned by a rectangular mesh. Let h = {K} be a rectangular mesh over Ω with size h. Let V h ⊂ H 1 0 (Ω) be a standard finite element space where Q m (K) denotes the space of polynomial functions of total degree lower or equal to m with respect to x and y on K .

Stability for fully discrete finite element scheme
In this subsection, we discuss the stability for the fully discrete finite element scheme. First, we consider the following two-dimensional variable coefficient time-fractional subdiffusion equation: where 0 < α < 1 and C 0 D α t denotes the Caputo fractional derivative, the symbol ∇· and ∇ denote the divergence and gradient operators, respectively, f (x, y, t) is the given smooth function. It is without loss of generality to assume that u(x, y, 0 y) and consider the problem with respect to v(x, y, t). b(x, y) is a smooth and bounded function, which satisfies following assumption: In the current work, we assume that the function u(x, y, t) can be extended by zero from the time bounded domain [0, T] to R. Noticing the equivalence between the Riemann- and using (2), we discretize the time-fractional derivative as follows: According to (7), we arrive at the following fully discrete finite element scheme for the problem (5) In order to obtain the stability and convergence of the fully discrete scheme (8), we prove firstly the following priori estimate.

Theorem 1 Let w k
h ∈ V h be the solution of the following scheme: Then we have Proof Taking v h = w k h in (9), using Cauchy-Schwarz inequality, we have Summing up the above inequality for k from 1 to n, we obtain Noticing w 0 h = 0, when τ is sufficiently small, with the help of Lemma 1 yields This completes the proof.
From Theorem 1, we can obtain the following stability conclusion.

Theorem 2
The fully discrete finite element scheme (8) is unconditionally stable with respect to the inhomogeneous term f .

Superconvergence estimate for the fully discrete finite element scheme
In this section, we discuss the error estimate for the fully discrete scheme. In order to analyze the spatial discretization error, we assume that the solution is sufficiently smooth. First, let us introduce the following two lemmas.
Proof Using (7) and (4) and recalling the analytical method and tools in the proof of [28,29], we similarly obtain the desired result.
By virtue of Lemmas 2 and 3, now we carry out the rigorous error analysis for the fully finite element scheme (8).
, let u k h and I h u k be the finite element solution and finite element interpolation of u(t k ), respectively. Then the following supercloseness estimate holds for 1 ≤ n ≤ N : where the constant c depends on Ω, T, α and b, but it is independent of h and τ .
Proof We split the error By comparing (5) and (8), we have the following error equation: where Taking into account v h = θ k in (11), then it follows that For the second term in the left-hand side of (12), noticing (6), we have b∇θ k , ∇θ k ≥ β 1 ∇θ k 2 .
Furthermore, we can obtain the global superconvergence result of the fully discrete scheme (8) by virtue of a proper postprocessing technique introduced by [32]. For the sake of completeness, we present the construction method of the operator Π 2h . First we combine four neighboring elements into a big elementK = 4 i=1 K i , which has the following properties (cf. [32,33]): where 2h consists of the four small elements K i (i = 1, 2, 3, 4) in h .
We can deduce the following global superconvergence result.
Theorem 4 Suppose u k h are the solution of the fully discrete finite element (8) and under the conditions of Theorem 3. Then we have the following result for 1 ≤ k ≤ N : This completes the proof. Remark 4.1 It should point that using the WSGD operator to approximate R -∞ D α t u(x, y, t k ), the temporal derivative provided that ∂ n u(x,y,t) ∂t n | t=0 = 0 (n = 0, 1, . . . , 4), which can be guaranteed in view of the work in [34] and [35], is sufficient but not necessary. We will further illustrate this point by numerical experiments (Example 2).
Remark 4.2 In the present work, we obtain the superclose and superconvergence results in L 1 (H 1 )-norm through approach for WSDG fully discrete finite element scheme, which improve the conclusions in [17,22,29]. It seems that these results have never be seen in the existing literature.

Numerical experiments
In this section, some test examples based on piecewise bilinear polynomials will be performed to illustrate the computational efficiency and convergence rate of the proposed schemes.
Example 1 In (5), set b(x, y) = 1 and let Ω = [0, π] × [0, π] and T = 1. In order to test the convergence rate of the proposed methods, we consider the exact solution of the problem as follows: u(x, y, t) = t 2+α sin x sin y. Then f (x, y, t) and u 0 (x, y) are chosen corresponding to the exact solution, respectively. Firstly, the numerical accuracies of the fully finite element scheme in time is tested. For different α, the numerical results are computed using the developed scheme with varying temporal stepsizes and fixed sufficiently small spatial stepsizes. The second-order convergence of both the schemes in time can be observed from the data of Table 1.
Secondly, the numerical accuracies of the fully finite element scheme in space is verified by the example. With fixed sufficiently small temporal stepsizes, the u nu n h , u nu n h 1 , I h u nu n h 1 and u n -Π 2h u n h 1 norm errors and spatial convergence orders of the scheme are illustrated in Tables 2 and 3 when α = 0.1, 0.3, respectively. As predicted by the theoretical estimates, the second-order supercloseness and superconvergence of the scheme in space variable for computing this example are verified. The figures of exact solution (top) and numerical solution (bottom) of both schemes are shown in Fig. 1.  We have pointed out that the condition ∂ n u(x,y,t) ∂t n | t=0 = 0 (n = 0, 1, . . . , 4) is sufficient but unnecessary for the numerical accuracy of the proposed schemes. In this example, we are concerned with the computational results of the fully finite element scheme for the cases β = 2, 3/2, 1, respectively. From Table 4, we can see that the second-order convergence in temporal direction can be achieved. Tables 5 and 6 show that the supercloseness and superconvergence results in spatial direction can still be achieved.
Tables 7 and 8 report the numerical errors and convergence orders in temporal direction for the problems with β = 3/2 and β = 1, respectively. The second-order convergence of the scheme in temporal direction cannot be achieved for the above two cases. Namely, certain conditions on the derivative values of the function u(x, y, t) with respect to t at t = 0 up to a necessary order are essential to ensure the second-order convergence of the in time variable, whereas the condition ∂ n u(x,y,t) ∂t n | t=0 = 0 (n = 0, 1, . . . , 4) is only sufficient but unnecessary. The figures of exact solution (top) and numerical solution (bottom) of both schemes are shown in Fig. 2.
We test the numerical scheme (8) for this Example. The results at the time T = 1 with h = π/300 are reported in Fig. 3. We show the temporal errors as the function of τ for different α. The slope is two in good agreement with the theoretical result of Theo-

Conclusions
In this work, we have developed a high-order approximation fully discrete finite element scheme in time for the variable coefficient time-fractional subdiffusion equation which involves Caputo fractional derivative in time. By using some of the same analytical techniques as [28,29], a supercloseness approximation between the interpolation of the exact solution and the finite element solution is derived. Then the global superconvergence O(τ 2 + h m+1 ) in the L 1 (H 1 )-norm is also obtained with the aid of a suitable postprocessing method. The global second-order convergence of the obtained schemes in time variable is attainable. In future work, the numerical solutions for the nonlinear fractional differential equations will be investigated.