A stabilizer free spatial weak Galerkin finite element methods for time-dependent convection-diffusion equations

Abstract

In this paper, we propose a stabilizer free spatial weak Galerkin (SFSWG) finite element method for solving time-dependent convection diffusion equations based on weak form Eq. (4). SFSWG method in spatial direction and Euler difference operator Eq. (37) in temporal direction are used. The main reason for using the SFSWG method is because of its simple formulation that makes this algorithm more efficient and its implementation easier. The optimal rates of convergence of $\mathcal{O}(h^{k})$ and $\mathcal{O}(h^{k+1})$ in a discrete $H^{1}$ and $L^{2}$ norms, respectively, are obtained under certain conditions if polynomial spaces $(P_{k}(K),P_{k}(e),[P_{j}(K)]^{2})$ are used in the SFSWG finite element method. Numerical experiments are performed to verify the effectiveness and accuracy of the SFSWG method.

Keywords

Stabilizer free weak Galerkin methods weak gradient time-dependent convection-diffusion problem Euler difference operator error estimates

1. Introduction

In this paper, we are concerned with the development of numerical methods for the following partial differential equation with initial-boundary conditions using a stabilizer free spatial weak Galerkin (SFSWG) finite element method

$\displaystyle u_{t}-\nabla\cdot(\alpha\nabla u)+\boldsymbol{\beta}\cdot\nabla u% +cu=f\quad\mbox{ in }\Omega\times(0,T],$ (1) $\displaystyle u=0\quad\mbox{ on }\partial\Omega\times(0,T],$ (2) $\displaystyle u(x,0)=u^{0}\quad\mbox{ in }\Omega,$ (3)

where $\Omega$ is a polygonal domain in $\mathbb{R}^{2}$ with Lipschitz-continuous boundary, $T>0$ , $\alpha=\alpha(x)$ is the diffusion coefficient matrix, $\boldsymbol{\beta}=\boldsymbol{\beta}(x)$ is the convection coefficient and $c=c(x)$ is the reaction coefficient in relevant applications. We suppose that $\alpha_{ij}(x)\in[W^{1,\infty}(\Omega)]^{2\times 2}$ , $0\leqslant c(x)\leqslant M$ , $\boldsymbol{\beta}\in[W^{1,\infty}(\Omega)]^{2}$ and $c-\frac{1}{2}\nabla\cdot\boldsymbol{\beta}>c_{0}>0$ for some constant $c_{0}$ and there exist positive constants $\alpha_{m}$ such that

$\displaystyle\alpha_{m}\xi^{T}\xi\leqslant\xi^{T}\alpha(x)\xi,\forall\xi\in% \mathbb{R}^{2},x\in\Omega.$

.

[1] Let $V$ be a Hillbert space with inner product $(\cdot,\cdot)$ and $T>0$ . Define

$\displaystyle L^{2}(0,T;V)=\left\{v|v(\cdot,t)\in V,\forall 0\leqslant t% \leqslant T,\int_{0}^{T}\|v(\cdot,t)\|^{2}_{V}dt<\infty\right\}.$

The standard weak form of the Eqs (1)–(3) is to find $u\in L^{2}(0,T;H_{0}^{1}(\Omega))$ such that

$\displaystyle(u_{t},v)+(\alpha\nabla u,\nabla v)+(\boldsymbol{\beta}\cdot% \nabla u,v)+(cu,v)=(f,v),\forall v\in H_{0}^{1}(\Omega),t\in(0,T].$ (4)

The convection-diffusion equation appears in broad engineering and science applications such as materials sciences, fluid flows, and image processing. There are many efforts devoted to detecting accurate and efficient approximation solutions for solving the convection-diffusion equation (see [2, 3]) and time-dependent convection-diffusion equation (see [4]). In [4], Shenglan Xie et al. proposed a fully discrete scheme, based on weak Galerkin method in space and Crank-Nicolson method in time, for time-dependent convection-dominated diffusion equation. The major idea of the weak Galerkin finite element method is to substitute the classical differential operators with weak differential operators on functions with discontinuity plus a stabilizer part that ensures a sufficient weak continuity for the approximating function. Recently, the weak Galerkin method has been developed to solve the elliptic equations [5, 6, 7], the nonlinear convection-diffusion problems in 1D and 2D [8, 9], singularly perturbed reaction-diffusion problems [10], the biharmonic problems [11], the Helmholtz equation [12], the coupled Burgers’ equations [13], p-Laplacian Problem [14] and the Maxwell equations [15].

The stabilizer free weak Galerkin (SFWG) finite element method, recently introduced by Ye and Zhang in [16], refers to the numerical techniques for solving Poisson equation on polytopal meshes in 2D or 3D, where there is a $j_{0}>0$ so that as long as the degree $j$ of the weak gradient satisfies $j\geqslant j_{0}$ , the SFWG scheme will work and the optimal rate of convergence can be obtained. In [17], Al-Taweel and Wang proved the optimal degree of weak gradient of the SFWG method to improve the efficiency of SFWG and to avoid the numerical difficulties associated with using high degree weak gradients. The benefits of utilizing the stabilizer free weak Galerkin method compared to the existing standard weak Galerkin method are two-part: firstly, it has a simple formulation which is closer to the weak form Eq. (4) and thus the implementation of the SFWG finite element method is easier than that of the standard weak Galerkin method; secondly and more importantly, it is more efficient than the standard WG method. The goal of this article is to study a stabilizer free spatial weak Galerkin finite element method (SFSWG) for solving convection-diffusion equations with time-dependent Eqs (1)–(3) on uniform triangular partitions and then establish the stability and error analysis in a discrete $H^{1}$ norm and $L^{2}$ norm.

This paper is organized as follows: In Section 2, the weak gradient, weak divergence, and a semi discretized SFSWG finite element scheme for the convection-diffusion equation are presented. In Section 3, we introduce a fully discretized SFSWG scheme for the Eqs (1)–(3). Numerical experiment results are presented in Section 4 to verify the theoretical results. Finally, concluding remarks are summarized in Section 5.

2. Weak Galerkin finite element schemes

Let ${\cal T}_{h}$ be a partition of the domain $\Omega$ consisting of convex polygons in 2D. Suppose that ${\cal T}_{h}$ is shaped regular in the sense defined by Eqs (15) and (16). Let ${\cal E}_{h}$ be the set of all edges in ${\cal T}_{h}$ , and let ${\cal E}_{h}^{0}={\cal E}_{h}\setminus\partial\Omega$ be the set of all interior edges. For each element $K\in{\cal T}_{h}$ , denote by $h_{K}$ the diameter of $K$ , and $h=\max_{K\in{\cal T}_{h}}h_{K}$ the mesh size of ${\cal T}_{h}$ .

On each $K$ , let $P_{k}(K)$ be the space of all polynomials with degree $k$ or less. Let $V_{h}$ be the weak Galerkin finite element space associated with ${\cal T}_{h}$ defined as follows:

$\displaystyle V_{h}=\{v=\{v_{0},v_{b}\}:v_{0}|_{K}\in P_{k}(K),v_{b}|_{e}\in P% _{k}(e),K\in{\cal T}_{h},e\in\partial K\},$ (5)

where $k\geqslant 1$ is a given integer. In this instance, the component $v_{0}$ symbolizes the interior value of $v$ , and the component $v_{b}$ symbolizes the edge value of $v$ on each $K$ and $e$ , respectively. Let $V^{0}_{h}$ be the subspace of $V_{h}$ defined as:

$\displaystyle V_{h}^{0}=\{v:v\in V_{h},v_{b}=0\mbox{ on }\partial\Omega\}.$ (6)

Now, we define the weak gradient and weak divergence operators as follows:

.

(Weak gradient) For any $v=\{v_{0},v_{b}\}\in V_{h}+H^{1}(\Omega)$ , the weak gradient $\nabla_{w}v\in[P_{j}(K)]^{2}$ where $j>k$ , is defined on $K$ as the unique polynomial satisfying

$\displaystyle(\nabla_{w}v,\mathbf{q})_{K}=-(v_{0},\mathbf{\nabla\cdot q})_{K}+% \langle v_{b},\mathbf{q\cdot\mathbf{n}}\rangle_{\partial K},\forall\mathbf{q}% \in[P_{j}(K)]^{2},$ (7)

where $\mathbf{n}$ is the unit outward normal vector of $\partial K$ .

.

(Weak divergence) For any $v=\{v_{0},v_{b}\}\in V_{h}$ , the weak divergence $\boldsymbol{\beta}\cdot\nabla_{w}v\in P_{k}(K)$ is defined on $K$ as the unique polynomial satisfying

$\displaystyle(\boldsymbol{\beta}\cdot\nabla_{w}v,w)_{K}=-(v_{0},\nabla\cdot(% \boldsymbol{\beta}w))_{K}+\langle v_{b},\boldsymbol{\beta}\cdot\mathbf{n}w% \rangle_{\partial K},\forall w\in P_{k}(K).$ (8)

Next, we define four global projections $Q_{0}$ , $Q_{b}$ , $Q_{h}$ , and $\mathbb{Q}_{h}$ as follows.

.

For each element $K\in{\cal T}_{h}$ ,

$\displaystyle Q_{0}:L^{2}(K)\longrightarrow P_{k}(K),$ $\displaystyle Q_{b}:L^{2}(e)\longrightarrow P_{k}(e),$ $\displaystyle\mathbb{Q}_{h}:[L^{2}(K)]^{2}\longrightarrow[P_{j}(K)]^{2},$

are the $L^{2}$ projections onto the associated local polynomial spaces. Finally, we define a projection operator $Q_{h}v=\{Q_{0}v$ , $Q_{b}v\}\in V_{h}$ for $v\in H^{1}(\Omega)$ .

For simplicity, we adopt the following notations,

$\displaystyle(v,w)_{{\cal T}_{h}}=\sum_{K\in{\cal T}_{h}}(v,w)_{K}=\sum_{K\in{% \cal T}_{h}}\int_{K}\textit{v w dx},$ $\displaystyle\langle v,w\rangle_{\partial{\cal T}_{h}}=\sum_{K\in{\cal T}_{h}}% \langle v,w\rangle_{\partial K}=\sum_{K\in{\cal T}_{h}}\int_{\partial K}% \textit{v w ds}.$

For any $v=\{v_{0},v_{b}\}$ and $w=\{w_{0},w_{b}\}$ in $V_{h}$ , we define the bilinear form as

$\displaystyle A(v,w)=(\alpha\nabla_{w}v,\nabla_{w}w)_{{\cal T}_{h}}+(% \boldsymbol{\beta}\cdot\nabla_{w}v,w_{0})_{{\cal T}_{h}}+(cv_{0},w_{0}).$ (9)

The semi discrete SFSWG finite element scheme for the Eqs (1)–(3) is to seek $u_{h}=\{u_{0,h},u_{b,h}\}\in L^{2}(0,T;V_{h}^{0})$ , such that the following equation holds

$\displaystyle\left\{\begin{array}[]{l}((u_{0,h})_{t},v_{0})+A(u_{h},v)=(f,v_{0% }),\forall v=\{v_{0},v_{b}\}\in V_{h}^{0},\\ u_{h}(0)=Q_{h}u^{0}.\end{array}\right.$ (10)

.

Let elements $T_{1}$ and $T_{2}$ have $e$ as a common edge, we define jump for a scalar function $v$ on $e$ as

$\displaystyle[v]_{e}=\left\{\begin{array}[]{ll}v|_{\partial T_{1}}\mathbf{n}_{% 1}+v|_{\partial T_{2}}\mathbf{n}_{2},&e\in{\cal E}_{h}^{0},\\ v\mathbf{n},&e\in\partial\Omega,\end{array}\right.$ (11)

where $\mathbf{n}_{1}$ and $\mathbf{n}_{2}$ are the unit outward normal vectors on common edge $e$ for elements $T_{1}$ and $T_{2}$ , respectively.

.

Discrete Poinccaré inequality [18]. Let $k\geqslant 0$ . There is $\sigma>0$ , independent of $h$ , such that $\forall v\in V_{h}$ ,

$\displaystyle|\!|v_{0}|\!|^{2}\leqslant\sigma\left(\sum_{K\in{\cal T}_{h}}|\!|% \nabla v_{0}|\!|^{2}_{K}+\sum_{e\in{\cal E}_{h}}\frac{1}{h_{e}}|\!|[v_{0}]|\!|% _{e}^{2}\right).$ (12)

.

If ${\cal T}_{h}$ satisfies the shape regular conditions of Lemma 4, then there exists $C>0$ , so that $\forall v\in V_{h}^{0}$ ,

$\displaystyle C\sum_{K\in{\cal T}_{h}}\frac{1}{h_{k}}|\!|v_{0}-v_{b}|\!|_{% \partial k}^{2}\geqslant\sum_{e\in{\cal E}_{h}}\frac{1}{h_{e}}|\!|[v_{0}]|\!|_% {e}^{2}.$

Proof..

By the assumptions of Lemma 4, there exists a constant $C>0$ so that, for each $e\in\partial K_{1}\cap\partial K_{2}$ , where $K_{i}\in{\cal T}_{h},i=1,2$ ,

$\displaystyle\frac{1}{h_{e}}\leqslant C\frac{1}{h_{K_{1}}}\mbox{ and }\frac{1}% {h_{e}}\leqslant C\frac{1}{h_{K_{2}}}.$

It is easy to see that

$\displaystyle|\!|[v_{0}]|\!|_{e}\leqslant|\!|v_{0}-v_{b}|\!|_{\partial K_{1}}+% |\!|v_{0}-v_{b}|\!|_{\partial K_{2}}.$

The proof then follows from the fact that if $e\in\partial K\cap\partial\Omega$ , then $|\!|[v_{0}]|\!|_{e}=|\!|v_{0}-v_{b}|\!|_{e}$ . ∎

Now, we define an energy norm ${|\hskip-1.4454pt|\hskip-1.4454pt|}\cdot{|\hskip-1.4454pt|\hskip-1.4454pt|}$ on $V_{h}+H^{1}(\Omega)$ as:

$\displaystyle{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.44% 54pt|}^{2}=(\alpha\nabla_{w}v,\nabla_{w}v)_{{\cal T}_{h}}+(cv,v)_{{\cal T}_{h}% }=|\!|\!|v|\!|\!|_{1}^{2}+(cv,v)_{{\cal T}_{h}},\forall v\in V_{h}+H^{1}(% \Omega).$ (13)

An $H^{1}$ semi-norm on $V_{h}+H^{1}(\Omega)$ is defined as:

$\displaystyle\|v\|_{1,h}^{2}=\sum_{K\in{\cal T}_{h}}(\|\nabla v_{0}\|^{2}_{K}+% h_{K}^{-1}\|v_{0}-v_{b}\|^{2}_{\partial K}).$

.

(Coercivity) Let ${\cal T}_{h}$ be shape regular as defined in Lemma 4. Then, There exists a constant $\gamma>0$ , such that

$\displaystyle A(v,v)\geqslant\gamma{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|% \hskip-1.4454pt|\hskip-1.4454pt|}^{2},\forall v\in V_{h}^{0}.$ (14)

Proof..

Let $\gamma=\frac{c_{0}}{M}$ . From the Eq. (9) and the definition of the $|\!|\!|v|\!|\!|$ in Eq. (13), we get

$\displaystyle A(v,v)=|\!|\!|v|\!|\!|_{1}^{2}-\frac{1}{2}\sum_{K\in{\cal T}_{h}% }(v_{0},\nabla\cdot\boldsymbol{\beta}v_{0})+(cv_{0},v_{0})\geqslant|\!|\!|v|\!% |\!|_{1}^{2}+c_{0}(v_{0},v_{0})\geqslant\frac{c_{0}}{M}|\!|\!|v|\!|\!|^{2}=% \gamma|\!|\!|v|\!|\!|^{2},$

which complets the proof. ∎

The following lemma shows that $\|\cdot\|_{1,h}$ is equivalent to ${|\hskip-1.4454pt|\hskip-1.4454pt|}\cdot{|\hskip-1.4454pt|\hskip-1.4454pt|}_{1}$ defined in Eq. (13).

.

(see [17, 1]) Suppose that $\forall K\in{\cal T}_{h}$ , $K$ is convex with at most $\mu$ edges and satisfies the following regularity conditions: for all edges $e_{t}$ and $e_{s}$ of $K$

$\displaystyle|e_{s}|<\alpha_{0}|e_{t}|;$ (15)

for any two adjacent edges $e_{t}$ and $e_{s}$ the angle $\theta$ between them satisfies

$\displaystyle\theta_{0}<\theta<\pi-\theta_{0},$ (16)

where $1\leqslant\alpha_{0}$ and $\theta_{0}>0$ are independent of $K$ and $h$ . Let $j_{0}=k+\mu-2$ or $j_{0}=k+\mu-3$ when each edge of $K$ is parallel to another edge of $K$ . Denote $\textit{deg}\nabla_{w}v$ be the degree of weak gradient when $\textit{deg}\nabla_{w}v=j\geqslant j_{0}$ , then there exist two constants $C_{1}$ , $C_{2}>0$ , such that for each $v=\{v_{0},v_{b}\}\in V_{h}$ , the following hold true

$\displaystyle C_{1}\|v\|_{1,h}\leqslant{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|% \hskip-1.4454pt|\hskip-1.4454pt|}_{1}\leqslant C_{2}\|v\|_{1,h},$

where $C_{1}$ and $C_{2}$ depend only on $\alpha_{0}$ and $\theta_{0}$ .

.

Let $j=k+\mu-2$ . Then there exists two constants $C_{1}$ , $C_{2}>0$ such that for any $v=\{v_{0},v_{b}\}\in V_{h}^{0}$ , we have

$\displaystyle C_{1}\|v\|_{1,h}\leqslant{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|% \hskip-1.4454pt|\hskip-1.4454pt|}\leqslant C_{2}\|v\|_{1,h}.$

where $\mu$ is the number of edges in each element.

Proof..

$|\!|\!|v|\!|\!|\geqslant C_{1}|\!|v|\!|_{1,h}$ follows from Lemma 4 and

$\displaystyle{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.44% 54pt|}\geqslant{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.% 4454pt|}_{1}.$

By Lemma 1 and Lemma 2

$\displaystyle|\!|\!|v|\!|\!|_{1}^{2}+(cv,v)<|\!|\!|v|\!|\!|_{1}^{2}+C_{3}\left% (\sum_{K\in{\cal T}_{h}}|\!|\nabla v|\!|^{2}_{K}+\sum_{K\in{\cal T}_{h}}|\!|[v% _{0}]|\!|_{\partial K}^{2}\right)<C_{2}|\!|v|\!|_{1,h}^{2}.$

Thus, $\forall v\in v_{h}^{0}$

$\displaystyle C_{1}|\!|v|\!|_{1,h}\leqslant|\!|\!|v|\!|\!|\leqslant C_{2}|\!|v% |\!|_{1,h}.$

This completes the proof of the lemma. ∎

.

When $\mu=3$ then $j_{0}=k+1$ . If all $K$ ’s are parallelograms, then $\mu=4$ and $j_{0}=k+1$ .

Now, we will introduce a stability inequality for our semi discretized Eq. (10) of the Eqs (1)–(3).

.

Let $f\in L^{2}(\Omega\times(0,T])$ and $u_{h}$ be the semi discretized approximation defined in Eq. (10). Then, there exists a constant $C>0$ , such that the following holds

$\displaystyle\|u_{0,h}(t)\|^{2}\leqslant C\left(\|u_{0,h}(0)\|^{2}+\int_{0}^{t% }\|f(s)\|^{2}ds\right),\forall t\in(0,T],$ (17)

where $C$ is independent of $t$ and $h$ .

Proof..

For each fixed $t$ taking $v=u_{h}$ in Eq. (10), we obtain

$\displaystyle\frac{1}{2}\frac{d}{dt}\|u_{0,h}\|^{2}+A(u_{h},v)=(f,u_{0,h}).$

From Eq. (13), we find that

$\displaystyle\frac{d}{dt}\|u_{0,h}\|^{2}\leqslant 2(f,u_{0,h})\leqslant\|f\|^{% 2}+\|u_{0,h}\|^{2}.$ (18)

Then Eq. (17) follows from Eq. (18) and Gronwall inequality. ∎

Next, we list important inequalities which will be used in error estimates.

.

(Trace inequality, see [19]) On each element $K\in{\cal T}_{h}$ , the following trace inequality holds true:

$\displaystyle\|\varphi\|_{e}^{2}\leqslant C(h_{K}^{-1}\|\varphi\|_{K}^{2}+h_{K% }\|\nabla\varphi\|_{1,K}^{2}),\varphi\in H^{1}(K),$ (19)

for some constant $C$ .

.

(Inverse Inequality, see [19]) There exists a constants $C$ such that for any piecewise polynomial $\varphi|_{K}\in P_{k}(K)$ ,

$\displaystyle|\!|\nabla\varphi|\!|_{K}\leqslant Ch^{-1}_{K}|\!|\varphi|\!|_{K}% ,\forall K\in{\cal T}_{h}.$ (20)

The following lemma presents estimates for the projection operator $Q_{0}$ and $\mathbb{Q}_{h}$ .

.

[20] Let $\varphi\in H^{k+1}(\Omega)$ and ${\cal T}_{h}$ be a finite element partition of $\Omega$ satisfying the shape regularity conditions Eqs (15) and (16). Then, the $L^{2}$ projections $Q_{0}$ and $\mathbb{Q}_{h}$ satisfy

$\displaystyle\sum_{K\in{\cal T}_{h}}(\|\varphi-Q_{0}\varphi\|_{K}^{2}+h_{K}^{2% }\|\nabla(\varphi-Q_{0}\varphi)\|_{K}^{2})\leqslant Ch^{2(s+1)}\|\varphi\|_{s+% 1}^{2},0\leqslant s\leqslant k,$ $\displaystyle\sum_{K\in{\cal T}_{h}}(\|\nabla\varphi-\mathbb{Q}_{h}\nabla% \varphi\|_{K}^{2}+h_{K}^{2}|\nabla\varphi-\mathbb{Q}_{h}\nabla\varphi)|_{1,K}^% {2})\leqslant Ch^{2s}\|\varphi\|_{s+1}^{2},0\leqslant s\leqslant k.$

The following lemma is given in [1] and will also be needed in error estimates. For the completeness we’ll give the proof here.

.

Let $Q_{h}$ and $\mathbb{Q}_{h}$ be the projection operators defined in definition (4) and $\phi\in H^{k+1}(\Omega)$ . Then for each element $K\in{\cal T}_{h}$ , we have

$\displaystyle\|\mathbb{Q}_{h}(\nabla\phi)-\nabla_{w}Q_{h}\phi\|_{K}\leqslant Ch% ^{k}_{K}|\phi|_{k+1,K}.$ (21)

Proof..

By definition (2) and integration by parts, we have

$\displaystyle(\mathbb{Q}_{h}(\nabla\phi),\mathbf{q})_{K}=-(\phi,\nabla\cdot% \mathbf{q})_{K}+\langle\phi,\mathbf{q}\cdot\mathbf{n}\rangle_{\partial K},$ (22) $\displaystyle(\nabla_{w}Q_{h}\phi,\mathbf{q})_{K}=-(Q_{0}\phi,\nabla\cdot% \mathbf{q})_{K}+\langle Q_{b}\phi,\mathbf{q}\cdot\mathbf{n}\rangle_{\partial K},$ (23)

for any $\mathbf{q}\in[P_{j}(K)]^{2}$ . Subtracting Eq. (23) from Eq. (22), using integration by parts, trace inequality Eq. (19), and the inverse inequality Eq. (20) yields

$\displaystyle(\mathbb{Q}_{h}(\nabla\phi)-\nabla_{w}Q_{h}\phi,\mathbf{q})_{K}=-% (\phi-Q_{0}\phi,\nabla\cdot\mathbf{q})_{K}+\langle\phi-Q_{b}\phi,\mathbf{q}% \cdot\mathbf{n}\rangle_{\partial K}=(\nabla(\phi-Q_{0}\phi),\mathbf{q})_{K}+% \langle Q_{0}\phi-Q_{b}\phi,\mathbf{q}\cdot\mathbf{n}\rangle_{\partial K}% \leqslant\|\nabla(\phi-Q_{0}\phi)\|_{K}\|\mathbf{q}\|_{K}+Ch^{-\frac{1}{2}}_{K% }\|\phi-Q_{0}\phi\|_{\partial K}\|\mathbf{q}\|_{K}\leqslant Ch^{k}_{K}|\phi|_{% k+1,K}\|\mathbf{q}\|_{K},$

which implies Eq. (21). ∎

.

Let $\phi\in H^{1}(\Omega)$ . Then for all $v\in V_{h}$ , we have

$\displaystyle(\alpha\nabla\phi,\nabla v_{0})_{K}=(\alpha\nabla_{w}(Q_{h}\phi),% \nabla_{w}v)_{K}+\langle(\mathbb{Q}_{h}(\alpha\nabla\phi)\cdot\mathbf{n},v_{0}% -v_{b}\rangle_{\partial K}{}+(\mathbb{Q}_{h}(\alpha\nabla\phi)-\alpha\nabla_{w% }Q_{h}\phi,\nabla_{w}v)_{K}.$ (24)

Proof..

Let $\mathbb{Q}_{h}(\alpha\nabla\phi)=\mathbf{q}$ and $Q_{h}\phi=P$ . Then

$\displaystyle(\mathbf{q},\nabla v_{0})_{K}=-(\nabla\cdot\mathbf{q},v_{0})_{K}+% \langle\mathbf{q}\cdot\mathbf{n},v_{0}\rangle_{\partial K},$ (25) $\displaystyle(\mathbf{q},\nabla_{w}v)_{K}=-(\nabla\cdot\mathbf{q},v_{0})_{K}+% \langle\mathbf{q}\cdot\mathbf{n},v_{b}\rangle_{\partial K}.$ (26)

Subtracting Eq. (26) from Eq. (25) yields

$\displaystyle(\mathbf{q},\nabla v_{0})_{K}=(\mathbf{q},\nabla_{w}v)_{K}+% \langle\mathbf{q}\cdot\mathbf{n},v_{0}-v_{b}\rangle_{\partial K}=(\alpha\nabla% _{w}P,\nabla_{w}v)_{K}+(\mathbf{q}-\alpha\nabla_{w}P,\nabla_{w}v)_{K}+\langle% \mathbf{q}\cdot\mathbf{n},v_{0}-v_{b}\rangle_{\partial K},$

which completes the proof. ∎

.

For all $v\in V_{h}$ , we have

$\displaystyle(\boldsymbol{\beta}\cdot\nabla u,v_{0})_{{\cal T}_{h}}=(% \boldsymbol{\beta}\cdot\nabla_{w}Q_{h}u,v_{0})_{{\cal T}_{h}}-\ell_{% \boldsymbol{\beta}}(u,v),$ (27)

where

$\displaystyle\ell_{\boldsymbol{\beta}}(u,v)=(u-Q_{0}u,\nabla\cdot(\boldsymbol{% \beta}v_{0}))_{{\cal T}_{h}}-\langle u-Q_{b}u,\boldsymbol{\beta}\cdot\mathbf{n% }(v_{0}-v_{b})\rangle_{\partial{\cal T}_{h}}$

Proof..

From integration by parts and the definition 3 we get

$\displaystyle(\boldsymbol{\beta}\cdot\nabla u,v_{0})_{{\cal T}_{h}}=-(u,\nabla% \cdot(\boldsymbol{\beta}v_{0}))_{{\cal T}_{h}}+\langle u,\boldsymbol{\beta}% \cdot\mathbf{n}v_{0}\rangle_{\partial{\cal T}_{h}}=-(Q_{0}u,\nabla\cdot(% \boldsymbol{\beta}v_{0}))_{{\cal T}_{h}}+(Q_{0}u-u,\nabla\cdot(\boldsymbol{% \beta}v_{0}))_{{\cal T}_{h}}+\langle Q_{b}u,\boldsymbol{\beta}\cdot\mathbf{n}v% _{0}\rangle_{\partial{\cal T}_{h}}+\langle u-Q_{b}u,\boldsymbol{\beta}\cdot% \mathbf{n}v_{0}\rangle_{\partial{\cal T}_{h}}=(\boldsymbol{\beta}\cdot\nabla_{% w}(Q_{h}u),v_{0})_{{\cal T}_{h}}-\ell_{\boldsymbol{\beta}}(u,v),$

where in the last equality we have used the fact that $\langle u-Q_{b}u,\boldsymbol{\beta}\cdot\mathbf{n}v_{b}\rangle_{\partial{\cal T% }_{h}}=0$ . This concludes the proof. ∎

.

Suppose that $u\in H^{k+1}(\Omega)$ . Then for any $v\in V_{h}^{0}$ , we have

$\displaystyle-(\nabla\cdot(\alpha\nabla u),v_{0})=(\alpha\nabla_{w}Q_{h}u,% \nabla_{w}v)_{{\cal T}_{h}}-\ell_{\alpha}(u,v),$ (28)

where

$\displaystyle\ell_{\alpha}(u,v)=\langle(\alpha\nabla u-\mathbb{Q}_{h}(\alpha% \nabla u))\cdot\mathbf{n},v_{0}-v_{b}\rangle_{\partial{{\cal T}_{h}}}+(\mathbb% {Q}_{h}(\alpha\nabla u)-\alpha\nabla_{w}Q_{h}u,\nabla_{w}v)_{{\cal T}_{h}}.$

Proof..

Using integration by parts and the fact that $\sum_{K\in{\cal T}_{h}}\langle\nabla u\cdot\mathbf{n},v_{b}\rangle_{\partial K% }=0$ , we get

$\displaystyle-(\nabla\cdot(\alpha\nabla u),v_{0})=\sum_{K\in{\cal T}_{h}}(% \alpha\nabla u,\nabla v_{0})_{K}-\sum_{K\in{\cal T}_{h}}\langle\alpha\nabla u% \cdot\mathbf{n},v_{0}\rangle_{\partial K}=\sum_{K\in{\cal T}_{h}}(\alpha\nabla u% ,\nabla v_{0})_{K}-\sum_{K\in{\cal T}_{h}}\langle\alpha\nabla u\cdot\mathbf{n}% ,v_{0}-v_{b}\rangle_{\partial K}.$ (29)

By letting $\phi=u$ in Eq. (24) and substituting it into Eq. (29), we get

$\displaystyle-(\nabla\cdot(\alpha\nabla u),v_{0})=\sum_{K\in{\cal T}_{h}}(% \alpha\nabla_{w}(Q_{h}u),\nabla_{w}v)_{K}-\sum_{K\in{\cal T}_{h}}\langle(% \alpha\nabla u-\mathbb{Q}_{h}(\alpha\nabla u))\cdot\mathbf{n},v_{0}-v_{b}% \rangle_{\partial K}+\sum_{K\in{\cal T}_{h}}(\mathbb{Q}_{h}(\alpha\nabla u)-% \alpha\nabla_{w}Q_{h}u,\nabla_{w}v)_{K}.$

This concludes the proof. ∎

Let $\ell_{c}(u,v)=-(cu-cQ_{0}u,v_{0})$ and $e_{h}(t)=Q_{h}u(t)-u_{h}(t)\in V_{h}^{0}$ . Then $e_{h}$ satisfies the following error equation.

.

For a fixed time $t\in(0,T]$ , $\forall v\in V_{h}^{0}$ , we have

$\displaystyle((e_{0,h})_{t},v_{0})+A(e_{h},v)_{{\cal T}_{h}}=\ell_{\alpha}(u,v% )+\ell_{\boldsymbol{\beta}}(u,v)+\ell_{c}(u,v).$ (30)

Proof..

By testing Eq. (1) against $v_{0}$ , we obtain

$\displaystyle(u_{t},v_{0})-(\nabla\cdot(\alpha\nabla u),v_{0})+(\boldsymbol{% \beta}\cdot\nabla u,v_{0})+(cu,v_{0})=(f,v_{0}).$

Using Lemmas 11, 12 and property of the projector $Q_{0}$ , we have

$\displaystyle(Q_{0}u_{t},v_{0})+(\alpha\nabla_{w}Q_{h}u,\nabla_{w}v)_{{\cal T}% _{h}}+(\boldsymbol{\beta}\cdot\nabla_{w}Q_{h}u,v_{0})_{{\cal T}_{h}}+(cQ_{0}u,% v_{0})=(f,v_{0})+\ell_{\alpha}(u,v)+\ell_{\boldsymbol{\beta}}(u,v)+\ell_{c}(u,% v).$

From the bilinear form Eq. (9), we get

$\displaystyle(Q_{0}u_{t},v_{0})+A(Q_{h}u,v)=(f,v_{0})+\ell_{\alpha}(u,v)+\ell_% {\boldsymbol{\beta}}(u,v)+\ell_{c}(u,v).$ (31)

Subtracting Eq. (10) from Eq. (31) generates Eq. (30), which completes the proof. ∎

.

Let $u\in L^{2}(0,T;H^{k+1}(\Omega))$ . If $\alpha$ is a piecewise constant matrix, then for any $v\in V_{h}^{0}$ and $0<t<T$ , the following estimates hold

$\displaystyle|\ell_{\alpha}(u,v)|\leqslant Ch^{k}|u|_{k+1}{|\hskip-1.4454pt|% \hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.4454pt|},$ (32) $\displaystyle|\ell_{\boldsymbol{\beta}}(u,v)|\leqslant Ch^{k}|u|_{k+1}{|\hskip% -1.4454pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.4454pt|},$ (33) $\displaystyle|\ell_{c}(u,v)|\leqslant Ch^{k}|u|_{k+1}{|\hskip-1.4454pt|\hskip-% 1.4454pt|}{v}{|\hskip-1.4454pt|\hskip-1.4454pt|}.$ (34)

Proof..

For the first estimate Eq. (32), it follows from Cauchy-Schwarz inequality and Lemma 9, that

$\displaystyle|\ell_{\alpha}(u,v)|\leqslant\sum_{K\in{\cal T}_{h}}|\langle(% \alpha\nabla u-\mathbb{Q}_{h}(\alpha\nabla u))\cdot\mathbf{n},v_{0}-v_{b}% \rangle_{\partial K}|+\sum_{K\in{\cal T}_{h}}|(\mathbb{Q}_{h}(\alpha\nabla u)-% \alpha\nabla_{w}Q_{h}u,\nabla_{w}v)_{K}|\leqslant C\sum_{K\in{\cal T}_{h}}\|% \nabla u-\mathbb{Q}_{h}\nabla u)\|_{\partial K}\|v_{0}-v_{b}\|_{\partial K}+Ch% ^{k}|u|_{k+1}{|\hskip-1.4454pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.44% 54pt|}\leqslant C\left(\sum_{K\in{\cal T}_{h}}h_{K}\|\nabla u-\mathbb{Q}_{h}% \nabla u\|_{\partial K}^{2}\right)^{1/2}\left(\sum_{K\in{\cal T}_{h}}h_{K}^{-1% }\|v_{0}-v_{b}\|_{\partial K}^{2}\right)^{1/2}+Ch^{k}|u|_{k+1}{|\hskip-1.4454% pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.4454pt|}.$

By using the trace inequality Eq. (19) and Lemma 8, we obtain

$\displaystyle\leqslant Ch^{k}|u|_{k+1}\left[\left(\sum_{K\in{\cal T}_{h}}h_{K}% ^{-1}\|v_{0}-v_{b}\|_{\partial K}^{2}\right)^{1/2}+{|\hskip-1.4454pt|\hskip-1.% 4454pt|}v{|\hskip-1.4454pt|\hskip-1.4454pt|}\right]\leqslant Ch^{k}|u|_{k+1}{|% \hskip-1.4454pt|\hskip-1.4454pt|}v{|\hskip-1.4454pt|\hskip-1.4454pt|}.$

The estimate for $\ell_{\boldsymbol{\beta}}(u,v)$ can be obtained from Lemma 5 in [4]. The last estimate Eq. (34) follows from the Cauchy-Schwarz inequality, and the Lemma 8

$\displaystyle|\ell_{c}(u,v)|=|(cu-cQ_{0}u,v_{0})|\leqslant c_{M}|(Q_{0}u-u,v_{% 0})|\leqslant Ch^{k}|u|_{k+1}|\!|\!|v|\!|\!|,$

which completes the proof. ∎

Now, we will provide an error estimate for the semi discretized SFSWG scheme Eq. (10).

.

Let $u_{h}=\{u_{0},u_{b}\}$ be the solution to the Eq. (10) and $u\in L^{2}(0,T;H^{k+1}(\Omega)\cap H_{0}^{1}(\Omega))$ be the solution to the Eqs (1)–(3). Assume that $e_{h}(t)=Q_{h}u(t)-u_{h}(t)\in V_{h}^{0}$ , $t\in(0,T]$ . If $\alpha$ is a piecewise constant matrix, then the following estimate holds

$\displaystyle\|e_{0,h}(t)\|^{2}+\int_{0}^{t}\gamma{|\hskip-1.4454pt|\hskip-1.4% 454pt|}e_{h}(s){|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}ds\leqslant Ch^{2k}\int_% {0}^{t}|u(s)|^{2}_{k+1}ds.$ (35)

Proof..

Taking $v=e_{h}$ in Eq. (30), we find

$\displaystyle((e_{0,h})_{t},e_{0,h})+A(e_{h},e_{h})_{{\cal T}_{h}}=\ell_{% \alpha}(u,e_{h})+\ell_{\boldsymbol{\beta}}(u,e_{h})+\ell_{c}(u,e_{h})\forall t% \in(0,T].$

From Eq. (14) and Lemma 14, we have

$\displaystyle\frac{1}{2}\frac{d}{dt}\|e_{0,h}\|^{2}+\gamma{|\hskip-1.4454pt|% \hskip-1.4454pt|}e_{h}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}\leqslant Ch^{k}|% u|_{k+1}{|\hskip-1.4454pt|\hskip-1.4454pt|}e_{h}{|\hskip-1.4454pt|\hskip-1.445% 4pt|}\leqslant Ch^{2k}|u|_{k+1}^{2}+\frac{1}{2}\gamma{|\hskip-1.4454pt|\hskip-% 1.4454pt|}e_{h}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}.$

Then it follows that

$\displaystyle\frac{d}{dt}\|e_{0,h}\|^{2}+\gamma{|\hskip-1.4454pt|\hskip-1.4454% pt|}e_{h}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}\leqslant Ch^{2k}|u|_{k+1}^{2}.$ (36)

Integrating Eq. (36) with respect to $t$ and noticing that $e_{h}(0)=0$ , yields Eq. (35). ∎

3. Fully discretized SFSWG scheme

Let $N$ be a positive integer and $\tau_{h}$ be the time-step size, we define $t_{n}=N\tau_{h}$ , $n=0,1,\ldots,N$ to get a fully discretization of Eq. (10). With regard to the temporal direction, we define the numerical derivative as follows:

$\displaystyle\partial_{t}v^{n}=(v^{n}-v^{n-1})/\tau_{h},\tau_{h}>0,$ (37)

at the $n$ th time step with step size $\tau_{h}$ by using the backward Euler difference operator, where $v^{n}=v(t_{n})$ . Suppose that $u_{h}^{n}:=u_{h}(t_{n})=\{u_{0,h}^{n},u_{b,h}^{n}\}$ and using the backward Euler difference quotient defined above to the semi discrtized Eq. (10), we have the fully discretized SFSWG finite element scheme for the Eqs (1)–(3): seek $u_{h}^{n}\in V_{h}^{0}$ such that

$\displaystyle\left\{\begin{array}[]{l}(\partial_{t}u_{0,h}^{n},v_{0})+A(u^{n}_% {h},v)_{{\cal T}_{h}}=(f^{n},v_{0}),\forall v\in V_{h}^{0},n=0,1,\ldots,N,\\ u_{h}^{0}=Q_{h}u^{0},\\ \end{array}\right.$ (38)

or, we can rewrite it as

$\displaystyle\left\{\begin{array}[]{l}(u_{0,h}^{n},v_{0})+\tau_{h}A(u^{n}_{h},% v)_{{\cal T}_{h}}=(u_{0,h}^{n-1},v_{0})+\tau_{h}(f^{n},v_{0}),\forall v\in V_{% h}^{0},\\ u_{h}^{0}=Q_{h}u^{0}.\\ \end{array}\right.$ (39)

Next, we prove an error estimate $Q_{h}u^{n}-u_{h}^{n}$ for any fixed $n$ . Let $e_{h}^{n}=Q_{h}u^{n}-u_{h}^{n}$ be an error, which will be utilized in the following error analysis.

.

For $v=\{v_{0},v_{b}\}\in V_{h}^{0}$ , we have

$\displaystyle(e_{0,h}^{n},v)+\tau_{h}A(e^{n}_{h},v)_{{\cal T}_{h}}=(e_{0,h}^{n% -1},v)-\tau_{h}(r^{n},v_{0})+\tau_{h}(\ell_{\alpha}(u^{n},v)+\ell_{\boldsymbol% {\beta}}(u^{n},v)+\ell_{c}(u^{n},v)),$ (40)

where

$\displaystyle r^{n}=u^{n}_{t}-\partial_{t}u^{n}=u^{n}_{t}-\frac{u^{n}-u^{n-1}}% {\tau_{h}}=-\int_{t_{n-1}}^{t_{n}}\frac{(t-t_{n-1})}{\tau_{h}}u_{tt}dt.$

Proof..

From Eq. (1), we have

$\displaystyle u^{n}_{t}-\nabla\cdot(\alpha\nabla u^{n})+\boldsymbol{\beta}% \cdot\nabla u^{n}+cu^{n}=f^{n}.$ (41)

Testing Eq. (41) by $v=\{v_{0},v_{b}\}\in V_{h}^{0}$ , we obtain

$\displaystyle(u^{n}_{t},v_{0})-(\nabla\cdot(\alpha\nabla u^{n}),v_{0})+(% \boldsymbol{\beta}\cdot\nabla u^{n},v_{0})+(cu^{n},v_{0})=(f^{n},v_{0}).$

Using Lemmas 11 and 12, we find

$\displaystyle(u^{n}_{t},v_{0})+(\alpha\nabla_{w}Q_{h}u^{n},\nabla_{w}v)_{{\cal T% }_{h}}+(\boldsymbol{\beta}\cdot\nabla_{w}Q_{h}u^{n},v_{0})_{{\cal T}_{h}}+(cQ_% {0}u^{n},v_{0})=(f^{n},v_{0})+\ell_{\alpha}(u^{n},v)+\ell_{\boldsymbol{\beta}}% (u^{n},v)+\ell_{c}(u^{n},v),$

which can be rewritten as

$\displaystyle\frac{1}{\tau_{h}}(Q_{0}u^{n}-Q_{0}u^{n-1},v_{0})+A(Q_{h}u^{n},v)% _{{\cal T}_{h}}=-(r^{n},v_{0})+(f^{n},v_{0})+\ell_{\alpha}(u^{n},v)+\ell_{% \boldsymbol{\beta}}(u^{n},v)+\ell_{c}(u^{n},v).$ (42)

Subtracting Eq. (39) from Eq. (42) yields Eq. (40). ∎

Now, we will present an error estimate for the fully discretized SFSWG finite element schemes Eq. (38) or Eq. (39).

.

Let $u_{h}\in L^{2}(0,T;V_{h}^{0})$ be the solution to the fully discretized formulation and $u\in L^{2}(0,T;H^{k+1}(\Omega)\cap H_{0}^{1}(\Omega))$ be the solution to the Eqs (1)–(3). For a fixed $0<n<N$ , let $e_{h}^{n}=\linebreak Q_{h}u(t_{n})-u_{h}^{n}\in V_{h}^{0}$ . If $\alpha$ is a piecewise constant matrix, then there exists a constant $C$ such that

for each $n=1,\ldots,N$ , where $|u|_{k+1}=\max_{1\leqslant m\leqslant n}|u(t_{m})|_{k+1}$ .

Proof..

Taking $v=e_{h}^{m}$ in Eq. (40) and by Lemma 3, we find

$\displaystyle\|e_{0,h}^{m}\|^{2}-\|e_{0,h}^{m-1}\|^{2}+2\tau_{h}\gamma{|\hskip% -1.4454pt|\hskip-1.4454pt|}e_{h}^{m}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}% \leqslant 2\tau_{h}(r^{m},e_{0,h}^{m})+2\tau_{h}(\ell_{\alpha}(u^{m},e_{h}^{m}% )+\ell_{\boldsymbol{\beta}}(u^{m},e_{h}^{m})+\ell_{c}(u^{m},e_{h}^{m})).$ (43)

By using Cauchy-Schwarz inequality, we have

$\displaystyle\|r^{m}\|^{2}=\int_{\Omega}\tau_{h}^{-2}\left(\int_{t_{m-1}}^{t_{% m}}(t-t_{m-1})u_{tt}dt\right)^{2}dx\leqslant\frac{1}{3}\int_{\Omega}\tau_{h}% \int_{t_{m-1}}^{t_{m}}u_{tt}^{2}\textit{dtdx}=\frac{\tau_{h}}{3}\int_{t_{m-1}}% ^{t_{m}}\int_{\Omega}u_{tt}^{2}\textit{dxdt}=\frac{\tau_{h}}{3}\int_{t_{m-1}}^% {t_{m}}\|u_{tt}\|^{2}dt.$ (44)

To bound the first term on the right-hand side of Eq. (43), we have

$\displaystyle 2\tau_{h}(r^{m},e_{0,h}^{m})\leqslant 2\tau_{h}\|r^{m}\|\|e_{0,h% }^{m}\|\leqslant C\tau_{h}\|r^{m}\|{|\hskip-1.4454pt|\hskip-1.4454pt|}e_{h}^{m% }{|\hskip-1.4454pt|\hskip-1.4454pt|}\leqslant C\frac{\tau_{h}}{\gamma}\|r^{m}% \|^{2}+\frac{1}{2}\tau_{h}\gamma{|\hskip-1.4454pt|\hskip-1.4454pt|}e_{h}^{m}{|% \hskip-1.4454pt|\hskip-1.4454pt|}^{2}\leqslant C\tau_{h}^{2}\int_{t_{m-1}}^{t_% {m}}\|u_{tt}\|^{2}dt+\frac{1}{2}\tau_{h}\gamma{|\hskip-1.4454pt|\hskip-1.4454% pt|}e_{h}^{m}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}.$ (45)

The second term on the right-hand side of Eq. (43) can be estimated by Lemma 14 and Young’s inequality as follows:

$\displaystyle 2\tau_{h}(\ell_{\alpha}(u^{m},e_{h}^{m})+\ell_{\boldsymbol{\beta% }}(u^{m},e_{h}^{m})+\ell_{c}(u^{m},e_{h}^{m}))\leqslant 2\tau_{h}Ch^{k}_{K}|u^% {m}|_{k+1}{|\hskip-1.4454pt|\hskip-1.4454pt|}e_{h}^{m}{|\hskip-1.4454pt|\hskip% -1.4454pt|}\leqslant C\frac{2\tau_{h}}{\gamma}h^{2k}_{K}|u^{m}|_{k+1}^{2}+% \frac{1}{2}\tau_{h}\gamma{|\hskip-1.4454pt|\hskip-1.4454pt|}e_{h}^{m}{|\hskip-% 1.4454pt|\hskip-1.4454pt|}^{2}.$ (46)

Combining Eqs (43), (45), and (46) yields

$\displaystyle\|e_{0,h}^{m}\|^{2}-\|e_{0,h}^{m-1}\|^{2}+\tau_{h}\gamma{|\hskip-% 1.4454pt|\hskip-1.4454pt|}e_{h}^{m}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}% \leqslant C\tau_{h}^{2}\int_{t_{m-1}}^{t_{m}}\|u_{tt}\|^{2}dt+C\tau_{h}h^{2k}_% {K}|u^{m}|_{k+1}^{2}.$ (47)

Since $e_{h}^{0}=0$ , for any fixed index $1\leqslant n\leqslant N$ , by summing up the inequalities in Eq. (47) from $m=1$ to $m=n$ , we find

$\displaystyle\|e_{0,h}^{n}\|^{2}+\tau_{h}\gamma\sum_{m=1}^{n}{|\hskip-1.4454pt% |\hskip-1.4454pt|}e_{h}^{m}{|\hskip-1.4454pt|\hskip-1.4454pt|}^{2}\leqslant C% \tau_{h}^{2}\int_{0}^{t_{n}}\|u_{tt}\|^{2}dt+C\sum_{m=1}^{n}\tau_{h}h^{2k}_{K}% |u|_{k+1}^{2}\leqslant C\tau_{h}^{2}\int_{0}^{t_{n}}\|u_{tt}\|^{2}dt+Ch^{2k}_{% K}|u|_{k+1}^{2},$

where $|u|_{k+1}=\max_{1\leqslant m\leqslant n}|u(t_{m})|_{k+1}$ . ∎

4. Numerical experiments

In this section, we will present the results of some numerical experiments to verify the theoretical results derived in previous sections and illustrate the efficiency of the SFSWG method by comparing it with the standard WG method on a square domain and an L-shaped domain with uniform triangular partitions. We will start with examples on square domain. The first two levels of meshes are plotted in Fig. 1.

.

In this example, we consider the Eqs (1)–(3) posed on the domain $\Omega\times(0,T]$ with the following data: $\Omega=(0,1)^{2},T=1,\alpha=I_{2},\boldsymbol{\beta}=(1,2)^{T},c=\sin(xy)$ , and $u(x,y,t)=\sin(\pi x)\sin(\pi y)\linebreak\exp(\sin(\pi t))$ . The initial boundary condition $u^{0}$ and the source term $f(x,y,t)$ are computed accordingly. Let $h=\sqrt{2}/N(N=2,4,8,16,32)$ be the mesh size for triangular grids. We use the polynomial pair $(P_{k}(K),P_{k}(e)),k=1,2$ for space Eq. (5) and set $j=k+1$ in the weak gradient Eq. (7) to find the SFSWG solution $u_{h}=\{u_{0},u_{b}\}$ in the spatial direction and backward Euler method in temporal direction. Table 1 lists errors and convergence rates in $|\!|\!|\cdot|\!|\!|$ -norm and $L^{2}$ -norm at the final time.

Table 1
Errors and convergence rates for $P_{k}(K),(k=1,2)$ elements with $[P_{j}(K)]^{2}$ weak gradient and $\tau_{h}=h^{k+1}$

$k$	$N$	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate
		When $\alpha=0.1I_{2}$				When $\alpha=I_{2}$
1	2	4.1772E-01	–	2.9577E-02	–	4.7949E-01	–	3.6141E-02	–
	4	3.5505E-01	0.23	1.6952E-02	0.80	4.0402E-01	0.25	2.0026E-02	0.85
	8	2.0487E-01	0.79	5.4297E-03	1.64	2.31251E-01	0.80	6.3011E-03	1.67
	16	1.0601E-01	0.95	1.4470E-03	1.91	1.1916E-01	0.96	1.6674E-03	1.92
	32	5.3513E-02	0.99	3.6766E-04	1.98	5.9935E-02	0.99	4.2237E-04	1.99
2	2	1.7926E-01	–	9.4441E-03	–	1.9825E-01	–	9.2873E-03	–
	4	6.2453E-02	1.52	1.6318E-03	2.53	6.9071E-02	1.52	1.6718E-03	2.47
	8	1.6175E-02	1.94	2.3153E-04	2.81	1.7801E-02	1.96	2.3902E-04	2.81
	16	4.0524E-03	2.00	3.0664E-04	2.92	4.4463E-03	2.00	3.1718E-05	2.91

Figure 1.

A triangulation meshes in computation.

Figure 2 shows the computational time (in seconds) comparison between SFSWG finite element method and weak Galerkin finite element method. As we can see in Fig. 2 that the SFSWG algorithm is running faster than the standard weak Galerkin algorithm. We can also see in Fig. 2 that the computation time with 2048 elements and $\mathcal{O}(\frac{1}{h^{k+1}})$ time steps by using the SFSWG is $3.1661e+03$ , which is much less than $3.4293e+03$ , needed by using the standard weak Galerkin algorithm. Therefore, when a large number of elements are used the computation time becomes a significant factor. The SFSWG method is more efficient in accuracy and computation time.

Table 2

Errors and convergence rates for $P_{k}(K),(k=1,2)$ elements with $[P_{j}(K)]^{2}$ weak gradient

$k$	$N$	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate
		When $\alpha=0.1I_{2}$				When $\alpha=I_{2}$
1	2	1.2535E-02	–	4.4207E-03	–	4.8481E-04	–	8.3798E-05	–
	4	6.2159E-03	1.01	1.9982E-03	1.15	2.3808E-04	1.03	2.1742E-05	1.94
	8	2.3337E-03	1.41	5.7518E-04	1.80	1.1645E-04	1.03	5.1005E-06	2.09
	16	1.0131E-03	1.20	1.5097E-04	1.93	5.8682E-05	0.99	1.2892E-06	1.98
	32	4.8384E-04	1.07	3.8295E-05	1.98	2.9471E-05	0.99	3.2471E-07	1.99
2	2	8.7835E-03	–	3.7630E-03	–	2.8773E-04	–	6.2074E-05	–
	4	1.4937E-03	2.56	6.3891E-04	2.56	3.3212E-05	3.11	7.1133E-06	3.13
	8	2.1831E-40	2.77	8.2225E-05	2.96	4.3021E-06	2.95	8.9704E-07	2.98
	16	3.7969E-05	2.52	1.0349E-05	2.99	5.9435E-07	2.85	1.1261E-07	2.99

Figure 2.

Comparison of computation times (in seconds) by using the $(P_{1}(k),P_{1}(e),[P_{2}(K)]^{2})$ elements on a uniform grids.

.

Rotating Gaussian pulse. This example is adopted from [4]. Let the Eqs (1)–(3) be posed on the domain $\Omega\times(0,T]$ with the following data: $\Omega=(-0.5,0.5)\times(-0.5,0.5),T=\pi/4,\boldsymbol{\beta}=(-4y,4x)$ , and $c=0$ . The initial boundary condition $u^{0}$ is given by

$\displaystyle u(x,y,0)=\exp\left(-\frac{(\bar{x}-x_{c})^{2}+(\bar{y}-y_{c})^{2% }}{2\sigma^{2}}\right).$

and the source term $f(x,y,t)=0$ . We set $x_{c}=-0.2,y_{c}=0$ , and $\sigma=0.1$ in the numerical test. The exact solution is given by

$\displaystyle u(x,y,t)=\frac{2\sigma^{2}}{2\sigma^{2}+4\alpha t}\exp\left(-% \frac{(\bar{x}-x_{c})^{2}+(\bar{y}-y_{c})^{2}}{2\sigma^{2}+4\alpha t}\right),$

where $\bar{x}=x\cos(4t)+y\sin(4t)$ and $\bar{y}=-x\sin(4t)+y\cos(4t)$ . The results obtained in Table 2 shows the errors of the SFSWG schemes and the convergence rates in the $L^{2}$ norm and ${|\hskip-1.4454pt|\hskip-1.4454pt|}\cdot{|\hskip-1.4454pt|\hskip-1.4454pt|}$ norm at final time for $\alpha=I_{2}$ and $\alpha=0.1I_{2}$ . The numerical calculations are done on triangular meshes with $k=1,2$ and time step size $\tau_{h}=h^{k+1}$ . The results show that the SFWG scheme with $P_{k}$ elements has convergence rate of $\mathcal{O}(h^{k+1})$ and $\mathcal{O}(h^{k})$ in $L^{2}$ -norm and $H^{1}$ norm, respectively. Let $h=\sqrt{2}/N(N=2,4,8,16,32)$ be the mesh size for triangular meshes.

Table 3

Error analysis and convergence rates for $P_{k}(K),(k=1,2)$ elements with $[P_{k+1}(K)]^{2}$ weak gradient and $\tau_{h}=h^{k+1}$

$k$	$N$	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate
		When $\alpha=0.1I_{2}$				When $\alpha=0.01I_{2}$
1	2	1.7071E-01	–	3.6779E-02	–	1.5094E-01	–	7.3556E-02	–
	4	1.0931E-01	0.64	1.0977E-02	1.74	5.8769E-02	1.36	1.8390E-02	1.99
	8	5.7971E-02	0.92	2.8701E-03	1.94	2.2973E-02	1.36	4.2138E-30	2.12
	16	2.9409E-02	0.98	7.2506E-04	1.98	9.9614E-03	1.20	9.8442E-04	2.01
2	2	3.5322E-02	–	1.6191E-02	–	7.5893E-02	–	3.8393E-02	–
	4	7.6857E-03	2.20	2.6338E-03	2.62	7.9017E-03	3.26	5.0591E-03	2.92
	8	1.6953E-03	2.18	3.4421E-04	2.94	9.7336E-04	3.02	6.3962E-04	2.98
	16	4.0601E-04	2.06	4.3631E-05	2.98	1.6125E-04	2.59	8.0811E-05	2.98

Table 4

Error analysis and convergence rates for the SFSWG scheme Eqs (10) and (38) on L-shaped domain $\Omega=[-1,1]^{2}\setminus(0,1)\times(-1,0)$

$k$	$N$	${\|\hskip-1.4454pt\|\hskip-1.4454pt\|}e_{h}{\|\hskip-1.4454pt\|\hskip-1.4454pt\|}$	Rate	$\\|e_{0}^{n}\\|$	Rate
1	2	2.6938E-01	–	1.5937E-02	–
	4	1.6357E-01	0.72	4.9222E-03	1.70
	8	8.5845E-02	0.93	1.2982E-03	1.92
	16	4.3447E-02	0.98	3.2896E-04	1.98
2	2	3.3114E-02	–	1.1655E-03	–
	4	3.3114E-02	1.85	1.5756E-04	2.89
	8	2.3323E-03	1.98	1.9805E-05	2.99
	16	5.8422E-04	2.00	2.4722E-06	3.00

Figure 3.

A triangulation of a L-shape domain.

.

L-shaped domain. In this example, we use a L-shaped domain $\Omega=[-1,1]^{2}\setminus(0,1)\times(-1,0)$ with the following data: $T=1,\boldsymbol{\beta}=(1,1)^{T},c=1$ , and the exact solution is

$\displaystyle u(x,y,t)=x(1-x)y(1-y)\exp(t(1-t)).$

The results reported in Table 3 shows the errors and the numerical convergence rates in the $L^{2}$ norm and ${|\hskip-1.4454pt|\hskip-1.4454pt|}\cdot{|\hskip-1.4454pt|\hskip-1.4454pt|}$ norm. The numerical calculations are done on basis of degrees $k=1,2$ and the polynomial degree $j=k+1$ for the weak gradient Eq. (7). The results show that for both $\alpha=0.1I_{2}$ and $\alpha=0.01I_{2}$ , the SFWG scheme with $P_{k}$ elements has convergence rate of $\mathcal{O}(h^{k+1})$ and $\mathcal{O}(h^{k})$ in $L^{2}$ -norm and $H^{1}$ norm, respectively. The first two levels of L-shaped domain are shown in Fig. 3.

.

As the final example, we use a L-shaped domain $\Omega=[-1,1]^{2}\setminus(0,1)\times(-1,0)$ with the following data: $T=1,\alpha=\begin{bmatrix}\par 2&0\\ 0&1\end{bmatrix},\boldsymbol{\beta}=(1,1)^{T},c=1$ , and the exact solution is

$\displaystyle u(x,y,t)=\exp(-t)x(1-x)y(1-y).$ (48)

Table 4 shows the performance of the SFSWG schemes for the Eq. (48) on polynomial of degrees $k=1,2$ and time step size $\tau_{h}=h^{k+1}$ . The results indicate that the SFSWG method with $P_{k}$ elements has convergence rate of $\mathcal{O}(h^{k+1})$ and $\mathcal{O}(h^{k})$ in $L^{2}$ -norm and $H^{1}$ norm, respectively. The numerical solution, the exact solution and their error are shown in Fig. 4.

Figure 4.

WG solution (Left), exact solutions (Middle) and the error (Right) at final time for $(P_{1}(K),P_{1}(e),[P_{2}(K)]^{2})$ element on L-shaped domain with $\tau_{h}=h^{2}$ .

5. Conclusion

A fully discretized scheme using SFSWG finite element method with backward Euler difference operator Eq. (37) in temporal direction is proposed in this paper. The error estimates and convergence of semi discretized and fully discretized SFSWG schemes are derived. Numerical results show that while the same rate of convergence can be obtained using standard weak Galerkin algorithm and SFSWG algorithm Eqs (10) and (38), the SFSWG algorithm is more efficient and easier to implement for time dependent convection-diffusion equation.

References

Al-Taweel

Hussain

Wang

. A stabilizer free weak Galerkin finite element method for parabolic equation. Journal of Computational and Applied Mathematics. 2021; 392: 113373.

Kim

Shin

. A high order discontinuous Galerkin method with skeletal multipliers for convection-diffusion-reaction problems. Computer Methods in Applied Mechanics and Engineering. 2019; 343: 207–233.

Stynes

. Finite volume methods for convection-diffusion problems. Journal of Computational and Applied Mathematics. 1995; 63(1-3): 83–90.

Xie

Zhu

Wang

. Error analysis of weak Galerkin finite element methods for time-dependent convection-diffusion equations. Applied Numerical Mathematics. 2019; 137: 19–33.

Hussain

Malluwawadu

Zhu

. A weak Galerkin finite element method for the second order elliptic problem with mixed boundary condition. J Appl Anal Comput. 2018; 8(5): 1452–1463.

Wang

. A computational study of the weak Galerkin method for second-order elliptic equations. Numerical Algorithms. 2013; 63(4): 753–777.

Malluwawadu

Hussain

. A weak Galerkin method for second order elliptic problems with polynomial reduction. Journal of Applied Analysis & Computation. 2019; 9(2): 655–670.

Cheichan

Kashkool

Gao

. A weak galerkin finite element method for solving nonlinear convection-diffusion problems in one dimension. International Journal of Applied and Computational Mathematics. 2019; 5(2): 1–15.

Cheichan

Kashkool

Gao

. A weak Galerkin finite element method for solving nonlinear convection-diffusion problems in two dimensions. Applied Mathematics and Computation. 2019; 354: 149–163.

10.

Al-Taweel

Hussain

Wang

Jones

. A P0–P0 weak Galerkin finite element method for solving singularly perturbed reaction-diffusion problems. Numerical Methods for Partial Differential Equations. 2020; 36(2): 213–227.

11.

Wang

. Weak Galerkin finite element methods for the biharmonic equation on polytopal meshes. Numerical Methods for Partial Differential Equations. 2014; 30(3): 1003–1029.

12.

Wang

. weak Galerkin finite element method for the Helmholtz equation with large wave number on polytopal meshes. Numerical Methods for Partial Differential Equations. 2015; 35(3): 1228–1255.

13.

Hussein

Kashkool

. Weak Galerkin finite element method for solving one-dimensional coupled Burgers’ equations. Journal of Applied Mathematics and Computing. 2020; 63(1): 265–293.

14.

Zhang

. A weak Galerkin finite element method for p-laplacian problem. East Asian Journal On Applied Mathematics. 2021; 11(2): 219–233.

15.

Wang

Zhang

. A weak Galerkin finite element method for the Maxwell equations. Journal of Scientific Computing. 2015; 65(1): 363–386.

16.

Zhang

. A stabilizer-free weak Galerkin finite element method on polytopal meshes. Journal of Computational and Applied Mathematics. 2020; 371: 112699.

17.

Al-Taweel

Wang

. A note on the optimal degree of the weak gradient of the stabilizer free weak Galerkin finite element method. Applied Numerical Mathematics. 2020; 150: 444–451.

18.

Di Pietro

Ern

. Mathematical aspects of discontinuous Galerkin methods. vol. 69. Springer Science & Business Media; 2011.

19.

Wang

. A weak Galerkin mixed finite element method for second order elliptic problems. Mathematics of Computation. 2014; 83(289): 2101–2126.

20.

Wang

. A weak Galerkin finite element method for second-order elliptic problems. Journal of Computational and Applied Mathematics. 2013; 241: 103–115.

A stabilizer free spatial weak Galerkin finite element methods for time-dependent convection-diffusion equations

Abstract

Keywords

1. Introduction

.

.

.

.

.

.

.

Proof..

.

Proof..

.

.

Proof..

.

.

Proof..

.

.

.

.

Proof..

.

Proof..

.

Proof..

.

Proof..

.

Proof..

.

Proof..

.

Proof..

.

Proof..

.

Proof..

.

Table 1 Errors and convergence rates for P k ( K ) , ( k = 1 , 2 ) elements with [ P j ⁢ ( K ) ] 2 weak gradient and τ h = h k + 1

.

.

.

References

Table 1
Errors and convergence rates for $P_{k}(K),(k=1,2)$ elements with $[P_{j}(K)]^{2}$ weak gradient and $\tau_{h}=h^{k+1}$