Skip to content
BY 4.0 license Open Access Published by De Gruyter October 21, 2020

Optimal Convergence Rates for Goal-Oriented FEM with Quadratic Goal Functional

  • Roland Becker , Michael Innerberger EMAIL logo and Dirk Praetorius

Abstract

We consider a linear elliptic PDE and a quadratic goal functional. The goal-oriented adaptive FEM algorithm (GOAFEM) solves the primal as well as a dual problem, where the goal functional is always linearized around the discrete primal solution at hand. We show that the marking strategy proposed in [M. Feischl, D. Praetorius and K. G. van der Zee, An abstract analysis of optimal goal-oriented adaptivity, SIAM J. Numer. Anal. 54 (2016), 3, 1423–1448] for a linear goal functional is also optimal for quadratic goal functionals, i.e., GOAFEM leads to linear convergence with optimal convergence rates.

MSC 2010: 65N30; 65N50; 65Y20; 41A25

1 Introduction

Let Ω R d , d 2 , be a bounded Lipschitz domain. For given f L 2 ( Ω ) and f [ L 2 ( Ω ) ] d , we consider a general linear elliptic partial differential equation (1.1)

- div A u + b u + c u = f + div f in Ω ,
u = 0 on Γ := Ω ,
where A ( x ) R sym d × d is a symmetric matrix, b ( x ) R d and c ( x ) R . As usual, we assume that A , b , c L ( Ω ) , that 𝑨 is uniformly positive definite and that the weak form (see (2.4) below) fits into the setting of the Lax–Milgram lemma.

While standard adaptivity aims to approximate the exact solution u H 0 1 ( Ω ) at optimal rate in the energy norm (see, e.g., [12, 23, 7, 25, 11] for some seminal contributions and [14] for the present model problem), goal-oriented adaptivity aims to approximate, at optimal rate, only the functional value G ( u ) R (also called quantity of interest in the literature). Usually, goal-oriented adaptivity is more important in practice than standard adaptivity and has therefore attracted much interest also in the mathematical literature; see, e.g., [3, 5, 13, 18] for some prominent contributions. Unlike standard adaptivity, there are only few works which aim for a mathematical understanding of optimal rates for goal-oriented adaptivity; see [20, 4, 15, 16]. While the latter works consider only linear goal functionals, the present work aims to address, for the first time, optimal convergence rates for goal-oriented adaptivity with a nonlinear goal functional. More precisely, we assume that our primary interest is not in the unknown solution u H 0 1 ( Ω ) , but only in the functional value

(1.2) G ( u ) := K u , u H - 1 × H 0 1

with a quadratic goal functional stemming from a bounded linear operator K : H 0 1 ( Ω ) H - 1 ( Ω ) . Here, , H - 1 × H 0 1 denotes the duality between H 0 1 ( Ω ) and its dual space H - 1 ( Ω ) := H 0 1 ( Ω ) with respect to the extended L 2 -scalar product. Possible examples for such goal functionals include, e.g., G ( u ) = Ω g ( x ) u ( x ) 2 d x for a given weight g L ( Ω ) or G ( u ) = Ω u ( x ) g ( x ) u ( x ) d x for a given weight g [ L ( Ω ) ] d .

In this work, we formulate a goal-oriented adaptive finite element method (GOAFEM), where the quantity of interest G ( u ) is approximated by G ( u ) for some FEM solution u u such that

(1.3) | G ( u ) - G ( u ) | 0 .

Moreover, if 𝒦 is compact, then the convergence (1.3) holds even with optimal algebraic rates. Compared to the available literature on convergence of goal-adaptive FEM for linear goal functionals [20, 4, 15, 16], we need to linearize the nonlinear goal functional in each step of the adaptive algorithm around the present discrete solution u . Put in explicit terms, the prior works consider dual problems which are independent of ℓ, while in the present work, the dual problems change in each step of the adaptive loop. It is our main contribution that the additional linearization error is thoroughly taken into account for the convergence analysis.

Outline

The paper is organized as follows. Section 2 formulates the finite element discretization together with two goal-oriented adaptive algorithms (Algorithms 1 and 2). Moreover, we state our main results. Proposition 1 and Proposition 4 guarantee convergence of Algorithm 1 and Algorithm 2, respectively. If the operator 𝒦 is even compact, then Theorem 2 yields linear convergence with optimal rates for Algorithm 1, while Theorem 5 yields convergence with almost optimal rates for the simpler Algorithm 2. Section 3 provides some numerical experiments which underline the theoretical predictions. Sections 47 are concerned with the proofs of our main results.

2 Adaptive Algorithm and Main Result

2.1 Variational Formulation

Define the bilinear form

(2.1) a ( u , v ) := Ω A u v d x + Ω b u v d x + Ω c u v d x .

We suppose that a ( , ) fits into the setting of the Lax–Milgram lemma, i.e., a ( , ) is continuous and elliptic on H 0 1 ( Ω ) . While continuity

(2.2) a ( u , v ) C cnt u H 1 ( Ω ) v H 1 ( Ω ) for all u , v H 0 1 ( Ω )

follows from the assumptions made with C cnt = A L ( Ω ) + b L ( Ω ) + c L ( Ω ) , the ellipticity

(2.3) a ( u , u ) C ell u H 1 ( Ω ) 2 for all u H 0 1 ( Ω )

requires additional assumptions on the coefficients, e.g.,

inf x Ω inf y R d \ { 0 } y A ( x ) y | y | 2 > 0    and    b H ( div ; Ω ) with inf x Ω ( 1 2 div b ( x ) + c ( x ) ) 0 .

The weak formulation of (1.1) reads

(2.4) a ( u , v ) = F ( v ) := Ω f v d x - Ω f v d x for all v H 0 1 ( Ω ) .

According to the Lax–Milgram lemma, (2.4) admits a unique solution u H 0 1 ( Ω ) . Given w H 0 1 ( Ω ) , the same argument applies and proves that the (linearized) dual problem

(2.5) a ( v , z [ w ] ) = b ( v , w ) + b ( w , v ) for all v H 0 1 ( Ω )

admits a unique solution z [ w ] H 0 1 ( Ω ) , where we abbreviate the notation by use of b ( v , w ) := K v , w H - 1 × H 0 1 . We note that b ( , ) is, in particular, a continuous bilinear form on H 0 1 ( Ω ) . Throughout, we denote by | | | v | | | 2 := Ω A v v d x the energy norm induced by the principal part of a ( , ) , which is an equivalent norm on H 0 1 ( Ω ) . Finally, we stress that all main results also apply to the case that a ( , ) satisfies only a Gårding inequality (instead of the strong ellipticity (2.3)) as long as the weak formulations (2.4) and (2.5) are well-posed; see Section 2.8 below.

2.2 Finite Element Method

For a conforming triangulation T H of Ω into compact simplices and a polynomial degree p 1 , we consider the conforming finite element space

X H := { v H H 0 1 ( Ω ) : for all T T H , v H | T is a polynomial of degree p } .

We approximate u u H X H and z [ w ] z H [ w ] X H . More precisely, the Lax–Milgram lemma yields the existence and uniqueness of discrete FEM solutions u H , z H [ w ] X H of

(2.6) a ( u H , v H ) = F ( v H ) and a ( v H , z H [ w ] ) = b ( v H , w ) + b ( w , v H ) for all v H X H .

2.3 Linearization of the Goal Functional

To control the goal error | G ( u ) - G ( u H ) | , we employ the dual problem. Note that

b ( u - u H , u - u H ) = b ( u , u ) - b ( u H , u ) - b ( u , u H ) + b ( u H , u H ) = [ G ( u ) - G ( u H ) ] - [ b ( u H , u ) + b ( u , u H ) - 2 b ( u H , u H ) ] = [ G ( u ) - G ( u H ) ] - [ b ( u H , u - u H ) + b ( u - u H , u H ) ] .

With the dual problem and the Galerkin orthogonality, we rewrite the second bracket as

b ( u - u H , u H ) + b ( u H , u - u H ) = (2.5) a ( u - u H , z [ u H ] ) = a ( u - u H , z [ u H ] - z H [ u H ] ) .

With continuity of the bilinear forms a ( , ) and b ( , ) , we thus obtain that

(2.7) | G ( u ) - G ( u H ) | = | a ( u - u H , z [ u H ] - z H [ u H ] ) + b ( u - u H , u - u H ) | | | | u - u H | | | [ | | | z [ u H ] - z H [ u H ] | | | + | | | u - u H | | | ] .

2.4 Mesh Refinement

Let T 0 be a given conforming triangulation of Ω. We suppose that the mesh refinement is a deterministic and fixed strategy, e.g., newest vertex bisection [26]. For each triangulation T H and marked elements M H T H , let T h := refine ( T H , M H ) be the coarsest triangulation, where all T M H have been refined, i.e., M H T H \ T h . We write T h T ( T H ) if T h results from T H by finitely many steps of refinement. To abbreviate notation, let T := T ( T 0 ) .

We further suppose that each refined element has at least two sons, i.e.,

(2.8) # ( T H \ T h ) + # T H # T h for all T H T and all T h T ( T H ) ,

and that the refinement rule satisfies the mesh-closure estimate

(2.9) # T - # T 0 C mesh j = 0 - 1 # M j for all N ,

where C mesh > 0 depends only on T 0 . This has first been proved for 2D newest vertex bisection in [7] and has later been generalized to arbitrary dimension d 2 in [26]. While both works require an additional admissibility assumption on T 0 , this has been proved unnecessary at least for 2D in [19]. Finally, it has been proved in [11, 25] that newest vertex bisection ensures the overlay estimate, i.e., for all triangulations T H , T h T , there exists a common refinement T H T h T ( T H ) T ( T h ) which satisfies that

(2.10) # ( T H T h ) # T H + # T h - # T 0 .

For meshes with first-order hanging nodes, (2.8)–(2.10) are analyzed in [8], while T-splines and hierarchical splines for isogeometric analysis are considered in [22, 21] and [9, 17], respectively.

2.5 Error Estimators

For T H T and v H X H , let η H ( T , v H ) 0 and ζ H ( T , v H ) 0 for all T T H be given refinement indicators. For U H T H , let

η H ( U H , v H ) := ( T U H η H ( T , v H ) 2 ) 1 / 2 and ζ H ( U H , v H ) := ( T U H ζ H ( T , v H ) 2 ) 1 / 2 .

To abbreviate notation, let η H ( v H ) := η H ( T H , v H ) and ζ H ( v H ) := ζ H ( T H , v H ) .

We suppose that the estimators η H and ζ H satisfy the following axioms of adaptivity from [10]. There exist constants C stab , C rel , C drel > 0 and 0 < q red < 1 such that, for all T H T and all T h T ( T H ) , the following assumptions are satisfied.

  1. Stability: for all v h X h , v H X H and U H T h T H , it holds that

    | η h ( U H , v h ) - η H ( U H , v H ) | + | ζ h ( U H , v h ) - ζ H ( U H , v H ) | C stab | | | v h - v H | | | .

  2. Reduction: for all v H X H , it holds that

    η h ( T h \ T H , v H ) q red η H ( T H \ T h , v H ) , ζ h ( T h \ T H , v H ) q red ζ H ( T H \ T h , v H ) .

  3. Reliability: for all w H 0 1 ( Ω ) , the Galerkin solutions u H , z H [ w ] X H to (2.6) satisfy that

    | | | u - u H | | | C rel η H ( u H ) , | | | z [ w ] - z H [ w ] | | | C rel ζ H ( z H [ w ] ) .

  4. Discrete reliability: for all w H 0 1 ( Ω ) , the Galerkin solutions u H , z H [ w ] X H and u h , z h [ w ] X h to (2.6) satisfy that

    | | | u h - u H | | | C drel η H ( T H \ T h , u H ) , | | | z h [ w ] - z H [ w ] | | | C drel ζ H ( T H \ T h , z H [ w ] ) .

We note that axioms (A1)(A4) are satisfied for, e.g., standard residual error estimators. Given w H 0 1 ( Ω ) , the mapping v b ( v , w ) + b ( w , v ) is linear and continuous by assumption. Hence, the Riesz theorem from functional analysis guarantees the existence (and uniqueness) of g [ w ] H 0 1 ( Ω ) such that

(2.11) ( g [ w ] , v ) H 1 := Ω g [ w ] v d x + Ω g [ w ] v d x = b ( v , w ) + b ( w , v ) for all v H 0 1 ( Ω ) .

With g [ w ] = - g [ w ] , we thus get that

(2.12) b ( v , w ) + b ( w , v ) = Ω g [ w ] v d x - Ω g [ w ] v d x for all v H 0 1 ( Ω ) ,

i.e., the right-hand sides of the primal problem (2.4) and the (linearized) dual problem (2.5) take the same form. With this[1], the residual error estimators read for v H X H as

η H ( T , v H ) 2 := h T 2 - div ( A v H + f ) + b v H + c v H - f L 2 ( T ) 2 + h T ( A v H + f ) n L 2 ( T Ω ) 2 ,
ζ H ( T , v H ) 2 := h T 2 - div ( A v H + g [ u H ] ) - b v H + ( c - div b ) v H - g [ u H ] L 2 ( T ) 2 + h T ( A v H + g [ u H ] ) n L 2 ( T Ω ) 2 ,
where denotes the jump across edges and 𝒏 is the outwards-facing unit normal vector. We stress that our experiments below directly provide g [ w ] L 2 ( Ω ) and g [ w ] [ L 2 ( Ω ) ] d satisfying representation (2.12) so that there is, in fact, no need to solve (2.11).

2.6 Adaptive Algorithm

We consider the following adaptive algorithm, which adapts the marking strategy proposed in [16].

Algorithm 1

Input: adaptivity parameters 0 < θ 1 and C mark 1 , initial mesh T 0 .

Loop: for all = 0 , 1 , 2 , , perform the following steps.

  1. Compute the discrete solutions u , z [ u ] X to (2.6).

  2. Compute the refinement indicators η ( T , u ) and ζ ( T , z [ u ] ) for all T T .

  3. Determine sets M ¯ u , M ¯ u z T of up to the multiplicative constant C mark minimal cardinality such that (2.13)

    (2.13a) θ η ( u ) 2 η ( M ¯ u , u ) 2 ,
    (2.13b) θ [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] [ η ( M ¯ u z , u ) 2 + ζ ( M ¯ u z , z [ u ] ) 2 ] .

  4. Let M u M ¯ u and M u z M ¯ u z with # M u = # M u z = min { # M ¯ u , # M ¯ u z } .

  5. Define M := M u M u z , and generate T + 1 := refine ( T , M ) .

Output: sequence of triangulations T with corresponding discrete solutions u and z [ u ] as well as error estimators η ( u ) and ζ ( z [ u ] ) .

With Algorithm 2 below, we give and examine an alternative adaptive algorithm that is seemingly cheaper in computational costs.

Our first result states that Algorithm 1 indeed leads to convergence.

Proposition 1

For any bounded linear operator K : H 0 1 ( Ω ) H - 1 ( Ω ) , the following statements hold.

  1. There exists a constant C rel > 0 such that

    (2.14) | G ( u ) - G ( u H ) | C rel η H ( u H ) [ η H ( u H ) 2 + ζ H ( z H [ u H ] ) 2 ] 1 / 2 for all T H T .

  2. For all 0 < θ 1 and 1 < C mark , Algorithm 1 leads to convergence

    (2.15) | G ( u ) - G ( u ) | C rel η ( u ) [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] 1 / 2 0 as .

    The constant C rel depends only on the constants from (A1)(A3), the bilinear form a ( , ) and the boundedness of 𝒦.

To formulate our main result on optimal convergence rates, we need some additional notation. For N N 0 , let T N := { T T : # T - # T 0 N } denote the (finite) set of all refinements of T 0 , which have at most 𝑁 elements more than T 0 . For s , t > 0 , we define

u A s := sup N N 0 ( ( N + 1 ) s min T H T N η H ( u H ) ) R 0 { } ,
z [ u ] A t := sup N N 0 ( ( N + 1 ) t min T H T N ζ H ( z H [ u ] ) ) R 0 { } .
In explicit terms, e.g., u A s < means that an algebraic convergence rate O ( N - s ) for the error estimator η is possible if the optimal triangulations are chosen.

The following theorem concludes the main results of the present work.

Theorem 2

For any compact operator K : H 0 1 ( Ω ) H - 1 ( Ω ) , even the following statements hold, which improve Proposition 1 (ii).

  1. For all 0 < θ 1 and C mark 1 , there exists 0 N 0 , C lin > 0 and 0 < q lin < 1 such that Algorithm 1 guarantees that, for all , n N 0 with n 0 ,

    (2.16) η n ( u n ) [ η n ( u n ) 2 + ζ n ( z n [ u n ] ) 2 ] 1 / 2 C lin q lin n - η ( u ) [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] 1 / 2 .

  2. There exist C opt > 0 and 0 N 0 such that Algorithm 1 guarantees that, for all

    0 < θ < θ opt := ( 1 + C stab 2 C drel 2 ) - 1 ,

    for all s , t > 0 with u A s + z [ u ] A t < , and all N 0 with 0 , it holds that

    (2.17) η ( u ) [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] 1 / 2 C opt u A s ( u A s + z A t ) ( # T - # T 0 ) - α ,

    where α := min { 2 s , s + t } .

The constants C lin , q lin and 0 depend only on 𝜃, q red , C stab , C rel , the bilinear form a ( , ) and the compact operator 𝒦. The constant C opt depends only on 𝜃, C mesh , C mark , C lin , q lin , 0 and (A1)(A4).

Remark 3

(i) We note that, according to the considered dual problem (2.5), the goal functional (1.2) is linearized around u in each step of the adaptive algorithm. Hence, we must enforce that the linearization error satisfies that | | | z [ u ] - z [ u ] | | | 0 as . This is guaranteed by Proposition 1 (ii) and Theorem 2 (i) since both factors of the product involve the primal error estimator η ( u ) .

(ii) For a linear goal functional and hence z [ u ] = z [ u ] , the work [16] considers plain ζ 2 (instead of η 2 + ζ 2 ) for the Dörfler marking (2.13b) and then proves a convergence behavior

| G ( u ) - G ( u ) | η ζ = O ( ( # T ) - α )

for the estimator product, where α = s + t with s > 0 being the optimal rate for the primal problem and t > 0 being the optimal rate for the dual problem. Instead, Algorithm 1 will only lead to O ( ( # T ) - α ) , where α = min { 2 s , s + t } ; see Theorem 2 (ii).

(iii) The marking strategy proposed in [4], where Dörfler marking is carried out for the weighted estimator

(2.18) ρ H ( T , u H , z H [ u H ] ) 2 := η H ( T , u H ) 2 ζ H ( z H [ u H ] ) 2 + η H ( u H ) 2 ζ H ( T , z H [ u H ] ) 2 ,

might be unable to ensure convergence of the linearization error | | | z [ u ] - z [ u ] | | | since, in every step, Dörfler marking is implied for either η H ( u H ) or ζ H ( z H [ u H ] ) ; cf. [16]. If one instead considers

(2.19) ϱ H ( T , u H , z H [ u H ] ) 2 := η H ( T , u H ) 2 [ η H ( u H ) 2 + ζ H ( z H [ u H ] ) 2 ] + η H ( u H ) 2 [ η H ( T , u H ) 2 + ζ H ( T , z H [ u H ] ) 2 ] ,

the present results and the analysis in [16] make it clear that this strategy implies convergence with rate min { 2 s , s + t } . Details are omitted.

2.7 Alternative Adaptive Algorithm

From the upper bound (2.14) in Proposition 1 (i), we can further estimate the goal error by

| G ( u ) - G ( u H ) | C rel [ η H ( u H ) 2 + ζ H ( z H [ u H ] ) 2 ] for all T H T .

This suggests the following algorithm, which marks elements solely based on the combined estimator.

Algorithm 2

Input: adaptivity parameters 0 < θ 1 and C mark 1 , initial mesh T 0 .

Loop: for all = 0 , 1 , 2 , , perform the following steps.

  1. Compute the discrete solutions u , z [ u ] X to (2.6).

  2. Compute the refinement indicators η ( T , u ) and ζ ( T , z [ u ] ) for all T T .

  3. Determine a set M T of up to the multiplicative constant C mark minimal cardinality such that

    (2.20) θ [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] [ η ( M , u ) 2 + ζ ( M , z [ u ] ) 2 ] .

  4. Generate T + 1 := refine ( T , M ) .

Output: sequence of triangulations T with corresponding discrete solutions u and z [ u ] as well as error estimators η ( u ) and ζ ( z [ u ] ) .

First, we note that Algorithm 2 also leads to convergence.

Proposition 4

For any bounded linear operator 𝒦, the following statements hold.

  1. There exists a constant C rel > 0 such that

    | G ( u ) - G ( u H ) | C rel [ η H ( u H ) 2 + ζ H ( z H [ u H ] ) 2 ] for all T H T .

  2. For all 0 < θ 1 and 1 < C mark , Algorithm 2 leads to convergence

    (2.21) | G ( u ) - G ( u ) | C rel [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] 0 as .

The constant C rel depends only on the constants from (A1)(A3), the bilinear form a ( , ) and the boundedness of 𝒦.

The following theorem proves linear convergence of Algorithm 2 with almost optimal convergence rate, where we note that β α for the rates in (2.17) and (2.22). By abuse of notation, we use the same constants as in Theorem 2.

Theorem 5

For any compact operator 𝒦, even the following statements hold, which improve Proposition 4 (ii).

  1. For all 0 < θ 1 and C mark 1 , there exists 0 N 0 , C lin > 0 and 0 < q lin < 1 such that Algorithm 2 guarantees that, for all , n N 0 with n 0 ,

    [ η n ( u n ) 2 + ζ n ( z n [ u n ] ) 2 ] C lin q lin n - [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] .

  2. There exist C opt > 0 and 0 N 0 such that Algorithm 1 guarantees that, for all

    0 < θ < θ opt := ( 1 + C stab 2 C drel 2 ) - 1 ,

    for all s , t > 0 with u A s + z [ u ] A t < , and all N 0 with 0 , it holds that

    (2.22) [ η ( u ) 2 + ζ ( z [ u ] ) 2 ] C opt ( u A s 2 + z A t 2 ) ( # T - # T 0 ) - β ,

    where β := min { 2 s , 2 t } .

The constants C lin , q lin and 0 depend only on 𝜃, q red , C stab , C rel , the bilinear form a ( , ) and the compact operator 𝒦. The constant C opt depends only on 𝜃, C mesh , C mark and (A1)(A4).

Note that Algorithm 2 has slightly lower computational costs than Algorithm 1, but achieves only a lower rate in general. However, if there holds s t , both algorithms achieve rate 2 s .

2.8 Extension of Analysis to Compactly Perturbed Elliptic Problems

For the ease of presentation, we have restricted ourselves to the case that the bilinear form a ( , ) from (2.1) is continuous (2.2) and elliptic (2.3). Actually, it suffices to assume that a ( , ) is continuous and that the energy norm | | | | | | induced by the principal part is an equivalent norm on H 0 1 ( Ω ) , e.g., by assuming that A L ( Ω ) is uniformly positive definite. Then a ( , ) is elliptic up to some compact perturbation (and hence satisfies a Gårding inequality). A prominent example for this problem class is the Helmholtz problem.

We have to assume that the primal formulation (2.4) is well-posed, i.e., for all w H 0 1 ( Ω ) , it holds that

[ a ( w , v ) = 0 for all v H 0 1 ( Ω ) ] w = 0 .

Then the Fredholm alternative and standard functional analysis imply that the primal formulation (2.4) as well as the dual formulation (2.5) admit unique solutions. Moreover, as soon as T H is sufficiently fine, also the FEM problems (2.6) admit unique solutions and, more importantly, the discrete inf-sup constants are uniformly bounded from below; see, e.g., [6, Section 2].

As noted in [6], such an analytical setting requires only two minor modifications of adaptive algorithms.

  1. Step (i) in Algorithm 1 or Algorithm 2: If the discrete solutions u and z [ u ] exist (and are hence also unique), then we proceed as before. If either u or z [ u ] does not exist, then the mesh T + 1 is obtained by uniform refinement of T , i.e., M := T .

  2. Step (iv) of Algorithm 1 or step (iii) of Algorithm 2: Having determined a set of marked elements M T , we select a superset M # M with # M # 2 # M as well as M # { T T : | T | | T | for all T T } and define the refined mesh T + 1 := refine ( T , M # ) via the extended set of marked elements.

It is observed in [6] that uniform refinement caused by modification (a) can only occur finitely many times. Moreover, modification (b) ensures that H 0 1 ( Ω ) = = 0 X ¯ so that the adaptive algorithm converges, indeed, to the right limit. For standard adaptive FEM, it is shown in [6] that this procedure still leads to optimal convergence rates. We note that the arguments from [6] obviously extend to the present goal-oriented adaptive FEM.

3 Numerical Experiments

In this section, we underline our theoretical findings by some numerical examples. As starting point of all examples, we use equation (1.1) with A = I , b = 0 and c = 0 on the unit square Ω = ( 0 , 1 ) 2 . The initial mesh T 0 on Ω is obtained from certain uniform refinements from the mesh shown in Figure 1. All examples are computed with conforming finite elements of order p = 1 and p = 2 , as outlined in Section 2.2.

In the following, we consider the marking strategies of Algorithm 1 and Algorithm 2 (denoted by A and B, respectively), as well as the marking strategies outlined in Remark 3 (iii), i.e., Dörfler marking for (2.18) and (2.19), which will be denoted by BET1 and BET2, respectively. If not stated otherwise, the marking parameter is θ = 0.5 for all experiments.

Figure 1 
               Initial mesh (left) and sets 
                     
                        
                           
                              
                                 U
                                 1
                              
                              ,
                              
                                 U
                                 2
                              
                              ,
                              
                                 U
                                 3
                              
                           
                        
                        
                        U_{1},U_{2},U_{3}
                     
                   (right) on the unit square 
                     
                        
                           
                              Ω
                              =
                              
                                 
                                    (
                                    0
                                    ,
                                    1
                                    )
                                 
                                 2
                              
                           
                        
                        
                        \Omega=(0,1)^{2}
                     
                  .
Figure 1

Initial mesh (left) and sets U 1 , U 2 , U 3 (right) on the unit square Ω = ( 0 , 1 ) 2 .

3.1 Weighted L 2 -Norm

Suppose some weight function λ L ( Ω ) with λ 0 a.e., whose regions of discontinuity are resolved by the initial mesh T 0 (i.e., 𝑔 is continuous in the interior of every element of T 0 ). Then we consider the weighted L 2 -norm

(3.1) G ( u ) = Ω λ ( x ) u ( x ) 2 d x = λ u , u H - 1 × H 0 1 = λ 1 / 2 u L 2 ( Ω ) 2

as goal functional. We note that b ( v , w ) = λ v , w H - 1 × H 0 1 , and hence (2.12) holds with g [ w ] = 2 λ w and g [ w ] = 0 . Moreover, we observe that K u = λ u L 2 ( Ω ) H - 1 ( Ω ) , where the embedding is compact, so that the goal functional from (3.1) fits in the setting of Theorem 2 and Theorem 5. We choose

λ ( x ) = { 1 , x U 1 , 0 , x U 1 ,

with U 1 = ( 0.25 , 0.75 ) 2 . This functional is evaluated at the solution of (1.1) with f = 2 x ( x - 1 ) + 2 y ( y - 1 ) and f = 0 . The solution of this equation, as well as the value of the goal functional, can be computed analytically to be u = x y ( 1 - x ) ( 1 - y ) and G ( u ) = U 1 u 2 d x = 41209 58982400 , respectively. The numerical results are visualized in Figure 2.

Figure 2 
                  Convergence rates of estimator product and goal error for the problem setting in Section 3.1 with 
                        
                           
                              
                                 p
                                 =
                                 1
                              
                           
                           
                           p=1
                        
                      (left) and 
                        
                           
                              
                                 p
                                 =
                                 2
                              
                           
                           
                           p=2
                        
                      (right).
Note that the lines marked with A, B and BET2 are almost identical.
Figure 2

Convergence rates of estimator product and goal error for the problem setting in Section 3.1 with p = 1 (left) and p = 2 (right). Note that the lines marked with A, B and BET2 are almost identical.

3.2 Nonlinear Convection

Suppose that λ [ L ( Ω ) ] 2 is some vector field, whose regions of discontinuity are resolved by the initial mesh T 0 . As goal functional, we consider the nonlinear convection term

(3.2) G ( u ) = Ω u ( x ) λ ( x ) u ( x ) d x = λ u , u H - 1 × H 0 1 .

We note that b ( v , w ) = λ v , w H - 1 × H 0 1 , and hence (2.12) holds with g [ w ] = λ w and g [ w ] = - w λ . Moreover, we observe that K u = λ u L 2 ( Ω ) H - 1 ( Ω ) , where the embedding is compact, so that the goal functional from (3.2) fits in the setting of Theorem 2 and Theorem 5.

We compute the solutions to the primal and the dual problem for f = 0 ,

f ( x ) = { 1 2 ( - 1 , 1 ) if x U 3 , 0 else ,    and    λ = σ 2 ( - 1 1 ) with σ = { 1 if x U 2 , - 1 else .

The sets U 3 := { x Ω : x 1 - x 2 0.25 } and U 2 := ( 0.5 , 1 ) × ( 0 , 0.5 ) are shown in Figure 1. The numerical results are visualized in Figure 3. Note that the primal problem in this case exhibits a singularity which is not induced by the geometry and thus is not present in the dual problem.

Figure 3 
                  Convergence rates of estimator product for the problem setting in Section 3.2 (left) and Section 3.3 (right).
Figure 3 
                  Convergence rates of estimator product for the problem setting in Section 3.2 (left) and Section 3.3 (right).
Figure 3

Convergence rates of estimator product for the problem setting in Section 3.2 (left) and Section 3.3 (right).

3.3 Force Evaluation

Let ε > 0 , and let 𝜓 be a cut-off function that satisfies

ψ ( x ) = { 1 if x U 1 , 0 if dist ( x , U 1 ) > ε .

For a given direction χ R 2 , consider a goal functional of the form

(3.3) G ( u ) := Ω ψ ( u u - 1 2 | u | 2 I ) χ d x .

This approximates the electrostatic force which is exerted by an electric potential 𝑢 on a charged body occupying the domain U 1 in direction χ (the part of the integrand in brackets is the so-called Maxwell stress tensor). We note that

b ( v , w ) = Ω ψ ( v w - 1 2 ( v w ) I ) χ d x ,

and hence (2.12) holds with g [ w ] = 0 and g [ w ] = ( ψ χ ) w - ( ψ w ) χ - ( χ w ) ψ . We stress that the goal functional from (3.3) does not fit in the setting of Theorem 2 and Theorem 5 since the corresponding operator 𝒦 is not compact. Hence, we cannot guarantee optimal rates for Algorithms 1 and 2. However, Proposition 1 and Proposition 4 still guarantee convergence of our algorithms.

For our experiments, we choose χ = 1 2 ( 1 , 1 ) , f = 1 and f = 0 . Furthermore, we choose 𝜓 to be in X 0 for p = 1 , i.e., 𝜓 is piecewise linear, and 𝜀 is chosen such that 𝜓 falls off to 0 exactly within one layer of elements around U 1 in T 0 . The results can be seen in Figure 3.

3.4 Discussion of Numerical Experiments

We clearly see from Figures 2 and 3 that Algorithm 1 and BET2 outperform BET1 and sometimes even Algorithm 2. From Figure 4, where we plot estimator product (and, if available, goal error) for different parameters θ = 0.1 , 0.2 , , 1.0 , we see that this behavior does not depend on the marking parameter 𝜃, generally speaking. It is striking that the strategy BET1 with θ < 1 fails to drive down the estimator product at the same speed as uniform refinement. This is likely due to the fact that the linearization error | | | z [ u ] - z [ u ] | | | is disregarded; see Remark 3.

In Figure 5, we plot the cumulative costs

(3.4) S [ τ ] # T , with S [ τ ] := { N : error τ } ,

where error is either the estimator product η ( u ) ζ ( z [ u ] ) or the goal error | G ( u ) - G ( u ) | in the ℓ-th step of the adaptive algorithm. We see that, for the setting from Section 3.1, where no singularity occurs, optimal costs are achieved by uniform refinement, as is expected. For the goal error, which is not known in general, the strategy BET1 performs better than Algorithms 1 and 2. However, for the estimator product, which is the relevant quantity in most applications (since the error is unknown), it is inferior. In the other settings, where there is a singularity, Algorithms 1 and 2 achieve their minimal cost around the value 0.7 for the marking parameter 𝜃.

Figure 4 
                  Variation of 𝜃 from 0.1 (light) to 1.0 (dark) in steps of 0.1.
Left: setting from Section 3.1 with 
                        
                           
                              
                                 p
                                 =
                                 2
                              
                           
                           
                           p=2
                        
                     , where the upper lines represent the estimator product and the lower ones the goal error.
Middle: estimator product for the setting from Section 3.2 with 
                        
                           
                              
                                 p
                                 =
                                 2
                              
                           
                           
                           p=2
                        
                     .
Right: estimator product for the setting from Section 3.3 with 
                        
                           
                              
                                 p
                                 =
                                 1
                              
                           
                           
                           p=1
                        
                     .
Figure 4

Variation of 𝜃 from 0.1 (light) to 1.0 (dark) in steps of 0.1. Left: setting from Section 3.1 with p = 2 , where the upper lines represent the estimator product and the lower ones the goal error. Middle: estimator product for the setting from Section 3.2 with p = 2 . Right: estimator product for the setting from Section 3.3 with p = 1 .

Figure 5 
                  Cumulative costs (3.4) for estimator product and goal error for the setting of Section 3.1 (top), for the estimator product for the setting of Section 3.2 (bottom left) and for the estimator product for the setting of Section 3.3 (bottom right).
The parameters 𝜃 are chosen uniformly in 
                        
                           
                              
                                 [
                                 0.1
                                 ,
                                 1
                                 ]
                              
                           
                           
                           [0.1,1]
                        
                      with step size 0.1.
Figure 5

Cumulative costs (3.4) for estimator product and goal error for the setting of Section 3.1 (top), for the estimator product for the setting of Section 3.2 (bottom left) and for the estimator product for the setting of Section 3.3 (bottom right). The parameters 𝜃 are chosen uniformly in [ 0.1 , 1 ] with step size 0.1.

4 Auxiliary Results

4.1 Axioms of Adaptivity

Clearly, z [ w ] and z H [ w ] depend linearly on 𝑤 (since 𝒦 is linear and hence b ( , ) is bilinear). Moreover, we have the following stability estimates.

Lemma 6

For all w H 0 1 ( Ω ) and all T H T ( T 0 ) , it holds that

(4.1) C 1 - 1 | | | z H [ w ] | | | | | | z [ w ] | | | C 2 | | | w | | | ,

where C 1 > 0 depends only on a ( , ) , while C 2 > 0 depends additionally on the boundedness of 𝒦.

Proof

The definition of the dual problem shows that

| | | z [ w ] | | | 2 a ( z [ w ] , z [ w ] ) = (2.5) b ( z [ w ] , w ) + b ( w , z [ w ] ) | | | w | | | | | | z [ w ] | | | ,

and hence | | | z [ w ] | | | | | | w | | | . Moreover, the stability of the Galerkin method yields that | | | z H [ w ] | | | | | | z [ w ] | | | . This concludes the proof. ∎

Next, we show that the combined estimator for the primal and dual problem satisfies assumptions (A1)(A4), where particular emphasis is put on (A3) and (A4). For the ease of presentation (and by abuse of notation), we use the same constants as for the original properties (A1)(A4), even though they now depend additionally on the bilinear form a ( , ) and the boundedness of 𝒦.

Proposition 7

Suppose (A1)(A4) for η H and ζ H . Let T H T and T h T ( T H ) . Then (A1)(A4) hold also for the combined estimator [ η H ( ) 2 + ζ H ( ) 2 ] 1 / 2 .

  1. For all v h , w h X h , v H , w H X H and U H T h T H , it holds that

    | [ η h ( U H , v h ) 2 + ζ h ( U H , w h ) 2 ] 1 / 2 - [ η H ( U H , v H ) 2 + ζ H ( U H , w H ) 2 ] 1 / 2 | C stab [ | | | v h - v H | | | + | | | w h - w H | | | ] .

  2. For all v H , w H X H , it holds that

    [ η h ( T h \ T H , v H ) 2 + ζ h ( T h \ T H , w H ) 2 ] 1 / 2 q red [ η H ( T H \ T h , v H ) 2 + ζ H ( T H \ T h , w H ) 2 ] 1 / 2 .

  3. The Galerkin solutions u H , z H [ u H ] X H to (2.6) satisfy that

    | | | u - u H | | | + | | | z [ u ] - z H [ u H ] | | | + | | | z [ u H ] - z H [ u H ] | | | C rel [ η H ( u H ) 2 + ζ H ( z H [ u H ] ) 2 ] 1 / 2 .

  4. The Galerkin solutions u H , z H [ u H ] X H and u h , z h [ u h ] X h to (2.6) satisfy that

    | | | u h - u H | | | + | | | z h [ u h ] - z H [ u H ] | | | C drel [ η H ( T H \ T h , u H ) 2 + ζ H ( T H \ T h , z H [ u H ] ) 2 ] 1 / 2 .

Proof

By the triangle inequality and [ a 2 + b 2 ] 1 / 2 a + b , (A1) follows from stability of η H and ζ H . Reduction (A2) follows directly from the corresponding properties of η H and ζ H . For (A3), we see with Lemma 6 that

| | | z [ u ] - z H [ u H ] | | | | | | z [ u ] - z [ u H ] | | | + | | | z [ u H ] - z H [ u H ] | | | (4.1) | | | u - u H | | | + | | | z [ u H ] - z H [ u H ] | | | .

Hence, (A3) follows from reliability of η H and ζ H . Discrete reliability (A4) follows from the same arguments. ∎

In the following, we recall some basic results of [10].

Lemma 8

Lemma 8 (Quasi-Monotonicity of Estimators [10, Lemma 3.6])

Let w H 0 1 ( Ω ) . Let T H T and T h T ( T H ) . Properties (A1)(A3) together with the Céa lemma guarantee that

(4.2) η h ( u h ) 2 C mon η H ( u H ) 2 as well as ζ h ( z h [ w ]