A weak Galerkin finite element approximation of two-dimensional sub-diffusion equation with time-fractional derivative

1.   Introduction

We consider the sub-diffusion initial-boundary value problem

where the spatial variable $\mathbf{x} = (x_1, x_2)$, $\Omega\subset { \mathbb{R}}^2$ denotes a bounded convex polygonal domain with boundary $\partial \Omega$, the variable diffusion coefficient $\mathcal{K}_\gamma(\mathbf{x})$ satisfies $\mathcal{K}_1\leq \mathcal{K}_\gamma(\mathbf{x}) \leq \mathcal{K}_2$ with positive constants $\mathcal{K}_1$ and $\mathcal{K}_2$, $f$ is the source term, the boundary value $\varphi$ and the initial function $u^0$ are given. Here $0 < \gamma < 1$ and $D^\gamma_{t}u(\mathbf{x}, t)$ is the fractional partial derivative of $u(\mathbf{x}, t)$ with respect to $t$ of order $\gamma$ in the Caputo form

with $\Gamma(\cdot)$ being the Gamma function.

Fractional differential models have been attracted considerable attention in recent years, to the authors knowledge, a novel technique for investigating the well-posedness of periodic solution for impulsive piecewise fractional functional differential equations was just developed in ^[1]. In this paper, we focus on the sub-diffusion equations with time-fractional derivative which usually models the solute transport in porous media with strong heterogeneity, in which the underlying particle movements have experience long waiting time and lead to a long tail in time direction. In addition, equation (1.1a) has also been applied in other fields, e.g., cytoplasmic crowding in living cells and diffusion of tracer microbeads in reconstituted ^[2,3]. In the last two decades, various effective numerical methods have been studied for sub-diffusion equations. Sun and Wu ^[4] derived an unconditionally stable finite difference method for one dimensional sub-diffusion equations by approximating the time fractional term with piecewise linear interpolation. This technique is well known as $L_1$ scheme with order $O(\tau^{2-\gamma})$ and has been widely used for solving the fractional differential equations with Caputo derivatives; for example, see ^{[5,6,7,8,9,10,11,12,13]}. Alikhanov studied a finite difference method for the sub-diffusion equation with multi-term variable-distributed order ^[5]. An unconditionally stable finite difference method for variable order sub-diffusion equation was presented by Chen et al. in ^[6]. Gao and Sun ^[7] developed a compact difference method for the sub-diffusion equation by applying a fourth order compact approximation for the space derivative. A space semidiscrete scheme and a fully discrete scheme based on the standard Galerkin finite element method using continuous piecewise linear functions were developed for a multi-term time-fractional diffusion equation by Jin et al ^[8]. Jin and Zhou also proposed an efficient Galerkin approximation scheme based on proper orthogonal decomposition for solving sub-diffusion equation ^[9]. In ^[11], Lin and Xu solved the sub-diffusion equation by a Legendre spectral collocation method in space. Ren et al. ^[12] proposed a compact difference method for the sub-diffusion equation with Neumann boundary conditions. Wang et al. considered a novel finite element method for the two-dimensional sub-diffusion equation with variable coefficients on anisotropic rectangular meshes in ^[13]. For the time-space fractional order nonlinear sub-diffusion equations ^[14], Li et al. proposed a semi-discrete and a fully discrete methods by using Galerkin finite element scheme for the space fractional operators and a finite difference scheme of $L_1$ type for the time Caputo derivative, respectively. Moreover, Galerkin finite element methods based on $L_1$ discretization for optimal control problems governed by time fractional diffusion equations were also studied ^[15,16,17].

Weak Galerkin finite element method was first introduced and analyzed by Wang and Ye in ^[18] for the second order elliptic equations. In general weak Galerkin approximations, differential operators (e.g., gradient, divergence, curl, Laplacian etc) are approximated by the weak forms as distributions for generalized functions.The local reconstruction of differential operators leads to a great flexibility in designing numerical schemes. Weak Galerkin technique has been successfully developed for linear parabolic equations ^[19], second order elliptic interface problems ^[20], biharmonic equations ^[21,22], Maxwell equations ^[23], Stokes equations ^[24], integro-differential equations ^[25,26] and Cahn-Hilliard equations ^[27] etc. Zhou et al. considered a weak Galerkin finite element method for multi-term time-fractional diffusion equations with one-dimensional space variable in ^[28]. However, weak Galerkin finite element methods for the two-dimensional sub-diffusion equation with time-fractional derivative are still limited. The goal of this paper is to present a fully weak Galerkin approximation scheme combining with $L_1$ discretization of Caputo time-fractional derivative to solve the sub-diffusion equation (1.1a) on triangular meshes.

The rest of this paper is organized as follows. In Section 2, we construct a fully discrete weak Galerkin finite element method for the sub-diffusion problem (1.1). In Section 3, we study the stability of the weak Galerkin scheme. The error estimates in $L^2$ and discrete $H^1$ norms are established in Section 4. Finally, we carry out numerical experiments to verify the convergence rate of the proposed scheme.

2.   Weak Galerkin finite element formulations

2.1.   Notation

For a multi-index $\alpha = (\alpha_1, \alpha_2)$, we denote its degree $|\alpha| = \alpha_1+\alpha_2$ and spatial weak derivative operator $\partial^\alpha = (\partial/\partial x_1)^{\alpha_1} (\partial/\partial x_2)^{\alpha_2}$. We use $L^2(\Omega)$ to denote the space of measurable functions whose square is Lebesgue integrable in domain $\Omega$ with the inner product and norm as

For a nonnegative integer $s$, we use

to denote the usual Sobolev space equipped with the norm

We also use the space $H^s_0(\Omega)$ as the closure of $C_0^\infty(\Omega)$ in $H^s(\Omega)$, where $C_0^\infty(\Omega)$ consists of functions in $C^\infty(\Omega)$ that have compact support in $\Omega$. Specially, $H^0_0(\Omega) = H^0(\Omega) = L^2(\Omega)$. For measurable function $w:[0, T]\rightarrow H^s(\Omega)$, let

and space $L^\infty(0, T;H^s(\Omega)) = \{w:\|w\|_{s, \infty} < +\infty\}$.

2.2.   A numerical discretization of the time-fractional derivative

We define a uniform partition on $[0, T]$ by $t_m = m\tau$ for $m = 0, 1, \cdots, M$ with $\tau = T/M$ and $M$ being a positive integer. By using the $L_1$ discretization of the Caputo time-fractional derivative at any time $t = t_m$ ^[11,29,30], we have

where $a_j = (j+1)^{1-\gamma}-j^{1-\gamma}$ for $j = 0, 1, \cdots, m-1$ and $L_t^{\gamma}$ is the corresponding difference operator. Then equation (1.1a) can be rewritten in a semidiscrete form

Occasionally an alternative form of operator $L_t^{\gamma}$ is employed as

for some $g(t)$ defined on $[0, T]$. Two properties of $L_1$ approximation are listed in the following lemma ^[4].

Lemma 2.1. 1. Given $0 < \gamma < 1$ and sequence $\{a_j\}_{j = 0}^{\infty}$ with $a_j = (j+1)^{1-\gamma}-j^{1-\gamma}$, we have that $\lim\limits_{j\rightarrow\infty}a_j = 0$, $a_0 = 1$ and

2. For $g(t)\in C^2([0, T])$ and $1\leq m \leq M$, we have

2.3.   A fully discrete weak Galerkin finite element scheme

In this part, we design a fully discrete weak Galerkin finite element scheme for the initial-boundary value problem (1.1). We consider the space of discrete weak functions and the discrete weak operator introduced in ^[18]. Let ${ \mathcal{T}}_h = \{K\}$ be a quasi-uniform triangulation partition of domain $\Omega$ with mesh size $h$. For each $K\in \mathcal{T}_h$, denote its interior and boundary by $K^0$ and $\partial K$, respectively. We choose the following kind of weak finite element space

where $l$ is a nonnegative integer, $P_l(K^0)$ and $P_l(\partial K)$ are the sets of polynomials with degree no more than $l$ on $K^0$ and each line segment of $\partial K$, respectively. Let $S_h^0(l)$ be the subspace of $S_h(l)$ with vanishing boundary values on $\partial \Omega$:

For each $v = \{v_0, v_b\}\in S_h(l)$, its discrete weak gradient $\nabla_d v\in RT_l(K)$ on each element $K$ is defined by the following local linear equation

where $RT_l(K)$ is the usual Raviart-Thomas element of order $l$ ^[31]. We consider a local $L^2$ projection $Q_h u = \{Q_0 u, Q_b u\}$ onto $P_l(K^0)\times P_l(\partial K)$, and a global elliptic projection $E_h u = \{E_0 u, E_b u\}$ onto the discrete weak space $S_h(l)$ defined by the following variational equation

with Dirichlet boundary condition $E_b u = Q_b \varphi(\mathbf{x}, t)$ for each $t\in[0, T]$. Here ${ a}(\cdot, \cdot)$ is the weak bilinear form defined by

Notice that ${ a}(v, v)\geq 0$ for any $v\in S_h(l)$ and the solvability of equation (2.9) has been proved in ^[18]. Then a fully discrete weak Galerkin finite element scheme based on the semidiscrete equation (2.6) and the discrete weak gradient operator $\nabla_d$ can be given as: find $u_h^m = \{u_0^m, u_b^m\}\in S_h(l)$ for $1\leq m\leq M$ satisfying $u_b^m = Q_b\varphi(\mathbf{x}, t_m)$ on $\partial \Omega$ and equation

with initial condition $u_h^0 = E_hu^0(\mathbf{x})$, where $f^m = f(\mathbf{x}, t_m)$ and $Q_b\varphi(\mathbf{x}, t_m)$ is the $L^2$ projection for each boundary segment.

3.   Stability of fully discrete weak Galerkin scheme

We discuss the stability of the fully discrete weak Galerkin finite element scheme (2.11) in the following theorem which implies the existence and uniqueness for the solution of scheme (2.11).

Theorem 3.1. Assume that $u_h^m = \{u_0^m, u_b^m\}$ ($1\leq m \leq M$) is the solution of the fully discrete weak Galerkin finite element scheme (2.11) with homogenous Dirichlet boundary condition, there exists a constant $C$ depends only on $\gamma$ and $T$ such that

and

for $m = 1, \cdots, M$, where the $L^2$ norm $\|\cdot\|$ is defined in (2.1).

Proof. First, for $m = 1$, let $v = u_h^1$ in (2.11), we have $(L_t^\gamma u_0^1, u_0^1) \leq (f^1, u_0^1)$. Thus

where $\rho = \Gamma(2-\gamma)\tau^\gamma$. Notice that $a_0 = 1$ in Lemma 2.1 and $\|u_h^1\| = \|u_0^1\|$, we have

by the Cauchy-Schwarz's inequality. For $m\geq 2$, by choosing $v = u_h^m$ in (2.11) we have

which leads to $(L^{\gamma}_{t}u_0^m, u_0^m) \leq (f^m, u_0^m)$, i.e.,

Using the Cauchy-Schwarz's inequality again, we arrive at

Considering the properties of $\{a_j\}_{j = 0}^{\infty}$, we have the following estimate by induction hypothesis

Notice that

Let $C = \Gamma(1-\gamma)T^\gamma$, then we have $Ca_{M-1} > \rho$ which implies that

From (3.6) we obtain

Next, Let $v = L_t^\gamma u_h^m = \{L_t^\gamma u_0^m, L_t^\gamma u_b^m\}$ in (2.11), we have

Using $\epsilon-$inequality for the term on the right hand side, we get

By exchanging the order of linear operators $L_t^\gamma$ and $\nabla_d$, we have

That is

Thus $\|\nabla_d u_h^m\|^2$ can be handled similarly to (3.6) as

This completes the proof.

4.   Error analysis

In this section, we establish the error estimates in $L^2$ and discrete $H^1$ norms for the weak Galerkin finite element scheme (2.11). According to the error estimates in ^[18], we first have the following lemma.

Lemma 4.1. Assume that $u \in L^\infty(0, T;H^{s+1}(\Omega))$ is the exact solution to the sub-diffusion problem (1.1), and $E_hu$ is defined by (2.9). There exists a positive constant $C$ independent of $h$ such that

and

provided that the mesh-size $h$ is sufficiently small, where the norms $\|\cdot\|_{0, \infty}$ and $\|\cdot\|_{s+1, \infty}$ are defined in (2.4).

For $1\leq m \leq M$, we denote $u^m = u(\mathbf{x}, t_m)$ with $u$ being the solution of equation (1.1a). Let

Our main goal here is to bound $u_h^m - Q_h u^m = \xi^m + \eta^m$. Since $\eta^m = E_hu^m - Q_hu^m$ can be handled by Lemma 4.1, we just focus on $\xi^m$.

Let $v = \{v_0, v_b\}\in S_h^0(l)$ be any test function. By testing the equation (1.1a) against the first component $v_0$ we get

Subtracting (4.1) from (2.11) with $t = t^m$, and using (2.9), we have the following error equation

The error estimates for the weak Galerkin finite element method in $L^2$ and discrete $H^1$ norms are provided in the following Theorem 4.2.

Theorem 4.2. Assume that $u$, $E_hu$ and $u_h^m(m = 1, \cdots, M)$ are the solutions of (1.1), (2.9) and (2.11), respectively. If $u, u_t, u_{tt}\in L^\infty(0, T;H^{s+1}(\Omega))$ with $0\leq s \leq l+1$, then there exists a constant $C$ independent of $h$ and $\tau$ such that

and

where the norms $\|\cdot\|_{0, \infty}$ and $\|\cdot\|_{s+1, \infty}$ are defined in (2.4).

Proof. We rewrite the error equation (4.2) as

where

and

Notice that $\xi^0 = 0$, $\nabla_d \xi^0 = 0$ and $\xi^m|_{\partial \Omega} = 0$ for $m = 0, \cdots, M$. By Theorem 3.1, we have

According to the definition of operator $L_t^\gamma$, the term $\| w_1^m \|$ can be bounded by

By Lemma 4.1 we have

Applying the error result of $L_1$ discretization in Lemma 2.1, we can bound $w_2^m$ by

We substitute (4.7)–(4.8) into (4.6) and then get

for $1\leq m \leq M$. Using Lemma 4.1 and the triangular inequality, we finally obtain that

for $1\leq m \leq M$. This completes the proof.

5.   Numerical experiments

In this section we carry out numerical experiments on three examples to demonstrate the convergence rate of the weak Galerkin finite element scheme (2.11) for the sub-diffusion problem (1.1). For each example, we compute the sub-diffusion equation on a square $\Omega = (0, \frac{\pi}{r})\times(0, \frac{\pi}{r})$ with optional $r\in\{1, \pi\}$ and the time interval $[0, T] = [0, 1]$. A combination of discrete weak spaces with $l = 0$ will be employed, i.e., $S_h(0)$ and $RT_0(K)$ for any $K\in \mathcal{T}_h$, where $\mathcal{T}_h$ is a uniform triangulation with mesh size $h$ as shown in Figure 1. We denote the error at time $T = 1$ by $e_{\tau, h} = u_h^M - Q_h u^M$ with $\tau$ being the time step.

Example 5.1. We set sub-diffusion equation (1.1a) with constant coefficient $\mathcal{K}_\gamma = 1$ and source term

The exact solution is

which satisfies homogeneous Dirichlet boundary condition, and the initial value $u^0(\mathbf{x}) = 0$. Actually, different $u(\mathbf{x}, t)$ will be approximated by weak Galerkin solutions when distinct combination $(r, \beta, \gamma)$ is chosen.

The numerical results for Example 5.1 are presented in Tables 1 and 2. The second and the fourth columns of Table 1 show the $L^2$ and discrete $H^1$ norms of $e_h$ with $(r, \beta, \gamma) = (\pi, 2, 0.8)$, respectively. Notice that the time step $\tau = \frac{1}{400}$ is fixed in Table 1. The convergence orders are given by $\log_{\frac{h_1}{h_2}}\dfrac{\|e_{\tau, h_1}\|}{\|e_{\tau, h_2}\|}$ and $\log_{\frac{h_1}{h_2}}\dfrac{\|\nabla_d e_{\tau, h_1}\|}{\|\nabla_d e_{\tau, h_2}\|}$. As the mesh size $h$ is decreased, it is observed that the numerical solution of weak Galerkin scheme (2.11) is convergent with optimal rate $O(h^2)$ in $L^2$ and $O(h)$ in discrete $H^1$ norms, which coincides the theoretical results in Theorem 4.2. Table 2 shows the error behavior by scheme (2.11) with fixed space mesh size $h = \frac{\pi}{128}$. The second column reports the $L^2$ norm of $e_h$ with $(r, \beta, \gamma) = (1, 6, 0.8)$. The fifth column corresponds to the discrete $H^1$ norm with $(r, \beta, \gamma) = (1, 10, 0.8)$. We employ $\log_{\frac{\tau_1}{\tau_2}}\dfrac{\|e_{\tau_1, h}\|}{\|e_{\tau_2, h}\|}$ and $\log_{\frac{\tau_1}{\tau_2}}\dfrac{\|\nabla_d e_{\tau_1, h}\|}{\|\nabla_d e_{\tau_2, h}\|}$ to reflect the convergence order of the weak Galerkin scheme. The results indicate that the convergence rates are both approximately $O(\tau^{2-\gamma})$, which also supports Theorem 4.2.

Example 5.2. In this example, we test the case with variable diffusion coefficient

and exact solution (5.1). The corresponding source term $f({ \mathbf{x}}, t)$ is given by

The numerical data for Example 5.2 is listed in Tables 3 and 4. We fixed a time step $\tau = \frac{1}{800}$ and show the convergence rate in Table 3 with setting $(r, \beta, \gamma) = (\pi, 2, 0.8)$. It is observed that the numerical solution is convergent with optimal rate $O(h^2)$ in $L^2$ and $O(h)$ in discrete $H^1$ norms. For $(r, \beta, \gamma) = (1, 10, 0.8)$ and a fixed mesh space size $h = \frac{\pi}{128}$, the second and third columns of Table 4 report the $L^2$ norm of error $e_{\tau, h}$ and the corresponding convergence rates. And for a fixed mesh space size $h = \frac{\pi}{256}$, we show the discrete $H^1$ norm of error and the convergence rates in the fifth and sixth columns, respectively. The results in Table 4 indicate that the convergence rates are both approximately $O(\tau^{2-\gamma})$.

Example 5.3. We also test the case with variable diffusion coefficient $\mathcal{K}_\gamma$ in a tensor form as

which is a variable symmetric positive definite matrix for two-dimensional diffusion problem in some anisotropic media. The exact solution is

The numerical results for Example 5.3 are reported in Table 5, which has a similiar format of Table 3. It is shown that the data also coincides with our theoretical analysis, which illustrates the effectiveness of the weak Galerkin finite element scheme (2.11) for fractional sub-diffusion models in porous media with heterogeneity and anisotropy.

6.   Conclusions

We have employed the weak Galerkin finite element method with $L_1$ discretization to solve the two-dimensional anomalous sub-diffusion equation with time-fractional derivative. A fully discrete weak Galerkin finite element scheme was presented and the optimal order error estimates in $L^2$ and discrete $H^1$ norms were established based on the stability of the numerical scheme. Fractional sub-diffusion equations with various settings have been tested to demonstrate the accuracy of our method.

Acknowledgments

The authors would like to express their sincere thanks to the referees for the valuable comments and suggestions which helped to improve the original paper. This work was supported by the fund of the Natural Science of Shandong Province (ZR2014AM033).

Conflict of interest

All authors declare no conflicts of interest in this paper.