Technische Universität Graz Space–time finite element methods for parabolic problems

We propose and analyze a space–time finite element method for the numerical solution of parabolic evolution equations. This approach allows the use of general and unstructured space–time finite elements which do not require any tensor product structure. The stability of the numerical scheme is based on a stability condition which holds for standard finite element spaces. We also provide related a priori error estimates which are confirmed by numerical experiments.


Introduction
Standard approaches to the numerical solution of parabolic evolution problems are based on semi-discretizations, either in space resulting in the method of lines [20], or in time yielding the Rothe method [11]. In both cases, the underlying structure is some tensor product of space and time variables. The sequential solution of the resulting linear systems complicates parallel solution algorithms, and adaptive refinement strategies which are done simultaneously in space and time. Due to the increasing computing capacities, an overall solution of a space-time discretized problem may become an alternative to think about. Applications of space-time finite element methods have already a long history [10]. However, equivalent parabolic operator equations were first analyzed in [17], see also [12,21] for a discussion of the involved stability condition. Note that in all cases, spacetime tensor product wavelet bases are used. Rather general space-time finite elements are already used in many applications, and also in the context of adaptive space-time mesh generation algorithms, see, e.g., [2,13,14,15,19], and the references given therein. Usually, discontinuous Galerkin methods are used, at least for time discretizations, which increase the number of degrees of freedom, and which involve certain parameters to be chosen to ensure stability of the numerical scheme. Therefore we are interested in the application of a standard finite element approach for the solution of evolution problems.
We assume the compatibility condition u 0 (0) = u 0 (1) = 0, and c H > 0 is a given heat capacity constant. It is well known that the solution of (1.1) allows the series representation u k e −(kπ) 2 t/c H sin kπx, u k = 2 The triangular space-time finite elements q ℓ are assumed to be shape-regular, see also Fig. 1. The related set of nodes is denoted by {(x k , t k )} M k=1 . With respect to the triangulation Q h we now define the piecewise linear nodal interpolation where ϕ k (x, t) are the usual two-dimensional continuous and piecewise linear basis functions. Since the solution (1.2) of the initial boundary value problem (1.1) turns out to be smooth, we can apply standard interpolation error estimates [4,5,18] to obtain (1.4) Since the solution (1.2) is known, we can compute the interpolation error (1.4) for a sequence of nested finite element spaces. Hence we can estimate the order of convergence (eoc) which confirms second order convergence, see Table 1. Starting from a coarse finite element mesh we can also drive an adaptive algorithm where we refine only those finite elements q ℓ where the local error satisfies To end up with an admissibile triangulation, we apply a red-green-blue refinement strategy [22]. To obtain a comparable error, the number of nodes, i.e., the number of degrees of freedom, is approximately a forth as for the case when using a uniform refinement strategy. The final mesh is given in Fig. 1. We observe that this adaptive mesh is completely unstructured, both in space and time. In Table 1 we present the interpolation errors for a sequence of space-time finite element meshes up to the seventh refinement level L. Obviously, second order convergence is confirmed. While in the case of a given function, i.e. of the given solution (1.2), we can easily determine the piecewise linear interpolations, we now may ask to find the optimal space-time finite element approximation by solving the initial boundary value problem (1.1) by using a space-time finite element method. In this paper we present a stability and a priori error analysis for a standard Galerkin finite element method for the numerical solution of parabolic evolution equations. In Sect. 2 we describe the model problem of a linear heat equation with Dirichlet boundary conditions, its Galerkin-Petrov variational formulation, and we prove a related stability condition which ensures unique solvability. The finite element discretisation is described and analysed in Sect. 3  stability of the proposed scheme, and in Theorem 3.3 on the a priori error estimate. For the particular case of a spatial one-dimensional problem we present a more detailed analysis on the interpolation error in anisotropic norms in Sect. 4. Numerical examples given in Sect. 5 confirm all theoretical results, and we also comment on a comparison with a space-time discontinuous Galerkin finite element approach as considered in [14]. We end this paper by some concluding remarks in Sect. 6, where we point out some open problems in the a priori error analysis, and we comment on challenging topics concerning a posteriori error estimation techniques and adaptive mesh refinement strategies, efficient preconditioned and parallel solution methods, space-time boundary element methods and their coupling with space-time finite element methods, and on applications including nonlinearities.
Since the Galerkin-Petrov variational formulation (2.2) is based on the use of different test and ansatz spaces, no concept of ellipticity is applicable to ensure unique solvability of (2.2). Instead, a suitable stability condition has to be established. For this we will make use of the trivial inclusion . (2.12) Proof. For any given u ∈ L 2 (0, T ; Due to the assumption u(x, 0) = 0 for x ∈ Ω we further conclude For the second term we obtain, by using (2.7), the Cauchy-Schwarz inequality, and (2.6), we finally conclude from which the stability condition (2.12) follows.
Remark 2.1 Instead of the energy norm (2.10) we may also use the alternative norm [21] u 2 Then we can prove the improved stability estimate However, to avoid the additional term, here we will use the energy norm (2.10) and the stability estimate (2.12).
Since the initial condition in (2.1) is seen as a Dirichlet condition in the space-time domain is satisfied for all v ∈ Y , where we have used the function spaces and Y := L 2 (0, T ; H 1 0 (Ω)). Obviously, we have X ⊂ Y . Now we are in a position to state unique solvability of the Galerkin-Petrov variational formulation (2.15), see, e.g., [3] or [18,Theorem 3.7].

Finite element discretisation
Let X h ⊂ X and Y h ⊂ Y be some finite dimensional spaces which will be specified later.
Before we discuss the unique solvability of the variational formulation (3.1) as well as related a priori error estimates, let us first consider the Galerkin discretization of the variational formulation (2.4) is satisfied for all v h ∈ Y h . By using standard arguments we conclude unique solvability of (3.2), in particular we may define the approximate Newton potential Instead of the energy norm (2.10) we will use the mesh dependent energy norm As in the continuous case we are now able to prove a discrete stability condition.
Then there holds the discrete stability condition Proof. The proof of Theorem 3.1 follows as in the continuous case, see Theorem 2.1.
Moreover, we have, by using u h (0) = 0, as well as which implies the discrete stability condition (3.6).
The discrete stability condition (3.6) implies unique solvability of the Galerkin-Petrov variational formulation (3.1). By combining the Galerkin variational formulation (3.1) with (2.15) and by using the inclusion Y h ⊂ Y we also conclude the Galerkin orthogonality With this we can derive the following quasi-optimal error estimate.
Theorem 3.2 Let u ∈ X and u h ∈ X h be the unique solutions of the variational formulations (2.15) and (3.1), respectively. Then there holds the a priori error estimate

8)
Proof. From the stability condition (3.6), by using the Galerkin orthogonality (3.7) and (2.9), we conclude, for arbitrary z h ∈ X h , and the assertion follows from the triangle inequality.
Now we are in a position to introduce the space-time finite element spaces X h ⊂ Y h . For the space-time domain Q = Ω × (0, T ) ⊂ R n+1 we consider a sequence of admissible decompositions Q h into shape regular simplicial finite elements q ℓ which are triangles (n = 1), tetrahedra (n = 2), or pentatops (n = 3, see, e.g., [15]), i.e.
we denote the related set of nodes (x k , t k ) ∈ R n+1 . As usual we introdue a reference element q ⊂ R n+1 which allows to describe a finite element q ℓ by using a local parametrization, where q ⊂ R n+1 are the usual reference elements. Hence we can introduce the volume ∆ ℓ of the finite element q ℓ , and the local mesh width for n = 1, 1 6 for n = 2, 1 24 for n = 3. With respect to the admissible decomposition Q h we define the space-time finite element space S 1 h (Q h ) = span{ϕ k } M k=1 of piecewise linear and continuous basis functions ϕ k . Now we are able to introduce the finite element spaces which obviously satisfy the assumptions of Theorem 3.1. Based on Theorem 3.2 we conclude the following a priori error estimate.
Theorem 3.3 Let u ∈ X and u h ∈ X h = S 1 h (Q h ) ∩ X be the unique solutions of the variational formulations (2.15) and (3.1), respectively. Assume u ∈ H 2 (Q). Then there holds the energy error estimate Proof. From the definition of the related norms and by using Theorem 3.1 we first have and it remains to define a suitable approximation z h of u to prove a related approximation property. In fact, let z h = P h u ∈ X h be the unique solution of the Galerkin variational formulation for all v h ∈ X h . When using the related energy norm, note that v(x, 0) = 0 for x ∈ Ω and v(x, t) = 0 for x ∈ Γ, i.e.
and by using standard arguments, i.e. Cea's lemma and standard approximation error estimates, we obtain Note that for n ≤ 2 we can use standard nodal interpolation error estimates, while for n = 3 we have to use projection type error estimates. Moreover, By taking into account the norm definition (2.10) the assertion follows from the previous estimates.
Corollary 3.4 Let u ∈ X and u h ∈ X h = S 1 h (Q h ) ∩ X be the unique solutions of the variational formulations (2.15) and (3.1), respectively. We now assume u ∈ H s (Q) for some s ∈ [1,2]. Then there holds the energy error estimate Proof. The assertion follows from the trivial error estimate, i.e. s = 1, , by using the energy error estimate (3.9), i.e. s = 2, and an interpolation argument for s ∈ (1, 2).

Space-time interpolation error estimates
As we will see in the numerical results, in particular of Example 5.1, the error estimates (3.9) and (3.11) are optimal in the case of a regular solution. However, in the particular case of solutions with some singular behavior, see Examples 5.2 and 5.3, the error estimate (3.10) turns out to be not optimal, or not applicable. This is mainly due to the unified handling of the spatial and temporal derivatives when proving the error estimate (3.9). In what follows we present a more detailed analysis of the interpolation error in the energy norm. To simplify the presentation, let us consider the spatial one-dimensional case only, where all space-time finite elements q ℓ are a right-angled isosceles triangles, as used in Fig. 1. Without loss of generality we consider the space-time finite elements q 1 h and q 2 h which are given by and The related linear interpolations I 1 h u and I 2 h u of a given function u then read and For the local interpolation error we then can prove the following result.  ((0, h) 2 ), and let I h u be the piecewise linear interpolation. Then there holds the error estimate Proof. To estimate the local error Hence, by using the Cauchy-Schwarz inequality twice, follows. In the same way we have and therefore follows. Summing up gives the desired result.
As an immediate consequence of Lemma 4.1 we can formulate the following interpolation error estimate. Although the interpolation error estimate (4.2) is not sufficient to establish an improved error estimate for the finite element solution u h of the Galerkin-Petrov variational formulation (3.1), it already indicates the expected order of convergence. However, since the second error term, i.e. ∂ t (u − z h ) L 2 (0,T ;H −1 (Ω)) , is measured in a Sobolev norm with a negative index, no interpolation error estimates are available. Instead, appropriate projection operators should be used to gain an additional order of convergence due to the Aubin-Nitsche trick. But since we have to use the same element z h ∈ X h within the a priori error estimate (3.9) we have to bound the projection error in different norms by the interpolation error.
Although we have shown the improved interpolation error estimate (4.2) only in the particular case of right-angled isosceles triangular elements, we claim that this error estimate remains true for more general situations, i.e. for unstructured meshes as used in Example 5.2, and for higher spatial dimensions. Since the related proofs are much more technical, we do not consider them here.
To close this section, let us recall that the error estimate (3.10) already covers the most general situation when assuming some regularity of the solution in the space-time domain which can be ensured in many cases.

Numerical results
In what follows we provide some numerical results to confirm the theoretical estimates as given in Sect. 3 and 4. As in [14] we consider the initial boundary value problem (2.1) for the spatial domain Ω = (0, 1), the final time T = 1, the heat capacity c H = 1, and A(x, t) = I. The discretization is done by considering a uniform and structured triangulation of the space-time domain Q = (0, 1) 2 as shown in Fig. 2 (left), and by using a more unstructured mesh as shown in Fig. 2 (right), and by using either first (p = 1) or second (p = 2) order continuous finite elements. Note that the structured mesh corresponds to the situation as considered in Lemma 4.1. The linear system is solved by a GMRES without preconditioning where it is sufficient to consider Y h = X h to build up the linear system, i.e. we neglect the basis functions with nonzero initial conditions.  u(x, t) = cos πt sin πx for (x, t) ∈ Q and we determine the initial datum u 0 and the given volume density f accordingly. Since the solution u is smooth we expect first and second order convergence when using piecewise linear (p = 1) and piecewise quadratic (p = 2) continuous finite element approximations, see the error estimate (3.11). This is confirmed by the numerical results as given in Table 2.  Note that the absolute errors and the estimated orders of convergence (eoc) almost coincide with the numerical results as given in [14, Table 2.1] in the case of a discontinuous Galerkin finite element approach where a related discrete energy norm is used which is not equivalent to the mesh dependent norm (3.4). It is worth to mention that there is a reduced order of convergence when using piecewise quadratic but continuous finite element approximations within the DG formulation [14, Table 2.7] which does not appear within our approach. This is probably due to the chosen energy norm of the DG approach which may not reflect the underlying function spaces in a proper way. While the results of Table 2 correspond to a globally uniform mesh, see Fig. 2 (left), we may also use a more unstructured mesh of the space-time domain as shown in Fig. 2 (right). The related results are shown in Table 3, and again they confirm the theoretical findings.
which has a line singularity at t = 1, i.e. u ∈ H 5/4−ε (Q) for any ε > 0. According to the error estimate (3.10) one would only expect an order of convergence of almost 0.25. However, due to the improved interpolation error estimate (4.2) we can expect linear convergence which is confirmed by the numerical results given in Table 4 in the case of the structured mesh, and in Table 5 when using the more unstructured mesh. Note that when using second order elements (p = 2) we observe less than linear convergence, but the order of convergence is much better than expected from the error estimate (3.11). The understanding of this behavior requires a more detailed analysis to be done.
For the numerical results in the case of a piecewise linear continuous approximation see Table 6. We observe an almost linear order of convergence. While in the case of discontinuous piecewise linear approximations a convergence order of 0.5 is observed, see [14, Table 6: Numerical results in the case of a singular solution (α = 0.5), structured mesh.

Concluding remarks
In this paper we have provided the stability and a priori error analysis of conforming spacetime finite element methods for the solution of parabolic evolution equations. In contrast to other space-time approaches, which rely on space-time tensor product ansatz spaces, see, e.g., [17,21], our approach is based on the use of arbitrary admissible finite element meshes of the space-time domain. Under rather mild assumptions on the regularity of the solution we have derived quasi-optimal error estimates which do not depend on the space dimension. In particular situations of non-smooth input data, less regular solutions may appear, where the above mentioned a priori error estimates may be either not optimal, or can not be applied. To get a better understanding of these phenomena we have presented a more detailed analysis of the piecewise linear interpolation error in the particular situation of right-angled isosceles triangular elements. It is obvious that additional work is required to improve these results, in particular to handle more general cases. To verify the theoretical results we have provided simple numerical examples for spatial one-dimensional domains. However, this paper extends previous work [14,15] where space-time discontinuous Galerkin finite element methods were used even for three-dimensional spatial domains. In fact, using the results of [15] it will be possible to realize the proposed approach using continuous space-time finite elements also for three-dimensional spatial domains, this will be one topic of future work. While in the present implementation of the continuous space-time finite element method we have used the GMRES algorithm without any preconditioner as an iterative solver, appropriate preconditioners and a parallel implementation are mandatory. Note that for the space-time discontinuous Galerkin approach related multigrid methods are analyzed and successfully implemented on state of the art parallel computer architectures showing optimal scaling properties [8,9]. Although the multigrid approach may be also applied to continuous space-time formulations, probably there are other preconditioning strategies such as geometric and algebraic multilevel techniques, or domain decomposition methods applicable. Obviously, this will be a second topic of future research.
The main advantage of the proposed space-time finite element approach is the possibility to handle general finite element meshes wich are adaptive simultaneously in space and time. To drive such an adaptive mesh refinement, appropriate a posteriori error estimators are required which are well established for elliptic boundary value problems, see, e.g., [22].
In the context of time-dependent problems, a posteriori error estimators are available as well, see, e.g., [1,7,16], but in most cases, adaptivity is considered either in space, or in time. So the design and analysis of appropriate space-time a posteriori error estimators is a third topic of future work.
The proposed strategy of rather general space-time approximations can be applied also to the numerical solution of boundary integral equations for parabolic partial differential equations. Then, the boundary Σ = Γ × (0, T ) of the space-time domain Q = Ω × (0, T ) has to be decomposed into general space-time boundary elements. Note that the theory of boundary integral equations in particular for the heat equation is well established [6]. Although space-time tensor products are used in [6] for the approximation, the theory already covers the general case of arbitrary space-time boundary elements. However, the implementation of general space-time boundary element methods is more involved due to the fundamental solution which is non-local. In any case, space-time boundary element methods and their coupling with space-time finite element methods will be interesting topics for future work as well.
Applications of the proposed space-time finite element approach are numerous, not only to handle evolution equations of parabolic type, but also hyperbolic equations with the wave equation as the most simple example. While space-time methods can be applied to nonlinear partial differential equations as well, see, e.g., [14] for a space-time DG approach for the Navier-Stokes equations, problems with moving interfaces such as in fluid-structure interaction problems seem to be more challenging. In particular, when the interface is unknown, even a coupled system of linear partial differential equations becomes nonlinear in the geometric description. It is obvious that there are much more challenging and interesting applications around, and we believe that space-time discretizations of partial differential equations can turn to be a method of choice.