Optimal complexity of goal-oriented adaptive FEM for nonsymmetric linear elliptic PDEs

Philipp Bringmann; Maximilian Brunner; Dirk Praetorius; Julian Streitberger

doi:10.1515/jnma-2023-0150

Enjoy 40% off

academic books on De Gruyter Brill *

Article Open Access

Optimal complexity of goal-oriented adaptive FEM for nonsymmetric linear elliptic PDEs

Philipp Bringmann , Maximilian Brunner , Dirk Praetorius and Julian Streitberger

Published/Copyright: November 4, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Journal of Numerical Mathematics Volume 33 Issue 2

Abstract

We analyze a goal-oriented adaptive algorithm that aims to efficiently compute the quantity of interest G(u ^⋆) with a linear goal functional G and the solution u ^⋆ to a general second-order nonsymmetric linear elliptic partial differential equation. The current state of the analysis of iterative algebraic solvers for nonsymmetric systems lacks the contraction property in the norms that are prescribed by the functional analytic setting. This seemingly prevents their application in the optimality analysis of goal-oriented adaptivity. As a remedy, this paper proposes a goal-oriented adaptive iteratively symmetrized finite element method (GOAISFEM). It employs a nested loop with a contractive symmetrization procedure, e.g., the Zarantonello iteration, and a contractive algebraic solver, e.g., an optimal multigrid solver. The various iterative procedures require well-designed stopping criteria such that the adaptive algorithm can effectively steer the local mesh refinement and the computation of the inexact discrete approximations. The main results consist of full linear convergence of the proposed adaptive algorithm and the proof of optimal convergence rates with respect to both degrees of freedom and total computational cost (i.e., optimal complexity). Numerical experiments confirm the theoretical results and investigate the selection of the parameters.

Keywords: goal-oriented adaptive finite element method; linear quantity of interest; iterative solver; nonsymmetric partial differential equations; optimal convergence rates; optimal complexity

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

1 Introduction

Adaptive finite element methods (AFEMs) are a cornerstone in the numerical solution of partial differential equations (PDEs). The abundant literature emphasizes significant progress and manifests a matured understanding of the topic; see, e.g. [1]–[9], for linear elliptic PDEs.

The variational formulation of a nonsymmetric second-order linear elliptic PDE with bilinear form b(⋅, ⋅) and right-hand side functional F on the Sobolev space X : = H 0 1 ( Ω ) seeks a weak solution u ^⋆ to

(1) b ( u ⋆ , v ) = F ( v ) ∀ v ∈ X .

While standard AFEM aims at an efficient approximation of the solution u ⋆ ∈ X , goal-oriented AFEM (GOAFEM) strives only to approximate a quantity of interest G(u ^⋆); see [10], [11], [12], [13] for early prominent contributions. However, to accurately approximate G(u ^⋆) for a continuous linear goal functional G : X → R , following the generic approach G(u _H) ≈ G(u ^⋆) leads to convergence rates determined by the error of the approximation u _H ≈ u ^⋆ to the primal problem (1). Instead, GOAFEM adopts a duality technique by additionally approximating z H ≈ z ⋆ ∈ X solving the dual problem

(2) b ( v , z ⋆ ) = G ( v ) ∀ v ∈ X .

Following [13], a discrete approximation G _H(u _H, z _H) ≈ G(u ^⋆) enables the control of the error for any u H , z H ∈ X by

(3) | G ( u ⋆ ) − G H ( u H , z H ) | ⩽ | b u ⋆ − u H , z ⋆ − z H | ⩽ L | | | u ⋆ − u H | | | | | | z ⋆ − z H | | | ,

where L > 0 is the continuity constant of b(⋅, ⋅) with respect to the energy norm  ||| ⋅ |||; see Section 2 for details. As seen in (3), this approach allows to add the convergence rates of the primal and dual problem. Moreover, it is not necessary – and may even lead to unnecessary computational expense – to compute approximations u _H ≈ u ^⋆ and z _H ≈ z ^⋆ across the entire domain with the same accuracy. Instead, a careful marking of elements for refinement enables a considerable reduction of the computational costs and makes GOAFEM highly relevant in both practical applications and mathematical research.

First rigorous convergence results of GOAFEM are found in [14]–[18], recent contributions in this context include [19], [20] and for a dual weighted-residual approach see, e.g., [21], [22], [23]. The works [14], [16], [17], [19], [20] focus on optimal convergence rates with respect to the degrees of freedom. However, the cumulative nature of adaptivity calls for optimal convergence rates with respect to the total computational effort, i.e., the overall computational time. Coined as optimal complexity initially for wavelet-based discretizations [24], [25], this notion was later adopted for AFEM with contributions including, e.g., [4], [26], [27], [28]. In the setting of GOAFEM, optimal complexity was established first in [14] for the Poisson problem and sufficiently small adaptivity parameters, and extended to a general second-order symmetric linear elliptic PDE with uniformly contractive algebraic solver in [29]. Since uniform contraction with respect to the PDE-related energy norm for nonsymmetric algebraic solvers such as GMRES is still open, as a remedy, the proof of the Lax–Milgram lemma motivates the application of an iterative symmetrization [28]. This results in a sequence of symmetric algebraic systems that allow the application of optimal algebraic solvers, e.g., [30], [31], [32]. Figure 1 illustrates the nested structure of the resulting goal-oriented adaptive iteratively symmetrized finite element method (GOAISFEM). The detailed Algorithm 1 is presented in Section 3 below. Table 1 displays the notation of the associated indices and quasi-error quantities, which are equivalent to the total error.

Figure 1:

Schematic overview of the GOAISFEM algorithm with nested symmetrization and inexact solver.

Table 1:

Iteration counters and quasi-errors for the GOAISFEM algorithm. We note that for the combination of the index sets, the quasi-errors are extended to the full index set by the last available quasi-error. We refer to Section 3 for details on the iteration counters and index sets and to the beginning of Section 5 for a detailed description of the quasi-errors and their extension to the full index set Q .

Iteration	Mesh refinement		Symmetrization		Algebraic solver		Index set	Quasi-error
	Running	Final	Running	Final	Running	Final
Primal	ℓ	ℓ ̲	m	m ̲	n	n ̲	Q u in (24a)	H ℓ m , n in (44a)
Dual	ℓ	ℓ ̲	μ	μ ̲	ν	ν ̲	Q z in (24b)	Z ℓ μ , ν in (44b)
Combined	ℓ	ℓ ̲	k	k ̲ = m a x { m ̲ , μ ̲ }	j	j ̲ = m a x { n ̲ , ν ̲ }	Q = Q u ∪ Q z	H ℓ k , j Z ℓ k , j in (45)

The first challenge in the analysis of the GOAISFEM algorithm consists of the nonlinear product structure attained by the combined quasi-error product as displayed in Table 1. The resulting nonlinear remainder term significantly complicates the proof compared to treating only the primal problem as in [28] and requires the application of a novel proof strategy from [33] that only utilizes summability of the remainder, denoted as tail-summability throughout. The second challenge arises from the combination of the primal and dual marking leading to a merged marked set. Thereby, either only the primal or only the dual estimator is guaranteed to satisfy the estimator reduction property. Since the estimator belongs to the quasi-error, this also leads to a failure of contraction for one of the two involved quasi-errors. While [29] solves this issue in the symmetric case, the additional symmetrization loop results in a more involved situation at hand. Adapting the novel approach of the tail-summability criterion from [33] enables the proof of full linear convergence and optimal complexity for the nonlinear quasi-error product in this paper. The analysis employs the generalized quasi-orthogonality from [34] to remedy the lack of a Pythagorean identity for nonsymmetric problems.

Our main result asserts full linear convergence of the quasi-error product H ℓ k , j Z ℓ k , j with respect to the total step counter |⋅, ⋅, ⋅| (measuring the total solver steps in the index set). Therein, we allow for an arbitrary symmetrization stopping parameter λ _sym and only require a small algebraic solver parameter λ _alg such that the product λ _sym λ _alg is sufficiently small. More precisely, Theorem 1 states that there exist constants C _lin > 0 and 0 < q _lin < 1 such that, for all ( ℓ , k , j ) , ( ℓ ′ , k ′ , j ′ ) ∈ Q with |ℓ′, k′, j′ |⩽| ℓ, k, j|,

H ℓ k , j Z ℓ k , j ⩽ C lin q lin | ℓ , k , j | − | ℓ ′ , k ′ , j ′ | H ℓ ′ k ′ , j ′ Z ℓ ′ k ′ , j ′ .

Note that, unlike [28], where full linear convergence is guaranteed only for sufficiently large ℓ ⩾ ℓ ₀, the current result is stronger in the sense that the result holds for ℓ ₀ = 0 owing to a generalized quasi-orthogonality from [34]. An immediate consequence of full linear convergence and the geometric series in Corollary 1 states that the rates with respect to the degrees of freedom coincide with the rates with respect to the cumulative computational work (i.e., computational time), i.e., for all r > 0, there holds

sup ( ℓ , k , j ) ∈ Q # T ℓ r H ℓ k , j Z ℓ k , j ⩽ sup ( ℓ , k , j ) ∈ Q ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k , j | # T ℓ ′ r H ℓ k , j Z ℓ k , j ⩽ C cost sup ( ℓ , k , j ) ∈ Q # T ℓ r H ℓ k , j Z ℓ k , j

along the sequence of meshes T ℓ generated by the GOAISFEM algorithm. The second main result of Theorem 2 proves that, for sufficiently small adaptivity parameters and any achievable rates s, t > 0 of the primal resp. dual problem (stated in terms of nonlinear approximation classes), the algorithm guarantees optimal complexity, i.e.,

sup ( ℓ , k , j ) ∈ Q ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k , j | # T ℓ ′ s + t H ℓ k , j Z ℓ k , j ⩽ C opt max ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t , H 0 0,0 Z 0 0,0 .

This means the convergence of the algorithm attains the optimal rate s + t with respect to the overall computational work, where ‖ u ⋆ ‖ A s < ∞ means that u ^⋆ can be approximated at rate s (along a sequence of unavailable optimal meshes) and likewise for z ^⋆.

The remaining parts of the paper are organized as follows. The preliminary Section 2 introduces the model problem, the assumptions on the solvers, and the axioms of adaptivity from [9], including the general quasi-orthogonality from [34]. Following the algorithm in Section 3 and its contraction properties in Sections 4 and 5 presents full linear convergence as the first main result of this paper. This allows to prove optimal complexity in Section 6 as the second main result, which is underlined by the numerical experiments in Section 7 including a thorough investigation of the adaptivity parameters. The paper concludes with a summary in Section 8.

2 Setting

In this section, we introduce the problem and explain the key components needed to design the adaptive algorithm in Section 3.

2.1 Continuous model problem

Let Ω ⊂ R d with d ⩾ 1 be a polygonal Lipschitz domain. Given right-hand sides f ∈ L ²(Ω) and f ∈ [ L 2 ( Ω ) ] d , we consider a general second-order linear elliptic PDE

(4) − div ( A ∇ u ⋆ ) + b ⋅ ∇ u ⋆ + c u ⋆ = f − div ( f ) in Ω subject to u ⋆ = 0 on ∂ Ω ,

with a pointwise symmetric and positive definite diffusion matrix A ∈ L ∞ ( Ω ) sym d × d , a convection coefficient b ∈ L ∞ ( Ω ) d , and a reaction coefficient c ∈ L ^∞(Ω). For well-definedness of the a posteriori error estimator in Section 2.6 below, we additionally require that A | T ∈ W 1 , ∞ ( T ) sym d × d and f | T ∈ H 1 ( T ) d for all T ∈ T 0 , where T 0 is an initial triangulation that subdivides Ω into compact simplices. Let ⟨ ⋅, ⋅ ⟩ denote the L ²(Ω)-scalar product. With the principal part a(u, v): = ⟨ A ∇u, ∇v⟩, the variational formulation of (4) seeks a solution u ⋆ ∈ X : = H 0 1 ( Ω ) to the so-called primal problem

(5) b ( u ⋆ , v ) : = a ( u ⋆ , v ) + ⟨ b ⋅ ∇ u ⋆ + c u ⋆ , v ⟩ = ⟨ f , v ⟩ + ⟨ f , ∇ v ⟩ = : F ( v ) ∀ v ∈ X .

We suppose that the bilinear form b(⋅, ⋅) from (5) is continuous and elliptic with respect to the norm  ‖ ⋅ ‖ X on X , i.e., there exist constants L′, α′ > 0 such that

(6) b ( u , v ) ⩽ L ′ ‖ u ‖ X ‖ v ‖ X , b ( v , v ) ⩾ α ′ ‖ v ‖ X 2 ∀ u , v ∈ X .

Then, the Lax–Milgram lemma proves existence and uniqueness of the solution u ^⋆ to (5). An elementary compactness argument shows that (6) implies ellipticity of the principal part a(⋅, ⋅) and thus a(⋅, ⋅) is a scalar product on X with induced energy norm  a ( ⋅ , ⋅ ) 1 / 2 = : | | | ⋅ | | | ≃ ‖ ⋅ ‖ X , cf. [35], Remark 3]. Therefore, b(⋅, ⋅) is also continuous and elliptic with respect to ||| ⋅ |||, i.e., there exist constants L, α > 0 such that

(7) b ( u , v ) ⩽ L | | | u | | | | | | v | | | , b ( v , v ) ⩾ α | | | v | | | 2 ∀ u , v ∈ X .

In the present paper, we suppose that the quantity of interest G is linear and reads for given data g ∈ L ²(Ω) and g ∈ L 2 ( Ω ) d ,

G ( v ) : = ∫ Ω g v + g ⋅ ∇ v d x .

In order to guarantee well-definedness of the error estimator in Section 2.6 below, we suppose g | T ∈ H 1 ( T ) d for all initial simplices T ∈ T 0 . In view of the continuity and coercivity of b(⋅, ⋅), the Lax–Milgram lemma yields existence and uniqueness of the solution z ⋆ ∈ X of the so-called dual problem: Find z ⋆ ∈ X such that

(8) b ( v , z ⋆ ) = G ( v ) ∀ v ∈ X .

2.2 Finite element discretization and discrete goal

For a polynomial degree p ∈ N and a conforming simplicial triangulation T H of Ω, the discrete ansatz space reads

(9) X H : = { v H ∈ X : ∀ T ∈ T H , v H | T is a polynomial of total degree ⩽ p } .

Since X H ⊂ X is conforming, the Lax–Milgram lemma ensures the existence and uniqueness of primal and dual discrete solutions u H ⋆ , z H ⋆ ∈ X H satisfying

(10) b u H ⋆ , v H = F ( v H ) , b v H , z H ⋆ = G ( v H ) ∀ v H ∈ X H .

It is well-known that conforming FEMs are quasi-optimal, i.e., there hold Céa-type estimates with constant C _Céa = L/α

(11) | | | u ⋆ − u H ⋆ | | | ⩽ C Céa min v H ∈ X H | | | u ⋆ − v H | | | , | | | z ⋆ − z H ⋆ | | | ⩽ C Céa min v H ∈ X H | | | z ⋆ − v H | | | .

For arbitrary approximations u H , z H , ∈ X H the linearity of the quantity of interest G as well as the primal and the dual problem (1) and (2) show that

G ( u ⋆ ) − G ( u H ) = G u ⋆ − u H = ( 2 ) b u ⋆ − u H , z ⋆ = ( 1 ) b u ⋆ − u H , z ⋆ − z H + F ( z H ) − b ( u H , z H ) .

The definition of the discrete goal quantity by G H ( u H , z H ) : = G ( u H ) + F ( z H ) − b ( u H , z H ) allows to control the goal error by continuity of b(⋅, ⋅)

(12) | G ( u ⋆ ) − G H ( u H , z H ) | ⩽ | b u ⋆ − u H , z ⋆ − z H | ⩽ L | | | u ⋆ − u H | | | | | | z ⋆ − z H | | | .

We emphasize that (12) holds for any u _H, z _H and, in particular, for those stemming from an iterative solution step. Moreover, if u H = u H ⋆ , then G ( u H , z H ) = G u H ⋆ as expected.

2.3 Zarantonello iteration

The discrete formulations (10) lead to positive definite, but nonsymmetric linear systems of equations. To reduce the formulation to symmetric and positive definite (SPD) problems, we follow previous own work [28] for the primal problem and employ the Zarantonello iteration [36]. Typically, the latter is used in the up-to-date proof of the Lax–Milgram lemma and also defines a linearization scheme for the treatment of a certain class of nonlinear elliptic PDEs (see, e.g., [33], [37], [38], [39]). In its core, it is a fixed-point method, thus also applicable in the nonsymmetric setting at hand. For a damping parameter δ > 0 and given u H , z H ∈ X H , the Zarantonello iterations Φ H u , Φ H z : ( 0 , ∞ ) × X H → X H compute the unique solutions Φ H u ( δ ; u H ) , Φ H z ( δ ; z H ) ∈ X H to the symmetric variational formulations

(13a) a ( Φ H u ( δ ; u H ) , v H ) = a ( u H , v H ) + δ F ( v H ) − b ( u H , v H ) ∀ v H ∈ X H ,

(13b) a ( v H , Φ H z ( δ ; z H ) ) = a ( v H , z H ) + δ G ( v H ) − b ( v H , z H ) ∀ v H ∈ X H .

The Riesz–Fischer theorem (and also the Lax–Milgram lemma) guarantees existence and uniqueness of Φ H u ( δ ; u H ) , Φ H z ( δ ; z H ) ∈ X H , i.e., the Zarantonello operators Φ H u ( δ ; ⋅ ) and Φ H z ( δ ; ⋅ ) are well-defined. In particular, the exact discrete solutions u H ⋆ = Φ H u δ ; u H ⋆ and z H ⋆ = Φ H z δ ; z H ⋆ are the unique fixed points for all δ > 0. Moreover, for a sufficiently small damping parameter δ, i.e., 0 < δ < δ ^⋆: = 2α/L ², the Banach fixed-point theorem [40], Section 25.4] guarantees that Φ H u ( δ ; ⋅ ) and Φ H z ( δ ; ⋅ ) are contractive with constant 0 < q sym ⋆ : = 1 − δ ( 2 α − δ L 2 ) 1 / 2 < 1 , i.e., for all functions v H , w H ∈ X H , it holds that

(14) max | | | Φ H u ( δ ; v H ) − Φ H u ( δ ; w H ) | | | , | | | Φ H z ( δ ; v H ) − Φ H z ( δ ; w H ) | | | ⩽ q sym ⋆ | | | v H − w H | | | .

The optimal value δ _opt = α/L ² yields the minimal contraction value q sym ⋆ = 1 − α 2 / L 2 .

2.4 Algebraic solver

A canonical candidate for solving (10) directly is a generalized minimal residual method [41], [42] with optimal preconditioner for the symmetric part. While this guarantees uniform contraction of the algebraic residuals in a discrete vector norm, the link between the algebraic residuals and the functional setting is still open [28]. Instead, after a symmetrization with the Zarantonello iteration, it remains to solve the SPD systems (13). Since large SPD problems are still computationally expensive and the exact solution cannot be computed in linear computational complexity, we employ an iterative algebraic solver whose iteration is expressed by the operator Ψ H : X ′ × X H → X H . More precisely, given a bounded linear functional ψ ∈ X ′ and an approximation w H ∈ X H of the exact solution w H ⋆ ∈ X H to a w H ⋆ , v H = ψ ( v H ) for all v H ∈ X H , the algebraic solver returns an improved approximation Ψ H ( ψ ; w H ) ∈ X H in the sense that there exists 0 < q _alg < 1 independent of ψ and X H such that

To simplify notation, we shall identify ψ with its Riesz representative w H ⋆ ∈ X H and write Ψ H w H ⋆ ; ⋅ instead of Ψ_H(ψ; ⋅), even though w H ⋆ is unknown in practice and will only be approximated by an optimal algebraic solver, e.g., [30], [31], [32]. In the following, we use the hp-robust multigrid method from [32] with localized lowest-order smoothing on intermediate levels and patchwise higher-order smoothing on the finest mesh as an innermost algebraic solver loop.

2.5 Mesh refinement

The mesh refinement employs newest-vertex bisection (NVB). We refer to [43] for NVB with admissible initial triangulation T 0 and d ⩾ 2, to [44], [45] for NVB with general T 0 for d ∈ {1, 2}, and to the recent work [46] for NVB with general T 0 in any dimension d ⩾ 2. For each triangulation T H and marked elements M H ⊆ T H , let T h : = r e f i n e ( T H , M H ) be the coarsest conforming refinement of T H such that at least 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 can be obtained from T H by finitely many steps of NVB, and T h ∈ T N ( T H ) if T h ∈ T ( T H ) with # T h − # T H ⩽ N with the number of additional elements N ∈ N 0 . To simplify notation, we write T : = T ( T 0 ) and T N : = T N ( T 0 ) . We note that the nestedness of meshes T h ∈ T ( T H ) implies nestedness of the corresponding finite element spaces X H ⊆ X h ⊂ X from (9).

2.6 A posteriori error estimation

For a triangle T ∈ T H ∈ T and v H ∈ X H , let n denote the outer unit normal vector and [ [ ⋅ ] ] the jump along inner edges of T H . We define the refinement indicators η _H(T; v _H) ⩾ 0 and ζ _H(T; v _H) ⩾ 0 for the primal and dual problem from (10), respectively, by

(16a) η H ( T ; v H ) 2 : = | T | 2 / d ‖ − div ( A ∇ v H − f ) + b ⋅ ∇ v H + c v H − f ‖ L 2 ( T ) 2 + | T | 1 / d ‖ [ [ A ∇ v H − f ⋅ n ] ] ‖ L 2 ( ∂ T ∩ Ω ) 2 , ζ H ( T ; v H ) 2 : = | T | 2 / d ‖ − div ( A ∇ v H − g ) − b ⋅ ∇ v H + c − div ( b ) v H − g ‖ L 2 ( T ) 2 + | T | 1 / d ‖ [ [ A ∇ v H − g ⋅ n ] ] ‖ L 2 ( ∂ T ∩ Ω ) 2 .

For any subset U H ⊆ T H , we abbreviate

(16b) η H ( U H ; v H ) 2 : = ∑ T ∈ U H η H ( T ; v H ) 2 , ζ H ( U H ; v H ) 2 : = ∑ T ∈ U H ζ H ( T ; v H ) 2

as well as η H ( v H ) : = η H ( T H ; v H ) and ζ H ( v H ) : = ζ H ( T H ; v H ) for all v H ∈ X H .

For details on residual-based error estimators, we refer to [47], [48]. Throughout the paper, the index of the estimators refer to the underlying mesh, e.g., η _h and ζ _h on the refinement T h ∈ T ( T H ) or η _ℓ and ζ _ℓ on a sequence of meshes T ℓ with ℓ ∈ N 0 . It is well-known that η _H, ζ _H satisfy the following axioms of adaptivity.

Lemma 1

(see [9], Section 6.1]). The error estimators η _H, ζ _H from (16) satisfy the following properties with constants C _stab, C _rel, C _drel, C _mon > 0 and 0 < q _red < 1 for any triangulation T H ∈ T and any conforming refinement T h ∈ T ( T H ) with the corresponding Galerkin solutions u H ⋆ , z H ⋆ ∈ X H , u h ⋆ , z h ⋆ ∈ X h to (10), any subset U H ⊆ T H ∩ T h , and arbitrary v H ∈ X H , v h ∈ X h :

(A1) stability: | η h ( U H ; v h ) − η H ( U H ; v H ) | + | ζ h ( U H ; v h ) − ζ H ( U H ; v H ) | ⩽ C s t a b | | | v h − v H | | | ;
(A2) reduction: η h ( T h \ T H ; v H ) ⩽ q r e d η H ( T H \ T h ; v H ) and ζ h ( T h \ T H ; v H ) ⩽ q r e d ζ H ( T H \ T h ; v H ) ;
(A3) reliability: | | | u ⋆ − u H ⋆ | | | ⩽ C r e l η H u H ⋆ and | | | z ⋆ − z H ⋆ | | | ⩽ C r e l ζ H z H ⋆ ;
(A3⁺) discrete reliability: | | | u h ⋆ − u H ⋆ | | | ⩽ C d r e l η H T H \ T h , u H ⋆ and | | | z h ⋆ − z H ⋆ | | | ⩽ C d r e l ζ H T H \ T h , z H ⋆ ;
(QM) quasi-monotonicity: η h u h ⋆ ⩽ C m o n η H u H ⋆ and ζ h z h ⋆ ⩽ C m o n ζ H z H ⋆ .

The constant C _rel depends only on the uniform γ-shape regularity of all T H ∈ T and on the space dimension d, while C _stab and C _drel additionally depend on the polynomial degree p. For NVB, reduction (A2) holds with q _red: = 2^−1/(2d). Moreover, the constant in quasi-monotonicity (QM) satisfies C _mon ⩽ min{1 + C _stab(1 + C _Céa)C _rel, 1 + C _stab C _drel}.

Reliability (A3) and stability (A1) verify

| | | u ⋆ − u H | | | ⩽ max { C rel , 1 + C stab C rel } η H ( u H ) + | | | u H ⋆ − u H | | | , | | | z ⋆ − z H | | | ⩽ max { C rel , 1 + C stab C rel } ζ H ( z H ) + | | | z H ⋆ − z H | | | .

In combination with the estimate (12), we finally conclude for C goal : = L max { C rel , 1 + C stab C rel } 2 the reliable goal-error estimate

which provides the core estimate of the proposed adaptive algorithm in Section 3 below.

The ellipticity of b(⋅, ⋅) from (7) ensures inf-sup stability of the elliptic problem at hand. Recall from [34] that inf-sup stability implies the generalized quasi-orthogonality, which will be an important tool in the subsequent analysis.

Proposition 1

(validity of quasi-orthogonality [34], Eq. (8)]). For any sequence X ℓ ⊆ X ℓ + 1 ⊂ X of nested discrete subspaces with ℓ ⩾ 0, there holds

(A4) quasi-orthogonality: There exist constants C _orth > 0 and 0 < δ < 1 such that the corresponding Galerkin solutions u ℓ ⋆ , z ℓ ⋆ ∈ X ℓ to (10) satisfy, for all ℓ , M ∈ N 0 ,

The constants C _orth and δ depend only on the dimension d, the elliptic bilinear form b(⋅, ⋅), and the chosen norm ||| ⋅ |||, but are independent of the spaces X ℓ .

3 Adaptive algorithm

In this section, we introduce our goal-oriented adaptive iteratively symmetrized algorithm. It utilizes specific stopping indices denoted by an underline, e.g., ℓ ̲ , m ̲ [ ℓ ] , n ̲ [ ℓ , k ] ∈ N 0 . For an overview, see Table 1 above. However, we may omit the dependence whenever it is apparent from the context, such as in the abbreviation n ̲ : = n ̲ [ ℓ , m ] for u ℓ m , n ̲ .

Algorithm 1

(GOAISFEM).

Input: Initial mesh T 0 , polynomial degree p ∈ N , marking parameters 0 < θ ⩽ 1, C _mark ⩾ 1, solver parameters λ _sym > 0, λ _alg > 0, Zarantonello damping parameter δ > 0, and initial guesses u 0 0,0 = u 0 0 , n ̲ , z 0 0,0 = z 0 0 , ν ̲ ∈ X 0 .

Adaptive loop: For all ℓ = 0, 1, 2, …, repeat the following steps (I)–(IV):

SOLVE & ESTIMATE (PRIMAL). For all m = 1, 2, 3, …, repeat (a)–(c):
1. Set u ℓ m , 0 : = u ℓ m − 1 , n ̲ and define for theoretical reasons u ℓ m , ⋆ ≔ Φ ℓ u δ ; u ℓ m − 1 , n ̲ .
2. For all n = 1, 2, 3, …, repeat the following steps (i)–(ii):
  1. Compute u ℓ m , n ≔ Ψ ℓ u ℓ m , ⋆ ; u ℓ m , n − 1 and corresponding refinement indicators η ℓ T ; u ℓ m , n for all T ∈ T ℓ .
  2. Terminate n-loop and define n ̲ [ ℓ , m ] ≔ n if
    (19) | | | u ℓ m , n − u ℓ m , n − 1 | | | ⩽ λ alg λ sym η ℓ u ℓ m , n + | | | u ℓ m , n − u ℓ m , 0 | | | .
3. Terminate m-loop and define m ̲ [ ℓ ] ≔ m if
  (20) | | | u ℓ m , n ̲ − u ℓ m , 0 | | | ⩽ λ sym η ℓ u ℓ m , n ̲ .
SOLVE & ESTIMATE (DUAL). For all μ = 1, 2, 3, …, repeat (a)–(c):
1. Set z ℓ μ , 0 ≔ z ℓ μ − 1 , ν ̲ and define for theoretical reasons z ℓ μ , ⋆ ≔ Φ ℓ z δ ; z ℓ μ − 1 , ν ̲ .
2. For all ν = 1, 2, 3, …, repeat the following steps (i)–(ii):
  1. Compute z ℓ μ , ν ≔ Ψ ℓ z ℓ μ , ⋆ ; z ℓ μ , ν − 1 and corresponding refinement indicators ζ ℓ T ; z ℓ μ , ν for all T ∈ T ℓ .
  2. Terminate ν-loop and define ν ̲ [ ℓ , μ ] ≔ ν if
    (21) | | | z ℓ μ , ν − z ℓ μ , ν − 1 | | | ⩽ λ alg λ sym ζ ℓ z ℓ μ , ν + | | | z ℓ μ , ν − z ℓ μ , 0 | | | .
3. Terminate μ-loop and define μ ̲ [ ℓ ] ≔ μ if
  (22) | | | z ℓ μ , ν ̲ − z ℓ μ , 0 | | | ⩽ λ sym ζ ℓ z ℓ μ , ν ̲ .
MARK. Determine sets
M ̄ ℓ u ∈ M ℓ u θ , u ℓ m ̲ , n ̲ : = U ℓ ⊆ T ℓ : θ η ℓ u ℓ m ̲ , n ̲ 2 ⩽ η ℓ U ℓ , u ℓ m ̲ , n ̲ 2 , M ̄ ℓ z ∈ M ℓ z θ , z ℓ μ ̲ , ν ̲ : = U ℓ ⊆ T ℓ : θ ζ ℓ z ℓ μ ̲ , ν ̲ 2 ⩽ ζ ℓ U ℓ , z ℓ μ ̲ , ν ̲ 2
satisfying the following Dörfler criterion [1] with quasi-minimal cardinality
(23) # M ̄ ℓ u ⩽ C mark min U ℓ ⋆ ∈ M ℓ u θ , u ℓ m ̲ , n ̲ # U ℓ ⋆ , # M ̄ ℓ z ⩽ C mark min U ℓ ⋆ ∈ M ℓ z θ , z ℓ μ ̲ , ν ̲ # U ℓ ⋆ .

As in [17], define the set of marked elements M ℓ ≔ M ℓ u ∪ M ℓ z , where M ℓ u ⊆ M ̄ ℓ u and M ℓ z ⊆ M ̄ ℓ z satisfy # M ℓ u = # M ℓ z = min # M ̄ ℓ u , # M ̄ ℓ z .
REFINE. Generate the new mesh T ℓ + 1 ≔ r e f i n e ( M ℓ , T ℓ ) by NVB and define u ℓ + 1 0,0 ≔ u ℓ + 1 0 , n ̲ ≔ u ℓ + 1 0 , ⋆ ≔ u ℓ m ̲ , n ̲ and z ℓ + 1 0,0 : = z ℓ + 1 0 , ν ̲ : = z ℓ + 1 0 , ⋆ : = z ℓ μ ̲ , ν ̲ (nested iteration).

Output: Sequences of successively refined triangulations T ℓ , successive discrete approximations u ℓ m , n , z ℓ μ , ν , and corresponding error estimators η ℓ u ℓ m , n , ζ z ℓ μ , ν .

Remark 1.

(i) Although the primal loop (I) and dual loop (II) in Algorithm 1 are displayed sequentially, they are independent of each other. Therefore, a practical implementation will realize these iterations simultaneously since the system matrix is the same (thanks to the symmetrization step).

(ii) In order to investigate the asymptotic behavior, it is reasonable to analyze Algorithm 1 in the present formulation with infinitely many steps. We note that a practical implementation will terminate with ℓ ̲ ≔ ℓ provided that the estimator product is smaller than a user-specified tolerance.

For the analysis of Algorithm 1, we define the index set Q ≔ Q u ∪ Q z with

(24a) Q u ≔ ( ℓ , m , n ) ∈ N 0 3 : u ℓ m , n is used in Algorithm 1 ,

(24b) Q z ≔ ( ℓ , μ , ν ) ∈ N 0 3 : z ℓ μ , ν is used in Algorithm 1 .

Furthermore, we require the following final indices and notice that these are consistent with those defined in Algorithm 1:

(25a) ℓ ̲ ≔ sup ℓ ∈ N 0 : ( ℓ , 0,0 ) ∈ Q u or ( ℓ , 0,0 ) ∈ Q z ∈ N 0 ∪ { ∞ } ,

(25b) m ̲ [ ℓ ] ≔ sup { m ∈ N : ( ℓ , m , 0 ) ∈ Q u } , μ ̲ [ ℓ ] ≔ sup { μ ∈ N : ( ℓ , μ , 0 ) ∈ Q z } ,

(25c) n ̲ [ ℓ , m ] ≔ sup { n ∈ N : ( ℓ , m , n ) ∈ Q u } , ν ̲ [ ℓ , μ ] ≔ sup { ν ∈ N : ( ℓ , μ , ν ) ∈ Q z } .

In addition, we set k ̲ [ ℓ ] ≔ max { m ̲ [ ℓ ] , μ ̲ [ ℓ ] } as well as j ̲ [ ℓ , k ] ≔ max { n ̲ [ ℓ , k ] , ν ̲ [ ℓ , k ] } .

Finally, we introduce the total step counter |⋅, ⋅, ⋅| defined for all ( ℓ , k , j ) ∈ Q by

| ℓ , k , j | = ∑ ℓ ′ = 0 ℓ − 1 ∑ k ′ = 0 k ̲ [ ℓ ′ ] ∑ j ′ = 0 j ̲ [ ℓ ′ , k ′ ] 1 + ∑ k ′ = 0 k − 1 ∑ j ′ = 0 j ̲ [ ℓ , k ′ ] 1 + ∑ j ′ = 0 j − 1 1 .

This definition indeed provides a lexicographic ordering on Q , if the solver steps Algorithm 1(I) for u ℓ m , n and Algorithm 1(II) for z ℓ μ , ν are done in parallel. We note that one solver step of an optimal geometric multigrid method on graded meshes can be performed in O ( # T ℓ ) operations; see, e.g., [30], [32]. For given u ℓ m , n , z ℓ μ , ν ∈ X ℓ , the simultaneous computation of the refinement indicators η ℓ T , u ℓ m , n and ζ ℓ T , z ℓ μ , ν requires O ( # T ℓ ) operations, hence the steps Algorithm 1(I)–(II) require O ( # T ℓ ) operations as well. Furthermore, Dörfler marking can be performed in O ( # T ℓ ) operations; see, e.g., [4], [49]. Therefore, the total work to compute u ℓ m , n and z ℓ μ , ν is (up to a constant) given by

(26) c o s t ( ℓ , k , j ) ≔ ∑ ( ℓ ′ , m ′ , n ′ ) ∈ Q u | ℓ ′ , m ′ , n ′ | ⩽ | ℓ , k , j | # T ℓ ′ + ∑ ( ℓ ′ , μ ′ , ν ′ ) ∈ Q z | ℓ ′ , μ ′ , ν ′ | ⩽ | ℓ , k , j | # T ℓ ′ ≃ ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k , j | # T ℓ ′ .

Since # Q = ∞ , we have either ℓ ̲ = ∞ , or k ̲ [ ℓ ̲ ] = ∞ , or j ̲ [ ℓ ̲ , k ̲ ] = ∞ . A further observation about Algorithm 1 is that the nested algebraic solver loop within the Zarantonello loop is guaranteed to terminate, and the latter case j ̲ [ ℓ ̲ , k ̲ ] = ∞ is therefore excluded.

Lemma 2

(finite termination of algebraic solver [28], Lemma 3.2]). Independently of the algorithmic parameters δ, θ, λ _sym, and λ _alg, the innermost n- and ν-loops of Algorithm 1 always terminate. In particular, j ̲ [ ℓ , k ] < ∞ for all ( ℓ , k , 0 ) ∈ Q .

4 A posteriori error analysis

Algorithm 1 does not provide the exact algebraic solutions u ℓ m , ⋆ and z ℓ μ , ⋆ to (13) but instead uses an inexact algebraic solver. However, the following result from [28] applies to the primal and the dual problem alike and shows that these inexact Zarantonello iterations remain contractions except for the final iterate on each mesh (see also [50] for an extended version).

Lemma 3

(contraction of inexact Zarantonello iteration [28], Lemma 5.1]). Choose any damping parameter 0 < δ < δ ^⋆ = 2α/L ² to ensure the contraction (14) of the Zarantonello iteration and

(27) 0 < λ alg ⋆ < ( 1 − q sym ⋆ ) ( 1 − q alg ) 4 q alg such that 0 < q sym ≔ q sym ⋆ + 2 q alg 1 − q alg λ alg ⋆ 1 − 2 q alg 1 − q alg λ alg ⋆ < 1 .

Then, for arbitrary λ _sym > 0 and any 0 < λ alg ⩽ λ alg ⋆ , we have for all ( ℓ , m , n ̲ ) ∈ Q u with 1 ⩽ m < m ̲ [ ℓ ] and all ( ℓ , μ , ν ̲ ) ∈ Q z with 1 ⩽ μ < μ ̲ [ ℓ ] that

Moreover, for m = m ̲ [ ℓ ] resp. μ = μ ̲ [ ℓ ] , it holds that

The subsequent lemma gathers a posteriori error estimates following directly from the corresponding contraction of the symmetrization, algebraic solver, and the inexact Zarantonello iteration. Further details of the elementary proof are omitted.

Lemma 4

(stability and a posteriori error control). For all ( ℓ , m , 0 ) ∈ Q u with m ⩾ 1, contraction (14) shows

Analogously, for all ( ℓ , m , n ) ∈ Q u with n ⩾ 1, the contraction (15) ensures

For all ( ℓ , m , n ̲ ) ∈ Q u with 1 ⩽ m < m ̲ [ ℓ ] , the contraction (28) leads to

The analogous estimates are also valid for the dual variable.

Finally, the following lemma shows that in the case of finitely many mesh-refinement steps, the Zarantonello iteration does not terminate and one of the two exact continuous solutions is already the discrete solution to (10).

Lemma 5

(case of finite mesh-refinement steps). Suppose that the inexact Zarantonello iteration satisfies contraction (28) and that η and ζ satisfy (A1)–(A3). If ℓ ̲ < ∞ , then k ̲ [ ℓ ̲ ] = ∞ and η ℓ ̲ u ℓ ̲ ⋆ = 0 (so that u ⋆ = u ℓ ̲ ⋆ ) or ζ ℓ ̲ z ℓ ̲ ⋆ = 0 (so that z ⋆ = z ℓ ̲ ⋆ ).

Proof.

By Lemma 2, we have j ̲ [ ℓ , k ] < ∞ . If ℓ ̲ < ∞ , then k ̲ [ ℓ ̲ ] = ∞ and, hence,

If (33) holds, then the inexact Zarantonello iterates u ℓ ̲ m , n ̲ are convergent with limit u ℓ ̲ ⋆ and we obtain by stability (A1) that

This proves that η ℓ ̲ u ℓ ̲ ⋆ = 0 , and we infer from reliability (A3) that u ℓ ̲ ⋆ = u ⋆ . The same arguments apply to z ℓ ̲ ⋆ in the case of (34).□

Due to the contraction of the inexact Zarantonello iteration (28), we have the following a posteriori error estimates for the final iterates.

Lemma 6

(stability of final iterates). Suppose that the inexact Zarantonello iteration satisfies (28). Then, for all ( ℓ + 1 , m ̲ , n ̲ ) ∈ Q u and ( ℓ + 1 , μ ̲ , ν ̲ ) ∈ Q z , there holds

(35) | | | u ℓ + 1 ⋆ − u ℓ + 1 m ̲ − 1 , n ̲ | | | ⩽ | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | , | | | z ℓ + 1 ⋆ − z ℓ + 1 μ ̲ − 1 , ν ̲ | | | ⩽ | | | z ℓ + 1 ⋆ − z ℓ μ ̲ , ν ̲ | | | ,

(36) | | | u ℓ + 1 m ̲ , n ̲ − u ℓ m ̲ , n ̲ | | | ⩽ 4 | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | , | | | z ℓ + 1 μ ̲ , ν ̲ − z ℓ μ ̲ , ν ̲ | | | ⩽ 4 | | | z ℓ + 1 ⋆ − z ℓ μ ̲ , ν ̲ | | | ,

(37) | | | u ℓ m ̲ , n ̲ − u ℓ m ̲ − 1 , n ̲ | | | ⩽ 4 | | | u ℓ ⋆ − u ℓ m ̲ − 1 , n ̲ | | | , | | | z ℓ μ ̲ , ν ̲ − z ℓ μ ̲ − 1 , ν ̲ | | | ⩽ 4 | | | z ℓ ⋆ − z ℓ μ ̲ − 1 , ν ̲ | | | .

Proof.

For ( ℓ + 1 , m ̲ , n ̲ ) ∈ Q u , nested iteration u ℓ + 1 0 , n ̲ = u ℓ m ̲ , n ̲ together with the contraction of the inexact Zarantonello iteration (28) and m ̲ [ ℓ + 1 ] ⩾ 1 prove (35) by

Let ( ℓ , m ̲ , n ̲ ) ∈ Q u . Contraction of the algebraic solver (15), the fact n ̲ [ ℓ , m ̲ ] ⩾ 1 , and nested iteration u ℓ m ̲ , 0 = u ℓ m ̲ − 1 , n ̲ show that

This and with the contraction of the exact Zarantonello iteration (14) result in

Consequently, the combination of (39) and (35) validates (36) via

| | | u ℓ + 1 m ̲ , n ̲ − u ℓ m ̲ , n ̲ | | | ⩽ | | | u ℓ + 1 ⋆ − u ℓ + 1 m ̲ , n ̲ | | | + | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | ⩽ ( 39 ) 3 | | | u ℓ + 1 ⋆ − u ℓ + 1 m ̲ − 1 , n ̲ | | | + | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | ⩽ ( 35 ) 4 | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | .

The estimate (39) also implies (37), because

| | | u ℓ m ̲ , n ̲ − u ℓ m ̲ − 1 , n ̲ | | | ⩽ | | | u ℓ ⋆ − u ℓ m ̲ , n ̲ | | | + | | | u ℓ ⋆ − u ℓ m ̲ − 1 , n ̲ | | | ⩽ ( 39 ) 4 | | | u ℓ ⋆ − u ℓ m ̲ − 1 , n ̲ | | | .

The same arguments prove the estimates for the dual variable and conclude the proof.□

The subsequent lemma states the estimator reduction for only one of the two error estimators. This poses a significant challenge in the proof of full linear convergence due to the required contraction of the nonlinear quasi-error product in Lemma 8 below.

Lemma 7

(estimator reduction and stability). Define the constant 0 < q ( θ ) : = 1 − ( 1 − q red 2 ) θ 1 / 2 < 1 and suppose that the estimators η and ζ satisfy (A1)–(A2). If the primal error estimator satisfies the Dörfler criterion, i.e., M ℓ u = M ̄ ℓ u ⊆ M ℓ in Algorithm 1(III), then

(40) η ℓ + 1 u ℓ + 1 m ̲ , n ̲ ⩽ q ( θ ) η ℓ u ℓ m ̲ , n ̲ + 4 C stab | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | ∀ ( ℓ + 1 , m ̲ , n ̲ ) ∈ Q u , ζ ℓ + 1 z ℓ + 1 μ ̲ , ν ̲ ⩽ ζ ℓ z ℓ μ ̲ , ν ̲ + 4 C stab | | | z ℓ + 1 ⋆ − z ℓ μ ̲ , ν ̲ | | | ∀ ( ℓ + 1 , μ ̲ , ν ̲ ) ∈ Q z .

If the dual error estimator satisfies the Dörfler criterion, i.e., M ℓ z = M ̄ ℓ z ⊆ M ℓ in Algorithm 1(III), then

(41) η ℓ + 1 u ℓ + 1 m ̲ , n ̲ ⩽ η ℓ u ℓ m ̲ , n ̲ + 4 C stab | | | u ℓ + 1 ⋆ − u ℓ m ̲ , n ̲ | | | ∀ ( ℓ + 1 , m ̲ , n ̲ ) ∈ Q u , ζ ℓ + 1 z ℓ + 1 μ ̲ , ν ̲ ⩽ q ( θ ) ζ ℓ z ℓ μ ̲ , ν ̲ + 4 C stab | | | z ℓ + 1 ⋆ − z ℓ μ ̲ , ν ̲ | | | ∀ ( ℓ + 1 , μ ̲ , ν ̲ ) ∈ Q z .

Proof.

For ( ℓ + 1,0,0 ) ∈ Q u , stability (A1) and reduction (A2) yield that

(42) η ℓ + 1 u ℓ m ̲ , n ̲ 2 = η ℓ + 1 T ℓ + 1 ∩ T ℓ ; u ℓ m ̲ , n ̲ 2 + η ℓ + 1 T ℓ + 1 \ T ℓ ; u ℓ m ̲ , n ̲ 2 ⩽ η ℓ T ℓ + 1 ∩ T ℓ ; u ℓ m ̲ , n ̲ 2 + q red 2 η ℓ T ℓ \ T ℓ + 1 ; u ℓ m ̲ , n ̲ 2 = η ℓ u ℓ m ̲ , n ̲ 2 − ( 1 − q red 2 ) η ℓ T ℓ \ T ℓ + 1 ; u ℓ m ̲ , n ̲ 2 .

The Dörfler marking in Algorithm 1(III) for the primal error estimator η and M ℓ ⊆ T ℓ \ T ℓ + 1 prove the contraction in (40)

(43) η ℓ + 1 u ℓ m ̲ , n ̲ 2 ⩽ η ℓ u ℓ m ̲ , n ̲ 2 − ( 1 − q red 2 ) η ℓ M ℓ ; u ℓ m ̲ , n ̲ 2 ⩽ q ( θ ) 2 η ℓ u ℓ m ̲ , n ̲ 2 .

For ( ℓ + 1 , m ̲ , n ̲ ) ∈ Q u , this and (36) lead to

For ( ℓ + 1 , μ ̲ , ν ̲ ) ∈ Q z , we argue analogously to (42) in order to obtain that ζ ℓ + 1 z ℓ μ ̲ , ν ̲ ⩽ ζ ℓ z ℓ μ ̲ , ν ̲ . Together with (36), it follows that

The proof holds verbatim in the case of Dörfler marking for the dual error estimator, albeit with reversed roles. This concludes the proof.□

5 Full linear convergence

This section presents full linear convergence of Algorithm 1 as the first main result of this work. Recall the goal-error estimate from (17) motivating the product structure of the respective primal and dual error components. Thus, we define the quasi-errors

The quasi-errors naturally extend to the full index set ( ℓ , k , j ) ∈ Q by

(45) H ℓ k , j : = H ℓ k , n ̲ if ( ℓ , k , 0 ) ∈ Q u but ( ℓ , k , j ) ∉ Q u , H ℓ m ̲ , n ̲ if ( ℓ , k , 0 ) ∉ Q u , Z ℓ k , j : = Z ℓ k , ν ̲ if ( ℓ , k , 0 ) ∈ Q z but ( ℓ , k , j ) ∉ Q z , Z ℓ μ ̲ , ν ̲ if ( ℓ , k , 0 ) ∉ Q z .

The following theorem asserts full linear convergence of the quasi-error product.

Theorem 1

(full linear convergence). Suppose that the estimators η and ζ satisfy (A1)–(A3) and (QM) and suppose (A4). Recall λ alg ⋆ and q _sym from Lemma 3. With the constant q(θ) from Lemma 7 and q ̄ : = max { q ( θ ) 1 / 2 , ( 1 + q sym ⋆ ) / 2 } < 1 , let

(46) 0 < λ ⋆ : = ( 1 − q alg ) ( q ̄ − q sym ⋆ ) ( 1 − q ̄ ) 10 q alg C stab .

Then, for arbitrary marking parameter 0 < θ ⩽ 1 and any solver parameters λ _sym > 0 and 0 < λ alg ⩽ λ alg ⋆ with λ _sym λ _alg ⩽ λ ^⋆, Algorithm 1 guarantees full linear convergence: There exist constants C _lin ⩾ 1 and 0 < q _lin < 1 such that the quasi-error product satisfies, for all ( ℓ , k , j ) , ( ℓ ′ , k ′ , j ′ ) ∈ Q with |ℓ′, k′, j′ |⩽| ℓ, k, j|

(47) H ℓ k , j Z ℓ k , j ⩽ C lin q lin | ℓ , k , j | − | ℓ ′ , k ′ , j ′ | H ℓ ′ k ′ , j ′ Z ℓ ′ k ′ , j ′ .

The constants C _lin and q _lin depend only on C _stab, C _rel, C _mon, C _orth, C _Céa, θ, q _red, q _sym, q sym ⋆ , q _alg, λ _sym, and λ _alg.

Three lemmas are required to prove Theorem 1. The characterization of R-linear convergence from [33], Lemma 5 and 10] is the primary tool for the proof of Theorem 1; see (70) below. The proof of Theorem 1 departs with the contraction of the quasi-error for the final iterates of the inexact Zarantonello loop up to a remainder on the mesh level ℓ. To this end, we define the simplified weighted quasi-error

where γ > 0 is a free parameter chosen in (51) below. This quasi-error quantity satisfies contraction up to a tail-summable remainder due to estimator reduction (40), (41).

Lemma 8

(contraction in mesh level up to tail-summable remainder). Under the assumptions of Theorem 1, there exists 0 < q < 1 such that the quasi-error product H_ℓ Z_ℓ from (48) satisfies contraction up to a remainder R _ℓ ⩾ 0,

(49) H ℓ + 1 Z ℓ + 1 ⩽ q H ℓ Z ℓ + q R ℓ ∀ ( ℓ + 1 , k ̲ , j ̲ ) ∈ Q .

The remainder R _ℓ satisfies

(50) R ℓ + M ≲ H ℓ Z ℓ and ∑ ℓ ′ = ℓ ℓ + M R ℓ ′ 2 ≲ ( M + 1 ) 1 − δ H ℓ 2 Z ℓ 2 ∀ ℓ , M ∈ N 0 with ℓ + M < ℓ ̲ .

Proof.

The proof consists of four steps.

Step 1 (choice of constants). Recall the constants 0 < q(θ) < 1 from Lemma 7 and λ ^⋆ > 0 and 0 < q ̄ < 1 defined in the statement of Theorem 1 and define the constants

C ( γ , λ ) ≔ 1 + 2 q alg 1 − q alg λ γ > 1 , 0 < q ctr ≔ max q sym ⋆ + 4 C stab C ( γ , λ ) γ , q ( θ ) C ( γ , λ ) .

Elementary calculations show that the choice of

(51) γ : = q ̄ ( q ̄ − q sym ⋆ ) 4 C stab < 1

ensures q sym ⋆ C ( γ , λ ) + 4 C stab γ C ( γ , λ ) 2 < 1 as well as, for all 0 < λ < λ ^⋆,

(52) C ( γ , λ ) = 1 + 2 q alg 1 − q alg λ γ < 1 + 1 − q ̄ q ̄ = 1 q ̄ ⩽ 1 q ( θ ) 1 / 2 .

Consequently, we have q(θ) C(γ,λ)² < 1 and thus 0 < q ctr ′ : = C ( γ , λ ) q ctr < 1 and q _ctr < 1.

Step 2 (contraction of H _ℓ and Z _ℓ ). Abbreviate λ: = λ _alg λ _sym. Recall that marking in Algorithm 1(III) ensures that the estimate (40) or (41) hold. If (40) is satisfied, the quasi-contraction of the inexact Zarantonello iteration (29) for the final iterate, the stability estimate (35), and the estimator reduction (40) lead, for all ( ℓ + 1 , k ̲ , j ̲ ) ∈ Q u , to

The same arguments yield, for all ( ℓ + 1 , μ ̲ , ν ̲ ) ∈ Q z ,

For 0 < q ctr < q ctr ′ = C ( γ , λ ) q ctr < 1 , the product of (53) and (54) reads

If (41) is satisfied, we obtain the same estimate with reversed roles in the derivation.

Step 3 (quasi-monotonicity of H _ℓ and Z _ℓ ). The Céa estimate (11), nestedness of the discrete spaces, reliability (A3), quasi-monotonicity (QM), stability (A1), and the definition (48) prove, for all ℓ ⩽ ℓ ′ ⩽ ℓ ″ ⩽ ℓ ̲ with ( ℓ , m ̲ , n ̲ ) ∈ Q u and ( ℓ , μ ̲ , ν ̲ ) ∈ Q z , that

(56a) | | | u ℓ '' ⋆ − u ℓ ' ⋆ | | | ≲ ( 11 ) | | | u ⋆ − u ℓ ' ⋆ | | | ≲ ( A 3 ) η ℓ ' u ℓ ' ⋆ ≲ ( Q M ) η ℓ u ℓ ⋆ ≲ ( A 1 ) η ℓ u ℓ m ̲ , n ̲ + | | | u ℓ ⋆ − u ℓ m ̲ , n ̲ | | | ≃ ( 48 ) H ℓ ,

(56b) | | | z ℓ ' ' ⋆ − z ℓ ' ⋆ | | | ≲ ( 11 ) | | | z ⋆ − z ℓ ' ⋆ | | | ≲ ( A 3 ) ζ ℓ ' z ℓ ' ⋆ ≲ ( Q M ) ζ ℓ z ℓ ⋆ ≲ ( A 1 ) ζ ℓ z ℓ m ̲ u , ν ̲ + | | | z ℓ ⋆ − z ℓ μ ̲ , ν ̲ | | | ≃ ( 48 ) Z ℓ ,

where the hidden constants depend only on γ ⁻¹, C _Céa, C _stab, C _rel, and C _mon.

Similarly to (53), the inexact Zarantonello contraction (29), stability (A1), and the stability estimate (36) yield for ℓ < ℓ ′ < ℓ ̲ and λ = λ _sym λ _alg,

(57) | | | u ℓ ' ⋆ − u ℓ ' m ̲ , n ̲ | | | ⩽ ( 29 ) q sym ⋆ | | | u ℓ ' ⋆ − u ℓ ' m ̲ − 1 , n ̲ | | | + 2 q alg 1 − q alg λ η ℓ ' u ℓ ' m ̲ , n ̲ ⩽ ( 28 ) q sym ⋆ q sym m ̲ [ ℓ ' ] − 1 | | | u ℓ ' ⋆ − u ℓ ' − 1 m ̲ , n ̲ | | | + 2 q alg 1 − q alg λ η ℓ ' u ℓ ' m ̲ , n ̲ ⩽ ( A 1 ) , ( 42 ) q sym ⋆ | | | u ℓ ' ⋆ − u ℓ ' − 1 m ̲ , n ̲ | | | + 2 q alg 1 − q alg λ η ℓ ' − 1 u ℓ ' − 1 m ̲ , n ̲ + 2 C stab q alg 1 − q alg λ | | | u ℓ ' m ̲ , n ̲ − u ℓ ' − 1 m ̲ , n ̲ | | | ⩽ ( 36 ) q sym ⋆ + 8 C stab q alg 1 − q alg λ | | | u ℓ ' ⋆ − u ℓ ' − 1 m ̲ , n ̲ | | | + 2 q alg 1 − q alg λ η ℓ ' − 1 u ℓ ' − 1 m ̲ , n ̲ ⩽ ( A 1 ) q sym ⋆ + 10 C stab q alg 1 − q alg λ | | | u ℓ ' − 1 ⋆ − u ℓ ' − 1 m ̲ , n ̲ | | | + q sym ⋆ + 8 C stab q alg 1 − q alg λ | | | u ℓ ' ⋆ − u ℓ ' − 1 ⋆ | | | + 2 q alg 1 − q alg λ η ℓ ' − 1 u ℓ ' − 1 ⋆ ⩽ ( 56 a ) q sym ⋆ + 10 C stab q alg 1 − q alg λ | | | u ℓ ' − 1 ⋆ − u ℓ ' − 1 m ̲ , n ̲ | | | + q sym ⋆ + 2 q alg 1 − q alg λ C mon 4 C stab ( 1 + C Céa ) C rel + 1 η ℓ u ℓ ⋆ .

The choice of λ ⩽ λ ^⋆ with λ ^⋆ from (46) ensures

(58) 0 < q : = q sym ⋆ + 10 C stab q alg 1 − q alg λ < 1 .

With C : = q sym ⋆ + 2 q alg 1 − q alg λ C mon 4 C stab ( 1 + C Céa ) C rel + 1 , a successive application of (57) and the geometric series shows

Hence, we have quasi-monotonicity of the quasi-error

The same argument proves

(60b) Z ℓ + M ≲ Z ℓ ∀ ℓ , M ∈ N 0 with ℓ + M < ℓ ̲ .

Step 4 (contraction of H _ℓ Z _ℓ up to tail-summable remainder). Define

The contraction (55) proves the quasi-contraction (49) via

The remainder term R _ℓ can be estimated via (56) and the Young inequality by

(61) R ℓ 2 ≲ ( 56 ) | | | u ℓ + 1 ⋆ − u ℓ ⋆ | | | Z ℓ + | | | z ℓ + 1 ⋆ − z ℓ ⋆ | | | H ℓ 2 ≲ | | | u ℓ + 1 ⋆ − u ℓ ⋆ | | | 2 Z ℓ 2 + | | | z ℓ + 1 ⋆ − z ℓ ⋆ | | | 2 H ℓ 2 .

Thus, the quasi-monotonicity (60) verifies

R ℓ + M ≲ H ℓ + M Z ℓ + M ≲ ( 60 ) H ℓ Z ℓ ∀ ℓ , M ∈ N with ℓ + M < ℓ ̲ .

Quasi-orthogonality (A4), reliability (A3), and the estimates (56) imply, for all ℓ , M ∈ N 0 with ℓ + M < ℓ ̲ ,

(62) ∑ ℓ ′ = ℓ ℓ + M | | | u ℓ ′ + 1 ⋆ − u ℓ ′ ⋆ | | | 2 ≲ ( A 4 ) ( M + 1 ) 1 − δ | | | u ⋆ − u ℓ ⋆ | | | 2 ≲ ( A 3 ) ( M + 1 ) 1 − δ η ℓ u ℓ ⋆ 2 ≲ ( 56 a ) ( M + 1 ) 1 − δ H ℓ 2 , ∑ ℓ ′ = ℓ ℓ + M | | | z ℓ ′ + 1 ⋆ − z ℓ ′ ⋆ | | | 2 ≲ ( A 4 ) ( M + 1 ) 1 − δ | | | z ⋆ − z ℓ ⋆ | | | 2 ≲ ( A 3 ) ( M + 1 ) 1 − δ ζ ℓ z ℓ ⋆ 2 ≲ ( 56 b ) ( M + 1 ) 1 − δ Z ℓ 2 .

Using (61), the quasi-monotonicity (60), and (62), we conclude the proof of (50), for all ℓ , M ∈ N 0 with ℓ + M < ℓ ̲ ,

□

The tail-summability in ℓ provides the basis for the proof of tail-summability on the mesh level ℓ together with the Zarantonello symmetrization index k for the final iterates of the algebraic solver. The main ingredients in the proof of tail-summability in (ℓ, k) are Lemma 8 and the following quasi-contraction in the symmetrization index k.

Lemma 9

(quasi-contraction of inexact Zarantonello symmetrization). There holds

(63) H ℓ k ' , j ̲ Z ℓ k ' , j ̲ ≲ q sym k ' − k H ℓ k , j ̲ Z ℓ k , j ̲ ∀ ( ℓ , k ' , j ̲ ) ∈ Q with 0 ⩽ k ⩽ k ' ⩽ k ̲ [ ℓ ] ,

(64) H ℓ 0 , j ̲ Z ℓ 0 , j ̲ ≲ H ℓ − 1 Z ℓ − 1 ∀ ( ℓ , 0 , 0 ) ∈ Q with ℓ ⩾ 1 .

Proof.

First, we note that the a posteriori error control (31) and the stopping criteria of the algebraic solver (19) and of the symmetrization (20) lead, for ( ℓ , m ̲ , n ̲ ) ∈ Q u , to

Since the two notions of quasi-errors H_ℓ and H ℓ k ̲ , j ̲ only differ by the middle term | | | u ℓ m ̲ , ⋆ − u ℓ m ̲ , n ̲ | | | and the fixed constant factor 0 < γ < 1, this and the analogous estimate for the dual variable show

(65) H ℓ ⩽ H ℓ k ̲ , j ̲ ≲ H ℓ , Z ℓ ⩽ Z ℓ k ̲ , j ̲ ≲ Z ℓ ∀ ( ℓ , k ̲ , j ̲ ) ∈ Q .

For 0 ⩽ k < k ′ < m ̲ [ ℓ ] < k ̲ [ ℓ ] (i.e., the primal iteration stops earlier than the dual iteration), the validity of the stopping criterion (19) for the algebraic solver and the failure of criterion (20) for the inexact Zarantonello symmetrization prove that

Moreover, for 0 ⩽ k < k ′ = m ̲ [ ℓ ] , stability (A1) and the estimate (37) verify

For 0 ⩽ k ⩽ m ̲ [ ℓ ] < k ′ ⩽ k ̲ [ ℓ ] , it follows H ℓ k ′ , n ̲ = H ℓ m ̲ , n ̲ ≲ q sym m ̲ [ ℓ ] − k H ℓ k , n ̲ . Finally, for m ̲ [ ℓ ] ⩽ k < k ′ ⩽ k ̲ [ ℓ ] , we have H ℓ k ′ , n ̲ = H ℓ m ̲ [ ℓ ] , n ̲ = H ℓ k , n ̲ . Notice that the same argumentation holds for the dual quasi-error Z ℓ k , ν ̲ in the remaining cases with μ ̲ [ ℓ ] < k ̲ [ ℓ ] (i.e., the dual iteration stops earlier than the primal iteration).

Since k ̲ [ ℓ ] = m ̲ [ ℓ ] or k ̲ [ ℓ ] = μ [ ℓ ] by definition, we obtain, for all ( ℓ , k ′ , j ̲ ) ∈ Q with 0 ⩽ k ⩽ k ′ ⩽ k ̲ [ ℓ ] ,

H ℓ k ′ , j ̲ ≲ q sym k ′ − k H ℓ k , j ̲ if k ̲ [ ℓ ] = m ̲ [ ℓ ] or Z ℓ k ′ , j ̲ ≲ q sym k ′ − k Z ℓ k , j ̲ if k ̲ [ ℓ ] = μ ̲ [ ℓ ] .

Furthermore, there holds H ℓ k ′ , j ̲ ≲ H ℓ k , j ̲ and Z ℓ k ′ , j ̲ ≲ Z ℓ k , j ̲ in any case. This yields (63) via

H ℓ k ′ , j ̲ Z ℓ k ′ , j ̲ ≲ q sym k ′ − k H ℓ k , j ̲ Z ℓ k , j ̲ ∀ ( ℓ , k ′ , j ̲ ) ∈ Q with 0 ⩽ k ⩽ k ′ ⩽ k ̲ [ ℓ ] ,

where the hidden constant depends only on C _stab, λ _sym, and q _sym.

Nested iteration u ℓ − 1 m ̲ , n ̲ = u ℓ 0 , n ̲ and z ℓ − 1 μ ̲ , ν ̲ = z ℓ 0 , ν ̲ and the estimates (56) yield, for all ( ℓ , 0,0 ) ∈ Q with ℓ > 0,

A multiplication of the two previous estimates proves (64).□

Finally, the quasi-contraction in (ℓ, k) from Lemma 9 together with a quasi-contraction in the algebraic solver index j leads to tail-summability in (ℓ, k, j).

Lemma 10

(quasi-contraction and stability by algebraic solver). There holds

(67) H ℓ k , j ′ Z ℓ k , j ′ ≲ q alg j ′ − j H ℓ k , j Z ℓ k , j ∀ ( ℓ , k , j ′ ) ∈ Q with 0 ⩽ j ⩽ j ′ ⩽ j ̲ [ ℓ , k ]

and, with the abbreviation (m − 1)₊: = max{m − 1, 0},

(68) H ℓ m , 0 ⩽ 3 H ℓ ( m − 1 ) + , n ̲ and Z ℓ μ , 0 ⩽ 3 Z ℓ ( μ − 1 ) + , ν ̲ ∀ ( ℓ , m , 0 ) ∈ Q u , ( ℓ , μ , 0 ) ∈ Q z .

Proof.

We recall that u ℓ 0,0 = u ℓ 0 , n ̲ = u ℓ 0 , ⋆ by definition and, hence, H ℓ 0,0 = H ℓ 0 , n ̲ = H ℓ 0 , j ̲ . Nested iteration u ℓ m , 0 = u ℓ m − 1 , n ̲ implies that

Therewith, we derive (68).

The combination of a posteriori error control (30) for the exact Zarantonello iteration, for the algebraic solver (31), and the failure of the stopping criterion (19) in Algorithm 1(I.b.ii) for the algebraic solver proves, for 0 ⩽ j < j ′ < n ̲ [ ℓ , m ] < j ̲ [ ℓ , m ] ,

For 0 ⩽ j < n ̲ [ ℓ , m ] ⩽ j ′ ⩽ j ̲ [ ℓ , m ] , stability (A1) and contraction of the algebraic solver (15) verify that

For n ̲ [ ℓ , m ] ⩽ j < j ′ ⩽ j ̲ [ ℓ , m ] , it holds that H ℓ m , j = H ℓ m , n ̲ = H ℓ m , j ′ . Since j ̲ [ ℓ , k ] = n ̲ [ ℓ , k ] or j ̲ [ ℓ , k ] = ν ̲ [ ℓ , k ] , we have, for all ( ℓ , k , j ′ ) ∈ Q with 0 ⩽ j ⩽ j ′ ⩽ j ̲ [ ℓ , k ] ,

H ℓ k , j ≲ q alg j − j ′ H ℓ k , j ′ if j ̲ [ ℓ , k ] = n ̲ [ ℓ , k ] or Z ℓ k , j ≲ q alg j − j ′ Z ℓ k , j ′ if j ̲ [ ℓ , k ] = ν ̲ [ ℓ , k ] .

Furthermore, we have H ℓ k , j ≲ H ℓ k , j ′ and Z ℓ k , j ≲ Z ℓ k , j ′ in any case. Hence, we obtain

H ℓ k , j Z ℓ k , j ≲ q alg j − j ′ H ℓ k , j ′ Z ℓ k , j ′ ∀ ( ℓ , k , j ) ∈ Q with 0 ⩽ j ′ ⩽ j ⩽ j ̲ [ ℓ , k ] ,

where the hidden constant depends only on q sym ⋆ , λ _sym, q _alg, λ _alg, and C _stab.□

Ultimately, synthesizing the preceding lemmas yields tail-summability of the quasi-error product and thus leads to the following proof of Theorem 1.

Proof of Theorem 1.

The proof consists of four steps.

Step 1 (tail-summability in mesh level ℓ ). We apply the tail-summability criterion from [33], Lemma 5] to the sequences a _ℓ := H_ℓ Z_ℓ and b ℓ : = q ctr ′ R ℓ . Therein, it is shown that R-linear convergence is equivalent to tail-summability and that, for tail-summability, it is sufficient to guarantee

(70) a ℓ + 1 ⩽ q a ℓ + b ℓ , b ℓ + M ⩽ C 1 a ℓ , ∑ ℓ ′ = ℓ ℓ + M b ℓ 2 ⩽ C 2 ( M + 1 ) 1 − δ a ℓ 2 ∀ ℓ , M ∈ N 0 .

Indeed, contraction up to a remainder from (49), the estimate of the remainder from (50), and the quasi-monotonicity of H_ℓ and Z_ℓ from (60) validate the assumptions of the tail-summability criterion (70) and lead to tail-summability

(71) ∑ ℓ ′ = ℓ + 1 ℓ ̲ − 1 H ℓ ′ Z ℓ ′ ≲ H ℓ Z ℓ ∀ ( ℓ , k ̲ , j ̲ ) ∈ Q .

Step 2 (tail-summability in (ℓ, k)). For ( ℓ , k , j ̲ ) ∈ Q , the estimates (63), (64) and the geometric series prove tail-summability

(72) ∑ ( ℓ ′ , k ′ , j ̲ ) ∈ Q | ℓ ′ , k ′ , j ̲ | > | ℓ , k , j ̲ | H ℓ k ′ , j ̲ Z ℓ k ′ , j ̲ = ∑ k ′ = k + 1 k ̲ [ ℓ ] H ℓ k ′ , j ̲ Z ℓ k ′ , j ̲ + ∑ ℓ ′ = ℓ + 1 ℓ ̲ ∑ k ′ = 0 k ̲ [ ℓ ′ ] H ℓ ′ k ′ , j ̲ Z ℓ ′ k ′ , j ̲ ≲ ( 63 ) H ℓ k , j ̲ Z ℓ k , j ̲ + ∑ ℓ ′ = ℓ + 1 ℓ ̲ H ℓ ′ 0 , j ̲ Z ℓ ′ 0 , j ̲ ≲ ( 64 ) H ℓ k , j ̲ Z ℓ k , j ̲ + ∑ ℓ ′ = ℓ ℓ ̲ − 1 H ℓ ′ Z ℓ ′ ≲ ( 71 ) H ℓ k , j ̲ Z ℓ k , j ̲ + H ℓ Z ℓ ≲ ( 65 ) H ℓ k , j ̲ Z ℓ k , j ̲ .

Step 3 (tail-summability in (ℓ, k, j) ). Finally, for all ( ℓ , k , j ) ∈ Q , we observe that

∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | > | ℓ , k , j | H ℓ ′ k ′ , j ′ Z ℓ ′ k ′ , j ′ = ∑ j ′ = j + 1 j ̲ [ ℓ , k ] H ℓ k , j ′ Z ℓ k , j ′ + ∑ k ′ = k + 1 k ̲ [ ℓ ] ∑ j ′ = 0 j ̲ [ ℓ , k ′ ] H ℓ k ′ , j ′ Z ℓ k ′ , j ′ + ∑ ℓ ′ = ℓ + 1 ℓ ̲ ∑ k ′ = 0 k ̲ [ ℓ ′ ] ∑ j ′ = 0 j ̲ [ ℓ ′ , k ′ ] H ℓ ′ k ′ , j ′ Z ℓ ′ k ′ , j ′ ≲ ( 67 ) H ℓ k , j Z ℓ k , j + ∑ k ′ = k + 1 k ̲ [ ℓ ] H ℓ k ′ , 0 Z ℓ k ′ , 0 + ∑ ℓ ′ = ℓ + 1 ℓ ̲ ∑ k ′ = 0 k ̲ [ ℓ ′ ] H ℓ ′ k ′ , 0 Z ℓ ′ k ′ , 0 ≲ ( 68 ) H ℓ k , j Z ℓ k , j + ∑ ( ℓ ′ , k ′ , j ̲ ) ∈ Q | ℓ ′ , k ′ , j ̲ | > | ℓ , k , j ̲ | H ℓ ′ k ′ , j ̲ Z ℓ ′ k ′ , j ̲ ≲ ( 72 ) H ℓ k , j Z ℓ k , j + H ℓ k , j ̲ Z ℓ k , j ̲ ≲ ( 67 ) H ℓ k , j Z ℓ k , j .

Step 4. Since the index set Q is linearly ordered with respect to the total step counter |⋅, ⋅, ⋅|, tail-summability in Step 3 and the equivalence of tail-summability and R-linear convergence from [33], Lemma 10] conclude the proof of (47) in Theorem 1.□

6 Optimal complexity of Algorithm 1

Full linear convergence (47) has a simple but crucial consequence. Using a geometric series argument, one can prove that the cumulative computational cost up to a given level is bounded by the cost of the said level; see [33], Corollary 14], where only the primal quasi-error H ℓ k , j has to be replaced by the quasi-error product H ℓ k , j Z ℓ k , j . As a consequence, the convergence rates with respect to the number of degrees of freedom (defined as M(r) in (73) below) and the rates with respect to the overall computational cost (cf. (26) and the discussion following the statement of Algorithm 1) coincide.

Corollary 1

(rates = complexity [33], Corollary 14]). Suppose the assumptions of Theorem 1. For all r > 0, the output ( T ℓ ) ℓ ∈ N 0 of Algorithm 1 satisfies

(73) M ( r ) : = sup ( ℓ , k , j ) ∈ Q # T ℓ r H ℓ k , j Z ℓ k , j ⩽ sup ( ℓ , k , j ) ∈ Q ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k , j | # T ℓ ′ r H ℓ k , j Z ℓ k , j ⩽ C cost ( r ) M ( r )

with the constant C cost ( r ) : = C lin / ( 1 − q lin 1 / r ) r > 0 .

While Theorem 1 only concerns R-linear convergence, a sufficiently small choice of the adaptivity parameters θ, λ _sym, and λ _alg even guarantees the optimal convergence rate r = s + t with respect to computational cost, i.e., the overall computational time. Here, we suppose that the primal solution u ^⋆ to (5) can be approximated at rate s and the dual solution z ^⋆ to (8) can be approximated at rate t. To formalize this idea, we introduce the notion of approximation classes [3], [4], [5], [9]. For s, t > 0, define

‖ u ⋆ ‖ A s : = sup N ∈ N 0 N + 1 s min T opt ∈ T N η opt u opt ⋆ , ‖ z ⋆ ‖ A t : = sup N ∈ N 0 N + 1 t min T opt ∈ T N ζ opt z opt ⋆ ,

where η _opt(⋅) and ζ _opt(⋅) denote the estimator values for the exact discrete solutions u opt ⋆ and z opt ⋆ on the unavailable optimal triangulations T opt ∈ T N ( T ) . We stress that ‖ u ⋆ ‖ A s and ‖ z ⋆ ‖ A t can equivalently be defined by energy error plus data oscillations [8], [9].

Theorem 2

(optimal complexity). Suppose that the estimators η and ζ satisfy (A1)–(A3⁺) and (QM) and suppose quasi-orthogonality (A4). Recall λ alg ⋆ from Lemma 3 and λ ^⋆ from (46) in Theorem 1. Define the constants

(74) λ sym ⋆ : = min 1 , C stab − 1 C alg − 1 ⩽ 1 with C alg : = 1 1 − q sym ⋆ 2 q alg 1 − q alg λ alg ⋆ + q sym ⋆ , θ ⋆ : = ( 1 + C stab 2 C rel 2 ) − 1 < 1 .

Suppose that θ, λ _sym, and λ _alg are sufficiently small in the sense of

(75) 0 < λ alg ⩽ λ alg ⋆ , 0 < λ sym < λ sym ⋆ , λ alg λ sym < λ ⋆ , 0 < θ mark : = ( θ 1 / 2 + λ sym / λ sym ⋆ ) 2 ( 1 − λ sym / λ sym ⋆ ) 2 < θ ⋆ < 1 .

Then, Algorithm 1 guarantees, for all s, t > 0, that

(76) sup ( ℓ , k , j ) ∈ Q ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k , j | # T ℓ ′ s + t H ℓ k , j Z ℓ k , j ⩽ C opt max ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t , H 0 0,0 Z 0 0,0 .

The constant C _opt depends only on C _stab, C _rel, C _drel, C _mark, C _mesh, C _lin, q _lin, # T 0 , and s + t. In particular, there holds optimal complexity of Algorithm 1.

The proof of Theorem 2 employs the following result from [50] providing estimator equivalence between the (unavailable) estimators for the exact discrete solutions u ℓ ⋆ , z ℓ ⋆ and the estimators at the computed approximations u ℓ m ̲ , n ̲ , z ℓ μ ̲ , ν ̲ .

Lemma 11

(estimator equivalence [50], Lemma 15]). Recall the constants λ sym ⋆ , C _alg > 0 from (74) and λ alg ⋆ > 0 from Lemma 3. Then, for all 0 < θ ⩽ 1, 0 < λ alg ⩽ λ alg ⋆ , 0 < λ sym < λ sym ⋆ , it holds that

(77) 1 − λ sym / λ sym ⋆ η ℓ u ℓ m ̲ , n ̲ ⩽ η ℓ u ℓ ⋆ ⩽ 1 + λ sym / λ sym ⋆ η ℓ u ℓ m ̲ , n ̲ ∀ ( ℓ , m ̲ , n ̲ ) ∈ Q u , 1 − λ sym / λ sym ⋆ ζ ℓ z ℓ μ ̲ , ν ̲ ⩽ ζ ℓ z ℓ ⋆ ⩽ 1 + λ sym / λ sym ⋆ ζ ℓ z ℓ μ ̲ , ν ̲ ∀ ( ℓ , μ ̲ , ν ̲ ) ∈ Q z .

Proof of Theorem 2.

By Corollary 1, it suffices to prove that, for any s, t > 0,

(78) sup ( ℓ , k , j ) ∈ Q # T ℓ s + t H ℓ k , j Z ℓ k , j ≲ max ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t , H 0 0,0 Z 0 0,0 .

Since the inequality becomes trivial if either ‖ u ⋆ ‖ A s = ∞ or ‖ z ⋆ ‖ A t = ∞ , we may assume ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t < ∞ . The proof consists of three steps.

Step 1. With 0 < θ mark : = ( θ 1 / 2 + λ sym / λ sym ⋆ ) 2 ( 1 − λ sym / λ sym ⋆ ) − 2 < θ ⋆ , the validity of (A3⁺) for both estimators and [16], Lemma 14] guarantee the existence of sets R ℓ ′ ⊆ T ℓ ′ with 0 ⩽ ℓ ′ < ℓ ̲ such that

(79a) # R ℓ ′ ≲ ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t 1 / ( s + t ) η ℓ ′ u ℓ ′ ⋆ ζ ℓ ′ z ℓ ′ ⋆ − 1 / ( s + t ) ,

(79b) θ mark η ℓ ′ u ℓ ′ ⋆ ⩽ η ℓ ′ R ℓ ′ , u ℓ ′ ⋆ or θ mark ζ ℓ ′ z ℓ ′ ⋆ ⩽ ζ ℓ ′ R ℓ ′ , z ℓ ′ ⋆ .

For 0 ⩽ ℓ ′ < ℓ ̲ , the estimator equivalence (77) in Lemma 11 leads to

1 − λ sym / λ sym ⋆ η ℓ ′ u ℓ ′ m ̲ , n ̲ ⩽ η ℓ ′ u ℓ ′ ⋆ , 1 − λ sym / λ sym ⋆ ζ ℓ ′ z ℓ ′ μ ̲ , ν ̲ ⩽ ζ ℓ ′ z ℓ ′ ⋆

and consequently with (79a) to

(80) # R ℓ ′ ≲ ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t 1 / ( s + t ) η ℓ ′ u ℓ ′ m ̲ , n ̲ ζ ℓ ′ z ℓ ′ μ ̲ , ν ̲ − 1 / ( s + t ) .

Note that the stopping criteria (20) and (22) lead to

and with (64) to

(81) H ℓ ′ + 1 0 , j ̲ Z ℓ ′ + 1 0 , j ̲ ≲ ( 64 ) H ℓ ′ Z ℓ ′ ≲ η ℓ ′ u ℓ ′ m ̲ , n ̲ ζ ℓ ′ z ℓ ′ m ̲ u , ν ̲ .

Hence, the combination of (80) and (81) reads

(82) # R ℓ ′ ≲ ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t 1 / ( s + t ) H ℓ ′ + 1 0 , j ̲ Z ℓ ′ + 1 0 , j ̲ − 1 / ( s + t ) .

Step 2. Recall from [29], Theorem 8] that the set R ℓ ′ satisfies the Dörfler criterion from Algorithm 1 with the same parameter θ. The quasi-minimality of M ℓ ′ implies

(83) # M ℓ ′ ⩽ C mark # R ℓ ′ ∀ 0 ⩽ ℓ ′ < ℓ ̲

with the constant C _mark ⩾ 1 from Algorithm 1.

Step 3. Let ( ℓ , k , j ) ∈ Q . Full linear convergence (47) from Theorem 1 yields that

(84) ∑ ( ℓ ' , k ' , j ' ) ∈ Q | ℓ ' , k ' , j ' | ⩽ | ℓ , k , j | H ℓ ' k ' , j ' Z ℓ ' k ' , j ' − 1 / ( s + t ) ≲ ( 47 ) H ℓ k , j Z ℓ k , j − 1 / ( s + t ) ∑ ( ℓ ' , k ' , j ' ) ∈ Q | ℓ ' , k ' , j ' | ⩽ | ℓ , k , j | ( q lin 1 / s ) | ℓ , k , j | − | ℓ ' , k ' , j ' | ≲ H ℓ k , j Z ℓ k , j − 1 / ( s + t ) .

NVB refinement satisfies the mesh-closure estimate [9], Eq. (2.9)] reading,

(85) # T ℓ − # T 0 ⩽ C mesh ∑ ℓ ' = 0 ℓ − 1 # M ℓ ' ∀ ℓ ⩾ 1 ,

where C _mesh > 1 depends only on T 0 . Thus, for ( ℓ , k , j ) ∈ Q , we have by the mesh-closure estimate (85), quasi-optimality of Dörfler marking (83), and the result (84) that

# T ℓ − # T 0 ≲ ( 85 ) ∑ ℓ ′ = 0 ℓ − 1 # M ℓ ′ ≲ ( 83 ) ∑ ℓ ′ = 0 ℓ − 1 # R ℓ ′ ≲ ( 82 ) ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t 1 / ( s + t ) ∑ ℓ ′ = 0 ℓ − 1 H ℓ ′ + 1 0 , j ̲ Z ℓ ′ + 1 0 , j ̲ − 1 / ( s + t ) ⩽ ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t 1 / ( s + t ) ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k , j | H ℓ ′ k ′ , j ′ Z ℓ ′ k ′ , j ′ − 1 / ( s + t ) ≲ ( 84 ) ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t 1 / ( s + t ) H ℓ k , j Z ℓ k , j − 1 / ( s + t ) .

Rearranging the terms and noting that 1 ⩽ # T ℓ − # T 0 implies # T ℓ − # T 0 + 1 ⩽ 2 ( # T ℓ − # T 0 ) , we obtain, for ℓ > 0, that

(86a) ( # T ℓ − # T 0 + 1 ) s + t H ℓ k , j Z ℓ k , j ≲ ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t .

Moreover, full linear convergence (47) proves that

(86b) ( # T 0 − # T 0 + 1 ) s + t H 0 k , j Z 0 k , j = H 0 k , j Z 0 k , j ≲ H 0 0,0 Z 0 0,0 .

We recall from [35], Lemma 22] that, for all T ℓ ∈ T , it holds

(87) # T ℓ − # T 0 + 1 ⩽ # T ℓ ⩽ # T 0 ( # T ℓ − # T 0 + 1 ) .

This shows, for all ( ℓ , k , j ) ∈ Q ,

( # T ℓ ) s + t H ℓ k , j Z ℓ k , j ≲ ( 87 ) ( # T ℓ − # T 0 + 1 ) s + t H ℓ k , j Z ℓ k , j ≲ ( 86 ) max ‖ u ⋆ ‖ A s ‖ z ⋆ ‖ A t , H 0 0,0 Z 0 0,0

and concludes the proof of (78).□

7 Numerical examples

In this section, we present numerical experiments using the open source software package MooAFEM [51].^[1] In the following, Steps (I) and (II) of Algorithm 1 employ the optimal hp-robust local multigrid method from [32] as an algebraic solver. If not explicitly stated otherwise, we choose the parameters θ = 0.5, δ = 0.5, λ _sym = λ _alg = 0.7 in Algorithm 1 throughout the numerical experiments.

7.1 Singularity in the goal functional

The first model problem is a nonsymmetric variant of the benchmark problem from [29], Section 4.1] with a singularity only in the goal functional. On the unit square Ω = ( 0,1 ) 2 ⊂ R 2 , we consider

(88) − Δ u ⋆ + x ⋅ ∇ u ⋆ + u ⋆ = f in Ω subject to u ⋆ = 0 on ∂ Ω ,

where the right-hand side is chosen such that the exact solution u ^⋆ reads

u ⋆ ( x ) = x 1 x 2 ( 1 − x 1 ) ( 1 − x 2 ) .

Consider g = 0 and g = χ _K (1,0)^⊤ in the quantity of interest

G ( u ⋆ ) : = ∫ K ∂ x 1 u ⋆ d x = 11 960 with K : = c o n v { ( 1 / 2,1 ) , ( 1,1 / 2 ) , ( 1,1 ) } .

Figure 2 (left) displays a mesh generated by Algorithm 1 and the support K of g . The error estimator captures and resolves the two point singularities induced by G.

$Figure 2: Left: Mesh T 15 ${\mathcal{T}}_{15}$ for the problem (88) generated by Algorithm 1 with # T 15 = 2315 $\#{\mathcal{T}}_{15}=2315$ . Right: Mesh T 18 ${\mathcal{T}}_{18}$ for the problem (89) with # T 18 = 2130 $\#{\mathcal{T}}_{18}=2130$ , where the Dirichlet boundary part Γ D is marked by red solid lines and the Neumann boundary part Γ N by green dashed lines.$

Figure 2:

Left: Mesh T 15 for the problem (88) generated by Algorithm 1 with # T 15 = 2315 . Right: Mesh T 18 for the problem (89) with # T 18 = 2130 , where the Dirichlet boundary part Γ_D is marked by red solid lines and the Neumann boundary part Γ_N by green dashed lines.

7.2 Geometric singularity and strong convection

The second benchmark problem investigates Ω = ( − 1,1 ) 2 \ c o n v { ( 0,0 ) , ( − 1,0 ) , ( − 1 , − 1 ) } ⊂ R 2 with the Dirichlet boundary Γ_D = conv{(−1, 0), (0, 0)} ∪ conv{(0, 0), (−1, − 1)} and Neumann boundary Γ_N = ∂Ω \Γ_D; see Figure 2 (right) for a visualization of the geometry. We consider

(89) − Δ u ⋆ + ( 5,5 ) ⊤ ⋅ ∇ u ⋆ = 1 in Ω subject to u ⋆ = 0 on Γ D and ∇ u ⋆ ⋅ n = 0 on Γ N .

Consider g = 0 and g = χ _S (1,1)^⊤ in the quantity of interest

G ( u ⋆ ) = ∫ S ∂ x 1 u ⋆ + ∂ x 2 u ⋆ d x with S : = ( − 1 / 2,1 / 2 ) 2 ∩ Ω .

The exact solution u ^⋆ is not known analytically in this case so that we do not have access to the exact goal error | G ( u ⋆ ) − G ℓ u ℓ m ̲ , n ̲ , z ℓ μ ̲ , ν ̲ | . Figure 2 (right) shows a mesh generated by Algorithm 1 as well as the configuration, i.e., the support S of g in blue, the Dirichlet boundary in red solid lines, and the Neumann boundary in green dashed lines.

7.2.1 Optimality of Algorithm 1

Figure 3 displays the estimator product η ℓ u ℓ m ̲ , n ̲ ζ ℓ ( z μ ̲ , ν ̲ ) and the goal error | G ( u ⋆ ) − G ℓ u ℓ m ̲ , n ̲ , z ℓ μ ̲ , ν ̲ | from (17) for the problem (88), due to higher-order approximations, we only show results prior to machine precision. For all investigated polynomial degrees p, the goal error and the estimator product are indeed equivalent. Algorithm 1 achieves the optimal rate −p with respect to the cumulative computational work and with respect to the cumulative computational time in Figure 3 for problem (88) and Figure 4 for problem (89). Figure 5 shows that the proposed algorithm indeed achieves linear complexity and is substantially faster than the Matlab built-in direct solver as the slightly larger slope of the latter indicates super-linear complexity. Table 2 displays the weighted costs

(90) η ℓ u ℓ m ̲ , n ̲ ζ ℓ z ℓ μ ̲ , ν ̲ ( ∑ ( ℓ ′ , k ′ , j ′ ) ∈ Q | ℓ ′ , k ′ , j ′ | ⩽ | ℓ , k ̲ , j ̲ | t i m e ( ℓ ′ , k ′ , j ′ ) ( p

of Algorithm 1 for polynomial degree p = 2 with t i m e ( ℓ ′ , k ′ , j ′ ) in seconds and highlights the corresponding optimal choices of the parameters. This justifies the selection of θ = 0.5 together with larger symmetrization parameter λ _sym = 0.7, and algebraic solver parameter λ _alg = 0.7. The table for the second benchmark problem from (89) leads to similar results and is therefore omitted. While the choice of the damping parameter 0 < δ < 2α/L ² in (13) is crucial in the case of large convection to guarantee the contraction property (14), the adaptivity parameters appear more robust with respect to other coefficients in (4).

$Figure 3: Convergence history plot of estimator product η ℓ u ℓ m ̲ , n ̲ ζ ℓ ( z μ ̲ , ν ̲ ) ${\eta }_{\ell }\left({u}_{\ell }^{\underline{m},\underline{n}}\right) {\zeta }_{\ell }\left({z}^{\underline{\mu },\underline{\nu }}\right)$ indicated by bullets and goal error from (17) indicated by diamonds with respect to the cumulative computational work (left) and with respect to the cumulative computational time (right) for the benchmark problem (88).$

Figure 3:

Convergence history plot of estimator product η ℓ u ℓ m ̲ , n ̲ ζ ℓ ( z μ ̲ , ν ̲ ) indicated by bullets and goal error from (17) indicated by diamonds with respect to the cumulative computational work (left) and with respect to the cumulative computational time (right) for the benchmark problem (88).

$Figure 4: Convergence history plot of estimator product η ℓ u ℓ m ̲ , n ̲ ζ ℓ ( z μ ̲ , ν ̲ ) ${\eta }_{\ell }\left({u}_{\ell }^{\underline{m},\underline{n}}\right) {\zeta }_{\ell }\left({z}^{\underline{\mu },\underline{\nu }}\right)$ with respect to the cumulative computational cost (left) and the cumulative computational time (right) for the benchmark problem (89).$

Figure 4:

Convergence history plot of estimator product η ℓ u ℓ m ̲ , n ̲ ζ ℓ ( z μ ̲ , ν ̲ ) with respect to the cumulative computational cost (left) and the cumulative computational time (right) for the benchmark problem (89).

Figure 5:

Comparison of cumulative time of the local multigrid solver with the Matlab built-in direct solver mldivide with respect to the cumulative computational cost for the benchmark problem (89).

Table 2:

Optimal selection of parameters with respect to the cumulative computational costs (overall computation time in seconds) for the experiment (88) with fixed polynomial degree p = 2 and δ = 0.5. For comparison, the table displays the value of the weighted costs from (90) (in 10⁻⁷) with overall stopping criterion η ℓ u ℓ m ̲ , n ̲ ζ ℓ u ℓ μ ̲ , ν ̲ < 5 ⋅ 1 0 − 10 for various choices of λ _sym, λ _alg, and θ. For each θ-block, we mark the row-wise optimal values in blue, the column-wise optimal values in yellow, and in green if both optimal values coincide.

×10⁻⁷	θ = 0.1					θ = 0.3					θ = 0.5
λ _alg/λ _sym	0.1	0.3	0.5	0.7	0.9	0.1	0.3	0.5	0.7	0.9	0.1	0.3	0.5	0.7	0.9
0.1	38.7	33.4	29.6	22.1	24.4	10.2	5.12	4.90	4.83	4.74	6.18	4.48	4.66	4.89	5.25
0.3	36.2	24.7	24.5	21.8	23.1	7.28	4.98	3.53	3.27	3.26	4.18	4.54	4.79	5.01	5.13
0.5	24.3	24.7	24.7	23.4	23.6	5.84	3.64	3.39	3.27	3.37	3.41	2.71	2.52	2.49	2.68
0.7	24.1	24.8	23.8	22.2	24.0	4.95	3.59	3.30	3.25	3.42	2.74	2.35	2.41	2.24	2.46
0.9	23.5	24.6	22.3	24.4	23.8	4.90	3.58	3.29	3.26	3.41	2.81	2.30	2.43	2.27	2.41
	θ = 0.7					θ = 0.8					θ = 0.9
0.1	5.82	5.18	5.43	5.40	5.93	8.53	6.10	7.31	6.67	7.77	11.6	8.86	9.12	9.87	9.97
0.3	4.65	4.86	5.35	5.98	6.67	6.27	5.92	7.20	7.46	7.57	8.62	8.40	9.27	10.6	11.5
0.5	3.69	2.89	2.88	2.95	3.13	5.09	3.61	3.66	3.63	3.66	7.27	5.32	4.84	4.93	5.12
0.7	2.99	2.56	2.64	2.62	2.89	3.75	3.12	3.23	3.03	3.11	4.58	3.95	4.04	4.43	4.79
0.9	2.89	2.49	2.65	2.66	2.89	3.79	3.11	3.19	3.13	3.27	4.67	4.06	4.16	4.35	4.61

Finally, in Figure 6, we display the number of total solver steps | ℓ , m ̲ , n ̲ | − | ℓ , 0,0 | resp. | ℓ , μ ̲ , ν ̲ | − | ℓ , 0,0 | on each mesh level for both benchmark problems (88) and (89). The plots show that the two iterations often stop after the same number of steps.

$Figure 6: Number of total solver steps | ℓ , m ̲ , n ̲ | − | ℓ , 0,0 | $\vert \ell ,\underline{m},\underline{n}\vert -\vert \ell ,0,0\vert $ resp. | ℓ , μ ̲ , ν ̲ | − | ℓ , 0,0 | $\vert \ell ,\underline{\mu },\underline{\nu }\vert -\vert \ell ,0,0\vert $ on each mesh level for the benchmark problems (88) (left) and (89) (right).$

Figure 6:

Number of total solver steps | ℓ , m ̲ , n ̲ | − | ℓ , 0,0 | resp. | ℓ , μ ̲ , ν ̲ | − | ℓ , 0,0 | on each mesh level for the benchmark problems (88) (left) and (89) (right).

8 Summary

In this work, we developed a cost-optimal goal-oriented adaptive finite element method for the efficient computation of the quantity of interest G(u ^⋆) with solution u ^⋆ to the general second-order linear elliptic partial differential equation (4). Since the current analysis of iterative algebraic solvers for nonsymmetric systems with optimal preconditioner only leads to contraction of the residual in a vector norm, we proposed a nested iterative solver for the primal and dual problem in parallel. The strategy consists of the Zarantonello iteration (13) as an outer solver loop and an optimal multigrid solver for the arising SPD system as an innermost solver loop. In recent own work [33], we have shown that the link between convergence rates with respect to the degrees of freedom and the total computational cost is full linear convergence of the quasi-error H ℓ k , j Z ℓ k , j . To this end, Theorem 1 shows that the proposed adaptive algorithm contracts (up to a multiplicative constant) the quasi-error product H ℓ k , j Z ℓ k , j in every step, independently of the algorithmic decision to employ mesh refinement, symmetrization, or the algebraic solver. A particular problem in the analysis is that the nested iterative solver procedure only guarantees contraction as long as 1 ⩽ k < k ̲ [ ℓ ] , whereas contraction for the final iterate is only guaranteed up to an estimator term (cf. (29)). Another difficulty arises from the nonsymmetric setting with a quasi-Pythagorean estimate (18) replacing the usual Pythagorean estimate. Therefore, the proof of Theorem 1 employs the equivalence of R-linear convergence and tail-summability of the quasi-error product H ℓ k , j Z ℓ k , j and leads to mild restriction on the product λ _sym λ _alg of the involved solver stopping parameters. The key ingredients to cost-optimality are an adaptive mesh-refinement algorithm with optimal convergence rate with respect to the number of degrees of freedom (under the assumption of exact solution) and an algebraic solver for the linear system of equations that is contractive with respect to the underlying Sobolev norm. In this regard, the analysis in this paper may guide the generalization to conforming discretizations of vector-valued elliptic problems. Finally, the numerical experiments in Section 7 suggest that the proposed strategy allows for large stopping parameter in practice and that a larger choice is beneficial in terms of total runtime. Admittedly, the development of an optimal solver for the nonsymmetric problem (10) would allow to prove full linear convergence with an arbitrary selection of the stopping parameter.

Corresponding author: Julian Streitberger, TU Wien, Institute of Analysis and Scientific Computing, Wien, Austria, E-mail: julian.streitberger@asc.tuwien.ac.at

Research ethics: Not applicable.
Author contributions: All authors have 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.
Competing interests: The authors state no conflict of interest.
Research funding: This research was funded in whole or in part by the Austrian Science Fund (FWF) through [10.55776/F65], [10.55776/I6802], and [10.55776/P33216]. For open access purposes, the author has applied a CC BY public copyright license to any author accepted manuscript version arising from this submission. Additionally, Maximilian Brunner and Julian Streitberger are supported by the Vienna School of Mathematics.
Data availability: Not applicable.

References

[1] W. Dörfler, “A convergent adaptive algorithm for Poisson’s equation,” SIAM J. Numer. Anal., vol. 33, no. 3, pp. 1106–1124, 1996. https://doi.org/10.1137/0733054.Search in Google Scholar

[2] P. Morin, R. H. Nochetto, and K. G. Siebert, “Data oscillation and convergence of adaptive FEM,” SIAM J. Numer. Anal., vol. 38, no. 2, pp. 466–488, 2000. https://doi.org/10.1137/s0036142999360044.Search in Google Scholar

[3] P. Binev, W. Dahmen, and R. DeVore, “Adaptive finite element methods with convergence rates,” Numer. Math., vol. 97, no. 2, pp. 219–268, 2004. https://doi.org/10.1007/s00211-003-0492-7.Search in Google Scholar

[4] R. Stevenson, “Optimality of a standard adaptive finite element method,” Found. Comput. Math., vol. 7, no. 2, pp. 245–269, 2007. https://doi.org/10.1007/s10208-005-0183-0.Search in Google Scholar

[5] J. M. Cascón, C. Kreuzer, R. H. Nochetto, and K. G. Siebert, “Quasi-optimal convergence rate for an adaptive finite element method,” SIAM J. Numer. Anal., vol. 46, no. 5, pp. 2524–2550, 2008. https://doi.org/10.1137/07069047x.Search in Google Scholar

[6] C. Kreuzer and K. G. Siebert, “Decay rates of adaptive finite elements with Dörfler marking,” Numer. Math., vol. 117, no. 4, pp. 679–716, 2011. https://doi.org/10.1007/s00211-010-0324-5.Search in Google Scholar

[7] J. M. Cascón and R. H. Nochetto, “Quasioptimal cardinality of AFEM driven by nonresidual estimators,” IMA J. Numer. Anal., vol. 32, no. 1, pp. 1–29, 2012. https://doi.org/10.1093/imanum/drr014.Search in Google Scholar

[8] M. Feischl, T. Führer, and D. Praetorius, “Adaptive FEM with optimal convergence rates for a certain class of nonsymmetric and possibly nonlinear problems,” SIAM J. Numer. Anal., vol. 52, no. 2, pp. 601–625, 2014. https://doi.org/10.1137/120897225.Search in Google Scholar

[9] C. Carstensen, M. Feischl, M. Page, and D. Praetorius, “Axioms of adaptivity,” Comput. Math. Appl., vol. 67, no. 6, pp. 1195–1253, 2014. https://doi.org/10.1016/j.camwa.2013.12.003.Search in Google Scholar PubMed PubMed Central

[10] R. Becker and R. Rannacher, “An optimal control approach to a posteriori error estimation in finite element methods,” Acta Numer., vol. 10, pp. 1–102, 2001. https://doi.org/10.1017/s0962492901000010.Search in Google Scholar

[11] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Basel, Springer Science & Business Media, 2003.10.1007/978-3-0348-7605-6Search in Google Scholar

[12] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, “Introduction to adaptive methods for differential equations,” Acta Numer., vol. 4, pp. 105–158, 1995. https://doi.org/10.1017/s0962492900002531.Search in Google Scholar

[13] M. B. Giles and E. Süli, “Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality,” Acta Numer., vol. 11, pp. 145–236, 2002. https://doi.org/10.1017/cbo9780511550140.003.Search in Google Scholar

[14] M. S. Mommer and R. Stevenson, “A goal-oriented adaptive finite element method with convergence rates,” SIAM J. Numer. Anal., vol. 47, no. 2, pp. 861–886, 2009. https://doi.org/10.1137/060675666.Search in Google Scholar

[15] R. Becker, E. Estecahandy, and D. Trujillo, “Weighted marking for goal-oriented adaptive finite element methods,” SIAM J. Numer. Anal., vol. 49, no. 6, pp. 2451–2469, 2011. https://doi.org/10.1137/100794298.Search in Google Scholar

[16] M. Feischl, G. Gantner, A. Haberl, D. Praetorius, and T. Führer, “Adaptive boundary element methods for optimal convergence of point errors,” Numer. Math., vol. 132, no. 3, pp. 541–567, 2016. https://doi.org/10.1007/s00211-015-0727-4.Search in Google Scholar

[17] M. Feischl, D. Praetorius, and K. G. van der Zee, “An abstract analysis of optimal goal-oriented adaptivity,” SIAM J. Numer. Anal., vol. 54, no. 3, pp. 1423–1448, 2016. https://doi.org/10.1137/15m1021982.Search in Google Scholar

[18] M. Holst and S. Pollock, “Convergence of goal-oriented adaptive finite element methods for nonsymmetric problems,” Numer. Methods Partial Differ. Equ., vol. 32, no. 2, pp. 479–509, 2016. https://doi.org/10.1002/num.22002.Search in Google Scholar

[19] R. Becker, M. Innerberger, and D. Praetorius, “Optimal convergence rates for goal-oriented FEM with quadratic goal functional,” Comput. Methods Appl. Math., vol. 21, no. 2, pp. 267–288, 2021. https://doi.org/10.1515/cmam-2020-0044.Search in Google Scholar

[20] R. Becker, M. Brunner, M. Innerberger, J. M. Melenk, and D. Praetorius, “Rate-optimal goal-oriented adaptive FEM for semilinear elliptic PDEs,” Comput. Math. Appl., vol. 118, pp. 18–35, 2022. https://doi.org/10.1016/j.camwa.2022.05.008.Search in Google Scholar

[21] B. Endtmayer, U. Langer, and T. Wick, “Multigoal-oriented error estimates for non-linear problems,” J. Numer. Math., vol. 27, no. 4, pp. 215–236, 2019. https://doi.org/10.1515/jnma-2018-0038.Search in Google Scholar

[22] B. Endtmayer, U. Langer, and T. Wick, “Two-side a posteriori error estimates for the dual-weighted residual method,” SIAM J. Sci. Comput., vol. 42, no. 1, pp. A371–A394, 2020. https://doi.org/10.1137/18m1227275.Search in Google Scholar

[23] V. Dolejší, O. Bartoš, and F. Roskovec, “Goal-oriented mesh adaptation method for nonlinear problems including algebraic errors,” Comput. Math. Appl., vol. 93, pp. 178–198, 2021. https://doi.org/10.1016/j.camwa.2021.04.004.Search in Google Scholar

[24] A. Cohen, W. Dahmen, and R. DeVore, “Adaptive wavelet methods for elliptic operator equations: convergence rates,” Math. Comput., vol. 70, no. 233, pp. 27–75, 2001. https://doi.org/10.1090/s0025-5718-00-01252-7.Search in Google Scholar

[25] A. Cohen, W. Dahmen, and R. DeVore, “Adaptive wavelet schemes for nonlinear variational problems,” SIAM J. Numer. Anal., vol. 41, no. 5, pp. 1785–1823, 2003. https://doi.org/10.1137/s0036142902412269.Search in Google Scholar

[26] C. Carstensen and J. Gedicke, “An adaptive finite element eigenvalue solver of asymptotic quasi-optimal computational complexity,” SIAM J. Numer. Anal., vol. 50, no. 3, pp. 1029–1057, 2012. https://doi.org/10.1137/090769430.Search in Google Scholar

[27] G. Gantner, A. Haberl, D. Praetorius, and S. Schimanko, “Rate optimality of adaptive finite element methods with respect to overall computational costs,” Math. Comput., vol. 90, no. 331, pp. 2011–2040, 2021. https://doi.org/10.1090/mcom/3654.Search in Google Scholar

[28] M. Brunner, M. Innerberger, A. Miraçi, D. Praetorius, J. Streitberger, and P. Heid, “Adaptive FEM with quasi-optimal overall cost for nonsymmetric linear elliptic PDEs,” IMA J. Numer. Anal., vol. 44, no. 3, pp. 1560–1596, 2024. https://doi.org/10.1093/imanum/drad039.Search in Google Scholar

[29] R. Becker, G. Gantner, M. Innerberger, and D. Praetorius, “Goal-oriented adaptive finite element methods with optimal computational complexity,” Numer. Math., vol. 153, no. 1, pp. 111–140, 2023. https://doi.org/10.1007/s00211-022-01334-8.Search in Google Scholar PubMed PubMed Central

[30] J. Wu and H. Zheng, “Uniform convergence of multigrid methods for adaptive meshes,” Appl. Numer. Math., vol. 113, pp. 109–123, 2017. https://doi.org/10.1016/j.apnum.2016.11.005.Search in Google Scholar

[31] L. Chen, R. H. Nochetto, and J. Xu, “Optimal multilevel methods for graded bisection grids,” Numer. Math., vol. 120, no. 1, pp. 1–34, 2012. https://doi.org/10.1007/s00211-011-0401-4.Search in Google Scholar

[32] M. Innerberger, A. Miraçi, D. Praetorius, and J. Streitberger, “hp-robust multigrid solver on locally refined meshes for FEM discretizations of symmetric elliptic PDEs,” ESAIM: Math. Modell. Numer. Anal., vol. 58, no. 1, pp. 247–272, 2024. https://doi.org/10.1051/m2an/2023104.Search in Google Scholar

[33] P. Bringmann, M. Feischl, A. Miraçi, D. Praetorius, and J. Streitberger, “On full linear convergence and optimal complexity of adaptive FEM with inexact solver,”arXiv:2311.15738, 2023. https://doi.org/10.48550/arXiv.2311.15738.Search in Google Scholar

[34] M. Feischl, “Inf-sup stability implies quasi-orthogonality,” Math. Comput., vol. 91, no. 337, pp. 2059–2094, 2022. https://doi.org/10.1090/mcom/3748.Search in Google Scholar

[35] A. Bespalov, A. Haberl, and D. Praetorius, “Adaptive FEM with coarse initial mesh guarantees optimal convergence rates for compactly perturbed elliptic problems,” Comput. Methods Appl. Mech. Eng., vol. 317, pp. 318–340, 2017. https://doi.org/10.1016/j.cma.2016.12.014.Search in Google Scholar

[36] E. Zarantonello, “Solving functional equations by contractive averaging, math,” Res. Cent. Rep., vol. 160, 1960.Search in Google Scholar

[37] S. Congreve and T. P. Wihler, “Iterative Galerkin discretizations for strongly monotone problems,” J. Comput. Appl. Math., vol. 311, pp. 457–472, 2017. https://doi.org/10.1016/j.cam.2016.08.014.Search in Google Scholar

[38] G. Gantner, A. Haberl, D. Praetorius, and B. Stiftner, “Rate optimal adaptive FEM with inexact solver for nonlinear operators,” IMA J. Numer. Anal., vol. 38, no. 4, pp. 1797–1831, 2018. https://doi.org/10.1093/imanum/drx050.Search in Google Scholar

[39] A. Haberl, D. Praetorius, S. Schimanko, and M. Vohralík, “Convergence and quasi-optimal cost of adaptive algorithms for nonlinear operators including iterative linearization and algebraic solver,” Numer. Math., vol. 147, no. 3, pp. 679–725, 2021. https://doi.org/10.1007/s00211-021-01176-w.Search in Google Scholar

[40] E. Zeidler, Nonlinear Functional Analysis and its Applications. Part II/A, New York, Springer-Verlag, 1990.10.1007/978-1-4612-0981-2Search in Google Scholar

[41] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. Philadelphia, PA, Society for Industrial and Applied Mathematics, 2003.10.1137/1.9780898718003Search in Google Scholar

[42] Y. Saad and M. H. Schultz, “GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM J. Sci. Comput., vol. 7, no. 3, pp. 856–869, 1986. https://doi.org/10.1137/0907058.Search in Google Scholar

[43] R. Stevenson, “The completion of locally refined simplicial partitions created by bisection,” Math. Comput., vol. 77, no. 261, pp. 227–241, 2008. https://doi.org/10.1090/s0025-5718-07-01959-x.Search in Google Scholar

[44] M. Aurada, M. Feischl, T. Führer, M. Karkulik, and D. Praetorius, “Energy norm based error estimators for adaptive BEM for hypersingular integral equations,” Appl. Numer. Math., vol. 95, pp. 15–35, 2015. https://doi.org/10.1016/j.apnum.2013.12.004.Search in Google Scholar

[45] M. Karkulik, D. Pavlicek, and D. Praetorius, “On 2D newest vertex bisection: optimality of mesh-closure and H1-stability of L2-projection,” Constr. Approx., vol. 38, no. 2, pp. 213–234, 2013. https://doi.org/10.1007/s00365-013-9192-4.Search in Google Scholar

[46] L. Diening, L. Gehring, and J. Storn, “Adaptive Mesh Refinement for Arbitrary Initial Triangulations,”arXiv.2306.02674, 2023.Search in Google Scholar

[47] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Ser. Pure and Applied Mathematics, New York, Wiley-Interscience, 2000.10.1002/9781118032824Search in Google Scholar

[48] R. Verfürth, “A posteriori error estimation and adaptive mesh-refinement techniques,” in Proceedings of the Fifth International Congress on Computational and Applied Mathematics (Leuven, 1992), vol. 50, pp. 67–83, 1994.10.1016/0377-0427(94)90290-9Search in Google Scholar

[49] C.-M. Pfeiler and D. Praetorius, “Dörfler marking with minimal cardinality is a linear complexity problem,” Math. Comput., vol. 89, no. 326, pp. 2735–2752, 2020. https://doi.org/10.1090/mcom/3553.Search in Google Scholar

[50] M. Brunner, M. Innerberger, A. Miraçi, D. Praetorius, J. Streitberger, and P. Heid, “Corrigendum to: adaptive FEM with quasi-optimal overall cost for nonsymmetric linear elliptic PDEs,” IMA J. Numer. Anal., vol. 44, no. 3, pp. 1903–1909, 2024. https://doi.org/10.1093/imanum/drad103.Search in Google Scholar

[51] M. Innerberger and D. Praetorius, “MooAFEM: an object oriented matlab code for higher-order adaptive FEM for (nonlinear) elliptic PDEs,” Appl. Math. Comput., vol. 442, Art. no. 127731, 2023.10.1016/j.amc.2022.127731Search in Google Scholar

Received: 2023-12-01

Accepted: 2024-07-06

Published Online: 2024-11-04

Published in Print: 2025-06-26

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

Articles in the same Issue

https://doi.org/10.1515/jnma-2023-0150

Keywords for this article

goal-oriented adaptive finite element method; linear quantity of interest; iterative solver; nonsymmetric partial differential equations; optimal convergence rates; optimal complexity

Creative Commons

BY 4.0