SUPG space-time scheme on anisotropic meshes for general parabolic equations

Ioannis Toulopoulos

doi:10.1515/jnma-2024-0053

Artikel Open Access

SUPG space-time scheme on anisotropic meshes for general parabolic equations

Ioannis Toulopoulos

Veröffentlicht/Copyright: 4. Juni 2025

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Manuskript einreichen Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Journal of Numerical Mathematics Band 33 Heft 4

Abstract

In this paper, a stabilized space-time finite element scheme on anisotropic quadrilateral meshes for general linear parabolic problems is considered. The scheme is devised on the basis of a unified space-time variational formulation, and uses continuous piece-wise polynomial spaces. The stabilization is achieved by incorporating Streamline-Upwind Petrov–Galerkin (SUPG) techniques. Defining appropriately the stability coefficients, we initially show anisotropic interpolation estimates and then a priori error estimates by following a classical finite element methodology. A series of numerical examples illustrates the theoretical findings.

Keywords: general linear parabolic equations; continuous space-time finite element discretizations; quadrilateral anisotropic meshes; a priori error estimates; asymptotic convergence rates

MSC 2010: 65N30; 65J15; 65N15; 65N12

1 Introduction

Many transport physical problems are often described by using general second order parabolic equations of the form u _t + Lu = f, where Lu ≔ −div(ɛ∇_x u) + β ⋅ ∇_x u + ru is the second order differential operator, and ∇_x u is the spatial gradient of u, β is a constant vector representing the convection velocity and the parameters ɛ > 0, r ⩾ 0 represent the diffusion and reaction coefficients, respectively, [1]. The numerical solution of these problems has been a subject of investigation by many authors in the past decades (cf. [2]). Standard Galerkin Finite Element Methods (FEMs) with continuous spaces may appear numerical instabilities, which are produced due to the advection character of the problem, (advection dominant case). Very often Streamline-Upwind (SU) stabilization terms are added to treat the problem numerically and to ensure the stability of the FE discretization, see, e.g., [3], [4], and in refs. [2], [5], [6] for an overview and computational results for SU methods. The full discretization of the problem is completed by applying a time-stepping scheme, e.g., Runge–Kutta, which results to sequential approximations of the solution in time (see, e.g., [7], [8], [9]). These approaches typically impose a restriction on the time step relative to the spatial mesh size, which can lead to additional difficulties when highly refined meshes are required.

In contrast to these methods, the last proposed space-time finite element methods (STFEMs) discretize time evolution problems by applying a unified and simultaneous finite element discretization in space and in time directions, [10]. The main idea is to see the time variable t as another spatial variable, let’s say, x d x + 1 , if x 1 , … , x d x , are the spatial variables, and the time derivative u _t as a convection in the direction x d x + 1 . In view of this, an associated global space-time variational formulation is derived and the time-dependent problem is considered as a stationary problem into a domain (i.e., the space-time cylinder) with one higher dimension. This idea is not new and has been used by many authors for developing different STFE discretizations for time evolution problems. Specifically, in the past few years the proposed STFEMs are mainly based on a corresponding space time variational formulation, see, e.g., [11], [12], [13], [14], and [15], [16] for alternative stabilization approaches, also [17] for nonlinear problems. In refs. [18], [19], space-time methods the Isogeometric Analysis framework are discussed and see the survey paper [20] for applications to engineering problems. The proposed STFEMs offer some further advantages compared to the traditional discretization methods, as for example, the direct use of existing FEM solvers developed for stationary problems, the generalization of adaptivity techniques which have been studied for stationary problems, the use of coarse/fine unstructured space-time meshes in order to compute highly accurate and efficient solutions without having time step limitations, and lastly are more flexible to be implemented in parallel environments and to perform computations in such architectures. We refer to refs. [21], [22] for a discussion on different parallelization approaches for STFEMs. However, many of these techniques do not directly carry over when solving problems in four-dimensional (4D) space-time domains. The construction of unstructured meshes in 4D space-time domains seems to be a challenging task, [23].

In this paper, inspired by the ideas presented in ref. [24], where SU stabilized finite element methods for steady advection–diffusion problems are analysed, we devise a stable STFEM for solving the general parabolic problem mentioned above. The method is considered on anisotropic quadrilateral meshes, which are aligned to the coordinate axes. The aim is to derive error estimates in the related energy norm uniformly with respect to the diffusion parameter ɛ, by taking into account the anisotropic mesh sizes of the triangulation of the space-time cylinder. Animated by the interpolation results given in ref. [25], anisotropic interpolation error estimates are derived here, where the associated interpolation constants depend on the directional stretching properties of the mesh. The additional SU stability terms appearing in the final space-time scheme are weighted by a numerical parameter, see (3.18), which is accordingly formed by the anisotropic character of the mesh. This parameter is determined through the stability of the resulting bilinear form and includes the local, i.e., per element, spatial mesh sizes. The numerical results confirm the theoretical findings. It is known that the solution of the aforementioned general parabolic problem can exhibit interior boundary layers. In practical applications, the resolution of this layer is of main interest and can be typically treated by applying an anisotropic meshing technique, [26], [27]. This work is the first step to devise STFEMs in this direction. Extensions of the proposed work to more general fluid flow problems with boundary layers are under preparation.

The outline of the paper can be stated as follows. In Section 2, some preliminaries together with the notation of the related Sobolev spaces are given. In Section 3, the general parabolic problem is given, and the weak space-time formulation is described. In the last part of Section 3 the ST-FE discretization in presented and the discretization error analysis is developed. Finally, in Section 4 we show a series of numerical examples for verifying the theoretical results. The paper closes with the conclusions.

2 Preliminaries

2.1 Notations

Let Q be a bounded Lipschitz domain in R d , d = 2, …, 4, with boundary Γ = ∂ Q . For any multi-index α = (α ₁, …, α _d) of non-negative integers α ₁, …, α _d, we use the following notations, (i) | a | = ∑ i = 1 d a i , (ii) set x a = x 1 a 1 x 2 a 2 … x d a d for x ∈ R d , (iii) moreover introduce the differential operator D α ≔ ∂ x a = ∂ x 1 α 1 … ∂ x d α d , with ∂ x j ( ⋅ ) = ∂ ( ⋅ ) / ∂ x j , j = 1, …, d. Let 1 ⩽ p ⩽ ∞ be fixed and ℓ be a non-negative integer. As usual, L p ( Q ) denotes the Lebesgue spaces for which ∫ Q | φ ( x ) | p d x < ∞ , endowed with the norm  ‖ φ ‖ L p ( Q ) = ∫ Q | φ ( x ) | p d x 1 / p , and W ℓ , p ( Q ) is the Sobolev space, which consists of the functions φ : Q → R such that their weak derivatives D ^α φ with | α | ⩽ ℓ exist and belong to L p ( Q ) . If φ ∈ W ℓ , p ( Q ) , then its norm is defined by

‖ φ ‖ W ℓ , p ( Q ) = ∑ 0 ⩽ | α | ⩽ ℓ ‖ D a φ ‖ L p ( Q ) p 1 / p , ‖ φ ‖ W ℓ , ∞ ( Q ) = max 0 ⩽ | α | ⩽ ℓ ‖ D a φ ‖ ∞

for 1 ⩽ p < ∞ and p = ∞, respectively. We further define the spaces

(2.1a) W 0 ℓ , p ( Q ) ≔ φ ∈ W ℓ , p ( Q ) such that φ | ∂ Q = 0 ,

(2.1b) W Σ * ℓ , p ( Q ) ≔ φ ∈ W ℓ , p ( Q ) such that φ | Σ * ⊂ ∂ Q = 0 .

Remark 2.1.

If p = 2 we usually use the notation H ℓ ( Q ) ≔ W ℓ , 2 ( Q ) , ℓ = 0,1,2 , …

Let the fixed integer ℓ ⩾ 1 and let the set A of all multi-indices having the form A ≔ {a:a = (…, m _i e _i, …)} where 0 ⩽ m _i ⩽ ℓ, i = 1, …, d, and e _i is the unit vector in the i-th direction. For φ ∈ W ℓ , p ( Q ) denote

(2.2) ∑ a ∈ A ‖ D a φ ‖ L p ( Q ) = ∑ i = 1 d ‖ ∂ x i m i φ ‖ L p ( Q ) .

We refer the reader to ref. [28] for more details about Sobolev spaces.

2.1.1 Differential operators on the space-time domain

Next, we define certain differential operators which are related to the time and the spatial variables. Let J = (0, T] be the time interval with some final time T > 0 and let Ω be a bounded domain in R d x , d x = 1,2 or 3. For later use, we consider the space-time cylinder Q ⊂ R d with d = d _x + 1, defined by Q = J × Ω, and its boundary parts Σ = ∂Ω × J, Σ_T = Ω × {T}, and Σ₀ = Ω × {0} such that ∂ Q = Σ ∪ Σ ̄ 0 ∪ Σ ̄ T . Accordingly to the definition of ∂ x α , we now define the operator ∂ x α d x and also define the spatial gradient ∇ x φ = ( ∂ x 1 φ , … , ∂ x d x φ ) , and the whole gradient ∇φ ≔ (∂_t φ, ∇_x φ).

2.2 Known inequalities and identities

The following inequalities are going to be used in several places in the text. Hölder’s and Young’s inequalities read: For any δ, 0 < δ < ∞, and 1 ⩽ p, q ⩽ ∞ such that 1 p + 1 q = 1 , for f ∈ L p ( Q ) and g ∈ L q ( Q ) , there holds

(2.3a) ∫ Q f g d x ⩽ ‖ f ‖ L p ( Q ) ‖ g ‖ L q ( Q ) ⩽ δ p ‖ f ‖ L p ( Q ) p + δ − q / p q ‖ g ‖ L q ( Q ) q .

Poincaré–Friedrichs inequality, see [28], [29]: Let Q ⊂ R d be a parallelepiped (cuboid) and let the face Σ * ⊂ ∂ Q vertical to the x _j, 1 ⩽ j ⩽ d, coordinate plane. Then for any f ∈ W Σ * 1 , p ( Q ) , it holds

(2.3b) ∫ Q | f | p d x ⩽ C ( Q ) ∫ Q | ∂ x i f | p d x , 1 ⩽ i ⩽ d .

Let the vector β = (β ₁, …, β _d), the function f ∈ W 1 , p ( Q ) and the outward normal vector n to ∂ Q . In several places we will use the identities:

(2.4a) ∇ ⋅ ( β f ) = β ⋅ ∇ f + ( ∇ ⋅ β ) f ,

(2.4b) 2 ∫ Q β ⋅ ∇ f f d x = − ∫ Q ∇ ⋅ β f 2 d x + 2 ∫ Q β ⋅ n f 2 d s .

3 The continuous problem

Let Ω be a bounded cuboid domain in R d x , with d _x = 1, 2, 3, with smooth boundary Γ = ∂Ω. We define the space-time cylinder Q ̄ T ≔ Ω ̄ × [ 0 , T ] , where T is the final time, and boundary ∂ Q T = Σ ∪ Σ ̄ 0 ∪ Σ ̄ T , where Σ ≔ Γ × (0, T) is the lateral boundary, Σ₀ ≔ Ω × {0} and Σ_T ≔ Ω × {T}. Consider the differential operator L in the form

(3.1) L u ≔ − d i v ( ε ∇ x u ) + β ⋅ ∇ x u + r u ,

where ∇_x u is the spatial gradient of u, and ɛ > 0 is the constant diffusion coefficient. The constant vector β ≔ (β _x, β _y, β _z) takes values in R d x and the reaction parameter r takes values in R + . We consider the following initial-boundary value problem: Find u ( x , t ) : Q ̄ T → R such that

(3.2a) u t + L u = f in Q T ,

(3.2b) u = u Σ = 0 on Σ ,

(3.2c) u ( x , 0 ) = u 0 ( x ) on Ω ,

where f, u ₀ are given functions. For simplicity, we only consider homogeneous Dirichlet boundary conditions on Σ. However, the analysis presented in this work can easily be generalized to other constellations of boundary conditions, cf., [10].

In literature, usually the main point of the concept for defining a weak formulation of (3.2) is to consider that u _t lives in the dual space (or in some sub-space) of the space where u lives, [1]. Anyway, as we mentioned before, in recent years appropriate space-time weak formulations for parabolic problems similar to (3.2) have been presented, where the regularity of the solution is considered uniformly in all the space time cylinder, i.e., u ∈ W ^2,p(Q _T), see, e.g., [10] and the reference therein.

3.1 Weak space-time form

Assume that u 0 ∈ H 0 1 ( Ω ) and f ∈ L ²(Q _T). Multiplying (3.2a) by a smooth function v(x, t) vanishing on Σ ∪ Σ_T, and after applying Green’s theorem and an integration by parts with respect to both t and x, and the usage of the boundary conditions, we can obtain

(3.3) ∫ Q T − u v t + ε ∇ x u ⋅ ∇ x v + β ⋅ ∇ x u v + r u v d x d t = ∫ Q T f v d x d t + ∫ Ω u 0 ( x ) v ( x , 0 ) d x .

Introducing the appropriate regularity properties for the data, global regularity properties can be shown for the generalized solution u of (3.3) in Q _T, i.e., u _t ∈ L ²(0, T; L ²(Ω)), [1], and furthermore, it can be inferred by using embeddings that u ∈ W ^1,p(Q _T), cf., [30], [31]. Applying a formal integration by parts with respect to time variable we can arrive at the space-time weak formulation: Find u ∈ H Σ 1 ( Q T ) with u(0, x) = u ₀(x), such that

(3.4) ∫ Q T u t v + ε ∇ x u ⋅ ∇ x v + β ⋅ ∇ x u v + r u v d x d t = ∫ Q T f v d x d t , v ∈ H Σ 1 ( Q T ) .

Assumption 3.1.

For the solution u of (3.4), assume that u ∈ V, with V = H Σ 1 ( Q T ) ∩ H ℓ ( Q T ) , ℓ ⩾ 2 .

Remark 3.1.

The space-time variational formulation (3.4) has a unique solution, see, e.g., analysis in refs. [30], [31], and also [12] for considerations in Gelfand triple spaces. In these works, beside existence and uniqueness results, one can also find useful a priori estimates and regularity results.

In view of (3.4), we define

(3.5a) B ( u , v ) ≔ ∫ Q T u t v + ε ∇ x u ⋅ ∇ x v + β ⋅ ∇ x u v + r u v d x d t ,

(3.5b) ℓ f ( v ) ≔ ∫ Q T f v d x d t

for v ∈ H Σ 1 ( Q T ) . Working in a formalistic way, and setting u = v in (3.4), and then using (2.4b) and the fact that ∇_x ⋅ β = 0, we can deduce that

(3.6) B ( u , u ) = ∫ Σ T 1 2 u 2 ( s ) d s + 2 ∫ Σ β ⋅ n u 2 d s + ∫ Q T ε | ∇ x u | 2 + r u 2 d x d t ≲ ∫ Q T | f u | d x d t + ∫ Ω u 0 2 ( x ) d x ≲ ( 2.31 ) ( 2.3 b ) c ( δ ) ∫ Q T | f | 2 d x d t + δ ∫ Q T | ∇ x u | p d x d t .

Recalling (3.2b), and choosing sufficiently small, i.e., δ = ɛ/4 in (3.6), we can have the a priory bound

(3.7) ∫ Σ T u 2 ( s ) d s + ϰ ( ε ) ∫ Q T | ∇ x u | 2 d x d t ≲ c ( δ ) ‖ f ‖ L 2 ( Q T ) 2 + ‖ u 0 ‖ L 2 ( Ω ) 2 .

Working in the same spirit for the case β = 0, make test with v ≔ u _t(x, t) in (3.4), and note that d d t ∫ Q T | ∇ x u | 2 = ∫ Q T ∇ x u ⋅ ∇ x u t , then

(3.8) ∫ Q T | u t | 2 d x d t + ε ∫ Σ T | ∇ x u | 2 d x + r ∫ Q T u u t d x d t ≲ c 2 ( δ ) ∫ Q T | f | 2 d x d t + δ ∫ Q T | u t | 2 d x d t + ∫ Ω ε | ∇ x u 0 | 2 d x

provided that ∇_x u ∈ C([0, T], L ²(Ω)). Choosing δ appropriate small, we obtain

(3.9) c 3 ( δ ) ∫ Q T | u t | 2 d x d t + ε ∫ Σ T r 2 u 2 + | ∇ x u | 2 d x ≲ c 2 ( δ ) ∫ Q T | f | 2 d x d t + ∫ Ω r 2 u 0 2 + ε | ∇ x u 0 | 2 d x .

Remark 3.2.

The estimate in (3.8) does not provide any bound for the u _t term. On the other hand the estimate in (3.9) gives a bound for u _t in L ²(Q _T). This aligns with the idea of employing streamline stabilization techniques, as described below, to develop a stable space-time discretization scheme for (3.4). Other ideas for producing stable space-time discretizations have been developed by means of Petrov–Galerkin techniques after an appropriate selection of trial and test space, see, e.g., [12], [13], [32].

3.2 The space-time finite element approximation

3.2.1 Basic concepts

We start by introducing the discrete setting. Let T h ≔ { E l } l = 1 , … , N be a conforming mesh partition, of the space-time cylinder Q _T into closed rectangular mesh elements, such that

(3.10) Q ̄ T = ⋃ l E l , E o , l 1 ∩ E o , l 2 = ∅ , 1 ⩽ l 1 ≠ l 2 ⩽ N ,

where E o , l is the interior of the mesh element. Denote h _E ≔ (h _1,E, …, h _i,E, …, h _d,E) to be the vector with the i-th directional mesh widths of the element E ∈ T h . For the analysis below some further mesh width quantities must be introduced. Define h = (h ₁, …, h _i, …, h _d) ≔ (max_E{h _1,E}, …, max_E{h _i,E}, …, max_E{h _d,E}) and h _min = (h _min,1, …, h _min,i, …, h _min,d) ≔ (min_E{h _1,E}, …, min_E{h _i,E}, …, min_E{h _d,E}). The diameter of every E l ∈ T h is denoted by h E l and we set h ≔ max E l h E l , h _t ≔ h ₁, h _x ≔ max_i⩾2{h _i}, h _E, min ≔ min_i⩾2{h _i,E}. In the sequel we write E ∈ T h instead of E l ∈ T h . Note that the mesh T h forms a Cartesian grid in Q _T and every mesh element E ∈ T h is affine equivalent to reference element E ̂ ≔ [ 0,1 ] d though the transformation

(3.11) T E : E ̂ → E , T E ( x ̂ ) ≔ x ( x ̂ ) = B E x ̂ + b E ,

with b E ∈ R d and B E ∈ R d × d to be a diagonal matrix of the form B E , i i = h i , E , i = 1 , … , d .

Assumption 3.2.

The partition T h is quasi-uniform, in the sense that there are positive constants σ _M, σ _m independent of h such that σ _M ⩾ h _i/h _min,i ⩾ σ _m for all i = 1, …, d.

Remark 3.3.

The quasi-uniformity properties Assumption 3.2 of T h allows the use of anisotropic mesh widths between the different axial directions, and furthermore to replace h _E by h in the analysis below.

Remark 3.4.

Consider a function f ∈ W ^1,2(E) and the function f ̂ = f ( T E ( x ̂ ) ) . It can be concluded by (3.11) that ∇ x ̂ f ̂ = B E ∇ x f . By applying the change of variables (3.11), it can be shown that ∫ E ̂ | ∂ x ̂ i f ̂ | p d x ̂ = | B E | − 1 | h i , E p ∫ E | ∂ x i f | p d x , i = 1, …, d, where p > 1, and | B E | is the determinant of B E .

On T h we define the finite dimensional space

(3.12) V h k = φ h ∈ C 0 ( Q ̄ ) : φ h | E ∈ Q k ( E ) ∀ E ∈ T h , φ h = 0 on Σ ∪ Σ 0 .

Here Q k ( E ) is the space of polynomials on E composed by tensor products of uni-variate Lagrange polynomials of degree at most k with respect to each variable, i.e., Q k ( E ) = ⊗ i = 1 d P k ( x i ) , where P k ( x i ) , i = 1 , … , d is the uni-variate Lagrange polynomials. In the analysis below, we consider the case of k = 1.

Proposition 3.1.

Consider a polynomial function φ h ∈ Q k ( E ) , E ∈ T h such that φ ̂ h = φ h ( T E ( x ̂ ) ) ∈ Q k ( E ̂ ) . There is a constant c _inv independent of h _E such that

(3.13) ∑ i = 1 d ‖ ∂ x i φ h ‖ L 2 ( E ) ⩽ c inv ∑ i = 1 d h i , E − 1 ‖ φ h ‖ L 2 ( E ) .

Proof.

Since all the norms of φ h ̂ ∈ Q k ( E ̂ ) are equivalent, there is a c > 0 such that

(3.14) ‖ ∂ x i φ ̂ h ‖ L 2 ( E ̂ ) ⩽ c ‖ φ ̂ h ‖ L 2 ( E ̂ ) .

Utilizing the form of the transformation (3.11) and applying the chain rule ∂ x i φ h = ∂ x ̂ i φ ̂ h ∂ x ̂ i ∂ x i , we can obtain

(3.15) ∫ E ( ∂ x i φ h ) 2 d x = h i , E − 2 | B E | − 1 ∫ E ̂ ( ∂ x ̂ i φ ̂ h ) 2 d x ̂ ( 3.14 ) ⩽ c h i E − 2 | B E | − 1 ∫ E ̂ ( φ ̂ h ) 2 d x ̂ ⩽ c h i , E − 2 ∫ E ( φ h ) 2 d x ,

from where we can deduce (3.13) by summing over i = 1, …, d. □

Corollary 3.1.

By inequality (3.13) we can infer

(3.16) ∑ i = 1 d ‖ D x i φ h ‖ L 2 ( E ) 2 ⩽ c inv ∑ i = 1 d h E , min − 2 ‖ φ h ‖ L 2 ( E ) 2 .

3.2.2 The unified space-time FE scheme

Assumption 3.3.

For simplicity in the discretization error analysis we suppose that u _0,h = u ₀ ≔ 0. For problems with u ₀ ≠ 0 we refer to refs. [10], [19].

Based on (3.4) and using (3.5) we consider the following discretization of problem (3.4): Find u h ∈ V h k such that u _h = u _0,h on Σ₀ and

(3.17) B ( u h , v h ) = ℓ f ( v h ) ∀ v h ∈ V h k .

In order to obtain stable solutions the scheme (3.17) is modified by adding an upwind stabilization term, and the final discrete problem is written: Find u h ∈ V h k such that u _h = u _0,h on Σ₀ and

(3.18a) B s ( u h , v h ) ≔ B ( u h , v h ) + S ( u h , v h ) + R s ( u h , v h ) = ℓ f v h + τ λ β Q T ⋅ ∇ v h ∀ v h ∈ V h k ,

where the vector β Q T ≔ ( 1 , β ) ∈ R d , and the streamline-upwind (SU) term S is

(3.18b) S ( u h , v h ) ≔ ∑ E ∈ T h ∫ E τ λ β Q T ⋅ ∇ u h β Q T ⋅ ∇ v h d x d t .

Here τ λ ≔ ϑ h E , min λ , with ϑ > 0, λ > 1 fixed parameters, which will be specified below. The residual term has the form

(3.18c) R s ( u h , v h ) ≔ ∑ E ∈ T h ∫ E − τ λ ε Δ u h − r u h β Q T ⋅ ∇ v h d x d t

and the linear form

(3.19) ℓ f v h + τ λ β Q T ⋅ ∇ v h ≔ ∑ E ∈ T h ∫ E f ( v h + τ λ β Q T ⋅ ∇ v h ) d x d t .

Remark 3.5.

Note that, in case of working with linear spaces, i.e., V h k = 1 , the residual terms take the form R s ( u h , v h ) = ∑ E ∈ T h ∫ E τ λ r u h β Q T ⋅ ∇ v h d x d t .

Remark 3.6

(consistency). Under the Assumption 3.1, the following localized variational form

(3.20) B s ( u , v h ) ≔ B ( u , v h ) + S ( u , v h ) + R s ( u , v h ) = ℓ f v h + τ λ β Q T ⋅ ∇ v h ∀ v h ∈ V h k

holds for the weak solution u.

In view of (3.20) and (3.18a), we have the following equation.

Corollary 3.2.

Let u be the solution of problem (3.4) and u _h the solution of problem (3.18a). Then the following error equation holds for v h ∈ V h k ,

(3.21) B s ( u − u h , v h ) = 0 .

Below the coercivity and boundedness properties of B _s(⋅, ⋅) are discussed.

Lemma 3.1.

Let v h ∈ V h k and the Assumption 3.3. Then

(3.22) B s ( v h , v h ) ⩾ 1 2 ‖ v h ‖ L 2 ( Σ T ) 2 + ∑ E ∈ T h ε 2 ‖ ∇ x v h ‖ L 2 ( E ) 2 + r 2 ‖ v h ‖ L 2 ( E ) 2 + τ λ 2 ‖ β Q T ⋅ ∇ v h ‖ L 2 ( E ) 2 .

Proof.

For E i , ∈ T h , let n E i = ( n t , E i , n x , E i ) be the unit normal vector on ∂E _i. Let n _e,ij = (n _t,ij, n _x,ij) be the unit normal vector on the common faces e _ij = ∂E _i ∩ ∂E _j for E i , E j ∈ T h . Using (2.4b) and the identity 2 v ( β Q T ⋅ ∇ v ) = β Q T ⋅ ∇ v 2 = ∇ ⋅ ( β Q T v 2 ) , it follows immediately that

(3.23) ∑ E i ∈ T h ∫ E i v h ( β Q T ⋅ ∇ v h ) d x d t = 1 2 ∑ E i ∈ T h ∫ E i ∇ ⋅ ( β Q T v h ) d x d t = 1 2 ∑ E i ∈ T h ∫ ∂ E i n E i ⋅ β Q T v h 2 d S = 1 2 ∑ e i j ∫ e i j β Q T ⋅ n e , i j v h 2 | E i − v h 2 | E j d S = 1 2 ∫ Σ T v h 2 d S .

Recalling the terms of B _s(⋅, ⋅) in (3.18) and following the same steps as in (3.6) and using (3.23), we can derive the bound

(3.24) B s ( v h , v h ) = ∑ E ∈ T h ∫ E β Q T ⋅ ∇ v h v h + ε ∇ x v h ⋅ ∇ x v h + r v h 2 d x d t + ∑ E ∈ T h ∫ E τ λ − ε Δ v h + β Q T ⋅ ∇ v h + r v h β Q T ⋅ ∇ v h d x d t ⩾ 1 2 ‖ v h ‖ L 2 ( Σ T ) 2 + ∑ E ∈ T h r ‖ v h ‖ L 2 ( E ) 2 + ε ‖ ∇ x v h ‖ L 2 ( E ) 2 + τ λ ‖ β Q T ⋅ ∇ v h ‖ L 2 ( E ) 2 + ∑ E ∈ T h ∫ E τ λ − ε Δ v h + r v h β Q T ⋅ ∇ v h d x d t .

Now, for the last sum in (3.24), we apply (2.3a) and obtain

(3.25) ∑ E ∈ T h ∫ E τ λ − ε Δ v h + r v h β Q T ⋅ ∇ v h d x d t ⩽ ∑ E ∈ T h ∫ E − τ λ / 2 ε Δ v h τ λ / 2 β Q T ⋅ ∇ v h d x d t + ∫ E τ λ / 2 r v h τ λ / 2 β Q T ⋅ ∇ v h d x d t ⩽ ∑ E ∈ T h τ λ ε 2 ‖ Δ v h ‖ L 2 ( E ) 2 + τ λ r 2 ‖ v h ‖ L 2 ( E ) 2 + τ λ 2 ‖ β Q T ⋅ ∇ v h ‖ L 2 ( E ) 2 .

Here, we make use of the discrete-inverse inequality, ‖ Δ v h ‖ L 2 ( E ) 2 ⩽ c inv 2 h E , min 2 ‖ ∇ x v h ‖ L 2 ( E ) 2 , see (3.13) and (3.16), and the identity a b ⩽ a 2 + 1 4 b 2 to find

(3.26) ∑ E ∈ T h ∫ E τ λ − ε Δ v h + r v h β Q T ⋅ ∇ v h d x d t ⩽ ∑ E ∈ T h ε 2 τ λ c inv h E , min 2 ‖ ∇ x v h ‖ L 2 ( E ) 2 + τ λ r 2 ‖ v h ‖ L 2 ( E ) 2 + τ λ 2 ‖ β Q T ⋅ ∇ v h ‖ L 2 ( E ) 2 .

Now, introducing the result (3.26) into (3.24), and using that 0 < ɛ ⩽ 1, we find

(3.27) B s ( v h , v h ) ⩾ 1 2 ‖ v h ‖ L 2 ( Σ T ) 2 + ∑ E ∈ T h r ‖ v h ‖ L 2 ( E ) 2 + ε ‖ ∇ x v h ‖ L 2 ( E ) 2 + τ λ ‖ ( β Q T ⋅ ∇ v h ) ‖ L 2 ( E ) 2 − ∑ E ∈ T h ε 2 τ λ c inv h E , min 2 ‖ ∇ v h ‖ L 2 ( E ) 2 + τ λ r 2 ‖ v h ‖ L 2 ( E ) 2 + τ λ 2 ‖ β Q T ⋅ ∇ v h ‖ L 2 ( E ) 2 = 1 2 ‖ v h ‖ L 2 ( Σ T ) 2 + ∑ E ∈ T h ε − ε 2 τ λ c inv h E , min 2 ‖ ∇ x v h ‖ L 2 ( E ) 2 + ( r − τ λ r 2 ) ‖ v h ‖ L 2 ( E ) 2 + τ λ − τ λ 2 ‖ β Q T ⋅ ∇ v h ‖ L 2 ( E ) 2 .

Making the appropriate choice for the parameter ϑ = 1 2 min 1 2 c inv 2 , 1 2 r , and setting λ = 2, we can infer the bound (3.22). □

In view of (3.22), we introduce the mesh-dependent norms

(3.28a) ‖ v ‖ s 2 = 1 2 ‖ v ‖ L 2 ( Σ T ) 2 + ∑ E ∈ T h ε 2 ‖ ∇ x v ‖ L 2 ( E ) 2 + r 2 ‖ v ‖ L 2 ( E ) 2 + τ λ 2 ‖ β Q T ⋅ ∇ v ‖ L 2 ( E ) 2 ,

(3.28b) ‖ v ‖ h , s 2 ≔ ‖ v ‖ s 2 + ∑ E ∈ T h τ − λ ‖ v ‖ L 2 ( E ) 2 + ∑ E ∈ T h τ λ ε 2 ‖ Δ v ‖ L 2 ( E ) 2 .

By (3.22) we can immediately write the estimate

(3.29) B s ( v h , v h ) ⩾ ‖ v h ‖ s 2 ∀ v h ∈ V h k .

Setting in (3.29) v h = z h − u h ∈ V h k and recalling the error equation (3.21), we get

(3.30) ‖ z h − u h ‖ s 2 ⩽ B s ( z h − u h , z h − u h ) + B s ( u h − u , z h − u h ) = B s ( z h − u , z h − u h ) .

Lemma 3.2.

Let u the weak solution of (3.5) under Assumption 3.1 and u ₀ = 0. Let u h ∈ V h k be the finite element solution in (3.17). The approximation error estimate

(3.31) ‖ u − u h ‖ s ⩽ C 0 , a p p r o x ‖ u − z h ‖ h , s

holds for all z h ∈ V h k , where the constant C _0,approx is independent of h.

Proof.

Let a z h ∈ V h k . Recalling the terms appearing in (3.18), we have

(3.32) B s ( u − z h , u h − z h ) ≔ B ( u − z h , u h − z h ) + S ( u − z h , u h − z h ) + R s ( u − z h , u h − z h )

with

(3.33) B ( u − z h , u h − z h ) ≔ ∫ Q T ∂ t ( u − z h ) ( u h − z h ) + ε ∇ x ( u − z h ) ⋅ ∇ x ( u h − z h )

+ β ⋅ ∇ x ( u − z h ) ( u h − z h ) + r ( u − z h ) ( u h − z h ) d x d t ,

(3.34) S ( u − z h , u h − z h ) ≔ ∑ E ∈ T h ∫ E τ λ β Q T ⋅ ∇ ( u − z h ) β Q T ⋅ ∇ ( u h − z h ) d x d t ,

(3.35) R s ( u − z h , u h − z h ) ≔ ∑ E ∈ T h ∫ E − τ λ ε Δ ( u − z h ) β Q T ⋅ ∇ ( u h − z h ) d x d t .

Using appropriately (2.3a), we derive the following bounds for the terms in (3.33):

(3.36) ∫ Q T ε ∇ x ( u − z h ) ⋅ ∇ x ( u h − z h ) d x d t ⩽ ε ‖ ∇ x ( u − z h ) ‖ L 2 ( Q T ) ‖ u h − z h ‖ s ,

(3.37) ∫ Q T ( u t − z h , t ) + β ⋅ ∇ x ( u − z h ) ( u h − z h ) d x d t = ∫ Q T β Q T ⋅ ∇ ( u − z h ) ( u h − z h ) d x d t = ( 2.4 ) − ∫ Q T β Q T ⋅ ∇ ( u h − z h ) ( u − z h ) d x d t + ∫ ∂ Q T β Q T ⋅ n ( u − z h ) ( u h − z h ) d s ⩽ τ λ ‖ β Q T ⋅ ∇ ( u h − z h ) ‖ L 2 ( Q T ) 2 1 / 2 τ − λ ‖ u − z h ‖ L 2 ( Q T ) 2 1 / 2 ( using BCs ) + ‖ u h − z h ‖ L 2 ( Σ T ) ‖ u − z h ‖ L 2 ( Σ T ) ⩽ 2 τ − λ ‖ u − z h ‖ L 2 ( Q T ) 2 1 / 2 + 2 ‖ u − z h ‖ L 2 ( Σ T ) ‖ u h − z h ‖ s ,

and also

(3.38) ∫ Q T r 2 1 / 2 ( u h − z h ) ( 2 r ) 1 / 2 ( u − z h ) d x d t ⩽ r 2 ‖ u h − z h ‖ L 2 ( Q T ) 2 1 / 2 2 r ‖ u − z h ‖ L 2 ( Q T ) 2 1 / 2 ⩽ C r ‖ u − z h ‖ L 2 ( Q T ) ‖ u h − z h ‖ s .

Treating together the terms in (3.34) and (3.35), we get

(3.39) ∑ E ∈ T h ∫ E τ λ / 2 − ε Δ ( u − z h ) + β Q T ⋅ ∇ ( u − z h ) + r ( u − z h ) τ λ / 2 β Q T ⋅ ∇ ( u h − z h ) d x d t ⩽ ∑ E ∈ T h τ λ / 2 ε ‖ Δ u − Δ z h ‖ L 2 ( E ) + ‖ β Q T ⋅ ∇ ( u − z h ) ‖ L 2 ( E ) + r ‖ u − z h ‖ L 2 ( E ) × τ λ / 2 ‖ β Q T ⋅ ∇ ( u h − z h ) ‖ L 2 ( E ) ≲ ∑ E ∈ T h τ λ ε 2 ‖ Δ u − Δ z h ‖ L 2 ( E ) 2 + ‖ β Q T ⋅ ∇ ( u − z h ) ‖ L 2 ( E ) 2 + r 2 ‖ u − z h ‖ L 2 ( E ) 2 1 / 2 ‖ u h − z h ‖ s .

Collecting the previous bounds, and using (3.32) we can obtain

(3.40) B s ( u − z h , u h − z h ) ⩽ ‖ u − z h ‖ h , s ‖ u h − z h ‖ s .

Combining (3.40) and (3.30) and adding ‖u − z _h‖_s on both sides of the inequality we can deduce

(3.41) ‖ u h − z h ‖ s + ‖ u − z h ‖ s ≲ 2 ‖ u − z h ‖ h , s .

The result (3.31) can be derived by applying triangle inequality. □

Note that the error estimate (3.31) includes bounds for the term ‖ ∂ t u − ∂ t u h ‖ L 2 , compare with Remarks 3.1 and 3.2.

3.2.3 Quasi-interpolation

Let V = W ^ℓ,p(Q _T), ℓ ⩾ 2, p > 1. We introduce the quasi-interpolant, [29], [33],

(3.42) Π h : V → V h k , Π h f = ∑ i = 1 dim V h k λ i ( f ) φ i , h ,

where λ _i are the linear dual functionals and φ _i,h the local shape functions of V h k . The multivariate quasi-interpolant given in (3.42) can be constructed by the corresponding uni-variate quasi-interpolants π ₁, …, π _d, i.e., Π_h f = (π ₁ ⊗ π ₂ ⋯ ⊗ π _d)(f), and each dual functional λ _i(f) can be similarly defined by the corresponding uni-variate dual functionals, see, e.g., [29], [33], [34], [35]. For convenience one can consider a dual basis satisfying λ _i(φ _j,h) = δ _ij, i , j = 1 , … , dim V h k , [33]. By this property it can be immediately inferred that Π_h(φ _h) = φ _h for all φ h ∈ V h k .

Remark 3.7.

A discussion for the dual functionals λ _i(f) when f belongs to tensor spaces is given in refs. [34], [35]. In ref. [25] an analysis is given for the case of low regularity functions f and anisotropic meshes.

The following stability bounds can be shown, [29], [33], [35].

Proposition 3.2.

Let f be a smooth function. The bounds

(3.43a) ‖ Π h f ‖ L p ( Q T ) ⩽ C 0 ‖ f ‖ L p ( Q T ) ,

(3.43b) ‖ Π h f ‖ W 1 , p ( Q T ) ⩽ C 1 ‖ f ‖ W 1 , p ( Q T )

hold, where C ₀, C ₁ are positive constants.

Suppose D ⊂ R d is a cuboid domain and h D = ( h 1 , … , h d ) is a vector containing the sizes h _i in each dimension x _i, i = 1, …, d. Define the set of multi-indices A = { m e i , i = 1 , … , d , m ∈ N } , where e _i are the orthonormal basis vectors, and denote with A ⁰ = {m ₀ = (m ₁, …, m _d):m _i < m} the set of multi-indices which corresponds to the tensor product polynomials of degree less than m. Note that A ⁰ can also be charactirized as A 0 = m 0 : D a x m 0 = 0 ∀ a ∈ A . For a smooth function f, we define the tensor Taylor polynomial of degree less than m evaluated at y:

(3.44a) T y m f ( x ) = ∑ m 0 ∈ A 0 1 m 0 ! D m 0 f ( y ) ( x − y ) m 0

and the averaged Taylor expansion over the ball B ⊂ D :

(3.44b) T y , φ A f ( x ) = ∫ B T y m f ( x ) φ ( y ) d y ,

where φ ∈ C 0 ∞ ( B ) is a cut-off function with ∫_B φ(y) dy = 1, see [29], Chapt. 4].

In the analysis below we focus mainly on the case m = 2. Let the multi-indices a ∈ A and β = (β ₁, …, β _d) such that | β | ⩽ 1.

Proposition 3.3.

Let f be a smooth function. The remainder R A f ≔ f − T y , φ A f has the form

(3.45a) R A f ( x ) = ∑ a ∈ A ∫ D k a ( x , y ) D a f ( y ) d y ,

where

(3.45b) k a ( x , y ) = C a , m ( x − y ) a k ( x , y )

with

(3.45c) k ( x , y ) ⩽ C | x − y | − d , | D x β k a ( x , y ) | ⩽ C | x − y | | a | − | β | − d .

Proof.

The bounds can be shown by combining the results given in refs. [29], [35], [36]. □

We have the commutativity result.

Proposition 3.4.

Let f be a smooth function then

(3.46) D β T y , φ A f = T y , φ A − β D β f .

Proof.

Note that for β > m ₀ holds that D x β x m 0 = 0 and thus

(3.47) D x β T y , φ A f ( x ) = ∑ m 0 ∈ A 0 ∫ B φ ( y ) D y m 0 f ( y ) ( x − y ) ( m 0 − β ) ( m 0 − β ) ! d y = ∑ γ ∈ ( A 0 − β ) ∫ B φ ( y ) D y β + γ f ( y ) ( x − y ) γ γ ! d y = ∑ γ ∈ ( A 0 − β ) ∫ B φ ( y ) D y γ D y β f ( y ) ( x − y ) γ γ ! d y = T y , φ A − β ( D β f ) ( x ) .

□

Let the multi-indices a ∈ A and β = (β ₁, …, β _d) such that | β | ⩽ 1. Following similar ideas as in refs. [25], [29], [35], we show the estimates.

Theorem 3.1.

Let f ∈ W ℓ , p ( D ) with ℓ = m. Then there exist a tensor Taylor polynomial T A f ≔ T y , φ A f such that

(3.48a) ‖ f − T A f ‖ L p ( D ) ⩽ C 0 ∑ a ∈ A h D a ‖ D a f ‖ L p ( D ) ,

(3.48b) ‖ D β ( f − T A f ) ‖ L p ( D ) ⩽ C 1 ∑ a ∈ A , γ = ( a − β ) h D γ ‖ D γ ( D β f ) ‖ L p ( D ) ,

where C ₀ > 0, C ₁ > 0 are independent of the mesh sizes h _i, i = 1, …, d.

Proof.

Assume initially that D ⊂ R d is a (reference) cuboid domain with d i a m ( D ) ∼ 1 , e.g., D ≔ D ̂ = [ − 1,1 ] d . By Proposition 3.3 we can deduce

(3.49) ∫ D ̂ k a ( x , z ) D a f ( z ) d z ⩽ C ‖ D a f ‖ L p ( D ̂ ) ∀ a ∈ A ,

which implies that

(3.50) ‖ f − T A f ‖ L p ( D ̂ ) ⩽ C 0 ∑ a ∈ A ‖ D a f ‖ L p ( D ̂ ) .

Further, combining the last bounds with Proposition 3.4, we can get

(3.51) ‖ D β ( f − T A f ) ‖ L p ( D ̂ ) = ‖ R A − β D β f ( x ) ‖ L p ( D ̂ ) ⩽ C ∑ a ∈ A , γ = ( a − β ) ‖ D γ ( D β f ( x ) ) ‖ L p ( D ̂ ) .

Utilizing now the tensor character of the mesh, we can apply the change of variables x = B D x ̂ + b , see (3.11), and to transform D ̂ to D . Then we can show that, see Remark 3.4,

(3.52) ‖ D β ( f − T A f ) ‖ L p ( D ̂ ) p = | B D | − 1 h D p β ‖ D β ( f − T A f ) ‖ L p ( D ) p , ‖ D a f ‖ L p ( D ̂ ) p = | B D | − 1 h D p a ‖ D a f ‖ L p ( D ) .

Inserting the results (3.52) into the estimates (3.50), (3.51), we can deduce the required estimates (3.48). □

Having (3.48) we can show bounds on how well the quasi-interpolant Π_h f approximates the function f. We follow standard ideas from the finite element methodology. Recall that a ∈ A and β = (β ₁, …, β _d) with | β | ⩽ 1.

Lemma 3.3.

Let a function f ∈ W ^ℓ,2(Q _T), ℓ = m, and Assumption 3.2 for the mesh T h . Then for the quasi-interpolant (3.42) the interpolation estimates

(3.53a) ‖ f − Π h f ‖ L 2 ( E ) 2 ⩽ c intp , 0 ∑ a ∈ A h E a ‖ D a f ‖ L 2 ( E ) , E ∈ T h ,

(3.53b) ‖ D β ( f − Π h f ) ‖ L 2 ( E ) ⩽ c intp , 1 ∑ a ∈ A , γ = ( a − β ) h E γ ‖ D γ ( D β f ) ‖ L 2 ( E ) , E ∈ T h ,

hold, where the positive constants c _intp,0, c _intp,1 are independent of h _E.

Proof.

By the results of Theorem 3.1 there exists a tensor polynomial p f ∈ P m 0 , m 0 ∈ A 0 such that

(3.54) ‖ f − p f ‖ L 2 ( E ) ≲ ∑ a ∈ A h E a ‖ D a f ‖ L 2 ( E ) , E ∈ T h .

Using that Π h ( v h ) = v h ∀ v h ∈ V h k and the stability properties (3.48a), (3.43a) we get

(3.55) ‖ f − Π h f ‖ L 2 ( E ) ⩽ ‖ f − p f ‖ L 2 ( E ) + ‖ Π h ( f − p f ) ‖ L 2 ( E ) ⩽ C ‖ f − p f ‖ L 2 ( E ) ⩽ ∑ a ∈ A h E a ‖ D a f ‖ L 2 ( E ) ⩽ ∑ a ∈ A h a ‖ D a f ‖ L 2 ( E ) .

Following the same steps as above and using the interpolation estimate (3.48b) we can prove (3.53b). □

Remark 3.8.

Note that the anisotropic character of the mesh appears in the interpolation estimates given in (3.53). Similar estimates have been shown in ref. [25] for low regularity functions on anisotropic triangular meshes.

Lemma 3.4.

Consider a mesh T h , with uniform directional mesh widths, i.e., h ₁ ∼ h ₂ ∼ ⋯ ∼ h _d and an element E ∈ T h with diam(E) ≔ h _E. Let a function f ∈ W ^ℓ,p(Ω), 1 ⩽ ℓ, 1 < p, and the local intepolation operator Π h ∈ Q k = l − 1 . Then the interpolation estimate

(3.56) | f − Π h f | W m , p ( E ) ≲ h E ℓ − m | f | W ℓ , p ( E )

holds for 0 ⩽ m ⩽ ℓ.

Proof.

See, cf. [29]. □

For isotropic meshes it is known (cf. [29]), that there is a constant C _trc > 0, such that

(3.57) ‖ v ‖ L p ( ∂ E ) 2 ⩽ C trc h − 1 ‖ v ‖ L 2 ( E ) + h ‖ ∇ v ‖ L 2 ( E ) 2 , E ∈ T h .

Proposition 3.5.

Let a face e ⊂ ∂E of a mesh element E ∈ T h such that |E| = |e| h _⊥,e, where h _⊥,e is the measure of the edge perpendicular to e. Let a function f ∈ H ¹(Q _T). Then there is a positive constant C _trc such that

(3.58) ‖ f ‖ L 2 ( e ) 2 ⩽ C trc h ⊥ , e − 1 ‖ f ‖ L 2 ( E ) 2 + ∑ i = 1 d h i , E 2 ‖ ∂ x i f ‖ L 2 ( E ) 2 .

Proof.

Let f ̂ ≔ f ◦ T E ∈ H 1 ( E ̂ ) . For the face e ̂ ⊂ ∂ E ̂ the trace ineguality has the form (cf. [29]),

(3.59) ‖ f ̂ ‖ L 2 ( e ̂ ) 2 ≲ ‖ f ̂ ‖ L 2 ( E ̂ ) 2 + ‖ ∇ f ̂ ‖ L 2 ( E ̂ ) 2 .

Recalling (3.11) and Remark 3.4, we can transform (3.59) onto E to deduce

(3.60) ‖ f ‖ L 2 ( e ) 2 | e | − 1 ≲ | E | − 1 ‖ f ‖ L 2 ( E ) 2 + ‖ B E ∇ f ‖ L 2 ( E ) 2 .

Owing to the form of B E we have ‖ B E ∇ f ‖ L 2 ( E ) 2 ≲ ∑ i = 1 d h i , E 2 ‖ ∂ x i f ‖ L 2 ( E ) 2 . Using that h _⊥,e|e| = |E| we can obtain (3.58). □

Remark 3.9.

For e ⊂ Σ_T inequality (3.58) is valid for h _⊥,e ∼ h _t.

Corollary 3.3.

Let the assumptions of Lemma 3.3 and assume f ∈ H ²(Q _T). Then the following interpolation estimate

(3.61) ‖ f − Π h f ‖ L 2 ( Σ T ) 2 ≲ ∑ E ∈ T h , ∂ E ⊂ Σ T h E − e 1 ∑ a h E 2 a + ∑ i = 1 d h E 3 e i + h i ∑ j ≠ i d h j 2 ‖ f ‖ H 2 ( E ) 2

holds for all multiindices a ∈ A and e _i to be the i-th unit normal vector.

Proof.

According to (3.58) and Remark 3.9,

(3.62) ‖ f − Π h f ‖ L 2 ( Σ T ) 2 = ∑ e ∈ E , e ⊂ Σ T ∫ e | f − Π h f | 2 d s ⩽ c ∑ E ∈ T h , ∂ E ⊂ Σ T h t − 1 ‖ f − Π h f ‖ L 2 ( E ) 2 + ∑ i = 1 d h i 2 ‖ ∂ x i ( f − Π h f ) ‖ L 2 ( E ) 2 ⩽ c ∑ E ∈ T h , ∂ E ⊂ Σ T h t − 1 ∑ a h E 2 a ‖ D a f ‖ L 2 ( E ) 2 + ∑ i = 1 d h i 2 ∑ γ = a − e i h E γ ‖ D γ ( ∂ x i f ) ‖ L 2 ( E ) 2 ⩽ c ∑ E ∈ T h , ∂ E ⊂ Σ T h E − e 1 ∑ a h E 2 a + ∑ i = 1 d ∑ γ = a − e i h E 2 e i h E γ ‖ f ‖ H 2 ( E ) 2 ⩽ c ∑ E ∈ T h , ∂ E ⊂ Σ T h E − e 1 ∑ a h E 2 a + ∑ i = 1 d ∑ a h E a + e i ‖ f ‖ H 2 ( E ) 2 .

Expanding the last sum in (3.62), we get

‖ f − Π h f ‖ L 2 ( Σ T ) 2 ⩽ c ∑ E ∈ T h , ∂ E ⊂ Σ T h E − e 1 ∑ a h E 2 a + ∑ i = 1 d h E 3 e i + h i ∑ j ≠ i d h j 2 f H 2 ( E ) 2

as reguired. □

Lemma 3.5.

Let f ∈ V satisfying the assumptions of Lemma 3.3, and let the associated interpolant Π_h f, see (3.53). Then there exist a constant independent of f and h _E such that the following quasi-interpolation estimate

(3.63) ‖ f − Π h f ‖ h , s 2 ≲ ∑ E ∈ T h ∑ a ∈ A h E 2 a − e 1 + ∑ i = 1 d ∑ a ∈ A h E a + e i − e 1 + ∑ a ∈ A τ − λ h E 2 a + ∑ j = 2 d ∑ a ∈ A τ λ ε 2 h E 2 ( a − 2 e j ) + ε 2 h E 2 ( a − e j ) + ∑ i = 1 d ∑ a ∈ A τ λ 2 h E 2 ( a − e i ) ‖ f ‖ H 2 ( E ) 2

holds true.

Proof.

Recall that

(3.64) ‖ f − Π h f ‖ h , s 2 = 1 2 ‖ f − Π h f ‖ L 2 ( Σ T ) 2 + ∑ E ∈ T h τ − λ ‖ f − Π h f ‖ L 2 ( E ) 2 + ∑ E ∈ T h τ λ ε 2 ‖ Δ f − Δ Π h f ‖ L 2 ( E ) 2 + ∑ E ∈ T h ε 2 ‖ ∇ x ( f − Π h f ) ‖ L 2 ( E ) 2 + r 2 ‖ f − Π h f ‖ L 2 ( E ) 2 + τ λ 2 ‖ β Q T ⋅ ∇ ( f − Π h f ) ‖ L 2 ( E ) 2 .

Under the regularity assumptions for f and utilizing (3.53) and (3.61), we can infer the following estimates

‖ f − Π h f ‖ L 2 ( Σ T ) 2 ≲ ∑ E ∈ T h , ∂ E ⊂ Σ T ∑ a ∈ A h E 2 a − e 1 + ∑ i = 1 d ∑ a ∈ A h E a + e i − e 1 ‖ f ‖ H 2 ( E ) 2

and also

∑ E ∈ T h τ − λ ‖ f − Π h f ‖ L 2 ( E ) 2 ≲ ∑ E ∈ T h ∑ a ∈ A τ − λ h E 2 a ‖ f ‖ H 2 ( E ) 2 , ∑ E ∈ T h τ λ ε 2 ‖ Δ ( f − Π h ) f ‖ L 2 ( E ) 2 ≲ ∑ E ∈ T h ∑ j = 2 d ∑ a ∈ A τ λ ε 2 h E 2 ( a − 2 e j ) ‖ f ‖ H 2 ( E ) 2 , ∑ E ∈ T h ε 2 ‖ ∇ x ( f − Π h f ) ‖ L 2 ( E ) 2 ≲ ∑ E ∈ T h ∑ j = 2 d ∑ a ∈ A ε 2 h E 2 ( a − e j ) ‖ f ‖ H 2 ( E ) 2 , ∑ E ∈ T h τ λ 2 ‖ β Q T ⋅ ∇ ( f − Π h f ) ‖ L 2 ( E ) 2 ≲ ∑ E ∈ T h ∑ j = 1 d ∑ a ∈ A τ λ 2 h E 2 ( a − e j ) ‖ f ‖ H 2 ( E ) 2 .

Inserting the previous estimates in (3.64) we derive (3.63). □

Example 1

(isotropic-uniform mesh). Assume Q T = [ 0 , L ] × [ 0 , L ] ⊂ R d = 2 , and assume a T h with uniform uni-directional grid sizes, that means h _t ∼ h _x ∼ h. Consider a function f ∈ H ²(Q _T) and the set of multi-indices A = {a = (a ₁, a ₂):a = 2e _i, i = 1, 2}. Then by (3.63) we can conclude that Π_h satisfies on isotropic meshes the following estimate

(3.65) ‖ f − Π h f ‖ h , s 2 ≲ h 2 + τ − λ h 4 + τ λ ε 2 + ε 2 h 2 + τ λ 2 h 2 ∑ E ∈ T h ‖ f ‖ H 2 ( E ) 2 .

Example 2

(independent directional mesh sizes for Q T ⊂ R d = 3 ). Assume that Q T = [ 0 , L ] 3 ⊂ R d = 3 , and for the uni-directional grid sizes assume h _E ≔ h = (h _t ≔ h ₁, h ₂ ≔ h _x, h ₃ ≔ h _y). Consider a function f ∈ H ²(Q _T) and the set of multi-indices A = {a = (a ₁, a ₂, a ₃):a = 2e _i, i = 1, 2, 3}. Then by the estimate given in (3.63) we can have

(3.66) ‖ f − Π h f ‖ h , s 2 ≲ ∑ a ∈ A h 1 2 a 1 − 1 h 2 2 a 2 h 3 2 a 3 + h 1 a 1 h 2 a 2 h 3 a 3 + h 1 a 1 − 1 h 2 a 2 + 1 h 3 a 3 + h 1 a 1 − 1 h 2 a 2 h 3 a 3 + 1 + τ − λ h 1 2 a 1 h 2 2 a 2 h 3 2 a 3 + τ λ ε 2 h 1 2 a 1 h 2 2 a 2 − 4 h 3 2 a 3 + h 1 2 a 1 h 2 2 a 2 h 3 2 a 3 − 4 + ε 2 h 1 2 a 1 h 2 2 a 2 − 2 h 3 2 a 3 + h 1 2 a 1 h 2 2 a 2 h 3 2 a 3 − 2 + τ λ 2 h 1 2 a 1 − 2 h 2 2 a 2 h 3 2 a 3 + h 1 2 a 1 h 2 2 a 2 − 2 h 3 2 a 3 + h 1 2 a 1 h 2 2 a 2 h 3 2 a 3 − 2 ∑ E ∈ T h ‖ f ‖ H 2 ( E ) 2 ≲ h 1 3 + h 1 − 1 h 2 4 + h 3 4 + h 1 2 + h 2 2 + h 3 2 + h 1 h 2 + h 1 − 1 h 2 3 + h 1 − 1 h 2 h 3 2 + h 1 h 3 + h 1 − 1 h 2 2 h 3 + h 1 − 1 h 3 3 + τ − λ h 1 4 + h 2 4 + h 3 4 + τ λ ε 2 2 + h 1 4 h 2 − 4 + h 2 − 4 h 3 4 + h 1 4 h 3 − 4 + h 2 4 h 3 − 4 + ε 2 h 1 4 h 2 − 2 + h 2 2 + h 2 − 2 h 3 4 + h 1 4 h 3 − 2 + h 2 4 h 3 − 2 + h 3 2 + τ λ 2 h 1 2 + h 1 − 2 h 2 4 + h 1 − 2 h 3 4 + h 1 4 h 2 − 2 + h 2 2 + h 2 − 2 h 3 4 + h 1 4 h 3 − 2 + h 2 4 h 3 − 2 + h 3 2 ∑ E ∈ T h ‖ f ‖ H 2 ( E ) 2 .

Theorem 3.2.

Let the solutions u and u _h satisfy the assumptions in Lemma 3.2 and let the Π_h satisfy the assumptions in Lemma 3.3. Then the following error convergence result holds

(3.67) ‖ u − u h ‖ h 2 ≲ ∑ E ∈ T h ∑ a ∈ A h E 2 a − e 1 + ∑ i = 1 d ∑ a ∈ A h E a + e i − e 1 + ∑ a ∈ A τ − λ h E 2 a + ∑ j = 2 d ∑ a ∈ A τ λ ε 2 h E 2 ( a − 2 e j ) + ε 2 h E 2 ( a − e j ) + ∑ i = 1 d ∑ a ∈ A τ λ 2 h E 2 ( a − e i ) ‖ u ‖ H 2 ( E ) 2 .

Proof.

We combine Lemmas 3.2 and 3.5 and the assertion follows. □

Corollary 3.4.

Let T h be a uniform mesh partition, i.e., h 1 ∼ ⋯ ∼ h d x + 1 ∼ h . Let further that u ∈ H ^ℓ(Q _T) with ℓ ⩾ 2 and u h ∈ V h k , k = 1 . Under the assumptions of Theorem 3.2 holds

(3.68) ‖ u − u h ‖ h ≲ c h ℓ − 1 ‖ u ‖ H ℓ ( Q T )

with c > 0 depending on the problem data.

Proof.

The convergence rate (3.68) follows immediately by applying the mesh uniformity properties in (3.63) and in (3.67). □

4 Numerical examples

In order to validate the estimates derived in the previous sections, a series of numerical tests are presented below choosing different values for the parameters of the problem. For the two-dimensional problems we set τ λ = ϑ h x λ with ϑ = 0.5, λ = 2, and the mesh sizes h _x, h _t are chosen independently depending on the purpose of the numerical computation. For the three dimensional problems we set τ λ = ϑ h E , min 2 . First, we start by considering the problem on a space time cylinder Q T ⊂ R 2 with a smooth solution and then with a less regular solution. Note that in this case β ≔ (1, β _x). Thereafter, we present computations considering the problem on Q T ⊂ R 3 . For all the test cases, we mainly use Q k with k = 1, local polynomial spaces, and in a few examples we provide results using k = 2 local polynomial spaces, see (3.12).

The examples have been solved on a series of mesh refinement levels, s = 0, 1, 2, … with h _s, h _s+1, …, where the asymptotic convergence behavior of the error ‖u − u _h‖_h is investigated. For the first problem the behavior of the error ‖ u − u h ‖ t , h ≔ ( ‖ u − u h ‖ L 2 2 + τ λ ‖ ∂ t u − ∂ t u h ‖ L 2 2 ) 1 / 2 is also studied. Note that a separation of the magnitudes of the two errors is not completely supported by the discretization analysis. We just provide this numerical test which investigates if ‖u − u _h‖_t,h and ‖u − u _s‖_s follow the same asymptotic behavior. Also, we mention that the “expected values” of the rates, which are written in the tables, are the values resulting from the discretization error analysis. In the numerical computations the initial condition u ₀ and the boundary data u _Σ are determined by the L ²-projection of the exact solution u onto polynomial space.

The linear system produced by the method (3.18a) is (in general) not symmetric and a direct LU method has been used to solve it. One can also apply GMRES iterative solvers for its solution. However, several efficient methods have been presented in the literature for speeding up the solution process of the system. We refer to refs. [19], [21], [22], [37] for discussions on developing algebraic multigrid methods and on different parallelization approaches for STFEMs.

The 2-d numerical examples have been performed using an in-house code which has been implemented on a Intel(R) Core(TM) i7-8700 CPU with Gentoo Linux optimized system.

The conclusion from the results presented below is that the proposed space-time FE scheme is stable, behaves well and the numerical convergence rates are in agreement with the theoretically predicted rates.

4.1 Two-dimensional space-time cylinders

4.1.1 Smooth solution, uniform meshes

In the first numerical example the domain is Q ̄ T = [ 0,1 ] × [ 0,1 ] and the exact solution is u(t, x) = sin(2πt) sin(2πx). The problem has been solved for the problem parameters {ɛ, β _x, r} = {0.1, 1, 1}. In Table 1, the results of the asymptotic convergence rates are displayed. The exact solution is smooth and optimal convergence rates are expected. We observe that the rates r _t and r _x have similar behavior and are a little higher than the expected rates in the first coarse meshes. As the meshes are progressively refined, both r _t and r _x tend to the expected values. In the last columns we can see the rates r _x,k=2, which correspond to Q k = 2 polynomial spaces. Moving to the finer meshes the rates have the expected values.

Table 1:

Example 1: smooth test case. The convergence rates r _x and r _t and r _x,k=2.

u smooth, with {ɛ, β _x, r} = {0.1, 1, 1}
Errors	‖u − u _h‖_s	‖u − u _h‖_t,h	‖u − u _h‖_s
Q k -space	k = 1	k = 1	k = 2
Expected rates	1	1	2
h ₀ = 0.2	Computed rates
h s = h 0 2 s	r _x	r _t	r _x,k=2
s = 1	3.20	2.84	4.2
s = 2	1.40	1.23	2.4
s = 3	1.40	1.08	2.14
s = 4	1.14	1.08	2.15
s = 5	1.17	1.02	2.01
s = 6	1.05	1.01	2.00
s = 7	1.00	1.04	2.02

4.1.2 Point singularity test case

We consider the problem on Q ̄ T = [ 0,1 ] × [ 0,1 ] with exact solution u ( x , t ) = ( ( x − 0.2 ) 2 + ( t − 0.2 ) 2 ) γ / 2 with γ ∈ {1.001, 0.51}. The singular point of the solution is located at the interior of the domain Q ̄ T . The problem has been solved for two sets of values for the parameters, i.e., {ɛ, β _x, r} = {0.1, 1, 1} and {ɛ, β _x, r} = {0.1, 10⁻¹⁰, 0.5}. In every case the associated solution u belongs to H ^ℓ=γ+1(Q _T) and the computational tests have been performed using Q k = 1 and Q k = 2 local polynomial spaces on uniform mesh partitions T h . In view of the interpolation results (3.56) and (3.63), see also Corollary 3.4, and the regularity of the solution u, the convergence rates expected to get values close to γ = ℓ − 1. We compute the rates on a sequence of meshes and present the results in Table 2. Note that the cases, where we use γ + 1 < 2 and Q k = 2 spaces, are not covered by the discretization analysis presented in the previous section. The associated results are added in Table 2 for comparison purposes only and for a more comprehensive computational study of the problem. The “expected values” were formed by a simple formalistic interpretation of the interpolation estimates.

Table 2:

Example 2: point singularity case. Convergence rates r _x and r _x,k=2.

u ∈ W ^ℓ,2(Q _T) with ℓ = γ + 1
Errors	‖u − u _h‖_s	‖u − u _h‖_s	‖u − u _h‖_s	‖u − u _h‖_s	‖u − u _h‖_s
Q k -space	k = 1	k = 2	k = 1	k = 2	k = 1
Parameters	γ = 1.001		γ = 1.001		γ = 0.6
	{ɛ, β _x, r} = {0.1, 0.5, 1}		{ɛ, β _x, r} = {0.1, 10⁻¹⁰, 0.5}		{ɛ, β _x, r} = {0.1, 0.5, 0.5}
Expected rates	1	1	1	1	0.5
h ₀ = 0.2	Computed rates
h s = h 0 2 s	r _x	r _x,k=2	r _x	r _x,k=2	r _x
s = 1	0.5	0.98	0.5	1.5	0.44
s = 2	0.95	1.04	0.92	1.82	0.63
s = 3	0.90	0.98	0.90	1.45	0.57
s = 4	0.89	0.94	0.89	1.34	0.55
s = 5	0.89	0.91	0.89	0.87	0.54
s = 6	0.91	0.90	0.90	1.25	0.52
s = 7	0.92	0.89	0.92	1.24	0.51

Looking at the table, we observe that for the tests with γ = 1.001, β _x = 0.5, see first columns, the rates r _x and r _x,k=2 have similar behavior and are close to the expected values, even for the first refinement steps. For the second test where we set β = 10⁻¹⁰, the values of r _x,k=2 are little higher during the first meshes but tend to get the expected values during the last meshes. In the last test we set γ = 0.6, see last column, the r _x rates are higher than the expected values on the first meshes. Moving to the last refinement steps the rates tend to get the expected values.

4.1.3 Anisotropic meshes

In this test, we investigate the behavior of the on anisotropic meshes. To illustrate this consider a problem where its solution has anisotropic behavior, i.e., its variations are quite different in each different axial direction. We solve the problem using anisotropic meshes, where the directional size of the mesh elements is small in the direction where the solution has strong variations, and the size of the elements is relatively larger in the direction where the solution has less variations. The domain is chosen to be Q ̄ T = [ 0,1 ] × [ 0,1 ] and the parameters take the values {ɛ, β _x, r} = {0.1, 1, 1}. The exact solution is u(x, t) = sin(πx) sin(4πt) and varies ’four times more’ in the direction of the time axis. The problem is solved in a series of successive refined mesh levels s = 1, 2, …, where for each refinement mesh level we set h _s,x = λh _s,t, here the parameter λ ∈ {1, 2, 4} forms the anisotropic character of the mesh. For each mesh level s, three different computations have been performed, where each of which corresponds to a different value of the parameter λ. We note that, for a better comparison of the results, the three refinement procedures start from a different initial mesh.

Our goal is to investigate the asymptotic behavior of the rates r _x, how they are affected by the anisotropic variational properties of u and the anisotropic character of the mesh. Also we investigate the computing resources needed to reduce the error to the same levels as those which result when isotropic meshes are used. Table 3 shows the results of the numerical convergence rates. For each λ-case, we can observe that the values of r _x rates are a little higher during the first meshes than the theoretically predicted values. However, moving to the last refinement steps, the rates reduce and get the expected values derived by the error analysis, see also Example 2.

Table 3:

Example 3: anisotropic meshes. The convergence rates r _x for all λ.

u(x, t) smooth, {ɛ, β _x, r} = {0.1, 1, 1}
Errors	‖u − u _h‖_s	‖u − u _h‖_s	‖u − u _h‖_s
Parameter	λ = 1	λ = 2	λ = 4
Expected rates	1	1	1
Ref. step s	Computed rates
	r _x	r _x	r _x
s = 1	0.5	2.01	1.8
s = 2	2.3	1.87	1.62
s = 3	1.7	1.61	1.85
s = 4	1.31	1.27	1.33
s = 5	1.11	1.09	1.11
s = 6	1.03	1.02	1.03
s = 7	1.01	1.00	1.2
s = 8	1.02	1.00	1.00
s = 9	1.01	1.02	1.02
s = 10	1.00	1.00	1.01

The left part of Figure 1 shows the reduction of the error of the numerical solution relative to the number of elements for each anisotropic mesh. Starting with a mesh with h _x = h _t, we move to next refinement steps, s = 2, 3, …, such that the directional mesh sizes have the fixed relation h _s,x = λh _s,t, λ = 1, 2, 4. During the first meshes, the values of the corresponding errors are very close for the same number of elements. However, looking the values in the next refinement levels, we observe that the λ = 4, λ = 2 meshes require a significantly reduced number of elements compared to the λ = 1 case, in order to reduce the error to the same level, (see the micro box in the graph). Thus, to reduce the error to a desired level, a uniform refinement with h _s,x ≈ h _s,t (isotropic) is not necessary. This also illustrates the usefulness of the anisotropic meshes when solving such problems. This can be seen more clearly in the right graph in Figure 1, where the error reduction against CPU time is shown. Note that in every case, the corresponding CPU time considers the whole performance-solution of the numerical example and not only the solution of the linear system. As expected the performance of the test with λ = 4 requires the least computing resources, when we refer to the same mesh refinement step, i.e., the same h _s,t. It is observed that the error, which corresponds to the numerical solution of λ = 4 case, decreases at the same level and at the same rate with what corresponds to the isotropic case, i.e., the λ = 1 case. Figure 1 somehow indicates that an appropriate use of anisotropic meshes can lead to high resolution numerical solutions at low computational cost.

Figure 1:

Example 3: anisotropic meshes, (a) the graph on the left shows the reduction of the error against the number of elements for each anisotropic mesh. (b) The discretization error against the CPU time for each λ test case.

4.2 Examples on three-dimensional space-time cylinders

4.2.1 Smooth solution, isotropic meshes

The purpose of this example is to investigate the convergence behavior of the discretization error ‖u − u _h‖_s for the case of having a three-dimensional space-time cylinder. We consider the problem on Q ̄ T = [ 0,1 ] 3 with exact solution u(x, y, t) = cos(2πx) cos(2πy) cos(2πt). We solve the problem using k = 1 and k = 2 polynomial spaces for two groups of parameter values, the first is {ɛ, β _x, β _y, r} = {0.1, 1, 1, 1} and the second {ɛ, β _x, β _y, r} = {0.1, 0, 0, 1}. We compare the numerical results with the theoretical findings given in Theorem 3.2, see also (3.66). The numerical results are reported in Table 4. For each k = 1-case, we observe that the convergence rates r _x are close to the expected rate value equal to one. The rates r _x,k=2 are shown in third and the fifth column for the two parameter test cases. The rates are a little lower in the first meshes but moving to the last meshes they have the expected value (with respect to the regularity of the solution) and follow the theoretical convergence rates, compare to the “Smooth solution, uniform meshes” above.

Table 4:

Example 4: smooth solution for Q T ⊂ R 3 . Convergence rates r _x and r _x,k=2.

u smooth, Q T ⊂ R 3
Errors	‖u − u _h‖_s	‖u − u _h‖_s	‖u − u _h‖_s	‖u − u _h‖_s
Parameters	{ɛ, β _x, β _y, r} = {0.1, 1, 1, 1}		{ɛ, β _x, β _y, r} = {0.1, 0, 0, 1}
Expected rates	1	2	1	2
Q k	Computed rates
	k = 1	k = 2	k = 1	k = 2
h s = 1 2 s	r _x	r _x,k=2	r _x	r _x,k=2
s = 1	0.554	1.825	0.490	1.873
s = 2	1.055	1.985	1.026	1.952
s = 3	1.08	1.984	1.049	1.967
s = 4	1.058	1.987	1.037	1.976
s = 5	1.041	1.989	1.027	1.983
s = 6	1.030	1.992	1.020	1.987
s = 7	1.023	1.993	1.016	1.990

4.2.2 An example using anisotropic meshes

As it has been illustrated by the examples above that an appropriate construction and use of anisotropic meshes can help on relaxing the high time requirements for the computation of the solution. Also, we have seen the use of the anisotropic meshes preserves the convergence properties of the numerical solution. Similar points are going to be investigated in this numerical example. The problem is considered in Q _T = [0,1]³ with exact solution u ( t , x , y ) = sin ( 2 π t ) e − ( x − 0.5 ) 2 / 0.01 y and the parameters are {ϵ, β _x, β _y, r} = {0.1, 1, 1, 1}. The solution has different variational behavior in each direction and this supports the use of anisotropic meshes. The problem is solved using k = 1-local spaces applying an anisotropic mesh refinement strategy, which is specified by the relations h _t = ϰh _y, ϰ ∈ {1, 1/2} and h _x = λh _y, λ ∈ {1, 1/4}. The first starting mesh contains equal directional mesh sizes, i.e., h _y = h _x = h _t, and for the next refinement levels the mesh sizes h _x and h _t in x and t directions respectively are determined by the values of ϰ and λ. The associated convergence rates r _x of the error are shown in Table 5. The first column includes the rates which are computed by applying global uniform refinement (isotropic case). The next two columns include the rates for the anisotropic meshes. It can be seen that the rates are very close to the (theoretical) predicted values for all cases. Figure 2(a) illustrates the error variation with respect to the number of mesh elements for the three different meshes. It can be seen that the errors corresponding to the anisotropic meshes have lower values compared to the errors computed using the isotropic meshes (with the same number of elements). Hence, it appears that the anisotropic mesh computations manage more appropriately the number of elements, in the sense that more elements are stacked in the direction where the solution has the most abrupt variation. On the right graph in Figure 2(b) the CPU times are plotted against the error reduction. We observe that the anisotropic meshes bring significant improvements in the time performance of the method. The anisotropic meshes require less computational resources than the isotropic meshes to keep the error at the same low level.

Table 5:

Example 5: anisotropic meshes Q T ⊂ R 3 . Convergence rates.

u = sin ( 2 π t ) e − ( x − 0.5 ) 2 / 0.01 y , {ϵ, β _x, β _y, r} = {0.1, 1, 1, 1}
Errors	‖u − u _h‖_s
Parameters	λ = 1, ϰ = 1	λ = 1/4, ϰ = 1	λ = 1/4, ϰ = 1/2
Expected rates	1	1	1
Ref. step s	Computed rates
	r _x	r _x	r _x
s = 1	0.927	0.917	0.943
s = 2	0.958	0.957	0.973
s = 3	0.968	0.967	0.983
s = 4	0.975	0.982	0.988
s = 5	0.979	0.986	0.991
s = 6	0.982	0.988	0.993
s = 7	0.984	0.990	0.994
s = 8	0.986	0.991	0.995
s = 9	0.988	0.992	0.996

$Figure 2: Example 5: anisotropic meshes for Q T ⊂ R 3 ${Q}_{T}\subset {\mathbb{R}}^{3}$ . (a) Left: the variation of the error with respect the number of elements for each mesh, (b) right: the discretization error against the CPU time for each mesh test case.$

Figure 2:

Example 5: anisotropic meshes for Q T ⊂ R 3 . (a) Left: the variation of the error with respect the number of elements for each mesh, (b) right: the discretization error against the CPU time for each mesh test case.

5 Conclusions

In this work space-time FE methods have been developed and analyzed with continuous spaces on anisotropic quadrilateral meshes for solving general linear parabolic problems. The scheme was stabilized following usual upwind streamline methodology. Discretization error estimates were shown in a suitable norm. The proposed method was applied to problems having regular and less regular solutions on isotropic and anisotropic meshes. The numerical convergence rates were in agreement with the theoretical rates.

The proposed scheme can be extended to the case of using discontinuous Galerkin discretizations in time. This can help on solving the problem in a sequential manner, i.e., one space-time slice at a time, [38]. This approach can be further combined with time-Domain Decomposition (DD) iterative solvers materialized in a parallel environment. A incorporation of anisotropic refinement strategies to the proposed method can lead to en efficient method for solving problems with solutions with anisotropic behavior, e.g., boundary layers in fluid problems, re-entrant edges, etc. The development of this type of numerical methods is the subject of a work in progress.

Corresponding author: Ioannis Toulopoulos, Department of Informatics, University of Western Macedonia, Fourka Area, 52100 Kastoria, Greece, E-mail: itoulopoulos@uowm.gr

Acknowledgments

The author would like to thank Oliver Koch from the Institute of Computational Mathematics of Johannes Kepler University Linz for his assistance in computing the numerical results.

Research ethics: Not applicable.
Informed consent: Not applicable.
Author contributions: The author has accepted responsibility for the entire content of this manuscript and approved its submission.
Use of Large Language Models, AI and Machine Learning Tools: None declared.
Conflict of interest: The author states no conflict of interest.
Research funding: None declared.
Data availability: Not applicable.

References

[1] L. C. Evans, Partial Differential Equations, ser. Graduate Studies in Mathematics, vol. 19, 2nd ed. USA, American Mathematical Society, 2010.Suche in Google Scholar

[2] H.-G. Roos, M. Stynes, and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations. Convection–Diffusion–Reaction and Flow Problems, ser. Springer Series in Computational Mathematics, vol. 24, Berlin, Heidelberg, Springer, 2008.Suche in Google Scholar

[3] L. P. Franca, G. Hauke, and A. Masud, “Revisiting stabilized finite element methods for the advective–diffusive equation,” Comput. Methods Appl. Mech. Eng., vol. 195, nos. 13–16, pp. 1560–1572, 2006. https://doi.org/10.1016/j.cma.2005.05.028.Suche in Google Scholar

[4] A. N. Brooks and T. J. Hughes, “Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations,” Comput. Methods Appl. Mech. Eng., vol. 32, no. 1, pp. 199–259, 1982. https://doi.org/10.1016/0045-7825(82)90071-8.Suche in Google Scholar

[5] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Computational Differential Equations, Cambridge, U.K., Cambridge University Press, 1996.Suche in Google Scholar

[6] M. Augustin, et al.., “An assessment of discretizations for convection-dominated convection–diffusion equations,” Comput. Methods Appl. Mech. Eng., vol. 200, no. 47, pp. 3395–3409, 2011. https://doi.org/10.1016/j.cma.2011.08.012.Suche in Google Scholar

[7] V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, ser. Springer Series in Computational Mathematics, 2nd ed. Berlin, Heidelberg, Springer-Verlag, vol. 25, 2006.Suche in Google Scholar

[8] W. Huang, L. Kamenski, and J. Lang, “Stability of explicit Runge–Kutta methods for high order finite element approximation of linear parabolic equations,” Numer. Math. Adv. Appl., vol. 103, pp. 165–173, 2015, https://doi.org/10.1007/978-3-319-10705-9_16.Suche in Google Scholar

[9] E. Burman and A. Ern, “Implicit-explicit Runge–Kutta schemes and finite elements with symmetric stabilization for advection-diffusion equations,” ESAIM Math. Model. Numer. Anal., vol. 46, no. 4, pp. 681–707, 2012. https://doi.org/10.1051/m2an/2011047.Suche in Google Scholar

[10] U. Langer and O. Steinbach, Space-Time Methods: Applications to Partial Differential Equations, ser. Radon Series on Computational and Applied Mathematics, vol. 25, U. Langer and O. Steinbach, Walter de Gruyter GmbH and Co KG, 2019.Suche in Google Scholar

[11] R. E. Bank, P. S. Vassilevski, and L. T. Zikatanov, “Arbitrary dimension convection–diffusion schemes for space-time discretizations,” J. Comput. Appl. Math., vol. 310, pp. 19–31, 2017, https://doi.org/10.1016/j.cam.2016.04.029.Suche in Google Scholar

[12] S. Larsson and M. Molteni, “Numerical solution of parabolic problems based on a weak space-time formulation,” Comput. Methods Appl. Math., vol. 17, no. 1, pp. 65–84, 2016. https://doi.org/10.1515/cmam-2016-0027.Suche in Google Scholar

[13] C. Mollet, “Stability of Petrov–Galerkin discretizations: application to the space-time weak formulation for parabolic evolution problems,” Comput. Methods Appl. Math., vol. 14, no. 2, pp. 231–255, 2014. https://doi.org/10.1515/cmam-2014-0001.Suche in Google Scholar

[14] C. Schwab and R. Stevenson, “Space-time adaptive wavelet methods for parabolic evolution problems,” Math. Comput., vol. 78, no. 267, pp. 1293–1318, 2009. https://doi.org/10.1090/s0025-5718-08-02205-9.Suche in Google Scholar

[15] R. Stevenson and J. Westerdiep, “Stability of Galerkin discretizations of a mixed space-time variational formulation of parabolic evolution equations,” IMA J. Numer. Anal., vol. 41, no. 1, pp. 28–47, 2020. https://doi.org/10.1093/imanum/drz069.Suche in Google Scholar

[16] I. Toulopoulos, “Space-time finite element methods stabilized using bubble function spaces,” Appl. Anal., vol. 99, no. 7, pp. 1153–1170, 2018. https://doi.org/10.1080/00036811.2018.1522630.Suche in Google Scholar PubMed PubMed Central

[17] I. Toulopoulos, “Numerical solutions of quasilinear parabolic problems by a continuous space-time finite element scheme,” SIAM J. Sci. Comput., vol. 44, no. 5, pp. A2944–A2973, 2022. https://doi.org/10.1137/21m1403722.Suche in Google Scholar

[18] A. Mantzaflaris, F. Scholz, and I. Toulopoulos, “Low-rank space-time decoupled isogeometric analysis for parabolic problems with varying coefficients,” Comput. Methods Appl. Math., vol. 19, no. 1, pp. 123–136, 2019. https://doi.org/10.1515/cmam-2018-0024.Suche in Google Scholar

[19] C. Hofer, U. Langer, M. Neumüller, and I. Toulopoulos, “Time-multipatch discontinuous Galerkin space-time isogeometric analysis of parabolic evolution problems,” Electron. Trans. Numer. Anal., vol. 49, pp. 126–150, 2018.Suche in Google Scholar

[20] T. E. Tezduyar and K. Takizawa, “Space-time computations in practical engineering applications: a summary of the 25-year history,” Comput. Mech., vol. 63, no. 4, pp. 747–753, 2019. https://doi.org/10.1007/s00466-018-1620-7.Suche in Google Scholar

[21] R. Dyja, B. Ganapathysubramanian, and K. G. van der Zee, “Parallel-in-space-time, adaptive finite element framework for nonlinear parabolic equations,” SIAM J. Sci. Comput., vol. 40, no. 3, pp. C283–C304, 2018. https://doi.org/10.1137/16m108985x.Suche in Google Scholar

[22] M. Neumüller and I. Smears, “Time-parallel iterative solvers for parabolic evolution equations,” SIAM J. Sci. Comput., vol. 41, no. 1, pp. C28–C51, 2018.Suche in Google Scholar

[23] C. V. Frontin, G. S. Walters, F. D. Witherden, C. W. Lee, D. M. Williams, and D. L. Darmofal, “Foundations of space-time finite element methods: polytopes, interpolation, and integration,” Appl. Numer. Math., vol. 166, pp. 92–113, 2021, https://doi.org/10.1016/j.apnum.2021.03.019.Suche in Google Scholar

[24] T. G. Lube, “Anisotropic mesh refinement in stabilized Galerkin methods,” Numer. Math., vol. 74, no. 3, pp. 261–282, 1996. https://doi.org/10.1007/s002110050216.Suche in Google Scholar

[25] T. Apel, “Interpolation of non-smooth functions on anisotropic finite element meshes,” Math. Model. Numer. Anal., vol. 33, no. 6, pp. 1149–1185, 1999. https://doi.org/10.1051/m2an:1999139.Suche in Google Scholar

[26] S. Micheletti, S. Perotto, and M. Picasso, “Stabilized finite elements on anisotropic meshes: a priori error estimates for the advection-diffusion and the Stokes problems,” SIAM J. Numer. Anal., vol. 41, no. 3, pp. 1131–1162, 2003. https://doi.org/10.1137/s0036142902403759.Suche in Google Scholar

[27] L. Formaggia, S. Micheletti, and S. Perotto, “Anisotropic mesh adaptation in computational fluid dynamics: application to the advection–diffusion–reaction and the Stokes problems,” Appl. Numer. Math., vol. 51, no. 4, pp. 511–533, 2004. https://doi.org/10.1016/j.apnum.2004.06.007.Suche in Google Scholar

[28] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, ser. Pure and Applied Mathematics, 2nd ed. Oxford U.K., Academic Press-Imprint Elsevier Science, vol. 140, 2003.Suche in Google Scholar

[29] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods, ser. Texts in Applied Mathematics, 2nd ed. New-York, Springer, 2008.Suche in Google Scholar

[30] O. Ladyzhenskaya and N. Uraltseva, Linear and Quasilinear Elliptic Equations, New York, London, Academic Press, 1968.Suche in Google Scholar

[31] A. Koshelev, “Regularity Problem for Quasilinear Elliptic and Parabolic Systems, ser. Lecture Notesin Mathematics, Berlin, Heidelberg, Springer-Verlag, 1995, p. 1614.Suche in Google Scholar

[32] R. Adreev, “Stability of sparse space-time finite element discretizations of linear parabolic evolution equations,” IMA J. Numer. Anal., vol. 33, no. 1, pp. 242–260, 2013. https://doi.org/10.1093/imanum/drs014.Suche in Google Scholar

[33] A. Ern and J. L. Guermond, One-Dimensional Finite Elements and Tensorization, ser. Texts in Applied Mathematics, Cham, Switzerland, Springer International Publishing, vol. 72, 2021, no. 1614.Suche in Google Scholar

[34] W. Hackbusch and B. N. Khoromskij, “Tensor-product approximation to operators and functions in high dimensions,” J. Complex., vol. 23, no. 4, pp. 697–714, 2007. https://doi.org/10.1016/j.jco.2007.03.007.Suche in Google Scholar

[35] L. L. Schumaker, Spline Functions: Basic Theory, 3rd ed. Cambridge, University Press, 2007.Suche in Google Scholar

[36] T. Dupont and R. Scott, “Polynomial approximation of functions in Sobolev spaces,” Math. Comput., vol. 34, no. 150, pp. 441–463, 1980. https://doi.org/10.2307/2006095.Suche in Google Scholar

[37] O. Steinbach and H. Yang, “Comparison of algebraic multigrid methods for an adaptive space-time finite-element discretization of the heat equation in 3D and 4D,” Numer. Lin. Algebra Appl., vol. 25, no. 3, 2018, https://doi.org/10.1002/nla.2143.Suche in Google Scholar

[38] I. Toulopoulos, “A unified time discontinuous Galerkin space-time finite element scheme for non-Newtonian power law models,” Int. J. Numer. Methods Fluid., vol. 95, no. 5, pp. 851–868, 2023. https://doi.org/10.1002/fld.5170.Suche in Google Scholar

Received: 2024-04-10

Accepted: 2024-12-06

Published Online: 2025-06-04

Published in Print: 2025-12-17

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

https://doi.org/10.1515/jnma-2024-0053

Schlagwörter für diesen Artikel

general linear parabolic equations; continuous space-time finite element discretizations; quadrilateral anisotropic meshes; a priori error estimates; asymptotic convergence rates

Creative Commons

BY 4.0