Superconvergent DPG Methods for Second-Order Elliptic Problems

1.1 Summary of Results

We seek approximations uh∈𝒫p⁢(𝒯), 𝝈h∈𝒫p⁢(𝒯)d of the field variables u, 𝝈, where 𝒯 is a mesh of simplices and 𝒫p⁢(𝒯) denotes the space of 𝒯-piecewise polynomials of degree less than or equal to p∈ℕ0, and approximations u^h, σ^h of the traces u^, σ^ in spaces that will be defined later on. For sufficient regular solutions basic a priori analysis arguments give the estimate

where ∥⋅∥ denotes the L2⁢(Ω) norm and ∥⋅∥𝒮 is some appropriate norm for the traces. This estimate is optimal, since we seek approximations of u and 𝝈 in polynomial spaces of the same order and their errors are measured in L2⁢(Ω) norms. Nevertheless, it is unsatisfactory to some extent. Consider 𝑪 the identity, 𝜷=0, 𝒇=0 in (1.1). Then 𝝈=∇⁡u and we seek approximations of u and its gradient 𝝈 in polynomial spaces of the same order, which seems to be suboptimal. Fortunately, there exist at least two possibilities to achieve higher convergence rates under some assumptions on the regularity of solutions of (1.1) and its adjoint problem:

1. Augmenting the trial space: Instead of seeking approximations uh∈𝒫p⁢(𝒯) we seek approximations uh+∈𝒫p+1⁢(𝒯) and show that

   ∥u-uh+∥=𝒪⁢(hp+2).

2. Postprocessing: We use a common local postprocessing technique (see, e.g., the early works [11, 20]) to obtain an approximation u~h∈𝒫p+1⁢(𝒯) and prove that

   ∥u-u~h∥=𝒪⁢(hp+2).

Based on similar techniques we also provide a proof of the following:

1. DPG for ultra-weak formulations delivers the L2⁢(Ω) best approximation up to a higher-order term, i.e., for the approximation uh∈𝒫p⁢(𝒯) it holds that

   ∥u-uh∥≤∥u-Πp⁢u∥+𝒪⁢(hp+2),

   where Πp denotes the L2⁢(Ω) projection to 𝒫p⁢(𝒯).

The latter observation is quite interesting, because it shows that even though we do not aim for higher convergence rates (by increasing the polynomial degree in the trial space or by postprocessing) we get highly accurate approximations. We stress that this result has been observed in various numerical experiments, particularly also for more complex model problems like Stokes [17], but up to now a rigorous proof has not been given.

If 𝜷=0, we show that these results hold true when using different test norms (one of them is the so-called quasi-optimal test norm or graph norm). Surprisingly (at this point), for 𝜷≠0 the results are only valid if the quasi-optimal test norm is used, although all test norms under consideration are equivalent. This is also observed in our numerical studies.

1.2 Basic Ideas

For the proofs of the main results, we develop duality arguments and show approximation results (Lemma 8 and Lemma 10). To get the essential idea, consider the abstract formulation: Find 𝒖∈U such that

where U denotes the trial space and V the test space. With the trial-to-test operator Θ:U→V,

the ideal DPG method reads: Find 𝒖h∈Uh⊂U such that

Then we solve a dual problem: For some given g∈L2⁢(Ω), we determine 𝒗∈V and 𝒘=Θ-1⁢𝒗∈U, both unique, and employ Galerkin orthogonality to obtain

for arbitrary 𝒘h∈Uh.

For the case where we want to show that the approximation uh∈𝒫p⁢(𝒯) is nearly the L2⁢(Ω) best approximation, we have g=Πp⁢(u-uh)=Πp⁢u-uh. Therefore,

The latter estimate is what we have to show. Suppose that it holds. With the estimate for ∥𝒖-𝒖h∥U from above, it is straightforward to see that

Let us come back to the essential estimate

It holds if we would know that the higher derivatives of 𝒘 exist (in some sense) and can be bounded by the norm of g, so that, formally,

by some standard arguments. In our case we have that 𝒗∈H01⁢(Ω)×𝑯⁢(div;Ω)⊂V is the solution to the adjoint problem of (1.1) and under some assumptions has the higher regularity 𝒗∈H2⁢(Ω)×𝑯1⁢(𝒯)∩𝑯⁢(div;Ω), where 𝑯1⁢(𝒯) denotes 𝒯-piecewise Sobolev functions. Recall that 𝒘=Θ-1⁢𝒗. One difficulty is that the inverse of the trial-to-test operator does not map regular functions back to regular functions. However, it turns out (Lemma 8) that 𝒘 can be written as

where components of ~⁢𝒘∈U are related to the dual solution 𝒗, which is sufficiently regular and 𝒘⋆ is the solution of the (primal) problem (1.1) with data f and 𝒇 depending on the dual solution 𝒗 so that 𝒘⋆ has sufficient regularity as well. Let us point out that this idea used in the proofs is new and allows to treat different test norms. In [10], which deals with a simple reaction-diffusion problem and one specific test norm only, the representation of 𝒘 is obtained by integration by parts using the dual solution 𝒗 and it is not clear if that approach can be generalized to the present setting. Here, in the general case we have to consider the regularity of the dual solution 𝒗 and the regularity of the solution 𝒘⋆ of the primal problem. For the proofs it is also necessary that g is a function in the finite element space, so that we can choose 𝒘h=(g,0,0,0)+¯⁢𝒘h, where ¯⁢𝒘h is the best approximation of ~⁢𝒘+𝒘⋆. Then we show that the above estimates hold true.

Let us note that Θ is defined through the inner product in the test space. Thus, the representation of 𝒘=Θ-1⁢𝒗 from above strongly depends on the choice of the test norm and has to be analyzed for each norm individually (this is done in Lemma 8). Earlier Keith, Vaziri Astaneh and Demkowicz [15] considered the optimal test norm. In our notation this would yield 𝒘=(g,0,0,0), i.e., ~⁢𝒘=0=𝒘⋆.

Moreover, the ideas so far dealt with the ideal DPG method. In this paper we work out all results for the practical DPG method under standard assumptions, i.e., the existence of Fortin operators. This implies that we have to deal with additional discretization errors.

Finally, we note that higher convergence rates for the dual variable 𝝈 can not be obtained with the same arguments since the components of 𝒗=(v,𝝉) with

are less regular than in the case described above.

1.3 Outline

The remainder of the paper is organized as follows: Section 2 introduces basic notations, states the assumptions, and presents the main results (Theorem 3–5). The proofs of these theorems are postponed to Section 3, which also includes an a priori convergence estimate (Theorem 6) and the important auxiliary results Lemma 8, 10. In Section 4 we present two numerical experiments. The final Section 5 concludes this work with some remarks.

2.1 Notation

We make use of the notation ≲, i.e., A≲B means that there exists a constant C>0, which is independent of relevant quantities, such that A≤C⁢B. Moreover, A≃B means that both directions hold, i.e., A≲B and B≲A.

2.2 Mesh

Let 𝒯 denote a regular mesh of Ω consisting of simplices T and let 𝒮:={∂⁡T:T∈𝒯} denote the skeleton. We suppose that 𝒯 is shape-regular, i.e., there exists a constant κ𝒯>0 such that

where |T| denotes the volume measure of T∈𝒯. As usual, h:=h𝒯:=maxT∈𝒯⁡diam⁢(T) denotes the mesh-size.

2.3 Ultra-Weak Formulation

Before we derive the ultra-weak formulation of (1.1) in this subsection, we introduce some notation. Let T∈𝒯. We denote by (⋅,⋅)T the L2⁢(T) scalar product and with ∥⋅∥T the induced norm. On boundaries ∂⁡T, the L2⁢(∂⁡T) scalar product is denoted by 〈⋅,⋅〉∂⁡T and extended to the duality between the spaces H1/2⁢(∂⁡T) and H-1/2⁢(∂⁡T). Furthermore, we define the piecewise trace operators

where 𝒏T denotes the normal on ∂⁡T pointing from T to its complement. With these operators we define the trace spaces

These Hilbert spaces are equipped with minimum energy extension norms

We use the broken test spaces

and define the piecewise differential operators ∇𝒯:H1⁢(𝒯)→𝑳2⁢(Ω), div𝒯:𝑯⁢(div;𝒯)→L2⁢(Ω) on each T∈𝒯 by

Moreover, we define the following dualities for all u^∈H01/2⁢(𝒮), 𝝉∈𝑯⁢(div;𝒯), σ^∈H-1/2⁢(𝒮), v∈H1⁢(𝒯):

These dualities measure the jumps of 𝒗=(v,𝝉)∈H1⁢(𝒯)×𝑯⁢(div;𝒯), i.e.,

(2.1)

see, e.g., [3, Theorem 2.3].

The ultra-weak formulation is then derived from (1.1) by testing (1.1a) with 𝝉∈𝑯⁢(div;𝒯), (1.1b) with v∈H1⁢(𝒯), and piecewise integration by parts, i.e.,

Here, (⋅,⋅):=(⋅,⋅)Ω is the L2⁢(Ω) scalar product with norm ∥⋅∥. Set

and define F:V→ℝ and b:U×V→ℝ by

for all 𝒖=(u,𝝈,u^,σ^)∈U, 𝒗=(v,𝝉)∈V. The ultra-weak formulation then reads: Find 𝒖∈U such that

2.4 DPG Method and Approximation

In U we use the canonical norm,

For the test space V we define the three different norms

(2.3)

for 𝒗=(v,𝝉)∈V and denote by (⋅,⋅)V,⋆ the corresponding scalar products. Note that all norms in (2.3) are equivalent with equivalence constants depending on the coefficients 𝑪, 𝜷, γ. However, our main results hold for the quasi-optimal test norm ∥⋅∥V,qopt under mild assumptions on the coefficient 𝜷, whereas they hold for ∥⋅∥V,1, ∥⋅∥V,2 only if 𝜷=0, i.e., for symmetric problems.

We stress that b:U×V→ℝ is a bounded bilinear form and satisfies the inf–sup conditions with mesh independent constant. This can be proved with the theory developed in [3]. For our model problem we explicitly refer to [3, Example 3.7] for the details. There it is assumed that div⁡(𝑪-1⁢𝜷)=0 and γ≥0. We note that their analysis can also be done with our more general assumption (1.2).

The DPG method seeks an approximation 𝒖h∈Uh⊂U of the solution 𝒖∈U using the optimal test space Θ⁢(Uh), where Θ:U→V is defined by

Then 𝒖h∈Uh is the solution of

An essential feature of DPG is that inf–sup stability directly transfers to the discrete problem. However, in practice we replace Θ by a discrete version Θh:Uh→Vh⊂V defined by

Then the practical DPG method reads: Find 𝒖h∈Uh such that

In this work we deal with the piecewise polynomial trial spaces

and the piecewise polynomial test spaces

Here, we set

where

is the space of Raviart–Thomas functions (here ~⁢𝒫p⁢(T) denotes the space of homogeneous polynomials of degree p).

We also use the space C1⁢(𝒯):={v∈L∞⁢(Ω):v|T∈C1⁢(T¯),T∈𝒯}.

2.5 Fortin Operators

It is well known, see, e.g., [12], that (2.5) satisfies inf–sup conditions (and therefore admits a unique solution) if there exists a Fortin operator ΠF:V→Vh such that

Throughout, we suppose that a Fortin operator exists for the discrete polynomial trial and test spaces under consideration and that CF depends only on 𝑪, 𝜷, γ, p∈ℕ0, and the shape-regularity of 𝒯. Let us note that for general coefficients 𝑪, 𝜷, γ the existence of such operators is not known, except for some special cases, i.e., the Poisson model problem where 𝑪 is the identity and 𝜷=0=γ. Fortin operators for the latter problem on simplicial meshes have been constructed and analyzed in [3, 12]. We refer also to [16] for the construction and analysis of Fortin operators for second-order problems.

Supposing the existence of an Fortin operator, i.e., (2.6), we have:

2.6 Adjoint Problem and Regularity Assumptions

We define the adjoint problem (in the sense of L2⁢(Ω)-adjoints) of (1.1) as

(2.7)

Supposing (1.2) this problem admits a unique solution (v,𝝉)∈H01⁢(Ω)×𝑯⁢(div;Ω) for g∈L2⁢(Ω), 𝒈∈𝑳2⁢(Ω).

For our results we make use of the following assumptions:

Here, ∥⋅∥Hs⁢(Ω) is the usual notation for norms in the Sobolev space Hs⁢(Ω) (s>0), and ∥⋅∥ is the L2⁢(Ω) norm and ∥⋅∥𝑯s⁢(𝒯) the broken Sobolev norm for vector valued functions. We note that the constant C>0 strongly depends on the material data, i.e., the coefficients 𝑪, 𝜷, γ. For example, a small constant diffusion implies a blow-up of C.

2.7 Assumptions on Coefficients and Test Norms

Besides the assumptions on the coefficients and the domain to ensure unique solvability of problems (1.1) and (2.7) and estimates (2.8), we also need some additional assumptions on the coefficients that are listed in Table 1. We emphasize that 𝜷=0 in Cases (b) and (c) is also necessary in general. In particular, in Section 4 we provide a simple example where 𝜷≠0 and the choice ∥𝒗∥V=∥𝒗∥V,1 or ∥𝒗∥V=∥𝒗∥V,2 does not lead to higher convergence rates, whereas ∥𝒗∥V=∥𝒗∥V,qopt does.

2.8 L2⁢(Ω) Projection

Our first main result shows that the DPG method with ultra-weak formulation delivers up to a higher-order term the L2⁢(Ω) best approximation for the scalar field variable. To that end let Πp:L2⁢(Ω)→𝒫p⁢(𝒯) denote the L2⁢(Ω) projector.

2.9 Higher Convergence Rate by Increasing Polynomial Degree

Our second main result shows that higher convergence rates for the scalar field variable are obtained by increasing the polynomial degree in the approximation space.

2.10 Higher Convergence Rate by Postprocessing

Our third and final main result shows that higher convergence rates for the scalar field variable are obtained by postprocessing the solution: Let 𝒖h=(uh,𝝈h,u^h,σ^h)∈Uh:=Uh⁢p be the solution of (2.5). We define u~h∈𝒫p+1⁢(𝒯) on each element T∈𝒯 as the solution of the local Neumann problem

(2.9)

Let us note that this type of postprocessing is common in literature and can already be found in the early works [11, 20].

3.1 Projection Operators and Approximation Results

Throughout let p∈ℕ0. Let Πp:L2⁢(Ω)→𝒫p⁢(𝒯) denote the L2⁢(Ω) projector. For 𝝉∈𝑳2⁢(Ω) the term Πp⁢𝝉 is understood as the application of Πp to each component. We have the (local) approximation properties

(3.1)

where |⋅|Hn⁢(𝒯):=∥D𝒯n⋅∥ with D𝒯n denoting the 𝒯-elementwise n-th derivative operator. Let

denote the Scott–Zhang projection operator [19] or any other operator with the property