Space-time finite element methods stabilized using bubble function spaces

ABSTRACT

In this paper, a stabilized space-time finite element method for solving linear parabolic evolution problems is analyzed. The proposed method is developed on a base of a space-time variational setting, that helps on the simultaneous and unified discretization in space and in time by finite element techniques. Stabilization terms are constructed by means of classical bubble spaces. Stability of the discrete problem with respect to an associated mesh dependent norm is proved, and a priori discretization error estimates are presented. Numerical examples confirm the theoretical estimates.

KEYWORDS:

• Parabolic initial-boundary value problems
• space-time finite element methods
• bubble stabilization
• optimal convergence rates

1. Introduction

Parabolic evolution equations are used to describe numerous physical phenomena, as for example heat transfer. The traditional methods for parabolic problems usually apply a separate method for the time discretization, e.g. implicit Runge–Kutta methods. During the last decades, efficient discontinuous Galerkin finite element methods have been presented for the time discretization of parabolic problems, see, e.g. an analysis for Galerkin time-stepping methods in [1–3], we also refer to the monograph [4]. Adaptive algorithms based on a posteriori error estimates have also been presented and successfully tested for linear and nonlinear problems, see e.g. [5,6] and the references therein. In [7,8], space-time adaptive wavelet methods for parabolic evolution problems have been studied. Also in the literature, p and hp finite element methods for parabolic problems have been presented, see [9,10].

Another approach that has been followed is the derivation of space-time finite element methods, based on appropriate space-time variational settings. The basic idea is to consider the time variable t as just another variable, lets say $x_{d+1}$, if we consider that $x=(x_1,\dots ,x_d)$ are the spatial variables. In that way, the time derivative, which appears in the parabolic PDE model, plays the role of a convection term in the time direction $x_{d+1}$. By multiplying the given parabolic problem by a space-time test function and applying integration by parts, we can derive the weak space-time formulation. The derived weak formulation helps on the unified space-time discretization by finite element techniques. This means that we discretize the problem in space and in time by using a common finite element space. In this spirit, in [11], space-time finite element methods have been developed for elastodynamics. In particular, the method uses discontinuous Galerkin techniques for the time discretization and incorporates Petrov–Galerkin techniques, see [12], to ensure stability. Stream-line diffusion techniques that are presented in [12], have been also used for developing space-time finite element methods for conservation laws and fluid flow problems, see e.g. [13,14] and the references there. In [15,16], the stability of Petrov–Galerkin discretizations of parabolic problems have been studied and stable space-time trial and test functions have been proposed. In [17], conforming space-time finite element approximations to parabolic problems have been investigated. In [18], upwind-stabilized single-patch space-time Isogeometric analysis (IgA) schemes for parabolic evolution problems are proposed. Recently in [19,20], the authors based on [18], analyzed a time discontinuous Galerkin multipatch IgA scheme and demonstrated the efficiency of a space-time solver implemented on a parallel environment.

In this paper, we focus on the model problem ${\mathrm{\partial }}_{t}u-\kappa \mathrm{\Delta }u=f$, with zero initial and boundary conditions, and the diffusivity parameter κ is positive and constant. Following standard procedures, e.g. Nitsche techniques for imposing weakly boundary conditions, see [21], the current analysis can be extended to more general initial boundary problems. We propose a new space-time finite element method, which is stabilized by introducing classical bubble spaces, see [22–24]. The bubble basis functions vanish on the edges of the mesh elements and, in addition, do not affect the continuity properties of the solution. By enriching in that way the initial finite element space, the numerical solution consists of two components, where the first lives in the initial finite element space, and the second lives in the bubble space. For developing our analysis, we are motivated and inspired by the subgrid scale stabilization techniques presented in [25,26], for solving linear first-order problems. There, the idea is to couple the initial finite element space with new subgrid scale spaces and to construct artificial diffusion terms in these new spaces. The artificial diffusion terms are added in the numerical scheme in order to ensure stability. The innovation in our approach is that, instead of using subgrid spaces on different meshes, we use bubble spaces and the artificial diffusion terms are formed in these spaces. We include a positive parameter θ in the corresponding bubble diffusion terms, that can control the magnitude of the artificial diffusion in the time direction. We prove stability of the discrete problem with respect to the produced norm, which is a mesh depended norm. Also, optimal error estimates for the full numerical solution containing the bubble component are shown. The latter is achieved by choosing the value of θ to be close to the mesh size but independent of κ. During the discretization error analysis, we analytically present the dependence of the constants on to the diffusivity parameter κ. In the end, this helps to have a clear idea, about the form of the constants that appear in the error bounds for the difference between the solution u and the discrete solution $u_h$. In Section 5, we perform tests for different values of κ. To author's knowledge, it is the first time that this type of stabilized space-time finite element methods are presented and analyzed.

The paper is structured as follows. In Section 2, the model parabolic problem is presented. In Section 3, we formulate the stabilized space-time finite element scheme. In Section 4, we present the error analysis and derive the error estimates. We discuss numerical examples in Section 5. The paper closes with the conclusions.

2. The model problem

2.1. Preliminaries

Let Ω be a bounded Lipschitz domain in ${\mathbb{R}}^d,d\ge 1$. Let $\mathit{\alpha }=({\alpha }_1,\dots ,{\alpha }_d)$ be a multi-index of non-negative integers ${\alpha }_1,\dots ,{\alpha }_d$ with degree $|\mathit{\alpha }|=\sum _{j=1}^d{\alpha }_j$. For any $\mathit{\alpha }$, we define the differential operator $D_x^{\mathit{\alpha }}=D_{x_1}^{{\alpha }_1}\dots D_{x_d}^{{\alpha }_d}$, with $D_{x_j}=\mathrm{\partial }/\mathrm{\partial }{x}_j$, $j=1,\dots ,d$. As usual, $L^2\left(\mathrm{\Omega }\right)$ denotes the Sobolev space for which ${\int }_{\mathrm{\Omega }}|\phi \left(x\right){|}^2\mathrm{d}x<\mathrm{\infty }$, endowed with the norm $\parallel \phi {\parallel }_{L^2\left(\mathrm{\Omega }\right)}=\left({\int }_{\mathrm{\Omega }}|\phi \right(x){|}^2\mathrm{d}x{)}^{1/2}$. Let ℓ be a non-negative integer, define $H^{\ell }\left(\mathrm{\Omega }\right)=\{\phi \in L^2(\mathrm{\Omega }):D_x^{\mathit{\alpha }}\phi \in L^2(\mathrm{\Omega }),\mathrm{for}\mathrm{all}|\mathit{\alpha }|\le \ell \},$ the standard Sobolev spaces endowed with the following norms $\parallel \phi {\parallel }_{H^{\ell }\left(\mathrm{\Omega }\right)}=(\sum _{0\le |\mathit{\alpha }|\le \ell }\parallel D_x^{\mathit{\alpha }}\phi {\parallel }_{L^2\left(\mathrm{\Omega }\right)}^2{)}^{1/2},$ and seminorms $|\phi {|}_{H^{\ell }\left(\mathrm{\Omega }\right)}=(\sum _{|\mathit{\alpha }|=\ell }\parallel D_x^{\mathit{\alpha }}\phi {\parallel }_{L^2\left(\mathrm{\Omega }\right)}^2{)}^{1/2}$. Also we define the subspace $H_0^1\left(\mathrm{\Omega }\right)$ of $H^1\left(\mathrm{\Omega }\right)$ $H_0^1\left(\mathrm{\Omega }\right)=\{\phi \in H^1(\mathrm{\Omega }):\phi =0\mathrm{on}\mathrm{\partial \Omega }\}.$ Let $I=[0,T]$ with T>0 be the time interval. For later use, we define the space-time cylinder $Q=\mathrm{\Omega }\times (0,T)$ and its boundary parts $\mathrm{\Sigma }=\mathrm{\partial \Omega }\times (0,T)$, ${\mathrm{\Sigma }}_T=\mathrm{\Omega }\times \left\{T\right\}$ and ${\mathrm{\Sigma }}_0=\mathrm{\Omega }\times \left\{0\right\},$ see an illustration Figure 1(a). We denote the gradient by $\mathrm{\nabla }\phi =({\mathrm{\nabla }}_x\phi ,{\mathrm{\partial }}_t\phi )$, where ${\mathrm{\nabla }}_x\phi$ is the gradient with respect to the spatial variables. Similarly, we denote by $n=(n_x,n_t)$ the normal component on $\mathrm{\partial }Q$, with $n_x$ the components related to space direction and $n_t$ the component related to time direction. Let $\ell ,m$ be positive integers, for functions defined in Q, we define the Sobolev spaces (1a) $H^{\ell ,m}\left(Q\right)=\{\phi \in L^2(Q):D_x^{\mathit{\alpha }}\phi \in L^2(Q)\mathrm{with}0\le |\mathit{\alpha }|\le \ell ,\mathrm{and}{\mathrm{\partial }}_t^i\phi \in L^2(Q),i=1,\dots ,m\}$(1a) and the subspaces (1b) $H_0^{1,0}\left(Q\right)=\{\phi \in L^2(Q):{\mathrm{\nabla }}_x\phi \in [L^2\left(Q\right){]}^d,\phi =0\mathrm{on}\mathrm{\Sigma }\},$(1b) (1c) $H_{0,\overline{0}}^{1,1}\left(Q\right)=\{\phi \in L^2(Q):{\mathrm{\nabla }}_x\phi \in [L^2\left(Q\right){]}^d,{\mathrm{\partial }}_t\phi \in L^2\left(Q\right),\phi =0\mathrm{on}\mathrm{\Sigma },\phi =0\mathrm{on}{\mathrm{\Sigma }}_T\},$(1c) (1d) $H_{0,\underset{\_}{0}}^{1,1}\left(Q\right)=\{\phi \in L^2(Q):{\mathrm{\nabla }}_x\phi \in [L^2\left(Q\right){]}^d,{\mathrm{\partial }}_t\phi \in L^2\left(Q\right),\phi =0\mathrm{on}\mathrm{\Sigma },\phi =0\mathrm{on}{\mathrm{\Sigma }}_0\}.$(1d) For a function $\phi \in H^{\ell ,m}\left(Q\right)$ with $\ell ,m\ge 1$, we define the norms and the seminorms (2a) $\parallel \phi {\parallel }_{H^{\ell ,m}\left(Q\right)}:={\left(\sum _{|\mathit{\alpha }|\le \ell }\parallel D^{\mathit{\alpha }}\phi {\parallel }_{L^2\left(Q\right)}^2+\sum _{i=0}^m\parallel {\mathrm{\partial }}_t^i\phi {\parallel }_{L^2\left(Q\right)}^2\right)}^{1/2},$(2a) (2b) $|\phi {|}_{H^{\ell ,m}\left(Q\right)}:={\left(\sum _{|\mathit{\alpha }|=\ell }\parallel D^{\mathit{\alpha }}\phi {\parallel }_{L^2\left(Q\right)}^2+\parallel {\mathrm{\partial }}_t^m\phi {\parallel }_{L^2\left(Q\right)}^2\right)}^{1/2}.$(2b)

Figure 1. (a) The space-time domain Q with the boundary parts and its mesh $T_h\left(Q\right)$. (b) The bubble function on the reference triangular mesh element. (c) The bubble function on the reference rectangular mesh element.

We recall Hölder's and Young's inequalities (3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) that hold for all $u\in L^2\left(Q\right)$ and $v\in L^2\left(Q\right)$ and for any fixed $ϵ\in (0,\mathrm{\infty })$.

We will use the following Poincare's inequality: Let $\mathrm{\Omega }\subset {\mathbb{R}}^d,d=1,2,\dots$ be a bounded rectangular domain and let $\mathrm{\Gamma }\subset \mathrm{\partial \Omega }$ with $|\mathrm{\Gamma }|>0$. For simplicity we assume that Γ lies on the plane with $x_1=0$. Let $v\in C^{\mathrm{\infty }}\left(\mathrm{\Omega }\right)$ and $v\left(x_{\mathrm{\Gamma }}\right)=0$ for all $x_{\mathrm{\Gamma }}\in \mathrm{\Gamma }$. For any interior point $x=(x_1,\dots ,x_d)$, we have (4) $v(x_1,\dots ,x_d)=v\left(x_{\mathrm{\Gamma }}\right)+{\int }_{x_{\mathrm{\Gamma },1}}^{x_1}\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_1}(\tau ,x_2,\dots ,x_d)\mathrm{d}\tau .$(4) The first inequality in (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ) yields (5) $\left(v\right(x_1,\dots ,x_d){)}^2={\left({\int }_{x_{\mathrm{\Gamma },1}}^{x_1}1\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_1}(\tau ,x_2,\dots ,x_d)\mathrm{d}\tau \right)}^2\le C_{\mathrm{\Omega }}{\int }_{x_{\mathrm{\Gamma },1}}^{x_1}{\left|\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_1}(\tau ,x_2,\dots ,x_d)\right|}^2\mathrm{d}\tau ,$(5) where the constant $C_{\mathrm{\Omega }}$ depends on the length of Ω. Integrating (5(5) $\left(v\right(x_1,\dots ,x_d){)}^2={\left({\int }_{x_{\mathrm{\Gamma },1}}^{x_1}1\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_1}(\tau ,x_2,\dots ,x_d)\mathrm{d}\tau \right)}^2\le C_{\mathrm{\Omega }}{\int }_{x_{\mathrm{\Gamma },1}}^{x_1}{\left|\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_1}(\tau ,x_2,\dots ,x_d)\right|}^2\mathrm{d}\tau ,$(5) ) over Ω, we can obtain (6) ${\int }_{\mathrm{\Omega }}{v}^2\left(x\right)\mathrm{d}x\le C_{\mathrm{\Omega }}^2{\int }_{\mathrm{\Omega }}({\mathrm{\partial }}_{x_1}v{)}^2\mathrm{d}x.$(6) We refer to [27] for a proof of (6(6) ${\int }_{\mathrm{\Omega }}{v}^2\left(x\right)\mathrm{d}x\le C_{\mathrm{\Omega }}^2{\int }_{\mathrm{\Omega }}({\mathrm{\partial }}_{x_1}v{)}^2\mathrm{d}x.$(6) ) for more general domains Ω.

In what follows, positive constants c and C appearing in inequalities are generic constants which do not depend on the mesh-size h and on u. In many cases, we will indicate on what may the constants depend for an easier understanding of the proofs. We will write $a\sim b$ meaning that $ca\le b\le Ca$, with $c,C$ generic constants.

2.2. The model parabolic problem

In the space-time cylinder $\overline{Q}=\overline{\mathrm{\Omega }}\times [0,T]$, with $\mathrm{\Omega }\subset {\mathbb{R}}^d,\mathrm{d}=1,2,3$, we consider the initial boundary value problem (7) $\begin{array}{rl}{u}_t-\kappa \mathrm{\Delta }u& =f\mathrm{in}Q\text{and}\\ u& =0\mathrm{on}\mathrm{\Sigma },u(\cdot ,0)=u_0\mathrm{on}{\mathrm{\Sigma }}_0,\end{array}$(7) as a model problem, where the diffusivity parameter $0<\kappa \le 1$ is taken to be constant, $f:Q\to \mathbb{R}$, with $f\in L^2\left(Q\right)$, and $u_0:\mathrm{\Omega }\to \mathbb{R}$, with $u_0\in L^2\left(\mathrm{\Omega }\right)$ are given functions, and $u:\overline{Q}\to \mathbb{R}$ is the unknown. Using the standard procedure and integration by parts with respect to both x and t, we can easily derive the following space-time variational formulation of (7(7) $\begin{array}{rl}{u}_t-\kappa \mathrm{\Delta }u& =f\mathrm{in}Q\text{and}\\ u& =0\mathrm{on}\mathrm{\Sigma },u(\cdot ,0)=u_0\mathrm{on}{\mathrm{\Sigma }}_0,\end{array}$(7) ): find $u\in H_0^{1,0}\left(Q\right)$ such that (8) $\overline{a}(u,v)=\overline{l}\left(v\right),\mathrm{for}\mathrm{all}v\in H_{0,\overline{0}}^{1,1}\left(Q\right)$(8) with (9) $\overline{a}(u,v)=-{\int }_{Q}u(x,t){\mathrm{\partial }}_{t}v(x,t)\mathrm{d}x\mathrm{d}t+\kappa {\int }_Q{\mathrm{\nabla }}_{x}u(x,t)\cdot {\mathrm{\nabla }}_{x}v(x,t)\mathrm{d}x\mathrm{d}t,$(9) (10) $\overline{l}\left(v\right)={\int }_{Q}f(x,t)v(x,t)\mathrm{d}x\mathrm{d}t+{\int }_{\mathrm{\Omega }}{u}_0\left(x\right)v(x,0)\mathrm{d}x.$(10)

The variational problem (8(8) $\overline{a}(u,v)=\overline{l}\left(v\right),\mathrm{for}\mathrm{all}v\in H_{0,\overline{0}}^{1,1}\left(Q\right)$(8) ) is known to have a unique weak solution, see [28] and in [29,30] for a more comprehensive analysis on existence and uniqueness results. For simplicity, we only consider homogeneous Dirichlet boundary conditions on Σ and $u_0=0$. However, the analysis presented in this paper can easily be generalized to other constellations of boundary conditions.

Assumption 2.1

We assume that the solution u of (8(8) $\overline{a}(u,v)=\overline{l}\left(v\right),\mathrm{for}\mathrm{all}v\in H_{0,\overline{0}}^{1,1}\left(Q\right)$(8) ) belongs to $V_{0,\underset{\_}{0}}:=H_{0,\underset{\_}{0}}^{1,1}\left(Q\right)\cap H^s\left(Q\right)$ with some $s\ge 2$.

Note that in Assumption 2.1 we require uniform regularity properties for the solution in x and t directions. In the following Sections we are going to present the discretization error analysis based on Assumption 2.1. In [19,20], error estimates have been given for Isogeometric Analysis methods considering anisotropic regularity properties for the solution u.

3. The discrete problem

3.1. A usual finite element scheme

Let $T_h\left(Q\right)$ be a partition of space-time domain Q into triangular (or quadrilateral elements), that is $\overline{Q}={\cup }_{E\in T_h}E$, see Figure 1(a). We denote by $h_E$ the diameter of $E\in T_h\left(Q\right)$ and the mesh size is defined as $h=\underset{E}{max}\left\{h_E\right\}$. We assume that $T_h\left(Q\right)$ is uniform, i.e. there exists a positive number $C_{un}$, such that ${\rho }_E\ge h/C_{un}$, where ${\rho }_E$ is the diameter of the circle inscribed in $E\in T_h\left(Q\right)$. Associated with $T_h\left(Q\right)$, we define the finite element subspace $V_{h0}$ of $H_{0,\underset{\_}{0}}^{1,1}\left(Q\right)$, consisting of continuous functions in space and in time, by (11) $V_{h0}=\{v_h\in H_{0,\underset{\_}{0}}^{1,1}(Q):v_h{|}_E\in {\mathbb{P}}^p(E),\mathrm{for}\mathrm{every}E\in T_h(Q\left)\right\},$(11) where ${\mathbb{P}}^p\left(E\right)$ is the polynomial space of total degree p, see e.g. [31–33]. In our analysis we focus on the case p=1.

The usual finite element approximation of (7(7) $\begin{array}{rl}{u}_t-\kappa \mathrm{\Delta }u& =f\mathrm{in}Q\text{and}\\ u& =0\mathrm{on}\mathrm{\Sigma },u(\cdot ,0)=u_0\mathrm{on}{\mathrm{\Sigma }}_0,\end{array}$(7) ) would read: find $u_h\in V_{h0}$ such that (12a) $a(u_h,v_h)=(f,v_h)\mathrm{\forall }{v}_h\in V_{h0},$(12a) where (12b) $a(u_h,v_h)={\int }_Q{\mathrm{\partial }}_t{u}_h{v}_h\mathrm{d}x\mathrm{d}t+\kappa {\int }_Q{\mathrm{\nabla }}_x{u}_h\cdot {\mathrm{\nabla }}_x{v}_h\mathrm{d}x\mathrm{d}t.$(12b) Setting $v_h=u_h$ in (12a(12a) $a(u_h,v_h)=(f,v_h)\mathrm{\forall }{v}_h\in V_{h0},$(12a) ) and using the identity $2{\int }_Q\mathrm{\partial }{u}_h{u}_h\mathrm{d}x\mathrm{d}t={\int }_Q{\mathrm{\partial }}_t\left(u_h^2\right)\mathrm{d}x\mathrm{d}t={\int }_{{\mathrm{\Sigma }}_T}{u}_h^2\mathrm{d}s$, we have $\parallel u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2\le \parallel f{\parallel }_{L^2\left(Q\right)}\parallel u_h{\parallel }_{L^2\left(Q\right)}$, and by inequalities (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ) and (6(6) ${\int }_{\mathrm{\Omega }}{v}^2\left(x\right)\mathrm{d}x\le C_{\mathrm{\Omega }}^2{\int }_{\mathrm{\Omega }}({\mathrm{\partial }}_{x_1}v{)}^2\mathrm{d}x.$(6) ), we can obtain the following estimate (13) $\frac{1}{2}\parallel u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\frac{\kappa }{2}\parallel {\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2\le \frac{2C_{\mathrm{\Omega }}^2}{\kappa }\parallel f{\parallel }_{L^2\left(Q\right)}^2,$(13) where $C_{\mathrm{\Omega }}$ is the constant that appears in (6(6) ${\int }_{\mathrm{\Omega }}{v}^2\left(x\right)\mathrm{d}x\le C_{\mathrm{\Omega }}^2{\int }_{\mathrm{\Omega }}({\mathrm{\partial }}_{x_1}v{)}^2\mathrm{d}x.$(6) ). In the stability estimate (13(13) $\frac{1}{2}\parallel u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\frac{\kappa }{2}\parallel {\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2\le \frac{2C_{\mathrm{\Omega }}^2}{\kappa }\parallel f{\parallel }_{L^2\left(Q\right)}^2,$(13) ) the control of ${\mathrm{\nabla }}_x{u}_h$ can be very poor if the parameter κ is very small. Furthermore (13(13) $\frac{1}{2}\parallel u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\frac{\kappa }{2}\parallel {\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2\le \frac{2C_{\mathrm{\Omega }}^2}{\kappa }\parallel f{\parallel }_{L^2\left(Q\right)}^2,$(13) ) does not provide a direct control of $\parallel {\mathrm{\partial }}_t{u}_h{\parallel }_{L^2\left(Q\right)}$. Thus in the case where κ is small, the finite element scheme (12(12a) $a(u_h,v_h)=(f,v_h)\mathrm{\forall }{v}_h\in V_{h0},$(12a) ) may not perform well. Therefore it is crucial to improve the stability properties of (12(12a) $a(u_h,v_h)=(f,v_h)\mathrm{\forall }{v}_h\in V_{h0},$(12a) ). Next we present a technique for stabilizing the finite element scheme (12(12a) $a(u_h,v_h)=(f,v_h)\mathrm{\forall }{v}_h\in V_{h0},$(12a) ) by adding artificial diffusion terms, which do not deteriorate the approximate properties of the method. The idea is to enrich the original finite dimensional space by adding a bubble function space, and then to construct the artificial diffusion terms in this space.

3.2. The stabilized scheme

We introduce the larger finite subspace $V_{h,b}$ of $H_{0,\underset{\_}{0}}^{1,1}\left(Q\right)$ that can be written as a direct sum as follows (14) $V_{h,b}=\{v_h\in H_{0,\underset{\_}{0}}^{1,1}(Q):v_h{|}_E\in {\mathbb{P}}^1(E)\oplus V_B(E),\mathrm{for}\mathrm{every}E\in T_h(Q\left)\right\},$(14) with $V_B\left(E\right):=V_B{|}_E$, where $V_B$ denotes the space of bubble functions, that vanish entirely on the boundary of the mesh elements and have exactly one degree of freedom in each $E\in T_h\left(Q\right)$. For example, in case of triangular elements, it is spanned by a cubic functions $V_B:=\{v^b\in H_0^1(Q):v^b{|}_E=C_E{\lambda }_1{\lambda }_2{\lambda }_3\}$, where ${\lambda }_i,i=1,2,3$ are linear polynomials vanishing on one side of $\mathrm{\partial }E$ and taking the value one at the opposite vertex. The constant $C_E$ is chosen such that $\underset{x\in E}{max}{v}^b\left(x\right)=1$. Thus, every bubble basis function $\phi \in V_B\left(E\right)$ satisfies, (i) $\phi \left(x\right)>0$ for $x\in E$, (ii) $\phi \left(x\right)=0$ for $x\in \mathrm{\partial }E$, and (iii) ${\int }_E{\phi }^2\left(x\right)\mathrm{d}x=C_E{h}_E^2$, with $C_E$ depending on the uniformity of $T_h$ but is independent of $h_E$. An illustration of bubble functions on two-dimensional elements is presented in Figure 1(b,c).

Based on $V_{h,b}$ defined in (14(14) $V_{h,b}=\{v_h\in H_{0,\underset{\_}{0}}^{1,1}(Q):v_h{|}_E\in {\mathbb{P}}^1(E)\oplus V_B(E),\mathrm{for}\mathrm{every}E\in T_h(Q\left)\right\},$(14) ), any $v_h\in V_{h,b}$ can be decomposed into two parts, i.e. $v_h=v_h^1+v_h^b$ with $v_h^1\in V_{h0}$ and $v_h^b\in V_B$. In view of this, we introduce the discrete problem: find $u_h\in V_{h,b}$ such that (15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) where $a(\cdot ,\cdot )$ as in (12b(12b) $a(u_h,v_h)={\int }_Q{\mathrm{\partial }}_t{u}_h{v}_h\mathrm{d}x\mathrm{d}t+\kappa {\int }_Q{\mathrm{\nabla }}_x{u}_h\cdot {\mathrm{\nabla }}_x{v}_h\mathrm{d}x\mathrm{d}t.$(12b) ) and (15b) $b_h(u_h^b,v_h^b)=\theta h{\int }_Q{\mathrm{\partial }}_t{u}_h^b{\mathrm{\partial }}_t{v}_h^b\mathrm{d}x\mathrm{d}t,$(15b) with $\theta >0$ to be a positive constant, which will be determined later. We recall the following inverse estimate and the scaled trace inequality, see proofs in [31] and in [21].

Lemma 3.1.

Let $v\in H^1\left(Q\right),v_h\in V_{h,b}$ and let a mesh element $E\in T_h\left(Q\right)$. Then there exist constants $c_{\mathrm{inv}},c_{\mathrm{trac}}>0$ independent of h such that (16) $\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(E\right)}\le c_{\mathrm{inv}}{h}^{-1}\parallel v_h{\parallel }_{L^2\left(E\right)},$(16) (17) $\parallel v{\parallel }_{L^2\left(\mathrm{\partial }E\right)}\le c_{\mathrm{trac}}{h}^{-1/2}(\parallel v{\parallel }_{L^2\left(E\right)}+h\parallel \mathrm{\nabla }v{\parallel }_{L^2\left(E\right)}).$(17)

For convenience, we introduce the discrete bilinear form (18) $a_h(u_h,v_h)=a(u_h,v_h)+b_h(u_h^b,v_h^b),$(18) and the associated mesh dependent norm (19) $\parallel v_h{\parallel }_h=(\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{v}_h^b{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}.$(19) Note that, the terms with the time derivatives in (19(19) $\parallel v_h{\parallel }_h=(\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{v}_h^b{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}.$(19) ) are related to the bubble function space.

Lemma 3.2.

The discrete form $a_h(\cdot ,\cdot ):V_{h,b}\times V_{h,b}\to \mathbb{R}$ defined in (18(18) $a_h(u_h,v_h)=a(u_h,v_h)+b_h(u_h^b,v_h^b),$(18) ), is $V_{h,b}$-coercive with respect to the norm $\parallel \cdot {\parallel }_h,$ i.e., (20) $a_h(v_h,v_h)\ge C_s\parallel v_h{\parallel }_h^2,\mathrm{\forall }{v}_h\in V_{h,b}.$(20)

Proof.

Let $v_h\in V_{h,b}$. Since $v_h(x,0)=0$ and $n_t{|}_{\mathrm{\Sigma }}=0$, it follows by Green's formula ${\int }_Q{\mathrm{\partial }}_t{v}_h{v}_h+v_h{\mathrm{\partial }}_t{v}_h\mathrm{d}x\mathrm{d}t={\int }_{\mathrm{\partial }Q}{n}_t{v}_h^2\mathrm{d}s,$ that (21) ${\int }_Q{\mathrm{\partial }}_t{v}_h{v}_h\mathrm{d}x\mathrm{d}t=\frac{1}{2}{\int }_Q{\mathrm{\partial }}_t{v}_h^2\mathrm{d}x\mathrm{d}t=\frac{1}{2}{\int }_{{\mathrm{\Sigma }}_T}{v}_h^2\mathrm{d}s-\frac{1}{2}{\int }_{{\mathrm{\Sigma }}_0}{v}_h^2\mathrm{d}s=\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2.$(21) The definition of $a_h$ and (21(21) ${\int }_Q{\mathrm{\partial }}_t{v}_h{v}_h\mathrm{d}x\mathrm{d}t=\frac{1}{2}{\int }_Q{\mathrm{\partial }}_t{v}_h^2\mathrm{d}x\mathrm{d}t=\frac{1}{2}{\int }_{{\mathrm{\Sigma }}_T}{v}_h^2\mathrm{d}s-\frac{1}{2}{\int }_{{\mathrm{\Sigma }}_0}{v}_h^2\mathrm{d}s=\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2.$(21) ) imply (22) $\begin{array}{rl}{a}_h(v_h,v_h)& ={\int }_Q\frac{1}{2}{\mathrm{\partial }}_t{v}_h^2+\theta h({\mathrm{\partial }}_t{v}_h^b{)}^2+\kappa |{\mathrm{\nabla }}_x{v}_h{|}^2\mathrm{d}x\mathrm{d}t\\ & =\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{v}_h^b{\parallel }_{L^2\left(Q\right)}^2,\end{array}$(22) which is (20(20) $a_h(v_h,v_h)\ge C_s\parallel v_h{\parallel }_h^2,\mathrm{\forall }{v}_h\in V_{h,b}.$(20) ) with $C_s=1$ and this completes the proof.

Proposition 3.1.

Let $u_h$ be the solution given by (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ). Then there is a $C_{\kappa ,\mathrm{\Omega }}>0$ such that the solution $u_h$ satisfies the following a priori estimate (23) $\parallel u_h{\parallel }_h\le C_{\kappa ,\mathrm{\Omega }}\parallel f{\parallel }_{L^2\left(Q\right)}.$(23)

Proof.

Using $u_h\in V_{h,b}$ as a test function in (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ), and utilizing (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ) and (20(20) $a_h(v_h,v_h)\ge C_s\parallel v_h{\parallel }_h^2,\mathrm{\forall }{v}_h\in V_{h,b}.$(20) )) together with the Poincare inequality (6(6) ${\int }_{\mathrm{\Omega }}{v}^2\left(x\right)\mathrm{d}x\le C_{\mathrm{\Omega }}^2{\int }_{\mathrm{\Omega }}({\mathrm{\partial }}_{x_1}v{)}^2\mathrm{d}x.$(6) ), we successively obtain (24) $\begin{array}{rl}{C}_s\parallel u_h{\parallel }_h^2& \le a_h(u_h,u_h)\le \frac{1}{{\kappa }^{1/2}}\left|{\int }_Q{\kappa }^{1/2}f{u}_h\mathrm{d}x\mathrm{d}t\right|\le C_{\mathrm{\Omega }}\frac{1}{{\kappa }^{1/2}}\parallel f{\parallel }_{L^2\left(Q\right)}\parallel {\kappa }^{1/2}{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}\\ & \le C_{\mathrm{\Omega }}\frac{1}{\kappa }\parallel f{\parallel }_{L^2\left(Q\right)}(\kappa \parallel {\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2},\end{array}$(24) where we have previously used that $C_s=1$ and $\kappa \le 1$. Setting $C_{\kappa ,\mathrm{\Omega }}=C_{\mathrm{\Omega }}\left(1/\kappa \right)$ we get (23(23) $\parallel u_h{\parallel }_h\le C_{\kappa ,\mathrm{\Omega }}\parallel f{\parallel }_{L^2\left(Q\right)}.$(23) ).

Note that the estimate in (23(23) $\parallel u_h{\parallel }_h\le C_{\kappa ,\mathrm{\Omega }}\parallel f{\parallel }_{L^2\left(Q\right)}.$(23) ) provides a direct control of $\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2$ due to the appearance of $b(\cdot ,\cdot )$ in the finite element scheme (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ), cf. (13(13) $\frac{1}{2}\parallel u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\frac{\kappa }{2}\parallel {\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2\le \frac{2C_{\mathrm{\Omega }}^2}{\kappa }\parallel f{\parallel }_{L^2\left(Q\right)}^2,$(13) ). A direct result of (23(23) $\parallel u_h{\parallel }_h\le C_{\kappa ,\mathrm{\Omega }}\parallel f{\parallel }_{L^2\left(Q\right)}.$(23) ) and (20(20) $a_h(v_h,v_h)\ge C_s\parallel v_h{\parallel }_h^2,\mathrm{\forall }{v}_h\in V_{h,b}.$(20) ) is the following corollary.

Corollary 3.1.

The discrete problem defined in (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ) is well posed, i.e. it has a unique solution which satisfies the stability estimate (23(23) $\parallel u_h{\parallel }_h\le C_{\kappa ,\mathrm{\Omega }}\parallel f{\parallel }_{L^2\left(Q\right)}.$(23) ).

Next, we show the boundedness of $a(\cdot ,\cdot )$ on $V_{0,\underset{\_}{0}}\times V_{h,b}$, where the space $V_{0,\underset{\_}{0}}$ is defined in Assumption 2.1. We define the norms (25a) $\parallel v{\parallel }_{h,\ast }=(\kappa \parallel {\mathrm{\nabla }}_{x}v{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_{t}v{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2},$(25a) (25b) $\parallel v{\parallel }_{h,V}=(\kappa \parallel {\mathrm{\nabla }}_{x}v{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_{t}v{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+(\theta h{)}^{-1}\parallel v{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}.$(25b) We point out that the $\parallel \cdot {\parallel }_{h,\ast }$ is defined for all $v\in V_{0,\underset{\_}{0}}$ and is similar to the $\parallel \cdot {\parallel }_h$ given in (19(19) $\parallel v_h{\parallel }_h=(\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{v}_h^b{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}.$(19) ), which is defined for all $v_h\in V_{h,b}$. The mesh-dependent norm $\parallel \cdot {\parallel }_{h,V}$ will be used later for deriving bounds for the discretization error.

Lemma 3.3.

There is a constant $C_b(\kappa ,\theta ,h)>0$ such that (26) $|a(u,v_h)|\le C_b(\kappa ,\theta ,h)\parallel u{\parallel }_{h,V}\parallel v_h{\parallel }_h,\mathrm{\forall }(u,v_h)\in (V_{0,\underset{\_}{0}}\times V_{h,b}).$(26)

Proof.

Let $v_h=v_h^1+v_h^b\in V_{h,b}$. We treat every term of the form $a(\cdot ,\cdot )$ separately. We apply integration by parts and (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ) to infer (27) $\begin{array}{rl}{\int }_Q{\mathrm{\partial }}_{t}u{v}_h\mathrm{d}x\mathrm{d}t& =-{\int }_{Q}u{\mathrm{\partial }}_t{v}_h\mathrm{d}x\mathrm{d}t+{\int }_{{\mathrm{\Sigma }}_T}u{v}_h\mathrm{d}\sigma \\ & \le \left(\right(\theta h{)}^{-1}\parallel u{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_t{v}_h{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}+\parallel u{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}\\ & \stackrel{\left(16\right)}{\le }\left(\right(\theta h{)}^{-1}\parallel u{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}{\left(\frac{c_1\theta h}{h^2\kappa }\kappa \parallel v_h{\parallel }_{L^2\left(Q\right)}^2\right)}^{1/2}+2{\left(\frac{1}{2}\parallel u{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\right)}^{1/2}{\left(\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\right)}^{1/2}\\ & \stackrel{\left(6\right)}{\le }{\left((\theta h{)}^{-1}\parallel u{\parallel }_{L^2\left(Q\right)}^2\right)}^{1/2}\left(c_2\theta \right(\kappa h{)}^{-1}{)}^{1/2}{\left(\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{v}_h^b{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\right)}^{1/2}\\ & +2\parallel u{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}{\left(\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{v}_h^b{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel v_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\right)}^{1/2}\\ & \le \left(c_2\theta \right(\kappa h{)}^{-1}{)}^{1/2}\parallel u{\parallel }_{h,V}\parallel v_h{\parallel }_h+2\parallel u{\parallel }_{h,V}\parallel v_h{\parallel }_h\le c_3\left(\theta \right(\kappa h{)}^{-1}+1{)}^{1/2}\parallel u{\parallel }_{h,V}\parallel v_h{\parallel }_h,\end{array}$(27)

where $c_3$ depends on the constants appearing in (16(16) $\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(E\right)}\le c_{\mathrm{inv}}{h}^{-1}\parallel v_h{\parallel }_{L^2\left(E\right)},$(16) ) and (6(6) ${\int }_{\mathrm{\Omega }}{v}^2\left(x\right)\mathrm{d}x\le C_{\mathrm{\Omega }}^2{\int }_{\mathrm{\Omega }}({\mathrm{\partial }}_{x_1}v{)}^2\mathrm{d}x.$(6) ). Similarly for the second term, applying (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ), we get (28) ${\int }_Q{\kappa }^{1/2}{\mathrm{\nabla }}_{x}u\cdot {\kappa }^{1/2}{\mathrm{\nabla }}_x{v}_h\mathrm{d}x\mathrm{d}t\le (\kappa \parallel {\mathrm{\nabla }}_{x}u{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}(\kappa \parallel {\mathrm{\nabla }}_x{v}_h{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\le \parallel u{\parallel }_{h,V}\parallel v_h{\parallel }_h.$(28) Combining all the above bounds and setting $C_b(\kappa ,\theta ,h)=2c_3\left(\theta \right(\kappa h{)}^{-1}+1{)}^{1/2}$, we can derive the desired result.

4. Error analysis

Proposition 4.1

weak consistency

Let Assumption 2.1 hold and let $u_h$ the solution given by (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ), and furthermore let $z_h\in V_{h,b}$ and $v_h^1\in V_{h0}$. Then (29a) $a_h(u_h,z_h)=a(u,z_h),\mathrm{and}$(29a) (29b) $a_h(v_h^1,z_h)=a(v_h^1,z_h).$(29b)

Proof.

Multiplying (7(7) $\begin{array}{rl}{u}_t-\kappa \mathrm{\Delta }u& =f\mathrm{in}Q\text{and}\\ u& =0\mathrm{on}\mathrm{\Sigma },u(\cdot ,0)=u_0\mathrm{on}{\mathrm{\Sigma }}_0,\end{array}$(7) ) by $z_h\in V_{h,b}$, integrating over Q, and then applying integration by parts, we arrive at the variational identity (30) $a(u,z_h)=(f,z_h).$(30) We recall the problem (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ) and directly have (31) $a_h(u_h,z_h)=(f,z_h)=a(u,z_h).$(31) Furthermore, we observe that $b_h(v_h^1,z_h)=0$ for any $v_h^1\in V_{h0}$, and (29b(29b) $a_h(v_h^1,z_h)=a(v_h^1,z_h).$(29b) ) directly follows.

Lemma 4.1.

Let $u_h$ solve (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ) and let $z_h^1\in V_{h0}$. Under Assumption 2.1, there exists a c, independent of h such that (32) $(\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2)\le \frac{c}{\kappa }(\parallel \mathrm{\nabla }u-\mathrm{\nabla }{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2).$(32)

Proof.

Let $z_h^1$ be a function in $V_{h0}$. By (29(29a) $a_h(u_h,z_h)=a(u,z_h),\mathrm{and}$(29a) ) and by subtracting similar terms, we have that (33) $\begin{array}{rl}& {\int }_Q({\mathrm{\partial }}_t{u}_h-{\mathrm{\partial }}_t{z}_h^1){\phi }_h+\kappa {\int }_Q{\mathrm{\nabla }}_x(u_h-z_h^1)\cdot {\mathrm{\nabla }}_x{\phi }_h+\theta h{\int }_Q{\mathrm{\partial }}_t{u}_h^b{\mathrm{\partial }}_t{\phi }_h^b\\ & ={\int }_Q({\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1){\phi }_h+\kappa {\int }_Q{\mathrm{\nabla }}_x(u-z_h^1)\cdot {\mathrm{\nabla }}_x{\phi }_h,\mathrm{\forall }{\phi }_h\in V_{h,b}.\end{array}$(33) Setting above ${\phi }_h=u_h^1+u_h^b-z_h^1$ and applying integration by parts on the first term on the left side, we have (34) $\begin{array}{rl}& {\int }_{{\mathrm{\Sigma }}_T}|{\mathrm{\partial }}_t{u}_h-{\mathrm{\partial }}_t{z}_h^1{|}^2\mathrm{d}\sigma +\kappa {\int }_Q|{\mathrm{\nabla }}_x(u_h-z_h^1){|}^2\mathrm{d}x\mathrm{d}t+\theta h{\int }_Q|{\mathrm{\partial }}_t{u}_h^b{|}^2\mathrm{d}x\mathrm{d}t\\ & \le {\left(\frac{1}{\kappa }{\int }_Q|{\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{|}^2\right)}^{1/2}{\left(\kappa {\int }_Q|u_h-z_h^1{|}^2\right)}^{1/2}\\ & +{\left({\int }_Q\kappa |{\mathrm{\nabla }}_x(u-z_h^1){|}^2\right)}^{1/2}{\left({\int }_Q\kappa |{\mathrm{\nabla }}_x(u_h-z_h^1){|}^2\right)}^{1/2},\end{array}$(34) and by applying (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ) and (4(4) $v(x_1,\dots ,x_d)=v\left(x_{\mathrm{\Gamma }}\right)+{\int }_{x_{\mathrm{\Gamma },1}}^{x_1}\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_1}(\tau ,x_2,\dots ,x_d)\mathrm{d}\tau .$(4) ) on the right hand side, yields (35) $\begin{array}{rl}& \parallel u_h-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2\le \frac{1}{c_{ϵ}\kappa }\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2\\ & +c_{ϵ}\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{\kappa }{c_{ϵ}}\parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+c_{ϵ}\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2.\end{array}$(35) Gathering the same terms and setting $0<c_{ϵ}<\frac{1}{2}$, we get (36) $\begin{array}{rl}& (1-2c_{ϵ})(\parallel u_h-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2)\\ & \le \frac{1}{c_{ϵ}\kappa }\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{\kappa }{c_{ϵ}}\parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2.\end{array}$(36) Using $0<\kappa \le 1$, applying triangle inequality and setting $c=\left(1/(1-2c_{ϵ})\right)\left(1/c_{ϵ}\right)$, the assertion follows.

Below we give the main error bound for the finite element solution $u_h\in V_{h,b}$.

Theorem 4.1.

Let $u_h=u_h^1+u_h^b$ solve (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ). Under Assumption 2.1 and choosing $\theta \ge h,$ there exist a constant $c_{\ast ,V},$ depending on $c_{\mathrm{inv}}$ in (16(16) $\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(E\right)}\le c_{\mathrm{inv}}{h}^{-1}\parallel v_h{\parallel }_{L^2\left(E\right)},$(16) ), such that (37) $\parallel u-u_h{\parallel }_{h,\ast }^2\le c_{\ast ,V}^2\left({\mu }_1\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\mu }_2(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2),for{z}_h^1\in V_{h0},$(37) where ${\mu }_1(\kappa ,\theta ,h)=(1+{\tilde{\gamma }}^2(\kappa ,\theta ,h)+{\tilde{\gamma }}^3(\kappa ,\theta ,h)+\tilde{\gamma }(\kappa ,\theta ,h\left)\right)$ and ${\mu }_2(\kappa ,\theta ,h)={\tilde{\gamma }}^2(\kappa ,\theta ,h)h^{-2}$, with $\tilde{\gamma }(\kappa ,\theta ,h)=\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)$ and $\gamma (\kappa ,\theta ,h)=(1/(\theta h{)}^{1/2}+c_{\mathrm{inv}}{\kappa }^{1/2}/h)$.

Proof.

Let us consider the functions $z_h^1\in V_{h0}$ and ${\sigma }_h=(u_h^1+u_h^b)-z_h^1$. Using triangle inequality, we decompose the error as (38) $\begin{array}{rl}& \frac{1}{2}\parallel u-u_h{\parallel }_{h,\ast }^2=\frac{1}{2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{u}_h{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2)\\ & \le \stackrel{T_1}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}+\stackrel{T_2}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2}}\\ & +\stackrel{T_3}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^b-0\cdot {\mathrm{\partial }}_t{z}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u_h-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}\\ & \le \parallel u-z_h^1{\parallel }_{h,V}^2+T_2+T_3,\end{array}$(38) where we previously set $T_1:=\parallel u-z_h^1{\parallel }_{h,\ast }^2$, $T_2:=\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2$, $T_3:=\parallel {\sigma }_h{\parallel }_h^2$, and in the last inequality we used that $T_1\le \parallel u-z_h^1{\parallel }_{h,V}^2.$ We will proceed by giving bounds for every term appearing in (38(38) $\begin{array}{rl}& \frac{1}{2}\parallel u-u_h{\parallel }_{h,\ast }^2=\frac{1}{2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{u}_h{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2)\\ & \le \stackrel{T_1}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}+\stackrel{T_2}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2}}\\ & +\stackrel{T_3}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^b-0\cdot {\mathrm{\partial }}_t{z}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u_h-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}\\ & \le \parallel u-z_h^1{\parallel }_{h,V}^2+T_2+T_3,\end{array}$(38) ). We first show few auxiliary results. Let $v_h\in V_{h,b}$ and ${\sigma }_h^1=u_h^1-z_h^1$, then using (16(16) $\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(E\right)}\le c_{\mathrm{inv}}{h}^{-1}\parallel v_h{\parallel }_{L^2\left(E\right)},$(16) ), (18(18) $a_h(u_h,v_h)=a(u_h,v_h)+b_h(u_h^b,v_h^b),$(18) ) and (29(29a) $a_h(u_h,z_h)=a(u,z_h),\mathrm{and}$(29a) ) and the fact that $0<h,\theta ,\kappa \le 1$, we can obtain that (39) $\begin{array}{rl}& \parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}=\parallel {\mathrm{\partial }}_t{\sigma }_h^1{\parallel }_{L^2\left(Q\right)}\le \underset{v_h\in V_{h,b}}{sup}\frac{|{\int }_Q{\mathrm{\partial }}_t\left(\right(u_h^1+u_h^b)-z_h^1-u_h^b)v_h\mathrm{d}x\mathrm{d}t|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \le \underset{v_h\in V_{h,b}}{sup}\frac{\left|-\left({\int }_Q{\mathrm{\partial }}_t{u}_h^b{v}_h+\kappa {\mathrm{\nabla }}_x{\sigma }_h\cdot {\mathrm{\nabla }}_x{v}_h+\theta h{\mathrm{\partial }}_t{u}_h^b{\mathrm{\partial }}_t{v}_h^b\mathrm{d}x\mathrm{d}t\right)\right|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & +\underset{v_h\in V_{h,b}}{sup}\frac{\left|{\int }_Q{\mathrm{\partial }}_t(u-{\mathrm{\partial }}_t{z}_h^1)v_h+\kappa {\mathrm{\nabla }}_x(u-z_h^1)\cdot {\mathrm{\nabla }}_x{v}_h\right|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \stackrel{\left(16\right)}{\le }\frac{(\theta h{)}^{1/2}}{(\theta h{)}^{1/2}}(\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}+c_{\mathrm{inv}}\theta \parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}))+\parallel {\mathrm{\partial }}_t(u-z_h^1){\parallel }_{L^2\left(Q\right)}\\ & +{\kappa }^{1/2}({\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x(u-z_h^1){\parallel }_{L^2\left(Q\right)}+{\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}))\underset{v_h\in V_{h,b}}{sup}\frac{\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(Q\right)}}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \le \frac{c_0}{(\theta h{)}^{1/2}}(\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}\\ & +\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}(\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}\\ & +\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}({\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x(u-z_h^1\left){\parallel }_{L^2\left(Q\right)}\right)+\frac{(\theta h{)}^{1/2}}{(\theta h{)}^{1/2}}\parallel {\mathrm{\partial }}_t(u-z_h^1){\parallel }_{L^2\left(Q\right)}\\ & \le c_0\left(\frac{1}{(\theta h{)}^{1/2}}+\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}\right)\parallel {\sigma }_h{\parallel }_h+\left(\frac{1}{(\theta h{)}^{1/2}}+\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}\right)\parallel u-z_h^1{\parallel }_{h,V}\\ & \le c_0\gamma (\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h),\end{array}$(39) where we previously set $c_0=(1+\theta c_{\mathrm{inv}})$. Using (39(39) $\begin{array}{rl}& \parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}=\parallel {\mathrm{\partial }}_t{\sigma }_h^1{\parallel }_{L^2\left(Q\right)}\le \underset{v_h\in V_{h,b}}{sup}\frac{|{\int }_Q{\mathrm{\partial }}_t\left(\right(u_h^1+u_h^b)-z_h^1-u_h^b)v_h\mathrm{d}x\mathrm{d}t|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \le \underset{v_h\in V_{h,b}}{sup}\frac{\left|-\left({\int }_Q{\mathrm{\partial }}_t{u}_h^b{v}_h+\kappa {\mathrm{\nabla }}_x{\sigma }_h\cdot {\mathrm{\nabla }}_x{v}_h+\theta h{\mathrm{\partial }}_t{u}_h^b{\mathrm{\partial }}_t{v}_h^b\mathrm{d}x\mathrm{d}t\right)\right|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & +\underset{v_h\in V_{h,b}}{sup}\frac{\left|{\int }_Q{\mathrm{\partial }}_t(u-{\mathrm{\partial }}_t{z}_h^1)v_h+\kappa {\mathrm{\nabla }}_x(u-z_h^1)\cdot {\mathrm{\nabla }}_x{v}_h\right|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \stackrel{\left(16\right)}{\le }\frac{(\theta h{)}^{1/2}}{(\theta h{)}^{1/2}}(\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}+c_{\mathrm{inv}}\theta \parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}))+\parallel {\mathrm{\partial }}_t(u-z_h^1){\parallel }_{L^2\left(Q\right)}\\ & +{\kappa }^{1/2}({\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x(u-z_h^1){\parallel }_{L^2\left(Q\right)}+{\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}))\underset{v_h\in V_{h,b}}{sup}\frac{\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(Q\right)}}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \le \frac{c_0}{(\theta h{)}^{1/2}}(\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}\\ & +\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}(\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}\\ & +\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}({\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x(u-z_h^1\left){\parallel }_{L^2\left(Q\right)}\right)+\frac{(\theta h{)}^{1/2}}{(\theta h{)}^{1/2}}\parallel {\mathrm{\partial }}_t(u-z_h^1){\parallel }_{L^2\left(Q\right)}\\ & \le c_0\left(\frac{1}{(\theta h{)}^{1/2}}+\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}\right)\parallel {\sigma }_h{\parallel }_h+\left(\frac{1}{(\theta h{)}^{1/2}}+\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}\right)\parallel u-z_h^1{\parallel }_{h,V}\\ & \le c_0\gamma (\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h),\end{array}$(39) ) and the relation (40) $(\theta h{)}^{1/2}\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}+(\theta h{)}^{1/2}\parallel {\mathrm{\partial }}_t{u}_h^b-0{\mathrm{\partial }}_t{z}_h^b{\parallel }_{L^2\left(Q\right)}\ge (\theta h{)}^{1/2}\parallel {\mathrm{\partial }}_t(u_h^1+u_h^b)-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)},$(40) and observing that $\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\ge 1$, it follows that (41) $(\theta h{)}^{1/2}\parallel {\mathrm{\partial }}_t{\sigma }_h{\parallel }_{L^2\left(Q\right)}\le 2c_0(\theta h{)}^{1/2}\gamma (\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h).$(41) Furthermore, working as in the proof of (26(26) $|a(u,v_h)|\le C_b(\kappa ,\theta ,h)\parallel u{\parallel }_{h,V}\parallel v_h{\parallel }_h,\mathrm{\forall }(u,v_h)\in (V_{0,\underset{\_}{0}}\times V_{h,b}).$(26) ) and using (41(41) $(\theta h{)}^{1/2}\parallel {\mathrm{\partial }}_t{\sigma }_h{\parallel }_{L^2\left(Q\right)}\le 2c_0(\theta h{)}^{1/2}\gamma (\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h).$(41) ) and (3(3) $\left|{\int }_{Q}uv\mathrm{d}x\right|\le \parallel u{\parallel }_{L^2\left(Q\right)}\parallel v{\parallel }_{L^2\left(Q\right)}\text{and}\left|{\int }_{Q}uv\mathrm{d}x\right|\le \frac{ϵ}{2}\parallel u{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2ϵ}\parallel v{\parallel }_{L^2\left(Q\right)}^2,$(3) ), we can find (42) $\begin{array}{rl}& a(u-z_h^1,{\sigma }_h)\\ & =-{\int }_Q(u-z_h^1){\mathrm{\partial }}_t{\sigma }_h\mathrm{d}x\mathrm{d}t+{\int }_{{\mathrm{\Sigma }}_T}(u-z_h^1){\sigma }_{h}d\sigma +{\int }_Q{\kappa }^{1/2}{\mathrm{\nabla }}_x(u-z_h^1)\cdot {\kappa }^{1/2}{\mathrm{\nabla }}_x{\sigma }_h\mathrm{d}x\mathrm{d}t\\ & \le \left(\right(\theta h{)}^{-1}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_t{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}+\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}\\ & +(\kappa \parallel {\mathrm{\nabla }}_x(u-z_h^1\left){\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\right(\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\\ & \le 2c_0\left(\right(\theta h{)}^{-1}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h\left)\right(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h)+2\parallel u-z_h^1{\parallel }_{h,V}\parallel {\sigma }_h{\parallel }_h\\ & \le \left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_0{h}^{-1}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}h\gamma (\kappa ,\theta ,h)\parallel {\sigma }_h{\parallel }_h+2\parallel u-z_h^1{\parallel }_{h,V}\parallel {\sigma }_h{\parallel }_h\\ & \le \left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2+\epsilon h^2{\gamma }^2(\kappa ,\theta ,h)\parallel {\sigma }_h{\parallel }_h^2\\ & +c_{\epsilon }\parallel u-z_h^1{\parallel }_{h,V}^2+\epsilon \parallel {\sigma }_h{\parallel }_h^2\\ & \le \left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2+\epsilon (h{\theta }^{-1}+c_{\mathrm{inv}}\kappa )\parallel {\sigma }_h{\parallel }_h^2\\ & +c_{\epsilon }\parallel u-z_h^1{\parallel }_{h,V}^2+\epsilon \parallel {\sigma }_h{\parallel }_h^2,\text{(using that }\theta \ge h\text{)}\\ & \le 2c_{\epsilon }\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2+2\epsilon C_{inv,\kappa }\parallel {\sigma }_h{\parallel }_h^2,\end{array}$(42) where we used that $\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)>1$, we chose $c_{\epsilon }>1$, and in the last step we set $C_{inv,k}=1+c_{\mathrm{inv}}\kappa$. Furthermore, using the properties of $a_h(\cdot ,\cdot )$ and (29(29a) $a_h(u_h,z_h)=a(u,z_h),\mathrm{and}$(29a) ), we have (43) $C_s\parallel {\sigma }_h{\parallel }_h^2\le a_h({\sigma }_h,{\sigma }_h)=a_h(u_h,{\sigma }_h)-a_h(z_h^1,{\sigma }_h)\stackrel{\left(29\right)}{=}a(u,{\sigma }_h)-a(z_h^1,{\sigma }_h)=a(u-z_h^1,{\sigma }_h).$(43) By replacing $C_s=1$ and by choosing $\epsilon \le 1/4C_{inv,\kappa }$ in (42(42) $\begin{array}{rl}& a(u-z_h^1,{\sigma }_h)\\ & =-{\int }_Q(u-z_h^1){\mathrm{\partial }}_t{\sigma }_h\mathrm{d}x\mathrm{d}t+{\int }_{{\mathrm{\Sigma }}_T}(u-z_h^1){\sigma }_{h}d\sigma +{\int }_Q{\kappa }^{1/2}{\mathrm{\nabla }}_x(u-z_h^1)\cdot {\kappa }^{1/2}{\mathrm{\nabla }}_x{\sigma }_h\mathrm{d}x\mathrm{d}t\\ & \le \left(\right(\theta h{)}^{-1}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_t{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}+\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}\\ & +(\kappa \parallel {\mathrm{\nabla }}_x(u-z_h^1\left){\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\right(\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\\ & \le 2c_0\left(\right(\theta h{)}^{-1}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h\left)\right(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h)+2\parallel u-z_h^1{\parallel }_{h,V}\parallel {\sigma }_h{\parallel }_h\\ & \le \left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_0{h}^{-1}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}h\gamma (\kappa ,\theta ,h)\parallel {\sigma }_h{\parallel }_h+2\parallel u-z_h^1{\parallel }_{h,V}\parallel {\sigma }_h{\parallel }_h\\ & \le \left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2+\epsilon h^2{\gamma }^2(\kappa ,\theta ,h)\parallel {\sigma }_h{\parallel }_h^2\\ & +c_{\epsilon }\parallel u-z_h^1{\parallel }_{h,V}^2+\epsilon \parallel {\sigma }_h{\parallel }_h^2\\ & \le \left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2+\epsilon (h{\theta }^{-1}+c_{\mathrm{inv}}\kappa )\parallel {\sigma }_h{\parallel }_h^2\\ & +c_{\epsilon }\parallel u-z_h^1{\parallel }_{h,V}^2+\epsilon \parallel {\sigma }_h{\parallel }_h^2,\text{(using that }\theta \ge h\text{)}\\ & \le 2c_{\epsilon }\left(\theta h{)}^{1/2}\gamma \right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+2c_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2+2\epsilon C_{inv,\kappa }\parallel {\sigma }_h{\parallel }_h^2,\end{array}$(42) ), we obtain that (44) $\frac{1}{2}\parallel {\sigma }_h{\parallel }_h^2\le 2c_{\epsilon }\left(\right(\theta h{)}^{1/2}\gamma (\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2).$(44) Now, we can bound the terms in (38(38) $\begin{array}{rl}& \frac{1}{2}\parallel u-u_h{\parallel }_{h,\ast }^2=\frac{1}{2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{u}_h{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2)\\ & \le \stackrel{T_1}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}+\stackrel{T_2}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2}}\\ & +\stackrel{T_3}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^b-0\cdot {\mathrm{\partial }}_t{z}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u_h-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}\\ & \le \parallel u-z_h^1{\parallel }_{h,V}^2+T_2+T_3,\end{array}$(38) ). Recalling the definition of $\tilde{\gamma }(\kappa ,\theta ,h)$, inequality (44(44) $\frac{1}{2}\parallel {\sigma }_h{\parallel }_h^2\le 2c_{\epsilon }\left(\right(\theta h{)}^{1/2}\gamma (\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2).$(44) ) immediately implies (45) $T_3=\parallel {\sigma }_h{\parallel }_h^2\le C_{\epsilon }\tilde{\gamma }(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+C_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2.$(45) Also, combining (39(39) $\begin{array}{rl}& \parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}=\parallel {\mathrm{\partial }}_t{\sigma }_h^1{\parallel }_{L^2\left(Q\right)}\le \underset{v_h\in V_{h,b}}{sup}\frac{|{\int }_Q{\mathrm{\partial }}_t\left(\right(u_h^1+u_h^b)-z_h^1-u_h^b)v_h\mathrm{d}x\mathrm{d}t|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \le \underset{v_h\in V_{h,b}}{sup}\frac{\left|-\left({\int }_Q{\mathrm{\partial }}_t{u}_h^b{v}_h+\kappa {\mathrm{\nabla }}_x{\sigma }_h\cdot {\mathrm{\nabla }}_x{v}_h+\theta h{\mathrm{\partial }}_t{u}_h^b{\mathrm{\partial }}_t{v}_h^b\mathrm{d}x\mathrm{d}t\right)\right|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & +\underset{v_h\in V_{h,b}}{sup}\frac{\left|{\int }_Q{\mathrm{\partial }}_t(u-{\mathrm{\partial }}_t{z}_h^1)v_h+\kappa {\mathrm{\nabla }}_x(u-z_h^1)\cdot {\mathrm{\nabla }}_x{v}_h\right|}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \stackrel{\left(16\right)}{\le }\frac{(\theta h{)}^{1/2}}{(\theta h{)}^{1/2}}(\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}+c_{\mathrm{inv}}\theta \parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}))+\parallel {\mathrm{\partial }}_t(u-z_h^1){\parallel }_{L^2\left(Q\right)}\\ & +{\kappa }^{1/2}({\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x(u-z_h^1){\parallel }_{L^2\left(Q\right)}+{\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}))\underset{v_h\in V_{h,b}}{sup}\frac{\parallel \mathrm{\nabla }{v}_h{\parallel }_{L^2\left(Q\right)}}{\parallel v_h{\parallel }_{L^2\left(Q\right)}}\\ & \le \frac{c_0}{(\theta h{)}^{1/2}}(\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}\\ & +\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}(\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{\sigma }_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel {\sigma }_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2{)}^{1/2}\\ & +\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}({\kappa }^{1/2}\parallel {\mathrm{\nabla }}_x(u-z_h^1\left){\parallel }_{L^2\left(Q\right)}\right)+\frac{(\theta h{)}^{1/2}}{(\theta h{)}^{1/2}}\parallel {\mathrm{\partial }}_t(u-z_h^1){\parallel }_{L^2\left(Q\right)}\\ & \le c_0\left(\frac{1}{(\theta h{)}^{1/2}}+\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}\right)\parallel {\sigma }_h{\parallel }_h+\left(\frac{1}{(\theta h{)}^{1/2}}+\frac{c_{\mathrm{inv}}{\kappa }^{1/2}}{h}\right)\parallel u-z_h^1{\parallel }_{h,V}\\ & \le c_0\gamma (\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}+\parallel {\sigma }_h{\parallel }_h),\end{array}$(39) ) and (45(45) $T_3=\parallel {\sigma }_h{\parallel }_h^2\le C_{\epsilon }\tilde{\gamma }(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+C_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2.$(45) ), we have that (46) $\begin{array}{rl}{T}_2& =\theta h\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2\le 2{\tilde{\gamma }}^2(\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}^2+\parallel {\sigma }_h{\parallel }_h^2)\\ & \le 2C_{\epsilon }\left({\tilde{\gamma }}^2\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\tilde{\gamma }}^2(\kappa ,\theta ,h\left)\tilde{\gamma }\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\tilde{\gamma }}^2(\kappa ,\theta ,h)h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2)\\ & \le 2C_{\epsilon }\left(\right({\tilde{\gamma }}^2(\kappa ,\theta ,h)+{\tilde{\gamma }}^2(\kappa ,\theta ,h)\tilde{\gamma }(\kappa ,\theta ,h))\parallel u-z_h^1{\parallel }_{h,V}^2+{\tilde{\gamma }}^2(\kappa ,\theta ,h)h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2).\end{array}$(46) Finally, inserting the bounds (46(46) $\begin{array}{rl}{T}_2& =\theta h\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2\le 2{\tilde{\gamma }}^2(\kappa ,\theta ,h)(\parallel u-z_h^1{\parallel }_{h,V}^2+\parallel {\sigma }_h{\parallel }_h^2)\\ & \le 2C_{\epsilon }\left({\tilde{\gamma }}^2\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\tilde{\gamma }}^2(\kappa ,\theta ,h\left)\tilde{\gamma }\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\tilde{\gamma }}^2(\kappa ,\theta ,h)h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2)\\ & \le 2C_{\epsilon }\left(\right({\tilde{\gamma }}^2(\kappa ,\theta ,h)+{\tilde{\gamma }}^2(\kappa ,\theta ,h)\tilde{\gamma }(\kappa ,\theta ,h))\parallel u-z_h^1{\parallel }_{h,V}^2+{\tilde{\gamma }}^2(\kappa ,\theta ,h)h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2).\end{array}$(46) ) and (45(45) $T_3=\parallel {\sigma }_h{\parallel }_h^2\le C_{\epsilon }\tilde{\gamma }(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+C_{\epsilon }{h}^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2.$(45) ) in (38(38) $\begin{array}{rl}& \frac{1}{2}\parallel u-u_h{\parallel }_{h,\ast }^2=\frac{1}{2}\left(\right(\theta h)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{u}_h{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2)\\ & \le \stackrel{T_1}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_{t}u-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}+\stackrel{T_2}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^1-{\mathrm{\partial }}_t{z}_h^1{\parallel }_{L^2\left(Q\right)}^2}}\\ & +\stackrel{T_3}{\overbrace{\left(\theta h\right)\parallel {\mathrm{\partial }}_t{u}_h^b-0\cdot {\mathrm{\partial }}_t{z}_h^b{\parallel }_{L^2\left(Q\right)}^2+\kappa \parallel {\mathrm{\nabla }}_x{u}_h-{\mathrm{\nabla }}_x{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\frac{1}{2}\parallel u_h-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2}}\\ & \le \parallel u-z_h^1{\parallel }_{h,V}^2+T_2+T_3,\end{array}$(38) ) and using that ${\tilde{\gamma }}^2\ge 1$, we obtain (47) $\begin{array}{rl}\frac{1}{2}\parallel u-u_h{\parallel }_{h,\ast }^2& \le 2C_{\epsilon }(1+{\tilde{\gamma }}^2(\kappa ,\theta ,h)+{\tilde{\gamma }}^2(\kappa ,\theta ,h\left)\tilde{\gamma }\right(\kappa ,\theta ,h)+\tilde{\gamma }(\kappa ,\theta ,h\left)\right)\parallel u-z_h{\parallel }_{h,V}^2\\ & +2C_{\epsilon }2{\tilde{\gamma }}^2(\kappa ,\theta ,h)h^{-2}\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2.\end{array}$(47) Setting $c_{\ast ,V}^2=4C_{\epsilon }$, we can derive estimate (37(37) $\parallel u-u_h{\parallel }_{h,\ast }^2\le c_{\ast ,V}^2\left({\mu }_1\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\mu }_2(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2),for{z}_h^1\in V_{h0},$(37) ).

Remark 4.1

Let us consider a fixed $T_h\left(Q\right)$ and let us recall the forms of ${\mu }_1(\kappa ,\theta ,h)$ and ${\mu }_2(\kappa ,\theta ,h)$ in (37(37) $\parallel u-u_h{\parallel }_{h,\ast }^2\le c_{\ast ,V}^2\left({\mu }_1\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\mu }_2(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2),for{z}_h^1\in V_{h0},$(37) ). Then (i) setting a fixed value $\theta ={\theta }_c$, we have ${\mu }_1(\kappa ,{\theta }_c,h)\sim {\theta }_c{h}^{-1}+{\theta }_c^{1+1/2}{h}^{-(1+1/2)}+{\theta }_c^{1/2}{h}^{1/2}$ and ${\mu }_2(\kappa ,{\theta }_c,h)\sim$ $h^{-2}+{\theta }_c{h}^{-3}$, (ii) setting $\theta \approx h$, we have that ${\mu }_1(\kappa ,\theta ,h)\sim c_{\mathrm{inv}}\kappa$ and ${\mu }_2(\kappa ,\theta ,h)\sim h^{-2}$. Note that in the previous cases (i) and (ii), the associated constants depend on $c_{\mathrm{inv}}$ and κ.

Below, we recall some approximation estimates of the finite element space. For the proof we refer to [31].

Lemma 4.2.

Let $s,m$ be integers such that $0\le m\le 1<s$ and let the space $V_{h0}$ defined in (11(11) $V_{h0}=\{v_h\in H_{0,\underset{\_}{0}}^{1,1}(Q):v_h{|}_E\in {\mathbb{P}}^p(E),\mathrm{for}\mathrm{every}E\in T_h(Q\left)\right\},$(11) ). Then for every $v\in V_{0,\underset{\_}{0}}$ with $V_{0,\underset{\_}{0}}=H_{0,\underset{\_}{0}}^{1,1}\left(Q\right)\cap H^s\left(Q\right),$ there exists a linear interpolation operator ${\pi }_{h}v:V\to V_{h0}$ such that (48) $|v-{\pi }_{h}v{|}_{H^m\left(Q\right)}\le c_{\mathrm{intp}}{h}^{min(p+1,s)-m}\parallel v{\parallel }_{H^s\left(Q\right)},$(48) where $c_{\mathrm{intp}}=c(m,s,Q)$ and p=1.

Lemma 4.3.

Let the space $V_{h0}$ defined in (11(11) $V_{h0}=\{v_h\in H_{0,\underset{\_}{0}}^{1,1}(Q):v_h{|}_E\in {\mathbb{P}}^p(E),\mathrm{for}\mathrm{every}E\in T_h(Q\left)\right\},$(11) ). Let $s\ge 2$ be an integer, and let a function $v\in V_{0,\underset{\_}{0}}$ with $V_{0,\underset{\_}{0}}=H_{0,\underset{\_}{0}}^{1,1}\left(Q\right)\cap H^s\left(Q\right)$. There exists a linear interpolation operator ${\pi }_{h}v:V\to V_{h0}$ such that (49a) $\parallel v-{\pi }_{h}v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\le c_1{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49a) (49b) $\parallel v-{\pi }_{h}v{\parallel }_{h,V}^2\le c_2{h}^{2r-1}(\kappa h^{-1}+\theta +{\theta }^{-1}+1)\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49b) where r=2 and $c_1,c_2$ depend on the constants appearing in (17(17) $\parallel v{\parallel }_{L^2\left(\mathrm{\partial }E\right)}\le c_{\mathrm{trac}}{h}^{-1/2}(\parallel v{\parallel }_{L^2\left(E\right)}+h\parallel \mathrm{\nabla }v{\parallel }_{L^2\left(E\right)}).$(17) ) and in (48(48) $|v-{\pi }_{h}v{|}_{H^m\left(Q\right)}\le c_{\mathrm{intp}}{h}^{min(p+1,s)-m}\parallel v{\parallel }_{H^s\left(Q\right)},$(48) ), but not on h and v.

Proof.

Introducing the operator ${\pi }_{h}v$ of Lemma 4.2 and by applying (17(17) $\parallel v{\parallel }_{L^2\left(\mathrm{\partial }E\right)}\le c_{\mathrm{trac}}{h}^{-1/2}(\parallel v{\parallel }_{L^2\left(E\right)}+h\parallel \mathrm{\nabla }v{\parallel }_{L^2\left(E\right)}).$(17) ) and (48(48) $|v-{\pi }_{h}v{|}_{H^m\left(Q\right)}\le c_{\mathrm{intp}}{h}^{min(p+1,s)-m}\parallel v{\parallel }_{H^s\left(Q\right)},$(48) ), we have (50) $\parallel u-{\pi }_{h}v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\le 2c_{\mathrm{trac}}^2{h}^{-1}(\parallel v-{\pi }_{h}v{\parallel }_{L^2\left(Q\right)}^2)+h^2\parallel \mathrm{\nabla }(v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2)\le c_1{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2.$(50) In the same way, we have (51) $\begin{array}{rl}\kappa \parallel {\mathrm{\nabla }}_x(v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2& \le c_{\mathrm{intp}}\kappa h^{2r-2}\parallel v{\parallel }_{H^s\left(Q\right)}^2,\\ \theta h\parallel {\mathrm{\partial }}_t(v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2& \le c_{\mathrm{intp}}\theta h^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,\\ (\theta h{)}^{-1}\parallel (v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2& \le c_{\mathrm{intp}}{\theta }^{-1}{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,\end{array}$(51) Collecting the estimates (50(50) $\parallel u-{\pi }_{h}v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\le 2c_{\mathrm{trac}}^2{h}^{-1}(\parallel v-{\pi }_{h}v{\parallel }_{L^2\left(Q\right)}^2)+h^2\parallel \mathrm{\nabla }(v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2)\le c_1{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2.$(50) ) and (51(51) $\begin{array}{rl}\kappa \parallel {\mathrm{\nabla }}_x(v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2& \le c_{\mathrm{intp}}\kappa h^{2r-2}\parallel v{\parallel }_{H^s\left(Q\right)}^2,\\ \theta h\parallel {\mathrm{\partial }}_t(v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2& \le c_{\mathrm{intp}}\theta h^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,\\ (\theta h{)}^{-1}\parallel (v-{\pi }_{h}v){\parallel }_{L^2\left(Q\right)}^2& \le c_{\mathrm{intp}}{\theta }^{-1}{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,\end{array}$(51) ), we easily obtain (52) $\parallel v-{\pi }_{h}v{\parallel }_{h,V}^2\le c_2(\kappa h^{2r-2}+\theta h^{2r-1}+{\theta }^{-1}{h}^{2r-1}+h^{2r-1})\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(52) which is (49b(49b) $\parallel v-{\pi }_{h}v{\parallel }_{h,V}^2\le c_2{h}^{2r-1}(\kappa h^{-1}+\theta +{\theta }^{-1}+1)\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49b) ).

Remark 4.2

The interpolation estimates given in (49(49a) $\parallel v-{\pi }_{h}v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\le c_1{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49a) ) have been derived for linear polynomial spaces, see (11(11) $V_{h0}=\{v_h\in H_{0,\underset{\_}{0}}^{1,1}(Q):v_h{|}_E\in {\mathbb{P}}^p(E),\mathrm{for}\mathrm{every}E\in T_h(Q\left)\right\},$(11) ). Analogous estimates can be derived for higher polynomial spaces. In that case we set $r=min(p+1,s)$.

In the proposed scheme (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ), the parameter θ can control the artificial diffusion terms. Its value can be adjusted according to the needs of the scheme for obtaining optimal approximation properties, see Remark 4.1. Below, we show that if θ is close to the grid size h, we obtain optimal order of convergence. The value of θ can be determined without having to tune it with respect to κ. However, in some finite element schemes the value of θ must be adjusted with respect to h and κ for getting optimal rates, e.g. see discussion for the Galerkin/least square methods in [24,33].

Theorem 4.2

error estimates

Let u be the solution of (7(7) $\begin{array}{rl}{u}_t-\kappa \mathrm{\Delta }u& =f\mathrm{in}Q\text{and}\\ u& =0\mathrm{on}\mathrm{\Sigma },u(\cdot ,0)=u_0\mathrm{on}{\mathrm{\Sigma }}_0,\end{array}$(7) ) and let $u_h$ be the solution of (15(15a) $a(u_h,v_h)+b_h(u_h^b,v_h^b)=(f,v_h),\mathrm{\forall }{v}_h\in V_{h,b},$(15a) ), respectively. Under the Assumption 2.1, the solution $u_h$ satisfies the estimate (53) $\parallel u-u_h{\parallel }_{h,\ast }\le c_{1}h\parallel u{\parallel }_{H^s\left(Q\right)},for\theta \approx h,$(53) and moreover for any $\theta >0$, (54) $(\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\le \frac{c_2}{\kappa }h\parallel u{\parallel }_{H^s\left(Q\right)},$(54) with $c_1$ depending on the constants in (37(37) $\parallel u-u_h{\parallel }_{h,\ast }^2\le c_{\ast ,V}^2\left({\mu }_1\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\mu }_2(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2),for{z}_h^1\in V_{h0},$(37) ) and in (49(49a) $\parallel v-{\pi }_{h}v{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2\le c_1{h}^{2r-1}\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49a) ) and $c_2$ on the constants in (32(32) $(\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2)\le \frac{c}{\kappa }(\parallel \mathrm{\nabla }u-\mathrm{\nabla }{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2).$(32) ) and (48(48) $|v-{\pi }_{h}v{|}_{H^m\left(Q\right)}\le c_{\mathrm{intp}}{h}^{min(p+1,s)-m}\parallel v{\parallel }_{H^s\left(Q\right)},$(48) ).

Proof.

The estimate (53(53) $\parallel u-u_h{\parallel }_{h,\ast }\le c_{1}h\parallel u{\parallel }_{H^s\left(Q\right)},for\theta \approx h,$(53) ) follows directly from (29(29a) $a_h(u_h,z_h)=a(u,z_h),\mathrm{and}$(29a) ), (37(37) $\parallel u-u_h{\parallel }_{h,\ast }^2\le c_{\ast ,V}^2\left({\mu }_1\right(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{h,V}^2+{\mu }_2(\kappa ,\theta ,h)\parallel u-z_h^1{\parallel }_{L^2\left(Q\right)}^2),for{z}_h^1\in V_{h0},$(37) ), Remark (4.1) and (49b(49b) $\parallel v-{\pi }_{h}v{\parallel }_{h,V}^2\le c_2{h}^{2r-1}(\kappa h^{-1}+\theta +{\theta }^{-1}+1)\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49b) ). The estimate (54(54) $(\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2{)}^{1/2}\le \frac{c_2}{\kappa }h\parallel u{\parallel }_{H^s\left(Q\right)},$(54) ) follows form (32(32) $(\parallel u-u_h{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2+\kappa \parallel {\mathrm{\nabla }}_{x}u-{\mathrm{\nabla }}_x{u}_h{\parallel }_{L^2\left(Q\right)}^2+\theta h\parallel {\mathrm{\partial }}_t{u}_h^b{\parallel }_{L^2\left(Q\right)}^2)\le \frac{c}{\kappa }(\parallel \mathrm{\nabla }u-\mathrm{\nabla }{z}_h^1{\parallel }_{L^2\left(Q\right)}^2+\parallel u-z_h^1{\parallel }_{L^2\left({\mathrm{\Sigma }}_T\right)}^2).$(32) ) and (48(48) $|v-{\pi }_{h}v{|}_{H^m\left(Q\right)}\le c_{\mathrm{intp}}{h}^{min(p+1,s)-m}\parallel v{\parallel }_{H^s\left(Q\right)},$(48) ).

Remark 4.3

In realistic cases, the solutions of parabolic evolution problems may present an anisotropic regularity behavior, for example different regularities properties with respect to time and to space direction. This may result in the case where the contribution of the discretization error in t direction can be different than the contribution of the discretzation error in x direction. In such cases, it is more appropriate to discretize the problem using anisotropic meshes, using small mesh size in the directions where the solution is less smooth and larger mesh size in the directions where the solution is smoother [34]. This is a topic that we will investigate in a forthcoming paper. A discretization of these type of problems using Isogeometric Analysis methodology is presented in [20].

Remark 4.4

For the sake of simplicity, we have presented the discretization analysis for spaces with degree p=1. However, assuming further regularity properties on u, we can follow the same analysis and to show optimal estimates for spaces with degree p>1, see Example 2 in Section 5.

5. Numerical examples

In this section, we present several numerical examples for validating the theoretical estimates. Although in the analysis, we used linear polynomial spaces, next we perform tests using both linear, (p=1), and second order, (p=2), polynomial spaces, which are combined with the associated cubic bubble space. We perform tests for $Q\subset {\mathbb{R}}^{d+1}$ with d=1 and d=2, using triangular and tetrahedral mesh elements correspondingly. One can apply the proposed method for solving problems on $Q\subset {\mathbb{R}}^{d+1}$ with d=3. The code that we have at our disposal does not support such calculations yet. Every example has been solved applying several mesh refinement steps with corresponding mesh size $h_s=\left(h_0/2^s\right),s=1,2,\dots$. Every $T_{h_s}\left(Q\right)$ satisfies the properties mentioned in Section 3. We present tables with the asymptotic convergence rates r of the error. The convergence rates r have been computed by the ratio $r=\ln \left(e_s/e_{s+1}\right)/\ln \left(h_s/h_{s+1}\right)$, where the error $e_s:=\parallel u-u_h{\parallel }_{h,\ast }$ is computed on $T_{h_s}\left(Q\right)$. We mention that we use highly smooth solutions, i.e. $min(p+1,s)=p+1$, and thus the expected values for the rates are r=p, see (49b(49b) $\parallel v-{\pi }_{h}v{\parallel }_{h,V}^2\le c_2{h}^{2r-1}(\kappa h^{-1}+\theta +{\theta }^{-1}+1)\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49b) ) and Remark 4.2. We consider test cases with $\kappa \in \{1,0.005\}$, and we study the behavior of the rates for $\theta =h_s$. Lastly, we point out that since the support of a bubble function is restricted to the interior of the element, we eliminate the associated variable from the produced linear system by static condensation.

Example 1: $Q\subset {\mathbb{R}}^2$, p=1. In the first example, the problem is considered in $Q=(0,1)\times (0,2)$. The exact solution is given by the formula (55) $u(x,t)=\sin \left(2\pi x\right)\sin \left(\pi t\right).$(55)

The source function f is determined by (55(55) $u(x,t)=\sin \left(2\pi x\right)\sin \left(\pi t\right).$(55) ). Note that u=0 on Σ and $u_0=0$, see (7(7) $\begin{array}{rl}{u}_t-\kappa \mathrm{\Delta }u& =f\mathrm{in}Q\text{and}\\ u& =0\mathrm{on}\mathrm{\Sigma },u(\cdot ,0)=u_0\mathrm{on}{\mathrm{\Sigma }}_0,\end{array}$(7) ). In Figure 2, we plot the exact solution on Q. We solve the problem using linear polynomials, i.e. p=1. We begin by first setting $\kappa =1$. The numerical convergence rates for the several levels of the mesh refinement are presented in the second column in Table 1. They are in good agreement with the theoretically predicted estimates given in (53(53) $\parallel u-u_h{\parallel }_{h,\ast }\le c_{1}h\parallel u{\parallel }_{H^s\left(Q\right)},for\theta \approx h,$(53) ). The numerical solution $u_h\in V_{h,b}$ gives optimal convergence rates, i.e. the values of r are very close to one for all the refinement steps. Next, we want to investigate the asymptotic behavior of the numerical convergence rates when the value of κ is small. We perform the same computations by setting $\kappa =0.005$. The associated convergence rates are presented in the last column in Table 1. We observe that for the first mesh levels, the rates r are a little higher than the expected value. However, as we move on to the last mesh levels, the values of r are close to one, and are in agreement with the values predicted by the theory.

Figure 2. Example 1: The solution u on Q.

Table 1. Example 1: The convergence rates r.

$h_s$	$p=1,\kappa =1$	$p=1,\kappa =0.005$
$h_0/2^s$	Convergence rates r
s=1	0.98	1.61
s=2	0.99	1.51
s=3	0.99	1.20
s=4	0.98	1.35
s=5	1.00	1.17
s=6	0.99	1.05
s=7	1.00	1.01

Example 2: $Q\subset {\mathbb{R}}^2,p=2$. In the second example, we consider the problem on $Q=(0,1)\times (0,1)$. The exact solution is given by the formula (56) $u(x,t)=\sin \left(2\pi t\right)\sin \left(2\pi x\right).$(56) The source function f is defined to match the solution in (56(56) $u(x,t)=\sin \left(2\pi t\right)\sin \left(2\pi x\right).$(56) ). In Figure 3, we plot the exact solution u on a relative coarse mesh with h=0.25. We solve the problem using second order, i.e. p=2, polynomial space. For the first group of computations we set $\kappa =1$ and $\theta =h_s$. In the second column in Table 2, we show the convergence rates r. The values of r are approaching the value two, and are the expected rates based (49b(49b) $\parallel v-{\pi }_{h}v{\parallel }_{h,V}^2\le c_2{h}^{2r-1}(\kappa h^{-1}+\theta +{\theta }^{-1}+1)\parallel v{\parallel }_{H^s\left(Q\right)}^2,$(49b) ) and Remark 4.2. We repeat the same computations setting $\kappa =0.005$ and keep the same values for θ. The produced rates r are shown in the last column in Table 2. We observe that, for the last mesh levels, the rates are approaching the expected value r=2. We can see again that the asymptotic convergence behavior of the error is the same for both values of κ. Example 3: $Q\subset {\mathbb{R}}^3$, p=1. In this example, the problem is considered on $Q=\mathrm{\Omega }\times (0,1)$ with $\mathrm{\Omega }=(0,1{)}^2$. The exact solution is given by the formula (57) $u(x,y,t)=(\cos (2\pi (x-y))-\cos (2\pi (x+y)\left)\right)\sin \left(2\pi t\right).$(57) Note that u=0 on Σ and $u_0=0$. The function f is determined by (57(57) $u(x,y,t)=(\cos (2\pi (x-y))-\cos (2\pi (x+y)\left)\right)\sin \left(2\pi t\right).$(57) ). In Figure 4, we plot the contours of u for t=0.8. The problem has been solved on a sequence of meshes, as in the previous tests, using linear polynomial space, i.e. p=1. We perform similar computations as before. In the second column in Table 3, we can see the convergence rates for $\kappa =1$. The last column in Table 3 shows the rates for $\kappa =0.005$. In both cases, the rates are optimal for linear polynomial spaces and are in agreement with the theoretical results in Theorem 4.2.

Figure 3. Example 2: The solution u on Q.

Figure 4. Example 3: The solution u on $Q\subset {\mathbb{R}}^3$.

Table 2. Example 2: The convergence rates r.

$h_s$	$p=2,\kappa =1$	$p=2,\kappa =0.005$
$h_0/2^s$	Convergence rates r
s=1	1.73	1.83
s=2	1.98	2.17
s=3	2.01	2.04
s=4	2.00	2.05
s=5	2.00	2.00
s=6	2.00	2.00
s=7	2.00	2.00

Table 3. Example 3: The convergence rates r.

$h_s$	$p=1,\kappa =1$	$p=1,\kappa =0.005$
$h_0/2^s$	Convergence rates r
s=1	1.09	0.60
s=2	0.90	1.06
s=3	0.92	1.16
s=4	0.97	1.15
s=5	1.05	1.10

Finally, we can conclude that the proposed bubble stabilization finite element scheme performs well for all the examples. The produced numerical solution gives optimal order of convergence in the $\parallel \cdot {\parallel }_{h,\ast }$-norm, when the problems have smooth solutions.

6. Conclusions

In this article, we have proposed and analyzed a bubble stabilized space-time finite element method for solving linear parabolic evolution problems. The construction of the method was based on a space-time variational formulation of the initial PDE problem, which allows the unified space-time discretization by finite element techniques. We presented a discretization error analysis and proved that the method has optimal convergence properties, when $\theta \approx h$ and the PDE problem has a smooth solution. The convergence properties are not affected by the choice of the value of the diffusion parameter κ, which appears in the PDE problem. The theoretical results have been verified by performing several numerical examples. A possible extension of our work here is to combine it with time or space-time mesh adaptivity techniques.

Acknowledgements

The author wish to thank Andreas Schafelner for the further validation of the results in Section 5.

Disclosure statement

No potential conflict of interest was reported by the author.

Additional information

Funding

This work was supported by the Austrian Science Fund (FWF) [grant number NFN S117-03].