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.


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][2][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 , . . . , 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 ∂ t u − κ 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][23][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.

Preliminaries
Let be a bounded Lipschitz domain in R d , d ≥ 1. Let α = (α 1 , . . . , α d ) be a multi-index of non-negative integers α 1 , . . . , α d with degree |α| = d j=1 α j . For any α, we define the differential operator D α As usual, L 2 ( ) denotes the Sobolev space for which |φ(x)| 2 dx < ∞, endowed with the norm φ L 2 ( ) = ( |φ(x)| 2 dx) 1/2 . Let be a non-negative integer, define the standard Sobolev spaces endowed with the following norms where ∇ x φ is the gradient with respect to the spatial variables. Similarly, we denote by n = (n x , n t ) the normal component on ∂Q, with n x the components related to space direction and n t the component related to time direction. Let , m be positive integers, for functions defined in Q, we define the Sobolev spaces and the subspaces For a function φ ∈ H ,m (Q) with , m ≥ 1, we define the norms and the seminorms We recall Hölder's and Young's inequalities that hold for all u ∈ L 2 (Q) and v ∈ L 2 (Q) and for any fixed ∈ (0, ∞).
We will use the following Poincare's inequality: Let ⊂ R d , d = 1, 2, . . . be a bounded rectangular domain and let ⊂ ∂ with | | > 0. For simplicity we assume that lies on the plane with x 1 = 0.
Let v ∈ C ∞ ( ) and v(x ) = 0 for all x ∈ . For any interior point The first inequality in (3) yields where the constant C depends on the length of . Integrating (5) over , we can obtain We refer to [27] for a proof of (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 ∼ b meaning that c a ≤ b ≤ C a, with c, C generic constants.

The model parabolic problem
In the space-time cylinderQ =¯ × [0, T], with ⊂ R d , d = 1, 2, 3, we consider the initial boundary value problem u t − κ u = f in Q and u = 0 on , u(·, 0) = u 0 on 0 , as a model problem, where the diffusivity parameter 0 < κ ≤ 1 is taken to be constant, f : Q → R, with f ∈ L 2 (Q), and u 0 : → R, with u 0 ∈ L 2 ( ) are given functions, and u :Q → 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): find u ∈ H 1,0 0 (Q) such that The variational problem (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) belongs to V 0,0 : 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.

A usual finite element scheme
Let T h (Q) be a partition of space-time domain Q into triangular (or quadrilateral elements), that is Q = ∪ E∈T h E, see Figure 1(a). We denote by h E the diameter of E ∈ T h (Q) and the mesh size is defined , consisting of continuous functions in space and in time, by where P p (E) is the polynomial space of total degree p, see e.g. [31][32][33]. In our analysis we focus on the case p = 1.
The usual finite element approximation of (7) would read: find u h ∈ V h0 such that Setting v h = u h in (12a) and using the identity , and by inequalities (3) and (6), we can obtain the following estimate where C is the constant that appears in (6). In the stability estimate (13) the control of ∇ x u h can be very poor if the parameter κ is very small. Furthermore (13) does not provide a direct control of ∂ t u h L 2 (Q) . Thus in the case where κ is small, the finite element scheme (12) may not perform well. Therefore it is crucial to improve the stability properties of (12). Next we present a technique for stabilizing the finite element scheme (12) 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.

The stabilized scheme
We introduce the larger finite subspace V h,b of H 1,1 0,0 (Q) that can be written as a direct sum as follows with 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 ∈ T h (Q). For example, in case of triangular elements, it is spanned by a cubic functions 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 In view of this, we introduce the discrete problem: where a(·, ·) as in (12b) and with θ > 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].
For convenience, we introduce the discrete bilinear form and the associated mesh dependent norm Note that, the terms with the time derivatives in (19) are related to the bubble function space.

Lemma 3.2:
The discrete form a h (·, ·) : The definition of a h and (21) imply which is (20) with C s = 1 and this completes the proof.

Proposition 3.1:
Let u h be the solution given by (15). Then there is a C κ, > 0 such that the solution u h satisfies the following a priori estimate Proof: Using u h ∈ V h,b as a test function in (15), and utilizing (3) and (20)) together with the Poincare inequality (6), we successively obtain where we have previously used that C s = 1 and κ ≤ 1. Setting C κ, = C (1/κ) we get (23).
Note that the estimate in (23) provides a direct control of ∂ t u b h 2 L 2 (Q) due to the appearance of b(·, ·) in the finite element scheme (15), cf. (13). A direct result of (23) and (20) is the following corollary.

Corollary 3.1:
The discrete problem defined in (15) is well posed, i.e. it has a unique solution which satisfies the stability estimate (23).
We point out that the · h, * is defined for all v ∈ V 0,0 and is similar to the · h given in (19), which is defined for all v h ∈ V h,b . The mesh-dependent norm · h,V will be used later for deriving bounds for the discretization error.

Lemma 3.3:
There is a constant C b (κ, θ, h) > 0 such that We treat every term of the form a(·, ·) separately. We apply integration by parts and (3) to infer where c 3 depends on the constants appearing in (16) and (6). Similarly for the second term, applying (3), we get Combining all the above bounds and setting C b (κ, θ, h) = 2c 3 (θ (κh) −1 + 1) 1/2 , we can derive the desired result.

Error analysis
Proof: Multiplying (7) by z h ∈ V h,b , integrating over Q, and then applying integration by parts, we arrive at the variational identity We recall the problem (15) and directly have Proof: Let z 1 h be a function in V h0 . By (29) and by subtracting similar terms, we have that Setting above φ h = u 1 h + u b h − z 1 h and applying integration by parts on the first term on the left side, we have and by applying (3) and (4) on the right hand side, yields Gathering the same terms and setting 0 < c < 1 2 , we get Using 0 < κ ≤ 1, applying triangle inequality and setting c = (1/(1 − 2 c ))(1/c ), the assertion follows.
Below we give the main error bound for the finite element solution u h ∈ V h,b . (15). Under Assumption 2.1 and choosing θ ≥ h, there exist a constant c * ,V , depending on c inv in (16), such that
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 ≤ m ≤ 1 < s and let the space V h0 defined in (11).

Lemma 4.3:
Let the space V h0 defined in (11). Let s ≥ 2 be an integer, and let a function v ∈ V 0,0 with where r = 2 and c 1 , c 2 depend on the constants appearing in (17) and in (48), but not on h and v.
Proof: Introducing the operator π h v of Lemma 4.2 and by applying (17) and (48), we have In the same way, we have Collecting the estimates (50) and (51), we easily obtain which is (49b).

Remark 4.2:
The interpolation estimates given in (49) have been derived for linear polynomial spaces, see (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), 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].
and moreover for any θ > 0, with c 1 depending on the constants in (37) and in (49) and c 2 on the constants in (32) and (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.

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 ⊂ 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 ⊂ 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 = (h 0 /2 s ), s = 1, 2, . . .. Every T h s (Q) 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(e s /e s+1 )/ln(h s /h s+1 ), where the error e s := u − u h h, * is computed on T h s (Q). 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) and Remark 4.2. We consider test cases with κ ∈ {1, 0.005}, and we study the behavior of the rates for θ = 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 ⊂ R 2 , p = 1. In the first example, the problem is considered in Q = (0, 1) × (0, 2). The exact solution is given by the formula u(x, t) = sin(2πx) sin(π t).
The source function f is determined by (55). Note that u = 0 on and u 0 = 0, see (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 κ = 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). The numerical solution u h ∈ 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 κ = 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. Example 2: Q ⊂ R 2 , p = 2. In the second example, we consider the problem on Q = (0, 1) × (0, 1). The exact solution is given by the formula u(x, t) = sin(2πt) sin(2πx).
The source function f is defined to match the solution in (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 κ = 1 and θ = 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) and Remark 4.2. We repeat the same computations setting κ = 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 ⊂ R 3 , p = 1. In this example, the problem is considered on Q = × (0, 1) with = (0, 1) 2 . The exact solution is given by the formula u(x, y, t) = (cos(2π(x − y)) − cos(2π(x + y))) sin(2πt).
Note that u = 0 on and u 0 = 0. The function f is determined by (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 κ = 1. The last column in Table 3 shows the rates for κ = 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.  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 · h, * -norm, when the problems have smooth solutions.

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 θ ≈ 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.