# Analysis of Backward Euler Primal DPG Methods

• Thomas Führer , Norbert Heuer and Michael Karkulik

## Abstract

We analyze backward Euler time stepping schemes for a primal DPG formulation of a class of parabolic problems. Optimal error estimates are shown in a natural norm and in the L 2 norm of the field variable. For the heat equation the solution of our primal DPG formulation equals the solution of a standard Galerkin scheme and, thus, optimal error bounds are found in the literature. In the presence of advection and reaction terms, however, the latter identity is not valid anymore and the analysis of optimal error bounds requires to resort to elliptic projection operators. It is essential that these operators be projections with respect to the spatial part of the PDE, as in standard Galerkin schemes, and not with respect to the full PDE at a time step, as done previously.

MSC 2010: 65N30; 65N12

## 1 Introduction

In this work we analyze a backward Euler primal DPG time stepping scheme for the parabolic problem

(1.1)

(1.1a) u ˙ - div 𝑨 u + 𝜷 u + γ u = f in  ( 0 , T ) × Ω ,
(1.1b) u = 0 in  ( 0 , T ) × Ω ,
(1.1c) u ( 0 , ) = u 0 in  Ω .

Here, Ω d ( d 2 , 3 ) denotes a bounded Lipschitz domain. We make the standing assumptions that the coefficients 𝑨 , 𝜷 , and γ are time-independent, that they are essentially bounded in space, that 𝑨 is a symmetric uniformly positive matrix almost everywhere in Ω, and that ( 𝜷 v + γ v , v ) 0 for all v H 0 1 ( Ω ) . The latter assumption plays a crucial role in our analysis (see Lemma 2 below). The L 2 ( Ω ) inner product is denoted with ( , ) . We assume that the data satisfy f ( t , ) L 2 ( Ω ) for all t [ 0 , T ] and u 0 L 2 ( Ω ) .

The discontinuous Petrov–Galerkin method with optimal test functions (DPG) pertains to the class of minimum residual methods and was introduced in a series of papers [5, 6, 10]. It has successfully been applied to elliptic problems, see, e.g., [7, 8] for the Poisson problem and [18, 16] for fourth-order problems. Through the use of optimal test functions, the discrete problem inherits the stability of the continuous problem. This comes in advantageous for problems where robustness is one of the main challenges, e.g., singularly perturbed problems, see [24] for reaction-dominated diffusion problems and [11, 3] for the convection-dominated case. Space-time DPG methods have been studied previously, see, e.g., [9, 13, 12]. For other space-time minimum residual methods we refer to [21, 1, 30, 28]. Approaches employing the DPG methodology for the time discretization of parabolic and hyperbolic initial value problems have recently been investigated, cf. [25, 26]. On the other hand, time-stepping methods for ODEs are frequently employed in combination with standard Galerkin finite element methods in space, cf. the monograph [29] for parabolic equations, but less so with DPG methods. To the best of our knowledge there exist only two works in this direction, dealing with time-stepping and spatial DPG methods for the heat equation, namely [19] and [27].

In [19], a backward Euler method is used to discretize in time, and then the DPG methodology is applied to the ultraweak variational formulation of the resulting equations. The a priori analysis given there employs the Galerkin projection with respect to these very equations, and hence a spatial discretization error has to be accounted for in every time step. This gives rise to a theoretical error bound of order 𝒪 ( h k ) + 𝒪 ( k ) for lowest-order discretizations with h and k being the spatial mesh-width and time step, respectively. As numerical experiments from [27] indicate, this asymptotic error bound is not optimal. The authors of [27] study general θ-schemes (including the backward Euler and the Crank–Nicolson time discretization) based on the primal DPG method [8] and the ultraweak DPG method [7], and provide an extensive numerical study and comparison of the different approaches, and it turns out that 𝒪 ( h + k ) is the optimal error bound which can be expected for the method from [19].

Our motivation for the present work is to give a theoretically sound explanation of the optimal convergence rates seen in the numerical experiments from [27] for a backward Euler primal DPG method. A DPG approximation is attractive, e.g., for singularly perturbed problems not considered here, cf. [11, 24]. An application for a primal DPG method can be found in the recent work [23]. We will consider general second order linear elliptic spatial differential operators. For the heat equation ( 𝑨 is the identity, 𝜷 = 0 , and γ = 0 ) optimal error estimates follow from the fact that the field solution component of the primal DPG method is identical to the solution of a standard Galerkin FEM (cf. Section 2.2.3 for details). Therefore, well-known results for time-stepping Galerkin FEMs, see e.g. [29], apply. However, for the general case 𝜷 0 this is not true anymore and a new analysis needs to be provided. The accumulation of the spatial discretization error in every time step, identified above as the reason for suboptimality of previous theoretical results, is usually avoided in standard FEMs by using the Galerkin projection with respect to the elliptic part of the parabolic equation only. This idea, nowadays being referred to as elliptic projection operators, was introduced in [31] and is by now one of the main tools in the analysis of time-stepping FEMs, as witnessed again in [29]. In the present work we prove optimal error estimates in the context of practical primal DPG methods by the use of an elliptic projection operator. We follow some ideas from [20] where an elliptic projection in the analysis of time-stepping first-order system least-squares finite element methods is used in the same spirit.

The remainder of this paper is organized as follows: In Section 2 we introduce the fully discrete method as well as the necessary notation. We prove stability of the method and provide quasi-optimality results for the elliptic projection operator. In Section 3 we use these results to show optimal error estimates in the H 1 ( Ω ) norm. Section 4 is devoted to optimal error estimates in weaker norms, particularly in L 2 ( Ω ) .

## 2 Time-Stepping DPG Formulation

### 2.1 Notation

The notation a b means that there exists a constant C > 0 that (possibly) depends on Ω, 𝑨 , 𝜷 , γ, and T, but is independent of involved functions. We write a b if a b and b a . Furthermore, a b stands for b a . If the dependence on T is not (at most) linear, we explicitly state the dependence. Moreover, we use the 𝒪 ( ) notation, e.g., g = 𝒪 ( h p ) with p > 0 means that | g | h p for 0 < h 0 .

We consider a time discretization 0 = t 0 < < t N = T , and for notational simplicity we use a uniform time step size k = t n - t n - 1 (but we stress that this is not necessary). By ( , ) and we denote the inner product and norm in L 2 ( Ω ) . We consider spatial discretizations based on a shape-regular conforming simplicial mesh 𝒯 of Ω. To any partition 𝒯 we associate its skeleton 𝒮 consisting of the boundaries of all elements, and the trace space

H - 1 2 ( 𝒮 ) := { σ ^ K 𝒯 H - 1 2 ( K ) | there exists  𝝈 𝑯 ( div ; Ω )  such that  σ ^ | K = 𝝈 𝐧 K | K  for all  K 𝒯 } ,

where 𝐧 K denotes the unit outward normal vector on K . This is a Hilbert space with norm

σ ^ - 1 2 , k := inf { ( 𝝈 2 + k div 𝝈 2 ) 1 2 | 𝝈 𝑯 ( div ; Ω ) , σ ^ | K = 𝝈 𝐧 K | K  for all  K 𝒯 } .

The trial space of our method will be

U := H 0 1 ( Ω ) × H - 1 2 ( 𝒮 )

equipped with the norm

( u , σ ^ ) U , k 2 := u 2 + σ ^ - 1 2 , k 2 ,

and the test space will be

V := H 1 ( 𝒯 ) := K 𝒯 H 1 ( K )

equipped with the norm

v V , k 2 := 1 k v 2 + 𝑨 1 2 𝒯 v 2 .

Here, 𝒯 v is the 𝒯 -piecewise gradient. For functions σ ^ H - 1 2 ( 𝒮 ) we define

σ ^ , v 𝒮 := ( div 𝝈 , v ) + ( 𝝈 , 𝒯 v ) for all  v H 1 ( 𝒯 ) ,

where 𝝈 𝑯 ( div ; Ω ) with 𝝈 𝒏 K | K = σ ^ | K for all K 𝒯 . We note that the definition is independent of the choice of 𝝈 and recall the following result.

### Lemma 1.

One has

σ ^ - 1 2 , k sup 0 v V σ ^ , v 𝒮 v V , k .

### Proof.

The proof follows by now well-established arguments, see [4] or [17] for problems where the norms depend on parameters. Note that if we use 𝑨 - 1 2 𝝈 instead of 𝝈 in the definition of - 1 2 , k , then equality holds. ∎

Throughout this work we will use that σ ^ , v 𝒮 = 0 for v H 0 1 ( Ω ) which follows from the definition of this duality and integration by parts, i.e.,

σ ^ , v 𝒮 = ( div 𝝈 , v ) + ( 𝝈 , 𝒯 v ) = ( div 𝝈 , v ) + ( 𝝈 , v ) = 𝝈 𝒏 Ω , v Ω = 0

for all σ ^ H - 1 2 ( 𝒮 ) , v H 0 1 ( Ω ) , where 𝝈 𝑯 ( div ; Ω ) is an extension of the trace σ ^ .

### 2.2 Primal DPG Formulation

For a simpler notation, we use superindices to indicate time evaluations, i.e., for a time-dependent function v = v ( t ; ) we use the notation v n ( ) := v ( t n ; ) . Moreover, we use similar notations for discrete quantities which however are not defined for all t [ 0 , T ] . For instance, later on u h n will denote an approximation of u n = u ( t n ; ) , where u is the exact solution of the parabolic equation. We approximate u ˙ n by a backward difference, i.e., u ˙ n ( u n - u n - 1 ) / k , this formally yields the elliptic PDE

1 k u n - div 𝑨 u n + 𝜷 u n + γ u n = f n + 1 k u n - 1 ,

which admits a unique solution u n H 0 1 ( Ω ) . By testing with v V and introducing σ ^ n as

σ ^ n | K := 𝑨 u n 𝐧 K | K for all  K 𝒯 ,

and integrating by parts one obtains the primal DPG formulation

1 k ( u n , v ) + ( 𝑨 u n , 𝒯 v ) + ( 𝜷 u n , v ) + ( γ u n , v ) - σ ^ n , v 𝒮 = ( f n , v ) + 1 k ( u n - 1 , v ) .

We introduce some bilinear forms and the right-hand side functional: For 𝑼 = ( u , σ ^ ) U , v V and w , g L 2 ( Ω ) set

b ( 𝑼 , v ) := ( 𝑨 u , 𝒯 v ) + ( 𝜷 u , v ) + ( γ u , v ) - σ ^ , v 𝒮 ,
a ( 𝑼 , v ) := 1 k ( u , v ) + b ( 𝑼 , v ) ,
F ( g , w ; v ) := ( g , v ) + 1 k ( w , v ) .

a ( 𝑼 n , v ) = F ( f n , u n - 1 ; v ) for all  v V .

We note that the bilinear form a ( , ) satisfies a continuous inf-sup condition using the norms U , k and V , k . This can be seen by extending the analysis of standard works like [4] to general elliptic PDEs. We do not pursue this further, as we ultimately want to analyze the practical DPG method. Furthermore, we note that a is not bounded independently of k using the above norms, but rather

(2.1) | a ( 𝑼 , v ) | ( k - 1 2 u + 𝑼 U , k ) v V , k ,

where 𝑼 = ( u , σ ^ ) .

#### 2.2.1 Fully Discrete Scheme

By 𝒫 p ( K ) we denote the space of polynomials of degree p 0 , and

𝒫 p ( 𝒯 ) := { v L 2 ( Ω ) : v | K 𝒫 p ( K )  for all  K 𝒯 }

For the discretization of traces we use the space of facewise polynomials denoted by 𝒫 p ( 𝒮 ) . In particular, we note that 𝒫 p ( 𝒮 ) is the (normal-)trace space of the p-th order Raviart–Thomas space. We consider the spaces

U h := 𝒫 p + 1 ( 𝒯 ) H 0 1 ( Ω ) × 𝒫 p ( 𝒮 ) ,
V h := 𝒫 p + d ( 𝒯 ) .

The (discrete) trial-to-test operator is given by

( Θ h 𝑼 h , v h ) V , k = a ( 𝑼 h , v h ) for all  v h V h .

We recall that the inner product is given by

( v , δ v ) V , k = 1 k ( v , δ v ) + ( 𝑨 𝒯 v , 𝒯 δ v ) .

The fully discrete scheme then reads: Given u h 0 L 2 ( Ω ) , solve

(2.2) a ( 𝑼 h n , Θ h 𝑾 h ) = F ( f n , u h n - 1 ; Θ h 𝑾 h ) for all  𝑾 h U h , n = 1 , 2 , .

The next lemma establishes coercivity of a ( , Θ h ) on U h .

#### Lemma 2.

Let 𝐔 h = ( u h , σ ^ h ) U h . Then

1 k u h 2 + 𝑨 1 2 u h 2 Θ h 𝑼 h V , k 2 .

If, additionally, h k - 1 2 C 1 for some constant C 1 > 0 , then

C 2 σ ^ h - 1 2 , k C 3 Θ h ( 0 , σ ^ h ) V , k Θ h 𝑼 h V , k .

Here, C 2 , C 3 > 0 denote generic constants independent of h , k and 𝐔 h .

#### Proof.

For the first part we use that 𝒫 p + 1 ( 𝒯 ) H 0 1 ( Ω ) V h and ( 𝜷 w + γ w , w ) 0 for w H 0 1 ( Ω ) : Let 𝑼 h = ( u h , σ ^ h ) U h be given. With the latter observation, the definition of the optimal test function and the fact that σ ^ h , w 𝒮 = 0 for w H 0 1 ( Ω ) we get that

1 k u h 2 + 𝑨 1 2 u h 2 1 k ( u h , u h ) + ( 𝑨 u h , u h ) + ( 𝜷 u h , u h ) + ( γ u h , u h )
= a ( ( u h , 0 ) , u h ) = ( Θ h 𝑼 h , u h ) V , k Θ h 𝑼 h V , k u h V , k
= Θ h 𝑼 h V , k ( 1 k u h 2 + 𝑨 1 2 u h 2 ) 1 2 .

This finishes the proof of the first estimate.

For the second one we use some results established in the literature on fully discrete DPG formulations (practical DPG), e.g., [22]. Let Π h : V V h be the Fortin operator defined in [22, Lemma 3.2] which in [22] is denoted by Π r grad with r = p + d . By Lemma 1 and the Fortin property we have that

(2.3) σ ^ h - 1 2 , k sup 0 v H 1 ( 𝒯 ) σ ^ h , v 𝒮 v V , k = sup 0 v H 1 ( 𝒯 ) σ ^ h , Π h v 𝒮 v V , k .

Moreover, from [22, proof of Lemma 3.2] we infer that

k - 1 2 Π h v + Π h v k - 1 2 v + ( 1 + h k - 1 2 ) v .

Thus, using the assumption h k - 1 2 1 , we see that

Π h v V , k v V , k for all  v V .

Combining this with (2.3) we conclude that

σ ^ h - 1 2 , k sup 0 v H 1 ( 𝒯 ) σ ^ h , Π h v 𝒮 v V , k sup 0 v H 1 ( 𝒯 ) σ ^ h , Π h v 𝒮 Π h v V , k
sup 0 v h V h σ ^ h , v h 𝒮 v h V , k = Θ h ( 0 , σ ^ h ) V , k .

Then Θ h ( 0 , σ ^ h ) V , k Θ h ( u h , σ ^ h ) V , k + Θ h ( u h , 0 ) V , k . It remains to estimate Θ h ( u h , 0 ) V , k : We use the notation Θ h ( u h , 0 ) = : v h ,

Θ h ( u h , 0 ) V , k 2 = a ( ( u h , 0 ) , Θ h ( u h , 0 ) ) = 1 k ( u h , v h ) + ( 𝑨 u h , 𝒯 v h ) + ( 𝜷 u h + γ u h , v h )
( 1 k u h 2 + 𝑨 1 2 u h 2 ) 1 2 v h V , k Θ h 𝑼 h V , k Θ h ( u h , 0 ) V , k .

Note that the involved constant only depends on the end time T and the coefficients, but is otherwise independent of k or h. To see the possible dependence on T we note that the triangle and Poincaré inequalities as well as boundedness of the coefficients imply that

( 𝜷 u h + γ u h , v h ) u h v h = k 1 2 u h k - 1 2 v h T 1 2 u h v h V , k .

This finishes the proof. ∎

It should be noted that the condition h k - 1 2 1 is only needed to get a k independent estimate for the trace norm. If this condition is not satisfied then the constant depends on k.

#### Theorem 3.

Problem (2.2) is well posed. In particular, the solutions are stable in the sense that

u h n ( u h n 2 + k 𝑨 1 2 u h n 2 ) 1 2 j = 1 n k f j + u h 0 .

#### Proof.

According to Lemma 2, Problem (2.2) admits unique solutions 𝑼 h n U h , n = 1 , 2 , . Then

Θ h 𝑼 h n V , k 2 = a ( 𝑼 h n , Θ h 𝑼 h n ) = F ( f n , u h n - 1 ; Θ h 𝑼 h n ) = ( f n + k - 1 u h n - 1 , Θ h 𝑼 h n )
k 1 2 f n + k - 1 u h n - 1 k - 1 2 Θ h 𝑼 h n k 1 2 f n + k - 1 u h n - 1 Θ h 𝑼 h n V , k .

This together with the estimate from Lemma 2 shows that

( u h n 2 + k 𝑨 1 2 u h n 2 ) 1 2 k 1 2 Θ h 𝑼 h n V , k k f n + u h n - 1 .

Iterating the arguments concludes the proof. ∎

#### 2.2.2 Remark: Known Results for a Bubnov–Galerkin Method for the Heat Equation

We recall results from [29, Section 1] for a standard Bubnov–Galerkin backward Euler method: Given an initial solution u h 0 L 2 ( Ω ) , find u h n 𝒫 p + 1 ( 𝒯 ) H 0 1 ( Ω ) ( n = 1 , 2 , ) such that

(2.4) 1 k ( u h n , v h ) + ( u h n , v h ) = ( f n , v h ) + 1 k ( u h n - 1 , v h ) for all  v h 𝒫 p + 1 ( 𝒯 ) H 0 1 ( Ω ) .

#### Proposition 4 ([29, p. 16 and Theorem 1.5]).

Let u denote the solution of (1.1) and let u h n P p + 1 ( T ) H 0 1 ( Ω ) denote the solution of (2.4). Suppose that u 0 | Ω = 0 . Then

( u n - u h n ) C 1 h p + 1 ( u 0 H p + 2 ( Ω ) + 0 t n u ˙ ( t ) H p + 2 ( Ω ) d t ) + C 2 k 0 t n u ¨ ( t ) d t .

If Ω is convex then,

u n - u h n C 3 h p + 2 ( u 0 H p + 2 ( Ω ) + 0 t n u ˙ ( t ) H p + 2 ( Ω ) d t ) + k 0 t n u ¨ ( t ) d t .

The constants C j , j = 1 , 2 , 3 , are independent of the solution u and T.

#### 2.2.3 Remark: Primal DPG for the Heat Equation — The Trivial Case

Let us consider the simplest model which is the heat equation where 𝑨 is the identity, 𝜷 = 0 , γ = 0 .

It is straightforward to see that for ( u h , 0 ) U h , Θ h ( u h , 0 ) = u h since

( u h , v h ) V , k = 1 k ( u h , v h ) + ( u h , 𝒯 v h ) = a ( ( u h , 0 ) , v h ) = ( Θ h ( u h , 0 ) , v h ) V , k for all  v h V h .

For any v h 𝒫 p + 1 ( 𝒯 ) H 0 1 ( Ω ) consider the test function Θ h ( v h , 0 ) in (2.2). Using the fact that σ ^ h , v h 𝒮 = 0 , we see that

a ( 𝑼 h n , Θ h ( v h , 0 ) ) = 1 k ( u h n , v h ) + ( u h n , v h ) = ( f n + k - 1 u h n - 1 , v h ) = F ( f n , u h n - 1 ; Θ h ( v h , 0 ) ) .

Let 𝑼 h n = ( u h n , σ ^ h n ) U h denote the solution to (2.2), then the solution component u h n satisfies

1 k ( u h n , v h ) + ( u h n , v h ) = ( f n , v h ) + 1 k ( u h n - 1 , v h ) for all  v h 𝒫 p + 1 ( 𝒯 ) H 0 1 ( Ω ) ,

which is the standard Galerkin FEM, see [29, Section 1]. Thus the primal DPG solution component u h is identical to the standard Galerkin FEM solution. In particular, optimal error estimates in L ( L 2 ) and L ( H 0 1 ) , see [29], are valid for the primal DPG solution. In [27] optimal error estimates L ( H 0 1 ) have been observed in numerical experiments, which can be explained with the observation above.

Note that this remark is true only if we consider 𝜷 = 0 , γ = 0 .

### 2.3 Elliptic Projection-Type Operator

To obtain optimal error estimates, we introduce an elliptic projection. The idea goes back to [31] to obtain optimal L 2 ( Ω ) a priori error estimates and is extensively used for Galerkin methods, but has not been studied for least-squares methods until only very recently in [20]. For DPG methods, an additional difficulty arises since test norms are mesh dependent. A main difference to Galerkin methods is that, although the elliptic part of the parabolic PDE might be symmetric, our elliptic projection operator always corresponds to a non-symmetric problem. The reason is that we need to use optimal test functions used in (2.2) in combination with the bilinear form b ( , ) , which models the elliptic part.

We define the elliptic projection operator h : U U h by

(2.5) b ( h 𝑼 , Θ h 𝑾 h ) = b ( 𝑼 , Θ h 𝑾 h ) for all  𝑾 h U h .

In the proofs below and in some results we will use the (semi-)norm

|| | 𝑼 | || k 2 := 𝑨 1 2 u 2 + Θ h ( 0 , σ ^ ) V , k 2 .

Note that from boundedness of 𝑨 and the definition of the discrete trial-to-test operator we infer that

|| | 𝑼 | || k 𝑨 1 2 u + Θ h ( 0 , σ ^ ) V , k = 𝑨 1 2 u + sup 0 v V h a ( ( 0 , σ ^ ) , v ) v V , k
u + sup 0 v V h σ ^ , v 𝒮 v V , k u + sup 0 v V σ ^ , v 𝒮 v V , k
u + σ ^ - 1 2 , k 𝑼 U , k ,

i.e., || | 𝑼 | || k 𝑼 U , k for 𝑼 U . Lemma 2 proves

σ ^ h - 1 2 , k Θ h ( 0 , σ ^ h ) V , k ,

hence, uniform positive definiteness of 𝑨 shows 𝑼 h U , k || | 𝑼 h | || k for 𝑼 h U h and the involved constants are independent of h , k if h k - 1 2 = 𝒪 ( 1 ) . In particular, the norm equivalence || | 𝑼 h | || k 𝑼 h U , k holds for 𝑼 h U h .

The next lemma establishes boundedness and an inf-sup condition of b ( , Θ h ) . Note that this is not as trivial as it seems since Θ h is calculated using the bilinear form a ( , ) and the inner product in V. Recall that a ( , ) is not bounded independently of k and that ( , ) V , k includes terms weighted with negative powers of the time step size k.

### Lemma 5.

The bilinear form b ( , Θ h ) is bounded,

| b ( 𝑼 , Θ h 𝑾 h ) | 𝑼 U , k 𝑾 h U , k ,
| b ( 𝑼 , Θ h 𝑾 h ) | 𝑼 U , k || | 𝑾 h | || k

for all 𝐔 U , 𝐖 h U h , and fulfills the inf-sup condition

|| | 𝑼 h | || k sup 0 𝑽 h U h b ( 𝑼 h , Θ h 𝑽 h ) || | 𝑽 h | || k

for all 𝐔 h U h . In particular, we conclude that Problem 2.5 admits a unique solution.

### Proof.

We show only the second boundedness estimate, as the first one follows from | | | | | | k U , k . To that end, let v := Θ h 𝑾 h = Θ h ( w h , χ ^ h ) . We use the splitting v = w h + w ~ h . From the definition of the optimal test function we see that

1 k ( w h , δ v ) + ( 𝑨 w h , 𝒯 δ v ) + 1 k ( w ~ h , δ v ) + ( 𝑨 𝒯 w ~ h , 𝒯 δ v )
= 1 k ( w h , δ v ) + ( 𝑨 w h , 𝒯 δ v ) + ( 𝜷 w h , δ v ) + ( γ w h , δ v ) - χ ^ h , δ v 𝒮 for all  δ v V h .

Therefore the component w ~ h V h satisfies

1 k ( w ~ h , δ v ) + ( 𝑨 𝒯 w ~ h , 𝒯 δ v ) = ( 𝜷 w h , δ v ) + ( γ w h , δ v ) - χ ^ h , δ v 𝒮 for all  δ v V h .

We note that a ( ( 0 , χ ^ h ) , δ v ) = - χ ^ h , δ v 𝒮 and therefore the definition of the optimal test function yields

| χ ^ h , δ v 𝒮 | Θ h ( 0 , χ ^ h ) V , k δ v V , k .

Standard estimates then show that

w ~ h V , k 2 k 1 2 𝜷 w h + γ w h k - 1 2 w ~ h + Θ h ( 0 , χ ^ h ) V , k w ~ h V , k .

Using that k T , boundedness of the coefficients and a Poincaré inequality, we infer that

(2.6) w ~ h V , k k 1 2 𝜷 w h + γ w h + Θ h ( 0 , χ ^ h ) V , k w h + Θ h ( 0 , χ ^ h ) V , k || | 𝑾 h | || k .

Observe that σ ^ h , w h 𝒮 = 0 since w h H 0 1 ( Ω ) . Finally, boundedness of the coefficients, a Poincaré inequality and the definition of | | | | | | k give

| b ( 𝑼 , v ) | | b ( 𝑼 , w h ) | + | b ( 𝑼 , w ~ h ) |
| ( 𝑨 u , w h ) + ( 𝜷 u + γ u , w h ) | + 𝑼 U , k w ~ h V , k
u w h + 𝑼 U , k w ~ h V , k 𝑼 U , k || | 𝑾 h | || k .

This shows boundedness. In order to show the inf-sup condition, we first establish the coercivity estimate

(2.7) 𝑨 1 2 u h 2 2 b ( 𝑼 h , Θ h 𝑼 h ) .

Recall that the optimal test function is characterized by

1 k ( Θ h 𝑼 h , δ v ) + ( 𝑨 𝒯 Θ h 𝑼 h , 𝒯 δ v ) = a ( 𝑼 h , δ v ) = 1 k ( u h , δ v ) + b ( 𝑼 h , δ v ) for all  δ v V h .

Setting δ v = v = Θ h 𝑼 h , the last identity together with Young’s inequality proves that

(2.8) b ( 𝑼 h , v ) = 1 k v 2 + 𝑨 1 2 𝒯 v 2 - 1 k ( u h , v ) 1 k v 2 + 𝑨 1 2 𝒯 v 2 - 1 2 k u h 2 - 1 2 k v 2 .

We are going to estimate k - 1 2 u h : From Lemma 2 we get that k - 1 u h 2 + 𝑨 1 2 u h 2 k - 1 v 2 + 𝑨 1 2 𝒯 v 2 , hence,

𝑨 1 2 u h 2 𝑨 1 2 𝒯 v 2 + 1 k v 2 - 1 k u h 2 .

Combining this with (2.8) yields

b ( 𝑼 h , v ) 1 k v 2 + 𝑨 1 2 𝒯 v 2 - 1 2 k u h 2 - 1 2 k v 2
= ( 1 2 𝑨 1 2 𝒯 v 2 + 1 2 k v 2 - 1 2 k u h 2 ) + 1 2 𝑨 1 2 𝒯 v 2
1 2 𝑨 1 2 𝒯 v 2 + 1 2 𝑨 1 2 u h 2 1 2 𝑨 1 2 u h 2 .

Now we are in a position to establish the inf-sup condition

|| | 𝑼 h | || k sup 0 𝑽 h U h b ( 𝑼 h , Θ h 𝑽 h ) || | 𝑽 h | || k .

With the boundedness estimates for b ( , ) we have that

Θ h ( 0 , σ ^ h ) V , k 2 = a ( ( 0 , σ ^ h ) , Θ h ( 0 , σ ^ h ) ) = b ( ( 0 , σ ^ h ) , Θ h ( 0 , σ ^ h ) )
= b ( 𝑼 h , Θ h ( 0 , σ ^ h ) ) - b ( ( u h , 0 ) , Θ h ( 0 , σ ^ h ) )
sup 0 𝑽 h U h b ( 𝑼 h , Θ h 𝑽 h ) || | 𝑽 h | || k Θ h ( 0 , σ ^ h ) V , k + 𝑨 1 2 u h Θ h ( 0 , σ ^ h ) V , k .

Applying the coercivity estimate (2.7) for the second term yields

𝑨 1 2 u h 2 b ( 𝑼 h , Θ h 𝑼 h ) = b ( 𝑼 h , Θ h 𝑼 h ) || | 𝑼 h | || k || | 𝑼 h | || k .

Combining the latter two estimates shows that

|| | 𝑼 h | || k 2 ( sup 0 𝑽 h U h b ( 𝑼 h , Θ h 𝑽 h ) || | 𝑽 h | || k ) 2 + sup 0 𝑽 h U h b ( 𝑼 h , Θ h 𝑽 h ) || | 𝑽 h | || k || | 𝑼 h | || k .

Young’s inequality finishes the proof of the inf sup condition. ∎

### Lemma 6.

Let E h 𝐔 := ( u h , σ ^ h ) U h be the solution of Problem (2.5). Then it holds

( u - u h ) || | 𝑼 - 𝑼 h | || k inf 𝑽 h U h 𝑼 - 𝑽 h U , k .

If, additionally, h k - 1 2 = O ( 1 ) (see Lemma 2), then

σ ^ - σ ^ h - 1 2 , k inf 𝑽 h U h 𝑼 - 𝑽 h U , k .

### Proof.

The best approximation properties follow from standard arguments (Babuška’s theorem). We give the details only for sake of completeness: Let 𝑾 h U h be arbitrary; then

|| | 𝑼 - 𝑼 h | || k || | 𝑼 - 𝑾 h | || k + || | 𝑼 h - 𝑾 h | || k ,

and using (2.5) and the preceding lemma,

|| | 𝑼 h - 𝑾 h | || k sup 0 𝑽 h U h b ( 𝑼 h - 𝑾 h , Θ h 𝑽 h ) || | 𝑽 h | || k sup 0 𝑽 h U h b ( 𝑼 - 𝑾 h , Θ h 𝑽 h ) || | 𝑽 h | || k 𝑼 - 𝑾 h U , k .

Since || | 𝑼 | || k u this shows the first of the claimed best-approximation estimates. The second follows by noting that σ ^ h - 1 2 , k Θ h ( 0 , σ ^ h ) V , k if h k - 1 2 = 𝒪 ( 1 ) (see Lemma 2) and considering

σ ^ - σ ^ h - 1 2 , k 𝑼 - 𝑼 h U , k 𝑼 - 𝑾 h U , k + 𝑼 h - 𝑾 h U , k 𝑼 - 𝑾 h U , k + || | 𝑼 h - 𝑾 h | || k .

The last term is handled as before which concludes the proof. ∎

Combining the latter quasi-best approximation result with standard approximation properties yields:

### Corollary 7.

Suppose that the components of 𝐔 = ( u , σ ^ ) are sufficiently smooth. Then

|| | 𝑼 - h 𝑼 | || k = 𝒪 ( h p + 1 ) .

If, additionally, h k - 1 2 = O ( 1 ) , then

𝑼 - h 𝑼 U , k = 𝒪 ( h p + 1 ) .

## 3 Optimal Error Estimate in Energy Norm

This section is devoted to proving optimal error estimates in the H 1 ( Ω ) norm of the n-th solution of the backward Euler method (2.2). We make the same assumptions on the regularity of solutions as in [29, Section 1], cf. Proposition 4.

## Theorem 8.

Let 𝐔 h n = ( u h n , σ ^ h n ) U h denote the solution of (2.2). Suppose that the components of 𝐔 ( t ; ) are sufficiently regular. Under these assumptions there exists k 0 (independent of h) such that for k < k 0 the solution satisfies

u n - u h n = 𝒪 ( h p + 1 + k ) + 𝒪 ( ( u 0 - u h 0 ) ) .

If, additionally, h k - 1 2 = O ( 1 ) , then

σ ^ n - σ ^ h n - 1 2 , k = 𝒪 ( h p + 1 + k ) + 𝒪 ( ( u 0 - u h 0 ) ) .

In both estimates, the dependence on T is exponential.

## Proof.

With the elliptic projection operator h (see (2.5)) we consider the splitting

𝑼 n - 𝑼 h n = ( 𝑼 n - h 𝑼 n ) + ( h 𝑼 n - 𝑼 h n ) .

We may also use the norm | | | | | | k defined in Section 2.3.

Step 1. By Corollary 7 we have

|| | 𝑼 n - h 𝑼 n | || k = 𝒪 ( h p + 1 ) .

Step 2. We derive error equations: First, write h 𝑼 n = ( h 1 𝑼 n , h 2 𝑼 n ) . Then, by (2.5),

a ( h 𝑼 n , v opt ) = 1 k ( h 1 𝑼 n , v opt ) + b ( h 𝑼 n , v opt ) = 1 k ( h 1 𝑼 n , v opt ) + b ( 𝑼 n , v opt )
= ( k - 1 h 1 𝑼 n + f n - u ˙ n , v opt ) for all  v opt Θ h ( U h ) .

Second, by (2.2),

a ( 𝑼 h n , v opt ) = ( f n + k - 1 u h n - 1 , v opt ) for all  v opt Θ h ( U h ) .

Third, combining both identities and writing 𝑾 n = ( w n , χ ^ n ) := h 𝑼 n - 𝑼 h n yields

(3.1) 1 k ( w n , v opt ) + b ( 𝑾 n , v opt ) = 1 k ( e h n + w n - 1 , v opt ) for all  v opt Θ h ( U h ) ,

where

(3.2) e h n := h 1 ( 𝑼 n - 𝑼 n - 1 ) - k u ˙ n .

Putting the term with w n - 1 on the left-hand side yields

(3.3) 1 k ( w n - w n - 1 , v opt ) + b ( 𝑾 n , v opt ) = 1 k ( e h n , v opt ) for all  v opt Θ h ( U h ) ,

Step 3. We use the test function v opt = Θ h ( v h , 0 ) = v h + v ~ h with v h = w n - w n - 1 . Reordering the terms in the error equations (3.3),

(3.4) 1 k w n - w n - 1 2 + 1 2 𝑨 1 2 w n 2 - 1 2 𝑨 1 2 w n - 1 2 + 1 2 𝑨 1 2 ( w n - w n - 1 ) 2 = 1 k ( e h n , v h + v ~ h ) - 1 k ( w n - w n - 1 , v ~ h ) - b ( 𝑾 n , v ~ h ) - ( 𝜷 w n , v h ) - ( γ w n , v h ) = r 1 + r 2 + r 3 + r 4 + r 5 .

Step 4. We estimate the contributions r j of the right-hand side of the error equation (3.4). Recall that v h = w n - w n - 1 . Throughout we use Young’s inequality with a parameter δ > 0 and the estimate from (2.6), i.e., k - 1 2 v ~ h v ~ h V , k k 1 2 ( w n - w n - 1 ) . First,

| r 1 | = 1 k | ( e h n , v h + v ~ h ) | k - 1 2 e h n k - 1 2 v h + k - 1 2 e h n k - 1 2 v ~ h
δ - 1 k - 1 e h n 2 + δ k - 1 w n - w n - 1 2 + δ k 𝑨 1 2 ( w n - w n - 1 ) 2 .

Second,

| r 2 | k - 1 2 w n - w n - 1 k - 1 2 v ~ h δ k - 1 w n - w n - 1 2 + δ - 1 k 𝑨 1 2 ( w n - w n - 1 ) 2 .

Third, recall the definition of the (discrete) optimal test function, i.e.,

| χ ^ n , v ~ h 𝒮 | = | a ( ( 0 , χ ^ n ) , v ~ h ) | = | ( Θ h ( 0 , χ ^ n ) , v ~ h ) V , k |

which together with the Cauchy–Schwarz inequality shows that

| r 3 | = | b ( 𝑾 n , v ~ h ) | ( w n + Θ h ( 0 , χ ^ n ) V , k ) k 1 2 v h
( w n - w n - 1 ) k 1 2 v h + k 1 2 w n - 1 v h + k 1 2 Θ h ( 0 , χ ^ n ) V , k v h .

We need to estimate Θ h ( 0 , χ ^ n ) V , k . To this end, using the error equation (3.1)

1 k ( w n , Θ h ( 0 , χ ^ n ) ) + b ( 𝑾 n , Θ h ( 0 , χ ^ n