An adjoint-based super-convergent Galerkin
approximation of functionals
A DISSERTATION
SUBMITTED TO THE FACULTY OF THE GRADUATE SCHOOL
OF THE UNIVERSITY OF MINNESOTA
BY
Shiqiang Xia
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF
Doctor of Philosophy
Advisor: Professor Bernardo Cockburn
May, 2022
© Shiqiang Xia 2022
ALL RIGHTS RESERVED
Acknowledgements
First and foremost, I would like to express my sincere gratitude to my advisor Prof.
Bernardo Cockburn not only for his continuous support of my Ph.D study and related
research, but also for his patience, motivation, enthusiasm and humor. His guidance
helped me a lot in my research and in the writing of this dissertation. The life wisdom
he passed to me and his way of critical thinking are also invaluable to me. I could not
have imagined having a better advisor and mentor for my Ph.D. work.
Besides my advisor, I would like to thank the rest of my dissertation committee:
Prof. Mitchell Luskin, Prof. Li Wang, and Prof. Yousef Saad, for their insightful
comments and encouragement.
Special thanks go to my parents, my girlfriend Xueyan, and my friends for their
constant support during all these years.
i
Dedication
To my parents
ii
Abstract
In this dissertation, we study an adjoint-based method for the computation of highly
accurate approximations of functionals of solutions of second-order elliptic problems.
First, we present a rigorous a priori error analysis of the method applied to linear
functionals. We also carry out comprehensive numerical experiments to display the
super-convergence properties of the method and to demonstrate the advantage of the
method over the standard method. We then extend the method to an important non-
linear problem, namely, that of the computation of eigenvalues. To handle the possible
lack of smoothness of the solution or that of the functional, we further develop a new
mesh adaptative algorithm which we incorporate into the adjoint-based method.
iii
Contents
Acknowledgements i
Dedication ii
Abstract iii
List of Tables vii
List of Figures x
1 Introduction 1
1.1 Background and motivation . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 The organization of the dissertation . . . . . . . . . . . . . . . . . . . . 3
2 The adjoint-based method and its a priori error analysis for linear
functionals 4
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 The adjoint-based method . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 The components of the method . . . . . . . . . . . . . . . . . . . 7
2.2.2 The definition of the adjoint-based method . . . . . . . . . . . . 13
2.3 The A Priori Error Estimates . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Proof of the main theorem . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.5.1 Illustration of the approximation error u− u∗h . . . . . . . . . . . 28
2.5.2 How close Ω0 can be to the boundary ∂Ω . . . . . . . . . . . . . 29
iv
2.5.3 Effect of the smoothness of the solution for variable Ω0 . . . . . 32
2.5.4 Effect of the smoothness of the functional . . . . . . . . . . . . . 35
2.6 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3 The adjoint-based method for the eigenvalue problem 41
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 The adjoint-based method for eigenvalue problem . . . . . . . . . . . . . 43
3.2.1 The components of the method . . . . . . . . . . . . . . . . . . . 44
3.2.2 The definition of the method . . . . . . . . . . . . . . . . . . . . 50
3.3 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.3.1 Proof of Theorem 3.1 . . . . . . . . . . . . . . . . . . . . . . . . 54
3.3.2 Proof of Proposition 3.1.1 . . . . . . . . . . . . . . . . . . . . . . 56
3.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.4.1 The one dimensional Laplace eigenproblem . . . . . . . . . . . . 58
3.4.2 The two-dimensional Laplace eigenproblem . . . . . . . . . . . . 60
3.4.3 The efficiency of the adjoint-based method . . . . . . . . . . . . 67
3.4.4 Application to the quantum harmonic oscillator . . . . . . . . . . 68
3.4.5 The relevance of the adjoint-correction term . . . . . . . . . . . . 70
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4 An adjoint-based adaptive error approximation of functionals by the
HDG method for second-order elliptic equations 77
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.2 The adjoint-based error approximation of functionals . . . . . . . . . . . 79
4.2.1 The adjoint-based method in an abstract form . . . . . . . . . . 80
4.2.2 The HDG method . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.2.3 Functional error approximation . . . . . . . . . . . . . . . . . . . 84
4.3 Examples of approximating functionals . . . . . . . . . . . . . . . . . . . 86
4.3.1 Volume integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.3.2 Normal flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.3.3 Eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.4 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
4.4.1 Uniform Refinement . . . . . . . . . . . . . . . . . . . . . . . . . 90
v
4.4.2 Adaptive Refinement . . . . . . . . . . . . . . . . . . . . . . . . . 97
4.5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5 An adjoint-based adaptive method with convolution filtering 107
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.2 The method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.3.1 Volume integral with corner singularity . . . . . . . . . . . . . . 110
5.3.2 Boundary integral with discontinuity . . . . . . . . . . . . . . . . 112
5.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6 Conclusion and Discussion 115
References 116
Appendix A. Supplements 124
A.1 Special projections and their properties . . . . . . . . . . . . . . . . . . 124
A.1.1 The M -decompositions . . . . . . . . . . . . . . . . . . . . . . . 124
A.1.2 The HDG Projection . . . . . . . . . . . . . . . . . . . . . . . . 125
A.1.3 The HDG projection for curved elements . . . . . . . . . . . . . 126
A.1.4 Estimate of the jump εuh − εûh . . . . . . . . . . . . . . . . . . . . 132
A.2 Auxiliary Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
vi
List of Tables
2.1 Local spaces V (K)×W (K) admitting an M(∂K)-decomposition. . . . 10
2.2 Coefficients of the convolution kernel K2k for k = 1, . . . , 5. Since the
kernel is symmetric, the coefficients are such that C−r = Cr. Polynomials
of degree 2k are reproduced by convolution with the kernel. . . . . . . . 13
2.3 History of convergence when the distance between Ω0 and ∂Ω is fixed
and equal to 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.4 History of convergence when the distance between Ω0 and ∂Ω is, for each
mesh, is 2k elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5 CPU time, error and efficiency comparison of the adjoint-based method
JFh (u
∗
h, v
∗
h) and the standard method J(uh) when u is a smooth function. 32
2.6 History of convergence when the solution u is smooth and k = 1, 2, 3.
The distance between Ω0 and ∂Ω is, for each mesh, 2k elements. . . . . 33
2.7 History of convergence when the solution u has a corner singularity and
k = 1, 2, 3. The distance between Ω0 and ∂Ω is, for each mesh, 2k elements. 34
2.8 CPU time, error and efficiency comparison of the adjoint-based method
JFh (u
∗
h, v
∗
h) and the standard method J(uh) when u has a corner singu-
larity. The set Ω0 varies with k and h. . . . . . . . . . . . . . . . . . . . 34
2.9 History of convergence for the boundary integral functional with v = π−θπ . 36
2.10 History of convergence of approximations to u = sin(πx) sin(πy) and to
v = π−θπ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.11 CPU time, error and efficiency comparison of the adjoint-based approxi-
mation JFh (u
∗
h, v
∗
h) and the standard approximation J(uh) for the bound-
ary integral when u = sin(πx) sin(πy) and to v = π−θπ . The set Ω0 varies
with k and h. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
vii
2.12 History of convergence of the boundary integral functional with u =
r
2
3 sin(23θ) and not smooth v. . . . . . . . . . . . . . . . . . . . . . . . . 38
2.13 CPU time, error and efficiency comparison of the adjoint-based approxi-
mation JFh (u
∗
h, v
∗
h) and the standard approximation J(uh) for the bound-
ary integral when u = r
2
3 sin(23θ) and is not smooth v. The set Ω0 varies
with k and h. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1 History of convergence of the first five eigenvalues provided by the HDG
method, λh(uh), (top) and by the adjoint-based method, λ
AC
h (u
∗
h), (bot-
tom) for the problem on (0, 1). . . . . . . . . . . . . . . . . . . . . . . . 59
3.2 History of convergence of the first five eigenvalues provided by the HDG
method, λh(uh), (top) and by the adjoint-based method, λ
AC
h (u
∗
h), (bot-
tom) for the problem on the unit square domain. . . . . . . . . . . . . . 61
3.3 History of convergence of the approximate (uh) and accuracy-enhanced
(u∗h) eigenfunctions on the unit square domain. . . . . . . . . . . . . . . 62
3.4 History of convergence for the L-shaped domain problem. . . . . . . . . 65
3.5 History of convergence for the harmonic oscillator. . . . . . . . . . . . . 69
3.6 History of convergence of λh(u
∗
h) on the unit square domain. . . . . . . . 71
3.7 The adjoint-correction term |ACh(u∗h)|/∥u∗h∥2L2(Ω) and the ratio ∥u∗h∥2L2(Ω)|λ−
λh(u
∗
h)|/ACh(u∗h) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
3.8 The adjoint-correction term |ACh(uh)|/∥uh∥2L2(Ω) and the ratio ∥uh∥2L2(Ω)|λ−
λh(uh)|/|ACh(uh)| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.9 History of convergence of λACh (uh) on the unit square domain. . . . . . . 73
3.10 The adjoint-correction term |ACh(uh)|/∥uh∥2L2(Ω) and the ratio ∥uh∥2L2(Ω)|λ−
λh(uh)|/ACh(uh) in an L-shaped domain. . . . . . . . . . . . . . . . . . 73
3.11 History of convergence of λACh (uh) on an L-shaped domain. . . . . . . . 73
4.1 History of convergence of approximating functional J(u) = (u, g)Ω uni-
formly on the unit square Ω by the HDG method. Here Eh is the approx-
imation of the error and Eh,|·| is the localized error bound. . . . . . . . . 91
4.2 History of convergence of approximating functional J(u) =< ∇u·n, ψ >∂Ω
on the unit square Ω by the HDG method. Here Eh is the approximation
of the error and Eh,|·| is the localized error bound. . . . . . . . . . . . . 93
viii
4.3 History of convergence of approximating the 1st eigenvalue λ1,1 of (4.3.3)
uniformly by the HDG method. Here Eh is the approximation of the error
and Eh,|·| is the localized error bound. . . . . . . . . . . . . . . . . . . . 95
5.1 Functional error comparison of the HDG method and adjoint-based adap-
tive method for the first functional. . . . . . . . . . . . . . . . . . . . . 111
5.2 The CPU time and DOFs comparison of the two methods for the first
functional. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.3 Functional error comparison of the HDG method and adjoint-based adap-
tive method for the second functional. . . . . . . . . . . . . . . . . . . . 113
5.4 The CPU time and DOFs comparison of the two methods for the second
functional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
ix
List of Figures
2.1 Illustration of star-shapedness (left, middle) and shape-regularity(right) 14
2.2 Approximation error u − u∗h when k = 1 and N = 40. The distance
between Ω0 and ∂Ω is 0.1. Note that the magnitude of the error in Ω\Ω0
is essentially the same as that in Ω0. . . . . . . . . . . . . . . . . . . . . 28
2.3 Slices of the approximation error u − u∗h at x = 0.50117 (left) and y =
0.20117 (right) when k = 1 and = 40. In each element, we evaluate u−u∗h
at 5 Gauss quadrature points. . . . . . . . . . . . . . . . . . . . . . . . 29
2.4 Log scale error comparison of the adjoint-based method JFh (u
∗
h, v
∗
h) (◦)
and the standard method J(uh) (□) when u is a smooth function and
k = 1. Note that to achieve the given error tolerance 9 × 10−5, the
adjoint-based method only takes 0.9s with a mesh N = 20, while the
standard method takes 49s with a mesh N = 90. . . . . . . . . . . . . . 32
3.1 An illustration of the definition of u∗h in the whole domain Ω. Here Ω0 is
the subdomain that is two elements away from the boundary ∂Ω. This
is how we have to define Ω0 to be able to carry out the postprocessing
with k = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.2 The errors of the first three eigenvalues in logarithmic scale for k = 1 with
N = 27, 28, 29, 210 elements on (0, 1). Here Ei = log10(|λi − λi,h(uh)|)
(dashed line) and E∗i = log10(|λi − λACi,h (u∗h)|) (solid line) for i = 1, 2, 3. . 60
3.3 Approximations to the first eigenfunction on a unit square domain with
k = 1 and h = 0.1: The HDG approximation uh (left column) and the
accuracy-enhanced approximation u∗h (right column.). The errors are
plotted at the bottom row. . . . . . . . . . . . . . . . . . . . . . . . . . 64
x
3.4 Approximations of the first eigenfunction in the L-shaped domain with
k = 1 and h = 0.1: The HDG approximation uh (left column) and
the accuracy-enhanced approximation u∗h (right column.). The errors
are plotted at the bottom row. They are concentrated at the corner
singularity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.5 The number of eigenvalues that have a relative error |λ−λh|/λ less than
the tolerance ϵ = 1%. The HDG method (♢) and the adjoint-based
method (◦) when k = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
3.6 The HDG approximation (left) and the accuracy-enhanced approxima-
tion (right) of the fourth eigenfunction of the harmonic oscillator equation. 70
3.7 The element-wise plot of the adjoint-correction term |ACh(uh)|K = ⟨uh−
ûh, (q̂h − qh) · n⟩∂K on an L-shaped domain when k = 1 and h = 0.1. It
captures well the location of the singularity of the eigenfucntion. . . . . 74
4.1 Plots of the error and the error approximation against the DOF in loga-
rithmic scale (left column) and plots of the effectivity index Ieff and I
loc
eff
(right column) when the functional J(u) = (u, g)Ω is approximated uni-
formly with k = 1 (top row), k = 2 (middle row), k = 3 (bottom row).
Here the solutions are smooth. . . . . . . . . . . . . . . . . . . . . . . . 92
4.2 Plots of the error and the error approximation against the DOF in log-
arithmic scale (left column) and plots of the effectivity index Ieff and
I loceff (right column) when the functional J(u) =< ∇u · n, ψ >∂Ω is ap-
proximated uniformly with k = 1 (top row), k = 2 (middle row), k = 3
(bottom row). Here the solutions are smooth. . . . . . . . . . . . . . . . 94
4.3 Plots of the error and the error approximation against the DOF in loga-
rithmic scale (left column) and plots of the effectivity index Ieff and I
loc
eff
(right column) when the 1st eigenvalue λ1,1 of (4.3.3) is approximated
uniformly with k = 1 (top row), k = 2 (middle row), k = 3 (bottom row).
Here the eigenfunctions are smooth. . . . . . . . . . . . . . . . . . . . . 96
4.4 Plots of the error and the error approximation against the DOF in loga-
rithmic scale for k = 1, 2, 3 (left column) and plots of the corresponding
effectivity index Ieff (right column) when approximating J1(u) = (u, g)Ω
and the solutions are smooth. . . . . . . . . . . . . . . . . . . . . . . . 98
xi
4.5 Example of adaptive meshes for the first functional when u is singular.
Left: k =2, 254 elements, the number of DOF is 5649 and error is 2.27e-
10; Right: k=3, 112 elements, the number of DOF is 3964 and error is
2.49e-10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.6 Plots of the error and the error approximation against the DOF in loga-
rithmic scale for k = 1, 2, 3 (left column) and plots of the corresponding
effectivity index Ieff (right column) when approximating J1(u) = (u, g)Ω
and u is singular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.7 Example of adaptive meshes for the second functional. Left: k =2, 2148
elements, the number of DOF is 48144, and error is 1.12e-10; Right: k=3,
310 elements, the number of DOF is 11076, and error is 1.97e-10. . . . . 101
4.8 Plots of the error and the error approximation against the DOF in log-
arithmic scale for k = 1, 2, 3 (left column) and plots of the correspond-
ing effectivity index Ieff(right column) when approximating J2(u) =<
∇u · n, ψ >∂Ω. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.9 Examples of adaptive meshes for the third functional. Top left: k =1,
2905 elements, the number of DOF is 34719, and error is 8.01e-04; Top
right: k=2, 244 elements, the number of DOF is 5427, and error is 9.01e-
04; Bottom: k=3, 170 elements, the number of DOF is 6044 and error is
1.00e-3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.10 Plots of the error and the error approximation against the DOF in loga-
rithmic scale for k = 1, 2, 3 (left column) and plots of the corresponding
effectivity index Ieff (right column) when approximating the 1st eigen-
value on the L-shaped domain. . . . . . . . . . . . . . . . . . . . . . . . 105
5.1 Example of adaptive meshes when u is singular with the initial mesh size
h0 = 0.05. Left: k = 1; Right: k = 2. . . . . . . . . . . . . . . . . . . . 111
5.2 Example of adaptive meshes when ψ is discontinuous with the initial
mesh size h0 = 0.05. Left: k = 1; Right: k = 2. . . . . . . . . . . . . . . 113
xii
Chapter 1
Introduction
In this chapter, we introduce the background and motivation of the adjoint-based
method we study and then give an overview of the dissertation.
1.1 Background and motivation
In many scientific and engineering applications, one is often interested in the ap-
proximations of certain quantities, such as the normal force on a surface, the mean
value of a solution, and the heat flux across certain boundary. These quantities can
be treated mathematically as functionals J(u), where u is typically one or several field
variables which are governed by partial differential equations (PDEs). For example, in
fluid mechanics, the lift and drag forces of an object in a viscous incompressible fluid
can be the quantities of interest. The fluid flow is governed by the Navier-Stokes equa-
tions and the lift and drag forces are surface integrals related to the stress tensor. Since
these quantities are often used to measure the efficiency and performance of engineering
designs, accurate approximation of functionals J(u) is of significant importance.
To approximate a functional J(u), the standard way consists in obtaining a numerical
approximation uh of the exact solution u and then using J(uh) as an approximation
of the functional. In 2000 [62], Pierce and Giles propose an adjoint-based recovery
method for approximating a functional J(u). Rather than simply using J(uh) as an
approximation, this method defines a new approximation Jh(uh) := J(uh) + ACh by
adding a carefully devised computable term, ACh. This term ACh is obtained by solving
1
2a dual (adjoint) PDE, hence the name adjoint-correction term. The new approximation
Jh(uh) gives a more accurate value for the functional and the method is general enough
to be applicable to various types of functionals, and not only to a wide variety of partial
differential equations but also to a wide class of numerical methods. Giles and Su¨li
further study this method in [41]. At around the same time, based on the pioneering
work of Eriksson, Estep, Hansbo and Johnson in [35], Becker and Rannacher develop the
dual-weighted-residual (DWR) method for a posteriori error estimate of J(u) − J(uh)
and mesh optimization within the context of finite element methods. See the thorough
review [5] and the reference therein. Since then, a considerable number of studies on
adjoint-based method are carried out. For example, Venditti and Darmofal apply the
adjoint-based error correction technique to viscous flows and inviscid flows in [73, 72];
Hartmann and Houston study the fluid dynamic applications focusing on discontinuous
Galerkin methods in [46]. See also [44, 45, 43, 52, 68, 34]. For general discretization,
Fidkowski and Darmofal provide a comprehensive literature review in [37].
In 2017, Cockburn and Wang introduce a novel adjoint-based method for approxi-
mating linear functionals in [27]. Extensive numerical experiments in [27] indicate that
this new method achieves a convergence rate of 4k when approximating smooth function-
als, where k is the polynomial degree used for the hybridizable discontinuous Galerkin
(HDG) method. In other words, the new method converges essentially with twice the
order of convergence of the classic methods while employing only twice the effort of com-
puting J(uh)! This is very attractive since for high dimensional problems, it typically
takes much more computational effort to double the convergence rate. The key idea
of this new method is to combine the adjoint-correction technique in [62] with another
accuracy-enhancing technique, called the convolution filtering technique. This technique
is first proposed by Bramble and Schatz for finite element methods for second-order el-
liptic equations in 1977 [9]. It is later extended to discontinuous Galerkin methods for
hyperbolic equations in [25]. More details about this technique will presented in the
following chapter.
In this dissertation, we continue the study of this new adjoint-based method pro-
posed by Cockburn and Wang in [27]. We rigorously prove the convergence property of
this method under suitable regularity assumptions and we extend this method to the
eigenvalue problem. To handle the lack of smoothness of the solution, we further develop
3an adjoint-based error approximation and mesh adaptation algorithm. The dissertation
is organized in four main chapters as explained below.
1.2 The organization of the dissertation
First, in Chapter 2, we introduce the components of the adjoint-based method.
To put the numerical experiments in [27] in firm mathematical ground, we present a
rigorous a priori error analysis of the adjoint-based method. We prove that under
suitable regularity assumptions, the error obtained by this new adjoint-based method
is guaranteed to converge with a rate of 4k when approximating linear functionals. The
material in this chapter is mainly from our publication [28].
Second, in Chapter 3, we extend the adjoint-based method to the eigenvalue prob-
lem. We treat the eigenvalue as a non-linear functional of its corresponding eigenfunction
and illustrate the method on a second-order elliptic eigenvalue problem. Our extensive
numerical results show that the approximate eigenvalues computed by our method con-
verge with a rate of 4k + 2 when tensor-product polynomials of degree k are used for
the HDG method. The material in this chapter is mainly from our publication [30].
Next, in Chapter 4, we present a novel, fully computable error approximation for
functionals. When the functional is not smooth enough, we utilize this error approxima-
tion to drive an adaptive algorithm. Unlike most adjoint-based error estimations in the
literature, the novelty of our method is that the error approximation is obtained without
requiring an auxiliary finer mesh or higher order approximation spaces for solving the
adjoint problem. This reduces the computational cost and eases the implementation of
the problem. The material in this chapter is mainly from our publication [29].
Then, in Chapter 5, we present an adaptive method that takes advantage of the
convolution filtering technique used in the adjoint-based method in Chapter 2 and 3, as
well as the error approximation technique in Chapter 4. We show that the new method
can achieve a significantly better approximation even when the exact solution of the
problem has corner singularities or when the given data is discontinuous.
Finally, in Chapter 6, we end with some concluding remarks and a sketch of future
work.
Chapter 2
The adjoint-based method and its
a priori error analysis for linear
functionals
2.1 Introduction
In this chapter, we introduce the adjoint-based Galerkin approximation method
initially proposed in [27] and we present the a priori error analysis of the method.
To approximate a functional J(u), the standard way is to obtain a numerical ap-
proximation uh of the exact solution u and then use J(uh) as an approximation of the
functional. The adjoint-based method we study here improves the accuracy of these
approximations in two ways. First, by increasing the accuracy of the numerical solution
uh by means of the convolution filtering technique of [9]. Second, by modifying the
formula to approximate the functional by using the adoint-correction approach of [62].
Let us briefly describe these two compnents of the method.
The first component is a convolution filtering technique which improves the accuracy
of the Galerkin approximation. This local post-processing technique was first proposed
in [9] for finite element methods for second-order elliptic equations. It takes advantage of
the well-known fact that the Galerkin solution must oscillate around the exact solution
in a certain pattern because of the Galerkin orthogonality property. Hence, convolving
4
5the Galerkin solution with a specific B-spine kernel filters out these oscillations and
provides a more accurate solution. [9] shows how this takes place in a subdomain
Ω0 ⊂⊂ Ω included in a set where the mesh is translation invariant. This approach was
later extended to discontinuous Galerkin methods for hyperbolic equations in [25].
Although this convolution filter has been proven to be effective only on a strict sub-
domain Ω0 of the original domain Ω, or with periodic boundary conditions, one can still
apply the convolution in the whole domain Ω and even when the meshes are not locally
translation invariant. Thus, Kirby, Ryan and their collaborators have further extended
this technique to various types of meshes including unstructured triangular meshes, see
[31, 55, 57, 56, 58, 50, 51]. In their work, this technique is referred to as smoothness-
increasing accuracy-conserving (SIAC) filtering technique. Also, to apply the filtering in
the whole domain, one-side kernel and position-dependent kernel approaches have been
studied in [64, 71, 49, 66]. The adjoint-based method considered here keeps the original
symmetric B-spline kernel employed by [9] to define a more accurate approximation in
Ω0. To extend it to the whole domain Ω, an auxiliary problem is solved on Ω \ Ω0, see
Section 2.2.2 for details.
The second component of the adjoint-based method for approximating functionals
more accurately is the adjoint-correction method of [62]; see the review in [41]. Roughly
speaking, this powerful technique consists in numerically solving the adjoint problem
for the functional J(·) so that, an extra, computable correction term can be added to
J(uh) which results in a much better approximation.
The idea of using the adjoint problem has been studied for decades. For example, the
dual-weighted residual method (DWR) uses the adjoint problems to devise a posteriori
error estimates and mesh adaptivity algorithms, see [5, 35] and the references therein.
Note that the DWR method and the adjoint-correction method we consider here are
different. To obtain the a posteriori error for the DWR method, the exact adjoint solu-
tion or a very high-order accurate approximation (obtained by a higher-order method
or a finer mesh) is needed. In contrast, the adjoint-correction method we consider here
does not require a more accurate approximation for the adjoint problem and is not used
within an adaptive algorithm driven by an a posteriori error estimate.
In [27], these two components are put together which resulted in the adjoint-based
6method we are analyzing here. Therein numerical experiments in one and two dimen-
sional spaces are carried out. The hybridizable discontinuous Galerkin (HDG) method
with polynomial degree k > 0 is used to obtain an approximations uh of the exact solu-
tion u. The results indicate that the approximation defined by this new adjoint-based
method converges to the functional J(u) with order h4k. Compared to the standard
approximation J(uh), which converges with order h
2k+1, this new method essentially
doubles the order of convergence by only doubling the computational effort for obtaining
J(uh).
In this chapter, we only study the convergence properties of the method when the
meshes are translation invariant in most of the interior of the domain Ω. The use
of the method with adaptive strategies will be considered in Chapter 4 and Chapter
5. Moreover, although the method we analyze here can be applied to quite general
functionals, as argued in [62, 41, 27], our results are for the following simple model
problem. The functional we are considering here is
J(u) :=
∫
Ω
g(x)u(x)dx, (2.1.1)
where u is the solution of a steady-state diffusion equation
−∆u = f in Ω, (2.1.2a)
u = uD on ∂Ω (2.1.2b)
on a bounded domain Ω with Lipschitz continuous boundary ∂Ω. For the error analysis,
we need to assume that the domain Ω is (k + 2)-regular. The precise definition of this
standard elliptic regularity property is given in Section 2.3.
The remainder of this chapter is organized as follows. In Section 2.2, we define the
adjoint-based method. We then state and briefly discuss our main result, the a priori
error estimates of Theorem 2.1 in Section 2.3. Its proof is detailed in Section 2.4. In
Section 2.5, we present numerical results designed to validate the theory and to explore
the convergence properties of the method in cases not covered by the theory. Finally,
in Section 2.6, we conclude and describe possible extensions of our results.
72.2 The adjoint-based method
This section is devoted to introduce the adjoint-based method.
2.2.1 The components of the method
We begin by defining the three components of the adjoint-based approximations,
namely, the Galerkin method (which we take it to be the HDG method), the adjoint-
correction method, and the convolution filtering technique. We would then be ready to
define the adjoint-based method.
The HDG Method
We use the HDG method here because it has a general structure that allows us
to extend the results to a large class of Galerkin methods, including the mixed finite
element method and the continuous Galerkin method.
To define the HDG method, we first partition the domain Ω into elements K forming
a conforming mesh Th. We set ∂Th := {∂K : K ∈ Th} and let Fh denote the set of faces
F of all the elements K ∈ Th. We also let F(K) denote the set of all faces F of the
element K ∈ Th.
We rewrite the model elliptic problem (2.1.2) as a system of first-order equations:
q +∇u = 0 in Ω, (2.2.1a)
∇ · q = f in Ω, (2.2.1b)
u = uD on Ω. (2.2.1c)
The HDG method seeks an approximation,(qh, uh, ûh), to the exact solution of the
model problem (2.2.1), (q|Ω, u|Ω, u|Fh), in the finite dimensional space V h ×Wh ×Mh,
where
8V h := {v ∈ L2(Th) :v|K ∈ V (K) ∀K ∈ Th},
Wh := {w ∈ L2(Th) :w|K ∈W(K) ∀K ∈ Th},
Mh := {µ ∈ L2(Fh) :µ|F ∈M(F ) ∀F ∈ Fh}.
If we use notation
(u, v)Th :=
∑
K∈Th
∫
K
uv dx and ⟨v , w⟩∂Th :=
∑
K∈Th
∫
∂K
vu ds,
then the HDG approximation is determined as the solution of the following weak for-
mulation:
(qh, r)Th − (uh,∇ · r)Th + ⟨ûh , r · n⟩∂Th = 0, (2.2.2a)
−(qh,∇w)Th + ⟨q̂h · n , w⟩∂Th = (f, w)Th , (2.2.2b)
⟨ûh , µ⟩∂Ω = ⟨uD , µ⟩∂Ω, (2.2.2c)
⟨q̂h · n , µ⟩∂Th\∂Ω = 0, (2.2.2d)
for all (r, w, µ) ∈ V h ×Wh ×Mh, where
q̂h · n := qh · n+ τ(uh − ûh) on ∂Th.
We obtain different methods by choosing different local spaces V (K), W (K) and
M(∂K) := {µ ∈ L2(∂K) : µ|F ∈M(F ) ∀F ∈ F(K)},
and stabilization functions τ . The hybridized version of mixed methods is obtained
by simply taking τ equal to zero, as pointed out in [21]. Since we consider the HDG
method in a fairly general setting, we make the following assumptions.
Assumption 1. When we say we use the HDG method with polynomial approximation
of degree k ≥ 0 for a collection of meshes {Th}h>0, we assume we choose proper local
spaces V (K)×W (K)×M(∂K) and τ such that the approximation errors in L2 norm
9satisfy
∥u− uh∥Th ≤ Chk+1, (2.2.3a)
∥q − qh∥Th ≤ Chk+1, (2.2.3b)
for smooth enough solution u and q, where C is a constant independent of h.
Assumption 2. For the HDG method with polynomial approximation of degree k, there
exists an element-wise auxiliary projection Πh(q, u) := (ΠV q, ΠWu) ∈ V (K)×W (K),
called the HDG projection, such that the so-called weak commutativity property
(∇ · q, w)K = (∇ ·ΠV q, w)K + ⟨τ(ΠWu− u) , w⟩∂K ∀w ∈W (K), (2.2.4a)
and the following property
(u,∇ · r)K = (ΠWu,∇ · r) ∀r ∈ V (K), (2.2.4b)
are satisfied and the HDG projection has the optimal approximation properties
∥u−ΠWu∥K ≤ Chk+1, (2.2.5a)
∥u−ΠWu∥∂K ≤ Chk+
1
2 , (2.2.5b)
∥q −ΠV q∥K ≤ Chk+1, (2.2.5c)
∥(q −ΠV q) · n∥∂K ≤ Chk+
1
2 , (2.2.5d)
for smooth enough solution u and q, where C is a constant independent of h.
These two assumptions do hold for local spaces that admit an M -decomposition,
see [19]. For the sake of completeness, we give in Appendix A.1.1 the definition of an
M -decomposition and in Appendix A.1.2 the definition of the HDG projection. For a
few simple shapes of elements, we give the local spaces that admit an M -decomposition
in Table 2.1.
10
Element V (K) W (K) M(∂K)
Triangle Pk(K) Pk(K) Pk(∂K)
Quadrilateral Pk(K)⊕ curl span{ξ4λk3 , ξ4λk4}
Square Qk(K)⊕ curl span{xk+1y, xyk+1} Qk(K) Pk(∂K)
Tetrahedron Pk(K) Pk(K) Pk(∂K)
Prisms1 Pk(K)⊕ curl span
{
zk+1(x∇y − y∇x),
zP˜k(x, y)(x∇y − y∇x)
}
Pk(K) Pk(∂K)
Cube Qk(K)⊕ curl span

xkyzk+1∇x,
xk+1z∇y,
xk+1ykz∇y
(1− x)x(1− z)zk∇y,
(1− x)x(1− y)yk∇z,
(1− x)x(1− y)ykzk∇z

Qk(K) Qk(∂K)
Table 2.1: Local spaces V (K)×W (K) admitting an M(∂K)-decomposition.
In Table 2.1, Pk(K) denotes the space of polynomials in element K with degree at
most k and Qk(K) denotes the space of tensor product polynomials of degree at most
k in each variable. Finally, P˜k(x, y) denotes the space of homogeneous polynomials
of degree k in x and y. The same notation is used for the three-dimensional case.
For a quadrilateral element K, let {vi}4i=1 be the set of vertices and let {ei}4i=1 be
the set of edges, where the edge ei connects vi and vi+1, where we set v5 = v1.
Then, for 1 ≤ i ≤ 4, we define λi to be the linear function that vanishes on edge ei
and reaches maximum value 1 in the closure of K, and let ξi be a rational function
such that ξi|ei ∈ P(ei) and ξi(vj) = δij , where δij is the Kronecker delta. More
details and examples about the construction of M -decomposition spaces in two and
three dimensions are provided in [19, 17, 18].
The adjoint-correction method
The adjoint-correction method is proposed in [62] for approximating functionals
J(u). Rather than simply using J(uh) as an approximation, this method obtains a
1The prism here is a unit prism, that is, the prism whose basis are unit triangles, and whose faces
on the xz− and yz−planes are unit squares. We also assume k ≥ 1 here. For other examples of prisms,
see [18].
11
new approximation Jh(uh) by adding a carefully devised computable term, ACh, called
the adjoint-correction term. Roughly speaking, this adjoint-correction term is zero or
significantly smaller than the error |J(u) − J(uh)| when the numerical solution uh is
obtained by a Galerkin method. However, for a solution u∗h that is not a Galerkin
solution, including this term results in a better approximation to J(u). Let us describe
the adjoint-correction method for the functional (2.1.1) . We follow [27].
Let (qh, uh, q̂h ·n, ûh) ∈ L2(Th)×H1(Th)×L2(∂Th)×L2(Fh) be any approximation
to (q, u, q ·n|∂Th , u|∂Th) in the model problem (2.2.1) such that ⟨ûh , µ⟩∂Ω = ⟨uD , µ⟩∂Ω
for any µ ∈ L2(Fh).
Similarly, consider the adjoint problem
p+∇v = 0 in Ω (2.2.6a)
∇ · p = g in Ω (2.2.6b)
v = 0 on ∂Ω (2.2.6c)
and let (ph, vh, p̂h ·n, v̂h) ∈ L2(Th)×H1(Th)×L2(∂Th)×L2(Fh) be any approximation
to (p, v,p · n|∂Th , v|∂Th) in the adjoint problem (2.2.6) such that v̂h = 0 on ∂Ω. Next
we define a new approximation of J(u) as
Jh(uh) = J(uh) +ACh, (2.2.7)
where the adjoint-correction term is as follows
ACh :=(f, vh)Th + (qh,∇vh)Th − ⟨q̂h · n , vh⟩∂Th
+ (qh +∇uh,ph)Th − ⟨uh − ûh , ph · n⟩∂Th
+ ⟨q̂h · n , v̂h⟩∂Th\∂Ω
+ ⟨uh − ûh , (ph − p̂h) · n⟩∂Th .
The derivation of ACh is presented in [27], Theorem 2.1. For general integral func-
tionals, including linear boundary integral functionals and non-linear domain integral
functionals, the derivation of ACh is discussed in [62, 41].
12
Filtering the oscillations of Galerkin Approximations
It is well-known that numerical solutions defined by a Galerkin method must oscillate
around the exact solutions due to the Galerkin orthogonality property. In [9], Bramble
and Schatz showed that it is possible to filter out these oscillations by convolving the
Galerkin solution with a B-spline kernel under the assumption that the test function
spaces are translation-invariant. As a result, a new approximation with a faster conver-
gence rate can be obtained. We use this convolution filtering technique to post-process
the approximate eigenfunction uh and obtain a more accurate approximation u
∗
h.
Let us first recall the definition of the B-spline kernel in one dimension. Let χ be
the function that is one on the interval (−12 , 12) and zero outside of it, and set ψ(0) = δ,
where δ is the Dirac delta function. The n-th order B-spline ψ(n) is defined as the
convolution of ψ(n−1) and χ, namely, ψ(n) = ψ(n−1) ∗ χ for n ≥ 1. The B-spline kernel
is defined as K2kh (x) :=
1
hK
2k(xh), where
K2k(x) :=
k∑
r=−k
Crψ
(k+1)(x− r). (2.2.8)
Here h is the size of the diameters of the locally-invariant mesh and k is a fixed integer,
usually the polynomial degree used in the Galerkin method. The coefficients of the
kernel Cr are determined by requiring that p ∗K2k = p for all polynomials p of degree
at most 2k. In Table 2.2, we display the coefficients Cr of the kernelK
2k for k = 1, . . . , 5.
They were computed according to the algorithm described in [70, Lemma 1] with the
parameters p = q = k + 1 by using a very short code in Mathematica.
13
k C0 C1 C2 C3 C4 C5
1 76 − 112
2 437320 − 97480 371920
3 122237560 − 9192520 3115040 − 417560
4 180179759289728 − 680345911612160 315395923224320 − 354111658880 15361792897280
5 1569217665280 −11796491330560 282511088 − 20813380160 307733991680 − 42017983360
Table 2.2: Coefficients of the convolution kernel K2k for k = 1, . . . , 5. Since the kernel
is symmetric, the coefficients are such that C−r = Cr. Polynomials of degree 2k are
reproduced by convolution with the kernel.
For the N-dimensional kernel, we define
K2kh (x) :=
1
hN
K2kN (
x
h
) where K2kN (x) :=
N∏
i=1
K2k(xi),
and the accuracy-enhanced approximation is then defined as
u∗h(x) := (K
2k
h ∗ uh)(x).
Note that the kernel K2kh has a finite support [−3k+12 h, 3k+12 h]N in N-dimensional space.
Hence, for a given point x, the convolution involves only a fixed number of neighboring
elements. This means that the convolution is local and it also implies that in general we
can only apply this technique to regions that are at least 3k+12 h away from the domain
boundary. In the next subsection, we present how we overcome this limitation and
obtain u∗h on the whole domain in the adjoint-based method. Details on the efficient
implementation of the convolution can be found in [57, 27].
2.2.2 The definition of the adjoint-based method
We are now ready to define the adjoint-based approximation of the functional J(u).
Let us start with two assumptions for the mesh.
14
Two assumptions for the meshes
Assumption 3. Let {Th}h>0 be a family of shape-regular meshes of Ω. We assume
that there exist subdomains Ω0 ⊂⊂ Ω′0 ⊂⊂ Ω such that
1. Th ∩ Ω′0 is a translation-invariant mesh with element size h0,
2. Any element K ∈ Th is fully contained in either Ω0 or Ω1 := Ω\Ω0,
3. Any boundary element is allowed to have one curved face.
We denote h := min{hK |K ∈ Th}, where hK is the diameter of an element K ∈ Th.
To get the convergence results in section 2.3, we also need the following assumptions on
the elements lying on the boundary; see Fig. 2.1 for an illustration of such elements.
With these assumptions, the classic trace and inverse estimates for elements with flat
faces can be extended to elements with curved faces; see [13, 12] for details.
Assumption 4. Let ∂Ω be Lipschitz continuous. Let K be a boundary element with
one possibly curved face E = K ∩ ∂Ω. We assume that:
1. (Star-shapedness) The element K is star-shaped with respect to all vertices opposite
to the face E. We assume K is also star-shaped with respect to all the midpoints
of the edges sharing a common vertex with the face E.
2. (Shape-regularity) There is a constant c such that m(x) · n(x) ≥ c|m(x)| uni-
formly across the mesh, for every vector m(x) = x − x0, with x ∈ E and x0
any vertex opposite E, and n(x) the unit outward normal vector of E at x. We
further assume |m(x)| = O(hK) uniformly for all boundary elements.
E E E
x0
m(x)
x
n(x)
Figure 2.1: Illustration of star-shapedness (left, middle) and shape-regularity(right)
15
The adjoint-based method
To define the adjoint-based method, we need the following notation. We denote by
Fh the set of all faces F of all elements K ∈ Th. Set, for i = 0, 1,
Tih :={K ∈ Th | K ⊂ Ωi},
∂Tih :={∂K ∈ Th | K ⊂ Ωi},
Fih :={F ∈ F(K) | K ∈ Tih}.
Definition 1. We define the adjoint-based method in four steps:
1. On the whole mesh Th, let (qkh, ukh, q̂ kh ·n, û kh ) and (pkh, vkh, p̂ kh ·n, v̂ kh ) be the approx-
imations of the model problem (2.1.2) and adjoint problem (2.2.6), respectively,
by the HDG method with polynomial degree k ≥ 1 .
2. Next, on T0,h, we use the filtering technique described in subsection 2.2.1 to obtain
u∗h := K
2k
h ∗ ukh and v∗h := K2kh ∗ vkh.
3. On T1,h, let (q2kh , u2kh , q̂ 2kh · n, û 2kh ) and (p2kh , v2kh , p̂ 2kh · n, v̂ 2kh ) be the approxima-
tions of the model problem and adjoint problem, respectively, by the HDG method
with polynomial degree 2k . To provide the boundary conditions on ∂Ω1\∂Ω, we
define û2kh := K
2k
h ∗ukh and v̂2kh := K2kh ∗vkh to be the Dirichlet boundary conditions
for the model and adjoint problems. Then we set the approximations u∗h := u
2k
h
and v∗h := v
2k
h on T1,h.
4. Finally, we compute Jh(u
∗
h, v
∗
h) = J(u
∗
h)+ACh as the adjoint-based approximation
to J(u), where
u∗h :=
K2kh ∗ ukh in K ∈ T0,h,u2kh in K ∈ T1,h,
and
v∗h :=
K2kh ∗ vkh in K ∈ T0,h,v2kh in K ∈ T1,h.
Note that the adjoint-correction term ACh defined in subsection 2.2.1 depends not
only on the scalar approximations u∗h, v
∗
h but also on their gradients and traces on each
16
face. To do that in the domain Ω0, we naturally use (∇u∗h, u∗h,∇u∗h|∂T0,h ·n, u∗h|F0h) and
(∇v∗h, v∗h,∇v∗h|∂T0,h ·n, v∗h|F0h). In the domain Ω1, however, we have two options, which
lead to two methods.
Method 1: Use the piecewise gradients in Ω1
The first method uses (∇u2kh , u2kh , q̂2kh · n, û 2kh ) and (∇v2kh , v2kh , p̂2kh · n, v̂ 2kh ) as the
approximation of the gradients and traces on Ω1. Since we use piecewise gradients, let
us denote this method with superscript “G”. We have
JGh (u
∗
h, v
∗
h) = J(u
∗
h) +AC
G
h , (2.2.9)
where
ACGh :=(f, v
∗
h)Th − (∇u∗h,∇v∗h)Th
+ ⟨q̂2kh · n , v̂ 2kh − v2kh ⟩∂T1,h
+ ⟨û 2kh − u2kh , p̂2kh · n⟩∂T1,h .
Method 2: Use the approximate fluxes in Ω1
For our second method, we use (q2kh , u
2k
h , q̂
2k
h ·n, û 2kh ) and (p2kh , v2kh , p̂2kh ·n, v̂ 2kh ) on
Ω1, which come from our Galerkin approximation. Let us denote it by the superscript
“F”. We have
JFh (u
∗
h, v
∗
h) = J(u
∗
h) +AC
F
h , (2.2.10)
where
ACFh :=(f, v
∗
h)T0,h − (∇u∗h,∇v∗h)T0,h
+ ⟨q̂2kh · n , v̂ 2kh ⟩∂Ω1\∂Ω
+ ⟨û 2kh − u2kh , p2kh − p̂2kh · n⟩∂T1,h .
For a detailed deduction of the above formulas, see [27].
17
2.3 The A Priori Error Estimates
We are now ready to present the a priori error estimates for the adjoint-based ap-
proximations. To do that, let us first introduce the notion of (s + 2)-regularity of the
domain Ω. For any integer s ≥ 0, we say Ω is (s + 2)-regular if the adjoint problem
(2.2.6) is uniquely solvable for g ∈ L2(Ω) and
∥p∥s+1,Ω + ∥v∥s+2,Ω ≤ C∥g∥s,Ω
for all g ∈ C∞0 (Ω), where C only depends on the domain Ω. Note that this inequality
does not necessarily hold for g ∈ Hs(Ω) but only for g ∈ C∞0 (Ω). For example, square
domains are (s+ 2)-regular according to this definition, see [60, Section 7, Example 3].
We can now state our main result.
Theorem 2.1. Suppose that k ≥ 1 and that Ω is (k+2)-regular. Let u , v ∈ H2k+3(Ω).
Let the HDG method satisfy Assumptions 1 and 2 where {Th}h>0 is a family of meshes
satisfying Assumptions 3 and 4. Finally, let JGh (u
∗
h, v
∗
h) and J
F
h (u
∗
h, v
∗
h) be the ap-
proximations of J(u) =
∫
Ω g(x)u(x)dx given by the adjoint-based methods (2.2.9) and
(2.2.10), respectively. Then, for h small enough, there exists a constant C such that
|J(u)− JGh (u∗h, v∗h)| ≤ C h4k,
|J(u)− JFh (u∗h, v∗h)| ≤ C h4k,
where C is independent of h, but depends on u, v, k and the subdomain Ω0.
Let us briefly discuss this result. First, we must note that, it is possible to show
that
|J(u)− J(uh)| ≤ C h2k+1,
just by using the smoothness of the functional J and by using that, for big enough
values of ℓ ∈ N, the H−ℓ(Ω)-norm of the error u−uh is of order h2k+1. There is no need
to assume that the meshes have to be translation invariant inside a fixed subdomain
of Ω. However, if we do assume this, the above theorem states that, by computing an
approximation vh to the adjoint solution v, and effectively doubling the computational
complexity, we can obtain an order of convergence of h4k.
18
It is worth noting that the numerical experiments in [27] indicate that the order of
convergence of J(uh), h
2k+1, and that of JGh (u
∗
h, v
∗
h) and J
G
h (u
∗
h, v
∗
h), h
4k, are actually
sharp. The numerical experiments carried out there used HDG methods with local
spaces admittingM -decompositions for triangular elements. Here, we use HDGmethods
with local spaces admitting M -decompositions for square elements. We verify that the
theoretical results do hold for this case, even though in some cases, the observed orders
of convergence are higher than the ones predicted.
2.4 Proof of the main theorem
In this section, we provide the proof of the error estimates of Theorem 2.1. We
proceed in several steps. First, we recall the formula for the error from [27] and use
it to get explicit expressions for the errors of the two methods we are considering. We
then estimate the terms of the error associated to the interior domain Ω0 and those
associated to the exterior domain Ω1. We conclude by putting those estimates together.
Step 1: A formula of the error
We begin by recalling a result that gives us an explicit formula of the approximation
error.
Theorem 2.2 ([27]). Let J(u) be the linear functional defined in (2.1.1). Let Jh(uh)
be the approximation defined in (2.2.7). Then we have that
J(u) = Jh(uh) + Eh,
where
Eh :=(q − qh,p− ph)Th
+ (q − qh,ph +∇vh)Th + (qh +∇uh,p− ph)Th
+ ⟨(q̂h − q) · n , vh − v̂h⟩∂Th + ⟨uh − ûh , (p̂h − p) · n⟩∂Th
If we apply this result to the method (2.2.9), for which Jh(uh) = J
G
h (u
∗
h, v
∗
h), we get
19
that
EGh :=J(u)− JGh (u∗h, v∗h)
=(q +∇K2kh ∗ ukh,p+∇K2kh ∗ vkh)T0,h
+ (q − q2kh ,p− p2kh )T1,h
+ ⟨(q̂2kh − q) · n , v2kh − v̂2kh ⟩∂T1,h + ⟨u2kh − û2kh , (p̂2kh − p) · n⟩∂T1,h ,
and if we apply it to the method (2.2.10), for which Jh(uh) = J
F
h (u
∗
h, v
∗
h), we obtain
EFh :=J(u)− JGh (u∗h, v∗h)
=(q +∇K2kh ∗ ukh,p+∇K2kh ∗ vkh)T0,h
+ (q − q2kh ,p− p2kh )T1,h
+ (q − q2kh ,p2kh +∇v2kh )T1,h + (q2kh +∇u2kh ,p− p2kh )T1,h
+ ⟨(q̂2kh − q) · n , v2kh − v̂2kh ⟩∂T1,h + ⟨u2kh − û2kh , (p̂2kh − p) · n⟩∂T1,h
Step 2: A basic result to get estimates involving convolutions
From the previous step, it is clear that we need to obtain estimates of functions of
the form w −K2kh ∗ wh. To state the basic result we want, we need to introduce some
standard notation.
For any integer m ≥ 0, open bounded set D ⊂ RN and sufficiently smooth function
u : D → R, we set
∥u∥m,D := (
∑
0≤ℓ≤m
|u|2ℓ,D)
1
2 , |u|ℓ,D := (
∑
|α|=ℓ
∫
D
|Dαu|2dx) 12 ,
∥u∥−m,D := sup
v∈C∞0 (D)
∫
D uvdx
∥v∥m,D ,
where C∞0 (D) denotes the space of infinitely differentiable functions on D with compact
support in D. If m = 0, we simply write ∥u∥0,D as ∥u∥D.
For any multi-index α = (α1, ..., αN ) and h > 0, let us define the difference quotient
as
∂αh u = ∂
α1
h,1 · · · ∂αNh,N u,
20
where
∂h,ju(x) =
1
h
(u(x+
h
2
ej)− u(x− h
2
ej)),
and ej is the unit vector with one in the j-th entry and zero otherwise.
We are now ready to recall the following key result.
Theorem 2.3 ([9, 25]). For a fixed integer k ≥ 1, define K2kh (x) := 1hNK2kN (x/h) as
described in 2.2.1. Let U be a function in L2(Ω), where Ω is an open set in RN , and
u ∈ H2k+1(Ω). Let Ω0 be an open subset of Ω and there exists h0 such that Ω0 +
2supp(K2kh ) ⊂⊂ Ω˜0 ⊂ Ω for any h < h0. Then for any h < h0, we have
∥u−K2kh ∗ U∥Ω0 ≤ C1h2k+1|u|2k+1,Ω˜0 + C2
∑
|α|≤k
∥∂αh (u− U)∥−k,Ω˜0 ,
where C1 and C2 are independent of u and h.
For this estimate to be useful, we need to obtain an estimate of negative-order norms
of the difference quotient ∂αh (u−U) for the case in which U is the component uh of the
approximation provided by the HDG method. The estimate we need is contained in the
following result. Its proof follows from the duality argument and the observation that
∂αh uh satisfies the local HDG equations. The proof will be omitted, as it is a variation
of that of the straight-faced elements.
Lemma 2.4. Suppose k ≥ 1 and Ω is (k+2)-regular. Let α be a fixed multi-index, and
assume the exact solution u ∈ Hk+2+|α|(Ω). Consider an HDG method with polynomial
degree k satisfying Assumptions 1 and 2, for a family of meshes of Ω, {Th}h>0, satisfying
Assumptions 3 and 4. Then for h small enough we have
∥∂αh (u− uh)∥−k,Ω˜0 ≤ Ch
k+l+1(∥u∥l+2+|α|,Ω + ∥u∥l+2,Ω),
where Ω0⊂⊂Ω˜0⊂⊂Ω′0⊂⊂Ω and 0 ≤ l ≤ k.
Step 3: Estimates of the terms defined in Ω0.
The following result contains the estimates of the terms associated to the subdomain
Ω0.
21
Lemma 2.5. We have
∥q +∇K2kh ∗ ukh∥T0,h ≤ Ch2k+1∥u∥2k+3,Ω,
where C only depends on k, Ω0, Ω
′
0 and Ω .
Proof. To prove the inequality, we begin by noting that,
∥q +∇K2kh ∗ ukh∥T0,h = ∥∇u−∇K2kh ∗ ukh∥T0,h ≤ I + II,
where I = ∥∇u−K2kh ∗ ∇u∥T0,h and II = ∥∇(K2kh ∗ (u− ukh))∥T0,h . To estimate I, we
apply Theorem 2.3 with u := U := ∂xiu, for each component i = 1, . . . , d, to get that
∥∇u−K2kh ∗ ∇u∥T0,h ≤ Ch2k+1∥∇u∥2k+1,Ω ≤ Ch2k+1∥u∥2k+2,Ω.
It remains to get estimate for IIi := ∥∂xi(K2kh ∗ (u− ukh))∥T0,h for i = 1, ..., d. We have
IIi ≤ C
∑
|α|≤k
∥Dα∂xi(K2kh ∗ (u− ukh))∥−k,Ω˜0 by [9, Lemma 2.2],
= C
∑
|α|≤k
∥Dα+eiK2kh ∗ (u− ukh)∥−k,Ω˜0 ,
≤ C
∑
|α|≤k
∥∂α+eih (u− ukh)∥−k,Ω˜′0 ,
by [9, Lemma 5.3]. Here, Ω0 ⊂⊂ Ω˜0 ⊂⊂ Ω˜′0 ⊂⊂ Ω′0. Now we use the negative-order
norm estimate of Lemma 2.4 for the HDG method with l := k and α := α+ ei to get
that, for h is small enough,
∥∂α+eih (u− ukh)∥−k,Ω˜′0 ≤ Ch
2k+1(∥u∥2k+3,Ω + ∥u∥k+2,Ω),
since |α| ≤ k + 1. This completes the proof.
22
Step 4: Estimates of the terms defined in Ω1.
Let us now estimate the terms of the errors defined in Ω1. In the following result,
we gather all the estimates we need. It is stated in terms of norms we define next:
∥w∥T1,h := (
∑
K∈T1,h
∥w∥2K)1/2 ∀ w ∈ L2(T1,h).
∥µ∥hα,∂T1,h := (
∑
K∈T1,h
hαK ∥µ∥2∂K)1/2 ∀ µ ∈ L2(∂T1,h),
∥µ∥hα0 ,∂Ω0 := (
∑
F∈∂Ω0
hα0 ∥µ∥2F )1/2 ∀ µ ∈ L2(∂Ω0),
where h0 is the diameter of the elements of the translation-invariant mesh T0,h in Ω0.
Lemma 2.6. Let (u2kh , q
2k
h , q̂
2k
h · n, û2kh ) be the HDG approximation with polynomial
degree 2k in the domain Ω1 as described in Definition 1. Let Πh(u, q) = (ΠWu,ΠV q)
be the corresponding HDG projection. Then for k ≥ 1, we have the estimates of the
errors
∥q − q2kh ∥T1,h ≤ C Θ1, (2.4.1a)
∥q · n− q̂2kh · n∥h,∂T1,h ≤ C(Θ1 +Θ2 +Θ3), (2.4.1b)
and the estimates of the residuals
∥q2kh +∇u2kh ∥T1,h ≤ C(Θ1 +Θ2), (2.4.1c)
∥u2kh − û2kh ∥h−1,∂T1,h ≤ C(Θ1 +Θ2), (2.4.1d)
where
Θ1 := ∥ΠV q − q∥T1,h + ∥PMu−K2kh ∗ ukh∥h−10 ,∂Ω0 ,
Θ2 := ∥u−ΠWu∥h−1,∂T1,h .
Θ3 := ∥(ΠV q − q) · n∥h,∂T1,h
Here PMu is the L
2 projection of u intoMh and the constant C depends on k, u,maxK∈T1h{τmaxK hK},Ω0
and Ω1.
23
The proof of this result entails carrying out a small modification of the standard a
priori error analysis of the HDG method. The modification is due to the fact that the
boundary condition on ∂Ω0 is given by the trace of K
2k
h ∗ ukh from Ω0, not by the exact
solution u. Since the proof is fairly long, we divide it in several steps.
Step i: The equations for the projection of the errors.
We begin by finding the equations satisfied by the projection of the errors, namely,
εqh :=ΠV q − q2kh , εuh := ΠWu− u2kh and εûh := u− û2kh .
Using the properties of the HDG projection (2.2.4), by [23, Lemma 3.1], we have
the following error equations,
(εqh, r)T1,h − (εuh,∇ · r)T1,h + ⟨εûh , r · n⟩∂T1,h = (ΠV q − q, r)Th , (2.4.2a)
−(εqh,∇w)T1,h + ⟨εq̂h · n , w⟩T1,h = 0, (2.4.2b)
⟨εûh , µ⟩∂Ω1 = ⟨PMu−K2kh ∗ ukh , µ⟩∂Ω0 , (2.4.2c)
⟨εq̂h · n , µ⟩∂T1,h\∂Ω1 = 0, (2.4.2d)
for all r ∈ V h, w ∈Wh and µ ∈Mh, where
εq̂h · n = εqh · n+ τ(εuh − εûh) on ∂T1h. (2.4.2e)
Note that the right-hand side of (2.4.2c) is not zero due to the fact that the Dirichlet
condition on ∂Ω0 is given by the trace of K
2k
h ∗ ukh from Ω0, not by the exact solution.
Step ii: Energy Argument.
To obtain the estimates of the errors, we use a standard energy argument. So, taking
r := εqh in (2.4.2a), w := ε
u
h in (2.4.2b), µ := −εq̂h ·n in (2.4.2c) and µ := −εûh in (2.4.2d)
and adding the resulting four equations, we obtain
Eh := (ε
q
h, ε
q
h)T1,h + ⟨τ(εuh − εûh) , (εuh − εûh)⟩∂T1,h = T1 + T2,
where T1 := (ΠV q − q, εqh)T1,h and T2 := ⟨PMu−K2kh ∗ ukh , −εq̂h · n⟩∂Ω0 .
24
Step iii: Estimates of the errors.
We estimate the first term as follows:
T1 ≤ ∥ΠV q − q∥T1,h ∥εqh∥T1,h ≤ ∥ΠV q − q∥T1,h E1/2h ,
since τ is non-negative.
Let us now estimate the second term. We have
T2 ≤ (
∑
F∈∂Ω0
h−10 ∥PMu−K2kh ∗ ukh∥2F )
1
2 (
∑
F∈∂Ω0
h0∥εq̂h · n∥2F )
1
2 .
≤ ∥PMu−K2kh ∗ ukh∥h−10 ,∂Ω0(
∑
K∈Ω1,∂K∩∂Ω0=F
hK∥εq̂h · n∥2F )
1
2 .
≤ ∥PMu−K2kh ∗ ukh∥h−10 ,∂Ω0∥ε
q̂
h · n∥h,∂T1,h .
By the definition of εq̂h · n, (2.4.2e), we have that
∥εq̂h · n∥h,∂T1,h = ∥εqh · n+ τ(εuh − εûh)∥h,∂T1,h
≤ ∥εqh∥h,∂T1,h + maxK∈T1,h(hKτ
max
K )
1
2 ⟨τ(εuh − εûh) , (εuh − εûh)⟩
1
2
T1,h
≤ C1,τ (∥εqh∥T1,h + ⟨τ(εuh − εûh) , (εuh − εûh)⟩
1
2
T1,h)= C1,τ E
1/2
h ,
where C1,τ = Cmax{1, (hKτmaxK )
1
2 , K ∈ T1,h}. Here we use the inverse inequality for
polynomials. As a consequence, we get
T2 ≤ C1,τ ∥PMu−K2kh ∗ ukh∥h−1,∂Ω0 E1/2h ,
and the first estimate (2.4.1a) immediately follows.
Let us obtain the second estimate. By the definition of εq̂h · n, (2.4.2e), we have
∥q · n− q̂2kh · n∥h,∂T1,h ≤ ≤ ∥q · n−ΠV q · n− τ(ΠWu− u)∥h,∂T1,h + ∥εq̂h · n∥h,∂T1,h
≤ ∥(q −ΠV q) · n∥h,∂T1,h + ∥τ(ΠWu− u)∥h,∂T1,h + ∥εq̂h · n∥h,∂T1,h
≤ ∥(q −ΠV q) · n∥h,∂T1,h + C2,τ∥ΠWu− u∥h−1,∂T1,h + ∥εq̂h · n∥h,∂T1,h ,
where C2,τ = max{1, (hKτmaxK ), K ∈ T1,h}. Then the estimate (2.4.1b) follows easily.
25
Step iv: Estimates of the residuals
It remains to prove the estimates of the residuals. Integrating by parts in the first
equation defining the HDG method in Ω1, (2.2.2a), we get
(q2kh +∇u2kh ,v)T1,h = ⟨u2kh − û2kh , v · n⟩∂T1,h .
Taking v := q2kh +∇u2kh on each element K ∈ T1h and zero elsewhere, we obtain that
∥q2kh +∇u2kh ∥2K ≤ ∥u2kh − û2kh ∥∂K∥(q2kh +∇u2kh ) · n∥∂K
≤ Ch−
1
2
K ∥u2kh − û2kh ∥∂K∥q2kh +∇u2kh ∥K ,
by an inverse inequality. Hence
∥q2kh +∇u2kh ∥K ≤ Ch
− 1
2
K ∥u2kh − û2kh ∥∂K .
And also note that
∥u2kh − û2kh ∥h−1,∂T1,h = ∥u2kh −ΠWu+ΠWu− u+ u− û2kh ∥h−1,∂T1,h
≤ ∥εuh − εûh∥h−1,∂T1,h + ∥ΠWu− u∥h−1,∂T1,h
≤ ∥q − qh∥T1,h + ∥ΠWu− u∥h−1,∂T1,h ,
where the estimate for ∥εuh − εûh∥h−1,∂T1,h in the last inequality is proven in detail in
Appendix A.1.4.
The estimates of the residuals easily follow. This completes the proof of Lemma 2.6.
Step 5: Estimates of the auxiliary quantities Θ1, Θ2, and Θ3
To be able to conclude, it only remains to estimate the three terms that define Θ1,
Θ2, and Θ3 in the lemma of the previous step. The estimates are displayed in the
following result.
26
Lemma 2.7. We have
∥ΠV q − q∥T1,h ≤ Ch2k+1∥u∥2k+2,Ω1 ,
∥(ΠV q − q) · n∥h,∂T1,h ≤ Ch2k+1∥u∥2k+2,Ω1 ,
∥u−ΠWu∥h−1,∂T1,h ≤ Ch2k∥u∥2k+2,Ω1 ,
∥PMu−K2kh ∗ ukh∥h−10 ,∂Ω0 ≤ Ch
2k
0 ∥u∥2k+2,Ω,
where C depends on k, Ω0, Ω
′
0 and Ω .
Proof. The first three inequalities follow from Assumption 1 and 2 about the HDG
method and the properties of the HDG projection; see Appendix A.1.2. Let us prove
the last estimate.
For each face F of the element K ∈ Th ⊂ Ω0, let UF,K be the function in P2k(K)
such that
UF,K = PMu on F,
(UF,K − u,w)K = 0 ∀w ∈ P2k(K) : w = 0 on F.
Then,
T := (
∑
F∈∂Ω0
h0 ∥PMu−K2kh ∗ ukh∥2F )
1
2 = (
∑
F∈∂Ω0
h0 ∥UF,K −K2kh ∗ ukh∥2F )
1
2
≤ C (
∑
K∈Ω0,∂K∩∂Ω0=F
∥UF,K −K2kh ∗ ukh∥2K)
1
2 ,
by an inverse inequality. Then T ≤ T1 + T2, where
T1 := C (
∑
K∈Ω0,∂K∩∂Ω0=F
∥UF,K − u∥2K)
1
2 ,
T2 := C (
∑
K∈Ω0
∥u−K2kh ∗ ukh∥2K)
1
2 .
A standard approximation theory gives us that
T1 ≤ Ch2k+1∥u∥2k+1,Ω0 ,
27
and, by Theorem 2.3 and by Lemma 2.4 with U := uh, the approximation of u given by
the HDG method, we get that
T2 = ∥u−K2kh ∗ ukh∥Ω0 ≤ Ch2k+1∥u∥2k+2,Ω.
This completes the proof.
Step 6: Conclusion
We are now ready to prove Theorem 2.1. We only have to recall the formulas of the
errors EGh and E
F
h in Step 1, apply the Cauchy-Schwarz inequality to each of the terms
of those formulas and use the estimates obtained in the previous Steps. Since each of
those terms is of the order of at least h2k, we get that both errors are of order h4k. This
completes the proof of Theorem 2.1.
2.5 Numerical results
In this section, we design several numerical experiments to explore the convergence
properties of the method under consideration. The numerical experiments in [27], car-
ried out with HDG methods using piecewise-polynomials of degree k ≥ 0 suggest that
the orders of convergence given by Theorem 2.1 are sharp. Here, we use HDG meth-
ods defined in squares and explore how close Ω0 can be to the boundary ∂Ω, how the
smoothness of the solution u affects the convergence rate, and how the smoothness of
the functional J affects the convergence rate.
We consider a unit square domain Ω = (0, 1)2. This domain has been shown to be
(k + 2)-regular for k ≥ 0 in [60, Section 7, Example 3]. For each natural number N ,
we obtain the mesh by dividing Ω into N2 uniform squares with h = 1N . Then for each
element K, we choose the local spaces
V (K) :=Qk(K)⊕ curlspan{xk+1y, xyk+1}, W (K) = Qk(K), M(∂K) = Qk(∂K),
so that Assumptions 1 and 2 are satisfied, see Table 2.1.
We implemented the numerical methods in MATLAB R2019a in a personal laptop
(MacBook Pro 2019). To show the high-order accuracy of the method, we used the
28
Multi-precision Computing Toolbox [1] and selected 32 digits of accuracy.
2.5.1 Illustration of the approximation error u− u∗h
We start by displaying the approximation error u−u∗h in the case in which the linear
functional is given by (2.1.1) with g(x, y) := 18π2 sin(3πx) sin(3πy). We choose the exact
solution to be u(x, y) := sin(3πx) sin(3πy) and determine the boundary conditions and
the source term f accordingly. We choose Ω0 to be region such that the distance between
Ω0 and the boundary ∂Ω is 0.1
In Fig 2.2, we display the approximation error plot u− u∗h when k = 1 and N = 40.
In Fig 2.3, we show 1D slices of the approximation error plot at x = 0.50117 and
y = 0.20117, respectively. From these figures, we can see that the size of the error in
region Ω\Ω0, using the HDG method with polynomial degree 2k, matches the size of
the error in Ω0. Similar error plots are obtained for the following cases so we omit them.
We refer the reader to [27] for more approximation error plots in 1D and 2D triangular
meshes.
Figure 2.2: Approximation error u−u∗h when k = 1 and N = 40. The distance between
Ω0 and ∂Ω is 0.1. Note that the magnitude of the error in Ω\Ω0 is essentially the same
as that in Ω0.
29
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1 10
-4
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1 10
-4
Figure 2.3: Slices of the approximation error u−u∗h at x = 0.50117 (left) and y = 0.20117
(right) when k = 1 and = 40. In each element, we evaluate u−u∗h at 5 Gauss quadrature
points.
2.5.2 How close Ω0 can be to the boundary ∂Ω
Next, we explore how close we can actually choose Ω0 to the boundary and how
the convergence rates might be effected. In our main theorem 2.1, we tacitly assume
Ω0 is independent of the mesh and of the polynomial degree k. In other words, our
theory holds if we fix the domain Ω0. In [27], however, numerical experiments in one-
and two-space dimensions were performed with Ω0 just (2k+1) elements away from the
boundary. The theoretical orders of convergence rate h4k, for k = 1, 2, 3 were actually
achieved when the solutions u and v are very smooth. From a practical standpoint,
the larger Ω0 ⊂⊂ Ω is, the less computational efforts are needed. This is because
the convolution filtering technique is fairly inexpensive when compared to solving the
boundary-value problem on Ω1 with the HDG method. Therefore, we want to choose
Ω0 to be as large as possible. On the other hand, based on Theorem 2.3, in order to use
the convolution filtering technique, Ω0 has to have at least 2k elements outside in each
direction. In the following numerical experiment, we consider the test problem of the
previous subsection. In this way, both solutions u and v are very smooth. In this ideal
case, we compare two cases: (1) Ω0 is a fixed domain such that the distance between
Ω0 and the boundary ∂Ω is 0.1, (2) Ω0 is defined by removing a boundary layer of 2k
elements from Ω.
30
Let us report the history of convergence for these two cases. The results of the first
case, where Ω0 is a fixed domain, is shown in Table 2.3. We see that the convergence rate
of u∗h is of order O(h2k+1), and that the accuracy of the two approximations JFh (u∗h, v∗h)
and JGh (u
∗
h, v
∗
h) are at least O(h4k) for k = 1, 2. These results are consistent with our
theoretical results.
k N ∥u− u∗h∥L2(Ω) order |J(u)− JFh (u∗h, v∗h)| order |J(u)− JGh (u∗h, v∗h)| order
1
40 4.82e-05 - 9.95e-07 - 8.50e-05 -
50 2.17e-05 3.59 2.58e-07 6.04 3.47e-05 4.02
60 1.14e-05 3.55 8.91e-08 5.84 1.67e-05 4.01
70 6.61e-06 3.51 3.70e-08 5.69 9.02e-06 4.00
2
40 5.4e-07 - 5.47e-11 - 6.69e-11 -
50 1.44e-07 5.93 3.85e-12 11.89 5.84e-12 10.93
60 4.87e-08 5.94 4.39e-13 11.91 8.91e-13 10.31
70 1.94e-08 5.95 6.99e-14 11.92 2.00e-13 9.71
Table 2.3: History of convergence when the distance between Ω0 and ∂Ω is fixed and
equal to 0.1.
We display a history of convergence when Ω0 varies with the mesh and k in Table
2.4. This case is not covered in the theory, but we can see that we still have super-
convergence. This implies that the constant C in our main theorem 2.1 can remain
bounded when we vary Ω0. Moreover, these results suggest that we should select Ω0
only O(h) away from the boundary since the computational effort is smaller and that
the results are better. The results also indicate that JFh (u
∗
h, v
∗
h) is more accurate than
JGh (u
∗
h, v
∗
h). This is most probably due to the better approximation to the gradient
used in the HDG method of degree 2k used for JFh (u
∗
h, v
∗
h).
31
k N ∥u− u∗h∥L2(Ω) order |J(u)− JFh (u∗h, v∗h)| order |J(u)− JGh (u∗h, v∗h)| order
1
40 4.79e-05 - 8.46e-07 - 4.57e-05 -
50 2.11e-05 3.67 1.94e-07 6.59 1.50e-05 5.00
60 1.08e-05 3.66 5.95e-08 6.50 6.03e-06 5.00
70 6.16e-06 3.65 2.21e-08 6.42 2.79e-06 5.00
2
40 5.61e-07 - 5.70e-11 - 6.73e-11 -
50 1.50e-07 5.91 4.03e-12 11.87 5.45e-12 11.26
60 5.07 e-08 5.95 4.59e-13 11.92 7.38e-13 10.96
70 2.02 e-08 5.97 7.28e-14 11.93 1.43e-13 10.63
Table 2.4: History of convergence when the distance between Ω0 and ∂Ω is, for each
mesh, is 2k elements.
Finally, in Table 2.5 we report a CPU time, error and efficiency comparison of
using the adjoint-based approximation JFh (u
∗
h, v
∗
h) and using the standard approximation
J(uh) with the HDG method. Here tAdj/tHDG denotes the ratio of the CPU time spent
on computing JFh (u
∗
h, v
∗
h) to that of computing J(uh), and eHDG/eAdj stands for the ratio
of the error |J(u) − J(uh)| to the error |J(u) − JFh (u∗h, v∗h)|. To measure the efficiency
of the method, we use the quantity ω := 1/(te). As we see in Table 2.5, the running
time of the adjoint-based method is about 3 times as much as that of the standard
method. However, the error of the adjoint-based method is several orders of magnitude
(from 103 to 105 times) smaller than that of the standard method. This is reflected
in the fact that the ratio of the efficiencies of the methods ωAdj/ωHDG is several orders
of magnitude bigger than one which shows the advantage of the adjoint-based method.
In Fig 2.4, we illustrate this advantage from another perspective. There we compare
the efficiencies of the adjoint-based method JFh (u
∗
h, v
∗
h) and the standard method J(uh)
when k = 1. We can see that to achieve the given error tolerance 9× 10−5, the adjoint-
based method only takes 0.9s with a mesh with N = 20, while the standard method
takes 49s with a mesh with N = 90. In other words, in this case, the adjoint-based
method is more than 50 times faster than the standard method.
32
k N
Ω0 fixed Ω0 varies with k and h
tAdj/tHDG eHDG/eAdj ωAdj/ωHDG tAdj/tHDG eHDG/eAdj ωAdj/ωHDG
1
40 3.23 1.07e+03 331 3.11 1.27e+03 408
50 3.41 2.15e+03 630 2.72 2.86e+03 1,051
60 3.43 3.58e+03 1,043 2.89 5.38e+03 1,861
70 3.43 5.38e+03 1,569 2.73 9.01e+03 3,300
2
40 4.16 1.27e+04 3,053 3.87 1.22e+04 3,152
50 3.62 5.56e+04 15,359 3.13 5.32e+04 16,997
60 2.94 1.87e+05 63,605 2.67 1.79e+05 67,041
70 3.07 5.26e+05 171,336 2.53 5.05e+05 199,604
Table 2.5: CPU time, error and efficiency comparison of the adjoint-based method
JFh (u
∗
h, v
∗
h) and the standard method J(uh) when u is a smooth function.
10 20 30 40 50 60 70 80 90 100
N
10-6
10-5
10-4
10-3
10-2
10-1
100
Fu
nc
tio
na
l e
rro
r
t=0.9s t=49s
| J(u) - J(uh) |
| J(u) - JF(u*h) |
Figure 2.4: Log scale error comparison of the adjoint-based method JFh (u
∗
h, v
∗
h) (◦) and
the standard method J(uh) (□) when u is a smooth function and k = 1. Note that to
achieve the given error tolerance 9 × 10−5, the adjoint-based method only takes 0.9s
with a mesh N = 20, while the standard method takes 49s with a mesh N = 90.
2.5.3 Effect of the smoothness of the solution for variable Ω0
In this section, we explore how the smoothness of the solution can affect the con-
vergence rate of the adjoint-based method. To do that, we compare the errors of J(uh)
and JFh (u
∗
h, v
∗
h) for solutions with different smoothness. Here we choose Ω0 to be the
subdomain of Ω which lies 2k elements away from the boundary for k=1,2,3.
33
Smooth solution u
As a first example, we consider the same problem as in the previous section. We
report the errors and the history of convergence in Table 2.6. As we can see, when the
exact solution is smooth, the HDG solution uh with polynomial degree k converges at
the rate of O(hk+1), the approximation J(uh) converges at the rate of O(h2k+1), the
post-processed solution u∗h converges at the rate of O(h2k+1), and the adjoint-based
method JFh (u
∗
h, v
∗
h) has the accuracy of at least O(h4k). We thus recover the orders of
convergence for the case in which the set Ω0 is fixed.
k N ∥u− uh∥L2(Ω) Order |J(u)− J(uh)| Order ∥u− u∗h∥L2(Ω) Order |J(u)− JFh (u∗h, v∗h)| Order
1
40 1.69e-03 - 1.07e-03 - 4.79e-05 - 8.46e-07 -
50 1.03e-03 2.21 5.54e-04 2.94 2.11e-05 3.67 1.94e-07 6.59
60 6.97e-04 2.15 3.19e-04 3.03 1.08e-05 3.66 5.95e-08 6.50
70 5.03e-04 2.11 1.99e-04 3.07 6.16e-06 3.65 2.21e-08 6.42
2
20 2.33e-04 - 2.88e-05 - 2.43e-05 - 1.38e-07 -
30 6.90e-05 3.00 3.23e-06 5.39 2.95e-06 5.21 1.65e-09 10.90
40 2.91e-05 3.00 6.97e-07 5.34 5.61e-07 5.77 5.70e-11 11.70
50 1.49e-05 3.00 2.14e-07 5.29 1.50e-07 5.91 4.03e-12 11.87
3
20 6.90e-06 - 1.22e-08 - 8.70e-07 - 1.74e-10 -
30 1.37e-06 4.00 4.89e-10 7.93 6.69e-08 6.32 1.04e-12 12.63
40 4.32e-07 4.00 4.94e-11 7.97 8.70e-09 7.09 1.49e-14 14.75
50 1.77e-07 4.00 8.46e-12 7.91 1.59e-09 7.63 4.67e-16 15.53
Table 2.6: History of convergence when the solution u is smooth and k = 1, 2, 3. The
distance between Ω0 and ∂Ω is, for each mesh, 2k elements.
Solution u with corner singularity
In this second example, we consider the same linear functional (2.1.1) and g(x, y) :=
ex+y, but we choose u(r, θ) := r
2
3 sin(23θ) in polar coordinates so that the exact solution
u is singular at the origin. We define the boundary conditions and the source term f
according to u. We use the adjoint-based method (2.2.10) to approximate the functional
and report the history of convergence in Table 2.7. As we can see, due to the lack of
smoothness of u, the convergence rate of the post-processed solution u∗h and the adjoint
approximation JFh (u
∗
h, v
∗
h) do not improve as the polynomial degree is increased. In
fact, in this case it can be shown that the HDG solution with polynomial degree k ≥ 1
converges in L2(Ω) with order O(h 53 ). Contrast this behavior against the case of smooth
solutions in the previous example in Table 2.6.
34
k N ∥u− uh∥L2(Ω) Order |J(u)− J(uh)| Order ∥u− u∗h∥L2(Ω) Order |J(u)− JFh (u∗h, v∗h)| Order
1
50 3.38e-05 - 4.68e-08 - 9.99e-06 - 8.60e-09 -
60 2.50e-05 1.66 2.88e-08 2.66 7.38e-06 1.66 5.32e-09 2.64
70 1.94e-05 1.66 1.91e-08 2.66 5.71e-06 1.66 3.54e-09 2.64
80 1.55e-05 1.66 1.34e-08 2.65 4.57e-06 1.66 2.49e-09 2.64
2
50 1.01e-05 - 8.42e-09 - 2.43e-06 - 1.01e-09 -
60 7.48e-06 1.66 5.18e-09 2.66 1.80e-06 1.66 6.20e-10 2.65
70 5.79e-06 1.66 3.44e-09 2.66 1.39e-06 1.66 4.12e-10 2.65
80 4.64e-06 1.66 2.41e-09 2.66 1.11e-06 1.66 2.89e-10 2.65
3
50 4.68e-06 - 3.04e-09 - 9.15e-07 - 1.52e-10 -
60 3.45e-06 1.66 1.87e-09 2.66 6.76e-07 1.66 9.38e-11 2.65
70 2.67e-06 1.66 1.24e-09 2.66 5.24e-07 1.66 6.22e-11 2.66
80 2.14e-06 1.66 8.70e-10 2.66 4.20e-07 1.66 4.37e-11 2.65
Table 2.7: History of convergence when the solution u has a corner singularity and
k = 1, 2, 3. The distance between Ω0 and ∂Ω is, for each mesh, 2k elements.
k N tAdj/tHDG eHDG/eAdj ωAdj/ωHDG
1
50 1.63 5.43 3.33
60 1.63 5.41 3.32
70 1.57 5.41 3.45
80 1.70 5.41 3.18
2
50 1.62 8.40 5.19
60 1.86 8.33 4.48
70 1.99 8.33 4.19
80 2.17 8.33 3.84
3
50 2.29 20.0 8.73
60 2.20 19.6 8.91
70 2.24 19.6 8.75
80 2.44 19.6 8.03
Table 2.8: CPU time, error and efficiency comparison of the adjoint-based method
JFh (u
∗
h, v
∗
h) and the standard method J(uh) when u has a corner singularity. The set Ω0
varies with k and h.
Let us also report the CPU time, error and efficiency comparison for this case
in Table 2.8. Even though the convergence rates of the adjoint-based approximation
JFh (u
∗
h, v
∗
h) is the same as the standard approximation J(uh), the error |J(u)−JFh (u∗h, v∗h)|
is still much smaller (from 5 to 20 times) than |J(u) − J(uh)| for k = 1, 2, 3 (see
eHDG/eAdj in Table 2.8). In addition, the running time of the adjoint-based approxi-
mation is only about twice that of the standard method (see tAdj/tHDG in Table 2.8).
Finally, the ratio of efficiencies ωAdj/ωHDG shows that the adjoint-based method is
around 3 times more efficient than the standard method for k = 1, 4 times for k = 2
and 8 times for k = 3. This indicates that even though the exact solution u displays a
35
corner singularity, we can still benefit from using the adjoint-based method.
2.5.4 Effect of the smoothness of the functional
The smoothness of a functional is reflected by the smoothness of the solution of the
adjoint problem. In this section, we consider a functional on the boundary
J(u) := ⟨∇u · n, ψ⟩∂Ω =< −q · n, ψ >∂Ω, (2.5.1)
where u is the solution of the model problem (2.1.2) and q = −∇u. If the dual problem
−∆v = 0 in Ω, v = ψ on ∂Ω,
has a weak solution v ∈ H1(Ω), then the functional can be written as
J(u, v) = (∇u,∇v)Ω − (f, v)Ω. (2.5.2)
Therefore, we can simply define an approximation J(uh, vh) = (qh,ph)− (f, vh), where
qh and ph are the approximations to −∇u and −∇v, respectively. We can now design
an adjoint-correction term ACh and define a new approximation
JIh = J(uh, vh) +ACh,
where
ACh =− (qh,ph) + (uh,∇ph)− ⟨ûh,ph · n⟩
− (qh,ph +∇vh)− ⟨v̂h − vh , ph · n⟩
+ ⟨(qh − q̂h) · n, v̂h − vh⟩.
After a simple calculation, we get that the error term is
Eh := J(u)− JIh
= (qh − q,ph +∇vh) + (u− uh,∇ · (ph − p))
+ ⟨(q̂h − q) · n , v − vh⟩ + ⟨(ph − p) · n , ûh − u⟩
36
In the following test cases, we compare two approximations of the functional defined by
(2.5.1), namely,
JBh :=< −q̂h · n, ψ >∂Ω,
J∗Ih := J(u
∗
h, v
∗
h) +ACh.
The first approximation only uses the HDG approximate solutions. The second uses
the approximations u∗h and v
∗
h, and the HDG flux approximation to compute ACh in
the same manner as (2.2.10).
Smooth solution u, but non-smooth solution v
In this first case, we choose the exact solution to be a smooth function u = sin(πx) sin(πy)
and let v = π−θπ , where θ is the angle component of the polar coordinates at (
1
8 , 0). We
define the boundary conditions accordingly. Note that ψ is a discontinuous function on
boundary since ψ(x, 0) = 1 when x > 18 and ψ(x, 0) = 0 when x <
1
8 .
k N |J(u)− JBh| Order |J(u)− J∗h| Order
1
14 2.85e-05 - 3.30e-07 -
22 7.82e-06 2.86 5.26e-08 4.06
30 3.18e-06 2.91 1.43e-08 4.20
38 1.59e-06 2.93 5.29e-09 4.20
2
14 3.34e-07 - 1.17e-10 -
22 4.89e-08 4.25 7.69e-12 6.02
30 1.33e-08 4.20 1.17e-12 6.07
38 4.98e-09 4.16 2.62e-13 6.32
3
14 4.42e-09 - 7.97e-15 -
22 4.37e-10 5.12 1.50e-16 8.79
30 9.05e-11 5.08 1.03e-17 8.63
38 2.72e-11 5.08 1.39e-18 8.46
Table 2.9: History of convergence for the boundary integral functional with v = π−θπ .
As shown in Table 2.9, the convergence rates of JBh is O(hk+2) and the convergence
rate of J∗Ih is O(h2k+2). Compared with the case in section 2.5.2, where the functional is
an integral over the whole domain Ω and the approximation accuracy is O(h2k+1) with
HDG method, and O(h4k) with the adjoint-based method, we see that the accuracy of
our approximation decreases. This result is reasonable because in this case the solution
of the adjoint problem v is not in H1(Ω) and the accuracy of our approximation solution
37
vh can only be O(h) as shown in Table 2.10.
k N ∥u− uh∥L2(Ω) order ∥u− u∗h∥L2(Ω) order ∥v − vh∥L2(Ω) order ∥v − v∗h∥L2(Ω) order
1
14 1.34e-03 - 5.97e-05 - 5.43e-03 - 3.48e-03
22 5.41e-04 2.02 1.24e-05 3.48 3.46e-03 1.00 2.20e-03 1.01
30 2.90e-04 2.01 4.29e-06 3.41 2.54e-03 1.00 1.63e-03 0.98
38 1.81e-04 2.01 1.94e-06 3.36 2.01e-03 1.00 1.29e-03 0.99
2
14 2.52e-05 - 4.13e-07 - 3.99e-03 - 2.44e-03 -
22 6.48e-06 3.00 2.94e-08 5.84 2.56e-03 0.98 1.56e-03 0.99
30 2.56e-06 3.00 4.68e-09 5.93 1.89e-03 0.99 1.15e-03 0.99
38 1.26e-06 3.00 1.15e-09 5.93 1.49e-03 0.99 9.09e-04 0.99
3
14 3.56e-07 - 4.44e-09 - 2.70e-03 - 1.39e-03 -
22 5.83e-08 4.00 1.74e-10 7.17 1.72e-03 0.99 8.89e-04 0.99
30 1.69e-08 4.00 1.51e-11 7.88 1.26e-03 1.00 6.53e-04 1.00
38 6.55e-09 4.00 2.30e-12 7.96 9.99e-04 1.00 5.16e-04 1.00
Table 2.10: History of convergence of approximations to u = sin(πx) sin(πy) and to
v = π−θπ .
We also report the CPU time, error and efficiency comparison in Table 2.11. We
can see that as we refine the mesh, the running time of the adjoint-based approximation
is only about 3 to 5 times that of the standard method (see tAdj/tHDG), however the
error of the adjoint-based method is significantly smaller (from 102 to 107 times) than
that of the standard method (see eHDG/eAdj). The values of the ratio of efficiencies
ωAdj/ωHDG are always bigger than one by, in most cases, several orders of magnitude,
showing that the adjoint-based method is clearly more efficient.
38
k N tAdj/tHDG eHDG/eAdj ωAdj/ωHDG
1
14 6.84 8.63e+01 13
22 5.25 1.49e+02 28
30 4.01 2.23e+02 56
38 3.61 3.00e+02 83
2
14 6.61 2.86e+03 433
22 7.23 6.36e+03 880
30 4.77 1.14e+04 1,333
38 4.38 1.90e+04 4,338
3
14 14.35 5.55e+05 38,676
22 11.57 2.91e+06 251,513
30 8.30 8.79e+06 1,059,036
38 5.98 1.96e+07 3,277,592
Table 2.11: CPU time, error and efficiency comparison of the adjoint-based approxima-
tion JFh (u
∗
h, v
∗
h) and the standard approximation J(uh) for the boundary integral when
u = sin(πx) sin(πy) and to v = π−θπ . The set Ω0 varies with k and h.
Solution u with corner singularity and non-smooth solution v
In this second case, we let the exact solution be a function with a singularity at the
origin, namely, u = r
2
3 sin(23θ) and set ψ = e
xχ[0,1]×[0]. A simple calculation gives us
that
J(u) = ⟨∇u · n, ψ⟩ = −2
3
∫ 1
0
x−
1
3 exdx
This integral is related to the Gamma function and its numerical value can be easily
computed to any precision (for example, using Mathematica).
k N |J(u)− JBh| Order |J(u)− J∗Ih| Order
1
20 1.24e-02 - 8.96e-03 -
30 9.51e-03 0.655 6.86e-03 0.661
40 7.87e-03 0.659 5.67e-03 0.663
50 6.79e-03 0.660 4.89e-03 0.663
2
20 9.12e-03 - 4.32e-03 -
30 6.97e-03 0.661 3.30e-03 0.665
40 5.76e-03 0.663 2.73e-03 0.666
50 4.97e-03 0.664 2.35e-03 0.666
3
20 6.09e-03 - 2.68e-03 -
30 4.65e-03 0.664 2.05e-03 0.666
40 3.84e-03 0.665 1.69e-03 0.666
50 3.31e-03 0.665 1.46e-03 0.666
Table 2.12: History of convergence of the boundary integral functional with u =
r
2
3 sin(23θ) and not smooth v.
39
k N tAdj/tHDG eHDG/eAdj ωAdj/ωHDG
1
20 3.89 1.38 0.35
30 3.49 1.39 0.40
40 3.42 1.39 0.41
50 2.87 1.39 0.48
2
20 4.49 2.11 0.47
30 4.33 2.11 0.49
40 3.28 2.11 0.64
50 2.98 2.11 0.71
3
20 7.58 2.27 0.30
30 7.23 2.27 0.31
40 4.97 2.27 0.46
50 3.68 2.27 0.62
Table 2.13: CPU time, error and efficiency comparison of the adjoint-based approxima-
tion JFh (u
∗
h, v
∗
h) and the standard approximation J(uh) for the boundary integral when
u = r
2
3 sin(23θ) and is not smooth v. The set Ω0 varies with k and h.
In Table 2.12, we see that we have less accuracy due to the lack of smoothness
of the solutions u and v. In Table 2.13, we see that the error of the adjoint-based
approximation is only about 2 times smaller than the standard approximation but the
running time of the adjoint-based approximation is at least 3 times longer than the
standard approximation. Since the values of the ratio of efficiencies ωAdj/ωHDG are less
than one, it is not advantageous to use the adjoint-based approximation. This suggests
the need of incorporating adaptive algorithms to deal with the lack of smoothness of
the solutions u and v.
2.6 Concluding Remarks
By combining the adjoint-correction method [62] and the technique of filtering by
convolution [9], two adjoint-based methods with super-convergent properties were ob-
tained in [27]. In this chapter, we provide the first a priori error analysis for these
methods when the solutions u and v are smooth enough. We also displayed extensive
numerical results showing the advantages of this method over the standard approach
and showed the need of incorporating adaptivity algorithms whenever the solutions u
and v lack sufficient smoothness.
Although we used the HDG method as the Galerkin method in our analysis, similar
proofs can be carried out for other Galerkin methods, like the mixed methods and the
40
continuous Galerkin method, as these methods typically provide a very small negative-
order norm estimate of the error.
The application of the adjoint-based method to other (linear or nonlinear) function-
als of solutions of (linear or nonlinear) partial differential equations, and the incorpora-
tion of adaptivity into the method, are studied in the next a few chapters.
Chapter 3
The adjoint-based method for the
eigenvalue problem
3.1 Introduction
In this chapter, we consider the eigenvalue problem arising from partial differen-
tial equations. We treat the eigenvalue as a non-linear functional of its corresponding
eigenfunction and we present an adjoint-based approach to approximate this non-linear
functional in a high-order accurate manner.
Eigenvalue problems play a fundamental role in many scientific and engineering
fields, such as structural dynamics, electromagnetism and quantum chemistry. For ex-
ample, eigenvalues of the time-independent Schro¨dinger equation describe the energy
levels of a molecular system in quantum chemistry. Another example is the Helmholtz
eigenvalue problem in acoustics. So, to achieve satisfactory accuracy with less compu-
tational cost, highly accurate methods for eigenvalue problems are of great practical
interest.
Among the numerous numerical methods proposed for eigenvalue problems, the
finite element methods constitute an important class. The mixed finite element method
for eigenvalue problems is analyzed in [54, 3] and further expanded upon in [8, 39].
The discontinuous Galerkin eigenvalue approximations are studied in [47, 2, 11]. More
recently, the hybridized mixed method and the HDG method have been considered for
eigenvalue problems in [42, 22]. Therein, the hybridization technique [20] is employed to
41
42
reduce the number of globally-coupled unknowns and render the method more efficient.
Moreover, it was proven that among these finite element methods, the approximate
eigenvalues obtained by the (hybridized) mixed finite element method exhibits the best
rate of convergence of O(h2k+2) if polynomial degree k is used. The rate of convergence
of the HDG eigenvalues is O(h2k+1) but numerical experiments in [42] indicate that,
with a local and inexpensive post-processing technique, the eigenvalues converge at rate
of O(h2k+2).
The method we propose builds upon our previous work [27, 28]. We use a Rayleigh
quotient-like formula and treat the eigenvalue as a non-linear functional of its eigen-
function. To achieve a highly accurate approximation of this non-linear functional, we
utilize two techniques. The first is the adjoint-correction method studied in [62, 41].
Thanks to this approach, we can devise an approximation of the non-linear functional
that has twice the order of accuracy of the numerical functions on which it is based.
The second technique is the accuracy-enhancing convolution proposed by Bramble and
Schatz in [9], see also [70]. This method was first applied to continuous Galerkin solu-
tions for second-order elliptic problems. It was shown that by using a local convolution,
the rate of convergence can be improved from O(hk+1) to O(h2k). This approach was
successfully applied to discontinuous Galerkin methods in [25, 65]. Ryan, Kirby and
their collaborators have further extended this method, see [71, 50, 49, 66, 51] and the
references therein.
Our main result, Theorem 3.1, contains the definition of the so-called adjoint-
corrected approximate eigenvalue and an identity for the difference with the exact
eigenvalue. Our adjoint-based method is an application of this result to the accuracy-
enhanced approximate eigenfunctions. Numerical results in Section 3.4 show that the
approximate eigenvalues by the adjoint-based method converge with a rate of O(h4k+2)
when the eigenfunctions are smooth enough. In addition, we show that, for a given
error tolerance, our method recovers more eigenvalues than using a Galerkin method.
We also apply our main result to the standard HDG method. We find that the
adjoint-corrected approximate eigenvalue converges faster than O(h2k+1). In fact, our
numerical results suggest a 2k + 2 convergence rate for k = 2 and 3. To the best of our
knowledge, this is the first time a convergence rate of 2k + 2 is obtained for the HDG
method without any post-processing, making it competitive with the mixed method.
43
Finally, we carry numerical experiments to explore the possibility of using the
adjoint-correction term as an a posteriori error estimate. The idea of using the adjoint-
correction term for error estimation was first considered by Pierce and Giles [63] in
2004. Therein, the adjoint-correction term was based on a cubic spline interpolation
of finite difference or finite volume solutions. In 2018, Sharbatdar and Ollivier-Gooch
[68] explored this approach for finite volume methods in unstructured meshes. Since
for the continuous Galerkin methods, the adjoint-correction term is zero, to obtain a
posteriori error estimates such as the one considered in [5, 35], an auxiliary approxi-
mation (obtained by a higher-order method or from a finer mesh) is usually needed.
Taking advantage that, for the particular HDG method we use in the experiments, the
adjoint-correction term is non-zero, we find experimentally that the adjoint-correction
term can be an asymptotically exact a posteriori error estimate. However, it does not
hold when the eigenvalues display singularities.
To illustrate the adjoint-based method, we consider a second order elliptic eigenvalue
problem
−∆u = λu in Ω and u = 0 on ∂Ω. (3.1.1)
Although this model problem is simple, the method is general enough to be adapted to
a wide variety of eigenvalue problems with known finite element methods.
The remainder of this chapter is organized as follows. In Section 3.2, we define
the adjoint-based method for the eigenvalue problem (3.1.1). We present and discuss
our main result, Theorem 3.1, and apply it to the standard HDG method and to the
new method we are proposing. Next, in Section 3.3, we present a detailed proof of our
main result. Then in Section 3.4, we provide extensive numerical experiments to show
the advantages of the adjoint-based method. Finally, in 3.5, we conclude with a short
discussion of possible extensions and a description of our forthcoming work.
3.2 The adjoint-based method for eigenvalue problem
This section is devoted to the definition of the adjoint-based method to approximate
eigenvalues.
44
3.2.1 The components of the method
We present the two of the three components of the adjoint-based method for the
eigenvalue problem. We start with the adjoint-correction technique in Section 3.2.1. It
contains the core idea about how we approximate eigenvalues. Next, in Section 3.2.1,
we introduce the HDG method for the eigenvalue problem. As to the third component,
the convolution filtering technique, it is the same as we present in Section 2.2. So we
omit it here.
The adjoint-correction technique
The adjoint-correction technique we consider here was proposed by Pierce and Giles
in [62] for approximating functionals; it was further explored in [41]. In our setting,
this method obtains a more accurate approximation of the eigenvalue λ(u) by adding
to the standard approximation, λh(uh), a carefully devised and computable term called
the adjoint-correction term, ACh. Let us describe the procedure.
First, let us rewrite the equations defining our eigenvalue problem, (3.1.1), in the
following mixed form:
q +∇u = 0 in Ω,
∇ · q = λ u in Ω,
u = 0 on ∂Ω.
(3.2.1)
Let Th be any conforming mesh of the domain Ω. We denote by ∂Th the set of all
boundaries ∂K of the elements K of Th. We denote by F(K) the set of faces F of the
element K ∈ Th and we denote by Fh the set of all faces F . Finally, Let us use the
notation
(u, v)Th :=
∑
K∈Th
∫
K
uv dx and ⟨u , v⟩∂Th :=
∑
K∈Th
∫
∂K
uv ds.
and ∥f∥2L2(Ω) = (f, f)
1/2
Th is the standard L
2 norm.
It is known that the Rayleigh quotient formula associates the eigenfunction u to the
45
corresponding eigenvalue λ(u) as a non-linear functional:
λ(u) =
(∇u,∇u)Ω
∥u∥2L2(Ω)
.
If we now use the second equation of our problem, we readily get that
λ(u) =
−(q,∇u)Ω + ⟨u, q · n⟩∂Ω
∥u∥2L2(Ω)
.
We use a discrete version of this formula to define the approximation to λ(u). For a
general approximation (qh, uh, q̂h ·n, ûh) of the exact solution (q|Th , u|Th , q ·n|∂Th , u|Fh)
of the model problem (3.2.1), we define what we can consider to be the standard ap-
proximation to λ(u) as
λh(uh) :=
−(qh,∇uh)Th + ⟨uh, q̂h · n⟩∂Th
∥uh∥2L2(Ω)
. (3.2.2)
Of course, λh(uh) depends also on qh and q̂h but we do not indicate such dependency
to alleviate the notation. This formula is also motivated by the fact that many finite
element methods, like the hybridized mixed methods, DG or HDG methods approx-
imate the second equation of the model problem (3.2.1) by using the following weak
formulation:
−(qh,∇w)Th + ⟨q̂h · n , w⟩∂Th = λh(uh, w)Th ∀w ∈Wh.
If we take w := uh the above equation becomes nothing but our definition of our
standard approximation, λh(uh) given by (3.2.2).
We are now ready to state our main result. It corresponds to [27, Theorem 2.1] for
the approximation of linear functionals. Its proof is presented in Section 3.3.
Theorem 3.1. Let (qh, uh, q̂h·n, ûh) be an approximation of the exact solution (q|Th , u|Th , q·
n|∂Th , u|Fh) of the model problem (3.2.1). Assume that ûh = 0 on ∂Ω. If we define
λh(uh) as in (3.2.2), we have that
λ(u) = λh(uh) +
ACh
∥uh∥2L2(Ω)
+
Eh(uh, u)
∥uh∥2L2(Ω)
,
46
where the adjoint-correction term ACh and error term Eh(uh, u) are given by
ACh(uh) :=− (qh, qh +∇uh)Th + ⟨uh − ûh, qh · n⟩∂Th
− ⟨ûh, q̂h · n⟩∂Th
−⟨uh − ûh, (q̂h − qh) · n⟩∂Th ,
Eh(uh, u) := λ(u) ∥uh − u∥2L2(Th) − ∥qh − q∥2L2(Th)
+ 2 (qh − q,∇uh + qh)Th − 2 ⟨uh − ûh, (qh − q) · n⟩.
Let us briefly discuss this result.
First, note that the identity for λ(u) in Theorem 3.1 holds for general approximations
(uh, qh, ûh, q̂h ·n). It is not limited to specific finite element approximations. It is even
valid for finite difference or finite volume methods, as these approximations can be first
created by interpolation through computed values at grid nodes.
A detailed proof of Theorem 3.1 can be found in Section 3.3, but we give here the
rough idea of the derivation of the adjoint-correction formula. Let us assume we have a
Rayleigh-quotient like formula R(·, ·) and we define
λh(uh) =
R(uh, uh)
∥uh∥L2(Ω)
as an approximation of λ(u). Let us also assume the eigenvalue problem can be written
in a variational form A(u, v) = λ(u, v) for v ∈ V and uh ∈ V , where V is some function
space. Then we have
(λ− λh)∥uh∥2L2(Ω) = λ∥uh∥2L2(Ω) −R(uh, uh)
= λ∥uh − u∥2L2(Ω) + 2λ(u, uh)− λ∥u∥2L2(Ω) −R(uh, uh)
= λ∥uh − u∥2L2(Ω) + 2A(u, uh)−A(u, u)−R(uh, uh)
= λ∥uh − u∥2L2(Ω) −A(u− uh, u− uh) +A(uh, uh)−R(uh, uh).
From here, we can gather the computable terms as the correction term.
Note also that the definition of ACh(uh) is computable and closely related to the
so-called Galerkin orthogonality property. If the approximation (uh, qhq̂h · n, ûh) is
obtained from Galerkin methods, the adjoint-correction term ACh(uh) is usually zero.
47
For example, for the continuous Galerkin method, qh = −∇uh, ûh = uh, and ACh(uh) =
0 if the normal component of the numerical flux q̂h is chosen to be single valued.
Moreover, the last two terms of the error Eh(uh, u) vanish and the identity of Theorem
3.1 becomes
λ(u)− λh(uh) =
λ(u) ∥uh − u∥2L2(Ω) − ∥∇uh −∇u∥2L2(Ω)
∥uh∥2L2(Ω)
This is why for the continuous Galerkin method the order of convergence of the error
λ(u)− λ(uh) is twice that of u− uh in the H1(Ω)-norm. We see then that Theorem 3.1
is an extension of the classic result for the continuous Galerkin method.
However, this is not the case for approximations not obtained by means of Galerkin
methods. As Pierce and Giles intended [62], to achieve the same order of convergence
typical of Galerkin methods, a new approximation to λ, called the adjoint-corrected
approximate eigenvalue, is defined as
λACh (uh) := λh(uh) +
ACh(uh)
∥uh∥2L2(Ω)
. (3.2.3)
Finally, note that since the model problem we consider here (3.1.1) is self-adjoint,
the adjoint eigenvalue problem is itself. So we do not actually solve an adjoint problem
like we did in [27, 28]. We use the term “adjoint-based” mainly because of the adjoint-
correction technique we used here. For a non self-adjoint problem, an adjoint equation
needs to be solved. See the discussion in Section 3.5.
The HDG method
In Section 2.2, we have introduced the HDG method for source problems . The
HDG method for a second-order elliptic eigenvalue problem is first studied in [42].
Although what we are presenting can be applied to any Galerkin method, we use the
HDG method here for several reasons. First, it has a general structure that allows us to
easily extend the results to a large class of Galerkin methods, including the mixed finite
element method and the continuous Galerkin method. Second, thanks to the static
condensation technique, the HDG method can reduce the system size by removing all
interior variables to yield a “condensed” system. Third, the HDG method has a flexible
48
stabilization mechanism and it can be easily adapted to different partial differential
equations.
The HDG method seeks an approximation, (λh, qh, uh, q̂h · n, ûh), of the exact so-
lution, (λ, q|Ω, u|Ω, q · n|∂Ω, u|Fh), of the model problem (3.2.1). Here λh ∈ R is the
discrete eigenvalue and (qh, uh, ûh) lies in the finite dimensional space V h ×Wh ×Mh,
where
V h := {v ∈ L2(Th) :v|K ∈ V (K) ∀K ∈ Th},
Wh := {w ∈ L2(Th) :w|K ∈W(K) ∀K ∈ Th},
Mh := {µ ∈ L2(Fh) :µ|F ∈M(F ) ∀F ∈ Fh}.
The HDG approximation is determined as the solution of the following weak formu-
lation:
(qh, r)Th − (uh,∇ · r)Th + ⟨ûh , r · n⟩∂Th = 0, (3.2.4a)
−(qh,∇w)Th + ⟨q̂h · n , w⟩∂Th = λh(uh, w)Th , (3.2.4b)
⟨ûh , µ⟩∂Ω = 0, (3.2.4c)
⟨q̂h · n , µ⟩∂Th\∂Ω = 0, (3.2.4d)
for all (r, w, µ) ∈ V h ×Wh ×Mh, where
q̂h · n := qh · n+ τ(uh − ûh) on ∂Th. (3.2.4e)
The different HDG methods are obtained by choosing different local spaces V (K),
W (K),
M(∂K) := {µ ∈ L2(∂K) : µ|F ∈M(F ) ∀F ∈ F(K)},
and stabilization functions τ . The hybridized version of mixed methods is obtained by
simply taking τ equal to zero, as pointed out in [21].
In [42], the HDG method defined on symplices and using polynomials of degree
k ≥ 0 was considered. It was shown that there is no pollution of the spectrum when
approximating (3.2.1). Moreover, under suitable regularity assumptions, the approxi-
mate eigenvalues and eigenfunctions were proven to converge at the rate of O(h2k+1)
49
and O(hk+1), respectively. In addition, it was shown that, static condensation can be
used to significantly reduce the size of the matrix system and it does not lose eigenvalues
up to size O(1/h). In other words, with static condensation, the physically important
lower range eigenvalues are preserved. See [42, Section 5] for details; in particular, see
the discussion after Theorem 5.3 therein.
Note that the definition of the HDG eigenvalue λh, (3.2.4), is consistent with the
definition of λh(uh) in (3.2.2), as the latter formula is obtained by taking w := uh in
(3.2.4b). We also note that, for the HDG method, the first two terms defining the
adjoint-correction term in Theorem 3.1 vanish since it is the first equation of the HDG
scheme (3.2.4a) by taking r := qh and integrating by parts. The third term of its
definition is also zero by the so-called transmission condition (3.2.4d) with µ := ûh.
Therefore, for the HDG method, the adjoint-correction term becomes
ACh = −⟨uh − ûh, (q̂h − qh) · n⟩∂Th ,
a non-zero term. This term contains the information about the stabilization of the HDG
method and it captures the contribution of the lack of conformity of the spaces.
To end, we gather some properties of general HDG methods in the following result.
We prove it in Section 3.3.
Proposition 3.1.1 (Properties of the HDG eigenvalues). We have that
(i) λh(uh) =
∥qh∥2L2(Th) + ⟨(uh − ûh), τ(uh − ûh)⟩∂Th
∥uh∥2L2(Ω)
≥ 0,
(ii) λACh (uh) =
∥qh∥2L2(Th)
∥uh∥2L2(Ω)
≤ λh(uh),
(iii) λ(u) = λACh (uh) +
Eh(uh, u)
∥uh∥2L2(Ω)
,
where
Eh(uh, u) = λ(u) ∥uh − u∥2L2(Th) − ∥qh − q∥2L2(Th)
− 2 ⟨uh − ûh, (P h(q)− q) · n⟩∂Th ,
50
where, on each element K ∈ Th, P h is the L2(K)−projection into V (K) +∇W (K).
Our numerical results show that, when we use HDG methods with polynomials of
degree k ≥ 1, λh(uh) converges with order h2k+1, in full agreement with the theoretical
results in [42]. They also suggest that λACh (uh) converges with an additional order. (For
another definition of an approximate HDG eigenvalue also converging with order h2k+2,
see [42].) This indicates that the ratio
∥uh∥2L2(Ω)|λ(u)− λh(uh)|/|ACh(uh)|,
remains extremely close to one, which we verify in our numerical experiments. This
shows that
ACh(uh)/∥uh∥2L2(Ω),
is a good candidate for an a posteriori error estimate for λ(u)− λh(uh).
3.2.2 The definition of the method
Let us define the adjoint-based method. We begin by describing our assumption on
the meshes.
Assumption 5. Let {Th}h>0 be a family of shape-regular meshes of Ω. There exist
subdomains Ω0 ⊂⊂ Ω′0 ⊂ Ω such that
(1) Th ∩ Ω′0 is a translation-invariant mesh with element size h,
(2) Any element K ∈ Th is fully contained in either Ω0 or Ω1 := Ω\Ω0
We require Th ∩ Ω′0 to be translation invariant so that we can use the accuracy-
enhancing convolution technique in the subdomain Ω0. Let us denote by T0,h := Th∩Ω0.
For the subdomain Ω1, we can simply take T1,h := Th ∩ Ω1 or we can define a different
mesh with size h on Ω1 as illustrated in Figure 3.1c.
51
(a) The HDG solution with
polynomial degree k in Ω
(b) The accuracy-enhanced so-
lution in Ω0
(c) The HDG solution with
polynomial degree 2k in Ω1
Figure 3.1: An illustration of the definition of u∗h in the whole domain Ω. Here Ω0 is
the subdomain that is two elements away from the boundary ∂Ω. This is how we have
to define Ω0 to be able to carry out the postprocessing with k = 1.
We are now ready to define the adjoint-based method in Algorithm 1 . See also the
schematic illustration in Figure 3.1. To compute the HDG eigenvalues in step (1) of
Algorithm 1, we follow the procedure proposed in [42]. To control the iterations of the
computation of λACh (u
∗
h), we use M as the maximum iteration steps and use an error
tolerance denoted by ϵ. In our numerical experiments, we take M = 50 and ϵ = 10−10
and we observe that λACh usually takes less than 5 iterations to converge.
52
Algorithm 1: The computation of one approximate eigenvalue λACh (u
∗
h)
Initialization
(1) compute, on the whole mesh Th, the approximate eigenfunctions,
(qkh, u
k
h, q̂
k
h · n, ûkh), of problem (3.2.1) by using the HDG method with local
spaces V (K), W (K) and M(F ) including the polynomials of degree k ≥ 1 for
all K ∈ Th and F ∈ Fh. (See Figure 3.1a.)
(2) compute u∗h := K
2k
h ∗ uh in Ω0 by applying the accuracy-enhancing
convolution technique introduced in Section 2.2. (See Figure 3.1b.)
(3) compute the initial eigenvalue approximations
λ0 :=
(∇u∗h,∇u∗h)T0,h − ⟨u∗h,∇u∗h · n⟩∂Ω0
∥u∗h∥2L2(Ω0)
Iterations
(4) for i=1:M do
(i) compute, on the mesh T1,h, the approximate eigenfunctions,
(q2kh , u
2k
h , q̂
2k
h · n, û2kh ), of the boundary value problem
q +∇u = 0 Ω1,
∇ · q − λ0 u = 0 Ω1,
u = u∗h ∂Ω1\∂Ω,
u = 0 ∂Ω,
(3.2.5)
by using the HDG method with local spaces V ∗(K), W ∗(K) and M∗(F )
including the polynomials of degree 2k for all K ∈ T1,h and F ∈ F1h. (See
Figure 3.1c.)
(ii) compute u∗h and q
∗
h in the whole domain Ω and û
∗
h and q̂
∗
h on
∂T0,h ∪ ∂T1,h as follows
u∗h :=
K2kh ∗ ukh in Ω0,u2kh in Ω1, q∗h :=
−∇u∗h in Ω0,q2kh in Ω1, (3.2.6a)
û∗h :=
u∗h on ∂T0,h,û2kh on ∂T1,h, q̂∗h :=
−∇u∗h on ∂T0,h,q̂2kh on ∂T1,h. (3.2.6b)
(iii) compute
λACh := λh(u
∗
h) +
ACh(u
∗
h)
∥u∗h∥2L2(Ω)
. (3.2.7)
as defined in (3.2.3).
(iv) Check
if |λ0 − λACh | > ϵ then
λ0 := λ
AC
h
else
Stop
end if
end for
(5) set λACh (u
∗
h) := λ
AC
h and stop.
53
To end, we gather some identities for the approximate eigenvalues obtained by the
adjoint-based method just described. We can prove them exactly as the identities of
Proposition 3.1.1. We only have to recall the definition of (u∗h, q
∗
h, û
∗
h, q̂
∗
h · n), (3.2.6),
and in particular, that q∗h = −∇u∗h and that u∗h lies in C1(Ω0).
Proposition 3.1.2 (Properties of the adjoint-based HDG eigenvalues). We have that
(i) λh(u
∗
h) =
(q∗h, q
∗
h)Th + ⟨(u∗h − û∗h), τ(u∗h − û∗h)⟩∂T1,h
∥u∗h∥2L2(Ω)
,
(ii) λACh (u
∗
h) =
(q∗h, q
∗
h)Th
∥u∗h∥2L2(Ω)
≤ λh(u∗h),
(iii) λ(u) = λACh (u
∗
h) +
Eh(u
∗
h, u)
∥uh∥2L2(Ω)
,
where
Eh(u
∗
h, u) = λ(u) ∥u∗h − u∥2L2(Ω) − ∥q∗h − q∥2L2(Ω)
− 2 ⟨u∗h − û∗h, (P ∗h(q)− q) · n⟩∂T1,h ,
where, on each element K ∈ Th, P ∗h is the L2(K)−projection into V ∗(K) +∇W ∗(K).
Our numerical results show that λACh (u
∗
h) converges to λ(u) with order h
4k+2 and
the error of λACh (u
∗
h) is several orders of magnitude smaller than that of λh(u
∗
h). This
indicates that the ratio
∥u∗h∥2L2(Ω)|λ(u)− λh(u∗h)|/|ACh(u∗h)|,
should remain extremely close to one, which we also verify in our numerical experiments.
In other words, the expression
ACh(u
∗
h)/∥u∗h∥2L2(Ω),
is a good candidate for an a posteriori error estimate for λ(u)− λh(u∗h).
54
3.3 Proofs
This section is devoted to providing detailed proofs of the expressions and corre-
sponding approximation errors for Theorem 3.1 and Proposition 3.1.1 .
3.3.1 Proof of Theorem 3.1
In this proof, to alleviate the notation, we write λ instead of λ(u), and λh instead of
λh(uh). Also, since there is no ambiguity, for the sake of clarity, we drop the subscripts
Th and ∂Th in (·, ·)Th and ⟨· , ·⟩∂Th .
The result follows if we show that ∥uh∥2L2(Ω)(λ − λh) = Eh(uh, u) + ACh. By the
definition of λh, we have that
∥uh∥2L2(Ω)(λ− λh) = λ(uh, uh) + (qh,∇uh)− ⟨uh, q̂h · n⟩
= λ(uh − u, uh − u) + 2λ (u, uh)− λ (u, u)
+ (qh,∇uh)− ⟨uh, q̂h · n⟩.
Since u is the solution of (3.1.1), λ(u, uh) = −(q,∇uh)+⟨uh , q ·n⟩ and λ(u, u) = (q, q).
This implies that
∥uh∥2L2(Ω)(λ− λh) = λ(uh − u, uh − u) + 2 (−(q,∇uh) + ⟨uh , q · n⟩)− (q, q)
+ (qh,∇uh)− ⟨uh, q̂h · n⟩.
Now, we only have to rearrange terms and carry out a few algebraic manipulations to
55
obtain the desired result. We have
∥uh∥2L2(Ω)(λ− λh) = λ(uh − u, uh − u)− (q, q)
− 2 (q,∇uh) + (qh,∇uh)
+ 2 ⟨uh , q · n⟩ − ⟨uh, q̂h · n⟩
= λ(uh − u, uh − u)− (q − qh, q − qh)− 2 (q, qh) + (qh, qh)
− 2 (q,∇uh) + (qh,∇uh)
+ 2 ⟨uh , q · n⟩ − ⟨uh, q̂h · n⟩
= λ(uh − u, uh − u)− (q − qh, q − qh)
+ 2 (qh − q,∇uh + qh)
− (qh,∇uh + qh) + 2 ⟨uh , q · n⟩ − ⟨uh, q̂h · n⟩.
= Eh(uh, u) + Θh
where, by the definition of Eh(uh, u), we have that
Θh = 2 ⟨uh − ûh , (qh − q) · n⟩ − (qh,∇uh + qh) + 2 ⟨uh , q · n⟩ − ⟨uh, q̂h · n⟩.
It only remains to show that Θh = ACh. But
Θh = − (qh,∇uh + qh) + 2 ⟨uh − ûh , qh · n⟩ − ⟨uh, q̂h · n⟩
because ⟨ûh , q·n⟩ = 0 due to the single-valuedness of ûh and the fact that, by hypothesis,
ûh is zero on the boundary. Finally, we get that
Θh = − (qh,∇uh + qh) + ⟨uh − ûh , qh · n⟩
− ⟨ûh, q̂h · n⟩
+ ⟨uh − ûh , (qh − q̂h) · n⟩
=ACh,
as wanted. This completes the proof of Theorem 3.1.
56
3.3.2 Proof of Proposition 3.1.1
The first identity can be proved as in the proof of [42, Lemma 2.2], but we give a
proof here for the sake of competeness. From the definition of λh(uh), (3.2.2),
∥uh∥2L2(Ω) λh(uh) = −(qh,∇uh)Th + ⟨uh, q̂h · n⟩∂Th
= (∇ · qh, uh)Th + ⟨uh, (q̂h − qh) · n⟩∂Th
= (qh, qh)Th + ⟨ûh, qh · n⟩∂Th + ⟨uh, (q̂h − qh) · n⟩∂Th ,
by the first equation defining the HDG method with r := qh. Completing squares, we
get
∥uh∥2L2(Ω) λh(uh) = (qh, qh)Th + ⟨(uh − ûh), (q̂h − qh) · n⟩∂Th + ⟨ûh, q̂h · n⟩∂Th ,
by the third and fourth equations defining the HDG method with µ := ûh. Inserting
the formula of the numerical trace q̂h, (3.2.4e), we get
λh(uh) =
(qh, qh)Th + ⟨(uh − ûh), τ(uh − ûh)⟩∂Th
∥uh∥2L2(Ω)
≥ 0,
Let us prove the second identity. We have
ACh =− (qh, qh +∇uh)Th + ⟨uh − ûh, qh · n⟩∂Th
− ⟨ûh, q̂h · n⟩∂Th
− ⟨uh − ûh, (q̂h − qh) · n⟩∂Th
=− ⟨ûh, q̂h · n⟩∂Th
− ⟨uh − ûh, (q̂h − qh) · n⟩∂Th
by the first equation defining the HDG method with r := qh. By the third and fourth
equations defining the HDG method,
ACh =− ⟨uh − ûh, (q̂h − qh) · n⟩∂Th = −⟨uh − ûh, τ(uh − ûh) · n⟩∂Th ,
57
by the definition of the numerical trace of the flux q̂h, (3.2.4e). This implies that
λACh (uh) =
(qh, qh)Th
∥uh∥2L2(Ω)
≤ λh(uh).
It remains to prove the third identity. By Theorem 3.1, we have that
λ(u) = λACh (uh) + Eh(uh, u),
where
Eh(uh, u) = λ(u) (uh − u, uh − u)Th−(qh − q, qh − q)Th
+ 2 (qh − q,∇uh + qh)Th − 2 ⟨uh − ûh, (qh − q) · n⟩
= λ(u) (uh − u, uh − u)Th−(qh − q, qh − q)Th
+ 2 (qh − P h(q),∇uh + qh)Th − 2 ⟨uh − ûh, (qh − q) · n⟩
and the result follows by using the first equation defining the HDG method with r :=
qh − P h(q).
This completes the proof of Proposition (3.1.1).
3.4 Numerical results
In this section, we present numerical experiments for the adjoint-based method. We
implemented all the experiments in MATLAB R2019a in a personal laptop (MacBook
Pro 2019). To show the high-order accuracy of the method, when needed, we used
the Multi-precision Computing Toolbox [1] and select 64 digits of accuracy. When the
accuracy-enhancing convolution technique is used, we choose Ω0 to be the region that
is 2k elements away from ∂Ω to make sure the support of the kernel K2kh is contained
in Ω, where k the polynomial degree used in the HDG method.
Our numerical experiments consist of two parts. In the first part, we compare the
eigenvalues obtained by the adjoint-based method and the standard HDG method in
Section 3.4.1 and 3.4.2. We show the efficiency of the adjoint-based method in Section
3.4.3 and we present one application in quantum chemistry in Section 3.4.4. In the
58
second part, we show the necessity of including the adjoint-correction term for the
adjoint-based method and we explore the possibility of using the adjoint-correction
term as an posteriori error estimate in Section 3.4.5.
3.4.1 The one dimensional Laplace eigenproblem
Let us first consider the model problem (3.1.1) on Ω := (0, 1). It is known that the
exact eigenvalues and eigenfunctions are λm = π
2m2 and um = sin(mπx) for m ∈ N+.
We report the history of the convergence of the first five eigenvalues obtained by
the HDG method and the adjoint-based method in Table 3.1. For the HDG method,
the eigenvalues λh(uh) converge with order h
2k+1, in full agreement with the theoretical
results in [42]. As for the adjoint-based method, we observe that the eigenvalues λACh (u
∗
h)
converge with order h4k+2 for k = 1, 2, 3. For the first three eigenvalues {λi}3i=1, a
logarithmic scale error comparison between the two methods is shown in Figure 3.2
when k = 1 . It clearly shows that the adjoint-based method converges twice as fast as
the HDG method.
59
First Second Third Fourth Fifth
k 1/h error order error order error order error order error order
Approximation given by λh(uh)
1
128 5.79e-06 – 4.18e-04 – 5.19e-03 – 3.26e-02 – 1.46e-01 –
256 7.20e-07 3.01 5.09e-05 3.03 6.09e-04 3.09 3.61e-03 3.17 1.47e-02 3.30
512 8.97e-08 3.00 6.29e-06 3.01 7.39e-05 3.04 4.27e-04 3.08 1.68e-03 3.13
1024 1.12e-08 3.00 7.81e-07 3.00 9.09e-06 3.02 5.19e-05 3.03 2.01e-04 3.06
2
128 3.47e-11 – 9.87e-09 – 2.67e-07 – 2.85e-06 – 1.85e-05 –
256 1.08e-12 5.00 3.04e-10 5.02 8.05e-09 5.05 8.32e-08 5.10 5.16e-07 5.17
512 3.37e-14 5.00 9.42e-12 5.01 2.47e-10 5.02 2.52e-09 5.05 1.53e-08 5.07
1024 1.05e-15 5.00 2.93e-13 5.00 7.66e-12 5.01 7.74e-11 5.02 4.66e-10 5.04
3
128 1.06e-16 – 1.20e-13 – 7.24e-12 – 1.34e-10 – 1.33e-09 –
256 5.04e-38 7.00 9.29e-16 7.01 5.51e-14 7.04 1.00e-12 7.07 9.62e-12 7.11
512 6.47e-21 7.00 7.22e-18 7.01 4.25e-16 7.02 7.66e-15 7.03 7.25e-14 7.05
1024 5.05e-23 7.00 5.62e-20 7.00 3.30e-18 7.01 5.92e-17 7.02 5.56e-16 7.03
Approximation given by λACh (u
∗
h)
1
128 8.64e-13 – 8.97e-10 – 5.88e-08 – 1.29e-06 – 1.64e-05 –
256 1.37e-14 5.97 1.37e-11 6.03 8.36e-10 6.14 1.63e-08 6.30 1.73e-07 6.56
512 2.16e-16 5.99 2.12e-13 6.02 1.25e-11 6.06 2.31e-10 6.14 2.29e-09 6.24
1024 3.39e-18 5.99 3.30e-15 6.01 1.91e-13 6.03 3.45e-12 6.07 3.31e-11 6.11
2
128 6.08E-23 – 1.02e-18 – 3.21e-16 – 2.04e-14 – 5.48e-13 –
256 6.30e-26 9.91 1.03e-21 9.95 3.12e-19 10.01 1.85e-17 10.10 4.53e-16 10.24
512 6.33e-29 9.96 1.02e-24 9.98 3.03e-22 10.00 1.75e-20 10.05 4.12e-19 10.10
1024 6.27e-32 9.98 1.01e-27 9.99 2.96e-25 10.00 1.68e-23 10.02 3.89e-22 10.05
3
128 7.47e-34 – 2.00e-28 – 3.15e-25 – 6.13e-23 – 3.85e-21 –
256 5.04e-38 13.86 1.33e-32 13.88 2.01e-29 13.94 3.71e-27 14.01 2.18e-25 14.11
512 3.22e-42 13.93 8.41e-37 13.94 1.26e-33 13.97 2.27e-31 14.00 1.29e-29 14.01
1024 2.01e-46 13.97 5.23e-41 13.97 7.76e-38 13.98 1.39e-35 14.00 7.81e-34 14.02
Table 3.1: History of convergence of the first five eigenvalues provided by the HDG
method, λh(uh), (top) and by the adjoint-based method, λ
AC
h (u
∗
h), (bottom) for the
problem on (0, 1).
60
2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3 3.1
log(1/h)
-18
-16
-14
-12
-10
-8
-6
-4
-2
lo
g(e
rro
r)
E1
E2
E3
E1*
E2*
E3*
1
0
3
1
0
6
Figure 3.2: The errors of the first three eigenvalues in logarithmic scale for k = 1 with
N = 27, 28, 29, 210 elements on (0, 1). Here Ei = log10(|λi − λi,h(uh)|) (dashed line) and
E∗i = log10(|λi − λACi,h (u∗h)|) (solid line) for i = 1, 2, 3.
3.4.2 The two-dimensional Laplace eigenproblem
In this subsection, we consider the model problem (3.1.1) in two dimensions. Nu-
merical results for a square domain and an L-shaped domain are presented. In the
implementation, we use rectangular meshes and choose the local spaces
V (K) :=Qk(K)⊕ curlspan{xk+1y, xyk+1},W (K) := Qk(K),M(∂K) := Qk(∂K),
for the HDG method with polynomial degree k. Here Qk denotes the space of tensor
product polynomials of degree at most k in each variable. The local spaces chosen here
admit an M-decomposition, see [19].
61
First Second Third Fourth Fifth
k 1/h error order error order error order error order error order
Approximation given by λh(uh)
1
20 3.44e-04 – 3.16e-03 – 3.16e-03 – 3.81e-03 – 1.70e-02 –
30 1.02e-04 3.00 9.25e-04 3.03 9.25e-04 3.03 1.50e-03 2.29 4.98e-03 3.03
40 4.29e-05 3.01 3.85e-04 3.05 3.85e-04 3.05 6.73e-04 2.80 2.04e-03 3.10
50 2.19e-05 3.01 1.95e-04 3.05 1.95e-04 3.05 3.51e-04 2.92 1.02e-03 3.11
2
c 20 8.58e-08 – 3.06e-06 – 3.06e-06 – 6.00e-06 – 3.95e-05 –
30 1.12e-08 5.02 3.88e-07 5.09 3.88e-07 5.09 7.64e-07 5.08 4.82e-06 5.19
40 2.65e-09 5.02 9.03e-08 5.07 9.03e-08 5.07 1.78e-07 5.07 1.10e-06 5.14
50 8.65e-10 5.01 2.92e-08 5.05 2.92e-08 5.05 5.76e-08 5.05 3.50e-07 5.12
3
20 1.07e-11 – 1.49e-09 – 1.49e-09 – 2.95e-09 – 4.24e-08 –
30 6.24e-13 7.02 8.44e-11 7.07 8.44e-11 7.07 1.68e-10 7.06 2.34e-09 7.14
40 8.30e-14 7.01 1.11e-11 7.05 1.11e-11 7.05 2.21e-11 7.05 3.01e-10 7.12
50 1.74e-14 7.01 2.30e-12 7.04 2.30e-12 7.04 4.59e-12 7.04 6.19e-11 7.09
Approximation given by λACh (u
∗
h)
1
20 1.02e-08 – 1.30e-06 – 1.30e-06 – 2.48e-06 – 3.11e-05 –
30 8.83e-10 6.04 1.15e-07 5.97 1.15e-07 5.97 2.21e-07 5.97 2.89e-06 5.86
40 1.56e-10 6.03 2.04e-08 6.01 2.04e-08 6.01 3.94e-08 6.00 5.19e-07 5.97
50 4.07e-11 6.02 5.33e-09 6.02 5.33e-09 6.02 1.03e-08 6.00 1.36e-07 6.00
2
20 4.24e-16 – 6.64e-13 – 6.64e-13 – 1.06e-12 – 6.62e-11 –
30 8.25e-18 9.71 1.56e-14 9.26 1.56e-14 9.26 2.88e-14 8.89 1.80e-12 8.90
40 4.77e-19 9.91 9.44e-16 9.74 9.44e-16 9.74 1.82e-15 9.60 1.16e-13 9.53
50 5.16e-20 9.96 1.04e-16 9.88 1.04e-16 9.88 2.04e-16 9.81 1.31e-14 9.77
3
20 3.16e-24 – 4.49e-20 – 4.49e-20 – 4.25e-20 – 3.34e-17 –
30 1.73e-26 12.85 4.34e-22 11.44 4.34e-22 11.44 6.80e-22 10.20 2.20e-19 12.39
40 3.42e-28 13.63 9.98e-24 13.12 9.98e-24 13.12 1.79e-23 12.65 5.63e-21 12.74
50 1.55e-29 13.86 4.82e-25 13.59 4.82e-25 13.59 9.08e-25 13.36 2.90e-22 13.30
60 1.22e-30 13.94 3.91e-26 13.78 3.91e-26 13.78 7.55e-26 13.64 2.43e-23 13.59
70 1.41e-31 13.98 4.61e-27 13.86 4.61e-27 13.86 9.02e-27 13.78 2.92e-24 13.74
Table 3.2: History of convergence of the first five eigenvalues provided by the HDG
method, λh(uh), (top) and by the adjoint-based method, λ
AC
h (u
∗
h), (bottom) for the
problem on the unit square domain.
62
First Fourth
∥u1 − u1,h∥L2(Ω) ∥u1 − u∗1,h∥L2(Ω) ∥u4 − u4,h∥L2(Ω) ∥u4 − u∗4,h∥L2(Ω)
k 1/h error order error order error order error order
1
10 5.43e-03 – 1.19e-04 – 6.47e-02 – 9.27e-04 –
20 1.31e-03 2.05 1.09e-05 3.45 6.55e-03 3.30 8.66e-05 3.42
30 5.80e-04 2.01 2.64e-06 3.49 2.53e-03 2.35 2.11e-05 3.48
40 3.26e-04 2.01 9.66e-07 3.50 1.36e-03 2.15 7.74e-06 3.49
50 2.08e-04 2.00 4.42e-07 3.50 8.56e-04 2.08 3.55e-06 3.50
2
10 1.38e-04 – 4.35e-08 – 1.11e-03 – 1.37e-06 –
20 1.73e-05 3.00 1.17e-09 5.22 1.38e-04 3.01 3.61e-08 5.24
30 5.11e-06 3.00 1.29e-10 5.44 4.09e-05 3.00 4.08e-09 5.38
40 2.16e-06 3.00 2.66e-11 5.48 1.73e-05 3.00 8.52e-10 5.45
50 1.10e-06 3.00 7.82e-12 5.49 8.84e-06 3.00 2.51e-10 5.48
3
20 1.71e-07 – 5.05e-14 – 2.73e-06 – 6.23e-12 –
25 7.00e-08 4.00 9.96e-15 7.28 1.12e-06 4.00 1.26e-12 7.15
30 3.37e-08 4.00 2.59e-15 7.38 5.40e-07 4.00 3.37e-13 7.25
35 1.82e-08 4.00 8.23e-16 7.44 2.91e-07 4.00 1.08e-13 7.35
Table 3.3: History of convergence of the approximate (uh) and accuracy-enhanced (u
∗
h)
eigenfunctions on the unit square domain.
A square domain
We consider the unit square domain Ω = (0, 1) × (0, 1). The exact eigenvalues are
λmn = (m
2 + n2)π2 and the exact eigenfunctions are umn = sin(mπx) sin(nπy) for
m,n ∈ N+. We list the numerical results of the first five eigenvalues in Table 3.2. As
we can see, the approximate eigenvalues λh(uh) obtained by the HDG method converge
with a rate of O(h2k+1), while the eigenvalues λACh (u∗h) provided by the adjoint-based
method converge with a rate of O(h4k+2).
We present the error and order of convergence of the HDG approximate and the
accuracy-enhanced eigenfunctions in Table 3.3. We observe that the HDG approximate
eigenfunctions uh converge with a rate of k+1, in agreement with the theoretical results
in [42]. The accuracy-enhanced approximate eigenfunctions u∗h converge with a rate
more than 2k + 1, in line with our previous studies in [28, 27].The extra observed rate
of convergence may result from a cancellation effect of the sin function. We also display
the plots of the approximations of the first eigenfunction and the errors in Figure 3.3.
In the first two rows, we can see that the HDG approximation uh is discontinuous while
the approximation u∗h obtained by the accuracy-enhancing technique is much smoother.
In the third row, note that the HDG error u−uh oscillates around zero in each element
63
while the error u−u∗h of the accuracy-enhanced approximation does not have oscillations
in the region Ω0 and that its magnitude is smaller than that of the errors of the HDG
method.
64
Figure 3.3: Approximations to the first eigenfunction on a unit square domain with
k = 1 and h = 0.1: The HDG approximation uh (left column) and the accuracy-
enhanced approximation u∗h (right column.). The errors are plotted at the bottom row.
65
An L-shaped domain
We consider an L-shaped domain Ω = Λ0\Λ1, where Λ0 := (0, 2)× (0, 2) and Λ1 :=
(1, 2)× (1, 2). Since Ω has a reentrant corner at the point (1, 1), the first eigenfunction
is singular and in [7] the first eigenvalue is calculated to 14 digits of accuracy as λ1 =
9.6397238440219. It is interesting to note that the third eigenfunction is actually smooth
and the corresponding exact eigenvalue is λ3 = 2π
2. Therefore, the L-shaped domain
has both singular and smooth eigenfunctions. It provides an interesting example to
study the changes in convergence rates due to the presence of singularities.
First Third
|λ1 − λ1,h(uh)| |λ1 − λAC1,h (u∗h)| |λ3 − λ3,h(uh)| |λ3 − λAC3,h (u∗h)|
k 1/h error order error order error order error order
1
10 1.38e-02 – 7.61e-03 – 2.65e-03 – 3.29e-06 –
20 5.90e-03 1.23 3.10e-03 1.29 3.44e-04 2.95 5.37e-08 5.94
30 3.50e-03 1.29 1.82e-03 1.31 1.02e-04 3.00 4.74e-09 5.99
40 2.40e-03 1.31 1.29e-03 1.32 4.29e-05 3.01 8.46e-10 5.99
2
10 7.86e-03 – 1.90e-03 – 2.81e-06 – 8.88e-13 –
20 3.16e-03 1.31 7.60e-04 1.32 8.58e-08 5.03 1.26e-15 9.46
30 1.85e-03 1.32 4.44e-04 1.33 1.12e-08 5.02 2.29e-17 9.88
40 1.26e-03 1.33 3.03e-04 1.33 2.65e-09 5.01 1.30e-18 9.96
Table 3.4: History of convergence for the L-shaped domain problem.
66
Figure 3.4: Approximations of the first eigenfunction in the L-shaped domain with
k = 1 and h = 0.1: The HDG approximation uh (left column) and the accuracy-
enhanced approximation u∗h (right column.). The errors are plotted at the bottom row.
They are concentrated at the corner singularity.
67
We report the errors and convergence of the first and third eigenvalues in Table 3.4.
We observe that for the first eigenvalue, both methods have a convergence rate of at
most O(h4/3). This is in line with the results in [22, 42]. Although no convergence
rate increase is observed for the adjoint-based method, we do see that the errors are
slightly reduced and the approximate first eigenfunction is improved visually even in this
singular case as shown in the first two rows of Figure 3.4. Since the exact eigenfunction
is unknown, to illustrate the errors of the two methods, we compute a HDG solution on
a 10 times finer mesh h = 0.01 with k = 2 and use it as the exact solution. The plots
are presented in the third row of Figure 3.4. We can clearly see that, compared with the
error of the HDG approximation, the error of the accuracy-enhanced approximation is
reduced in most parts of the domain and is concentrated near the singularity. For the
third eigenvalue, since the eigenfunction is smooth, we obtain expected results that the
HDG method converges at the rate of O(h2k+1) and the adjoint-based method converges
at the rate of O(h4k+2).
3.4.3 The efficiency of the adjoint-based method
Here, we show the efficiency of the adjoint-based method. We consider the eigenvalue
problem on a unit square as in Section 3.4.2. For a given mesh Th and error tolerance
ϵ, we count the number of eigenvalues that have relative errors |λ− λh|/λ less than ϵ.
As shown in Figure 3.5a, we observe that the adjoint-based method recovers signif-
icantly more eigenvalues than the HDG method when k = 1 and the error tolerance
ϵ = 1%. For example, when h = 0.05, the HDG method only recovers 60 eigenvalues
while the adjoint-based method obtains 242 eigenvalues, which is about four times more.
Let us also display the plot of the number of eigenvalues and the square root of the CPU
time needed to recover these eigenvalues in Figure 3.5b. We see that to recover the same
number of eigenvalues, the adjoint-based method takes much less time.
68
10 15 20 25 30 35
1/h
0
50
100
150
200
250
n
u
m
be
r o
f e
ig
en
va
lu
es
242
60
(a)
0 5 10 15 20 25 30
0
50
100
150
200
250
n
u
m
be
r o
f e
ig
en
va
lu
es
(b)
Figure 3.5: The number of eigenvalues that have a relative error |λ − λh|/λ less than
the tolerance ϵ = 1%. The HDG method (♢) and the adjoint-based method (◦) when
k = 1.
3.4.4 Application to the quantum harmonic oscillator
In this subsection, we consider the application of the adjoint-based method to the
Schro¨dinger’s equation for the quantum harmonic oscillator, which is an important
model in quantum chemistry. The harmonic oscillator equation is
−1
2
∆u+
1
2
|x|2u = λu, x ∈ Rd.
We consider the two dimensional case d = 2. We know the first three eigenvalues are 1,
2, 3 with multiplicity 1, 2, 3, respectively. Since the eigenfunctions decay exponentially,
we can consider a finite computational domain and impose zero boundary condition. In
our calculation, we use Ω = [−10, 10]2 and we define
λh(uh) =
−(qh,∇uh)Th + ⟨uh, q̂h · n⟩∂Th + (|x|2uh, uh)Th
2∥uh∥2L2(Ω)
.
Following the proof in Section 3.3, we can see the adjoint-correction term ACh for this
equation is the same as the one defined in Theorem 3.1. Therefore, we can define
λACh and the adjoint-based method accordingly. Let us use the same HDG method as
considered in Section 3.4.2. We report the errors and the rates of convergence of the
first, second and fourth eigenvalues in Table 3.5 when k = 1, 2, 3. We see that rates of
69
the convergence match with what we observe in previous sections. Eigenvalues obtained
from the HDG method converge with a rate of O(h2k+1), while the eigenvalues from the
adjoint-based method converge with a rate of O(h4k+2). Let us also display the plot
of the fourth eigenfunction in Figure 3.6. We observe that the approximation obtained
by the adjoint-based method (right) is smoother than the one obtained by the HDG
method (left), a clear consequence of the smoothing effect of the convolution with the
kernel K2kh .
First Second Fourth
k 1/h error order error order error order
Approximation given by λh(uh)
1
20 1.77e-02 – 5.33e-02 – 6.49e-02 –
30 5.89e-03 2.72 1.76e-02 2.73 2.57e-02 2.28
40 2.56e-03 2.90 7.76e-03 2.85 1.21e-02 2.62
50 1.33e-03 2.95 4.03e-03 2.93 6.47e-03 2.81
2
20 5.32e-04 – 1.70e-03 – 8.84e-03 –
30 6.99e-05 5.01 2.83e-04 4.42 4.81e-04 7.18
40 1.68e-05 4.96 6.82e-05 4.95 1.18e-04 4.88
50 5.50e-06 4.99 2.24e-05 5.00 3.89e-05 4.97
3
30 5.64e-07 – 2.85e-06 – 1.23e-05 –
40 7.59e-08 6.97 3.84e-07 6.96 1.65e-06 6.97
50 1.59e-08 7.00 8.05e-08 7.00 3.45e-07 7.02
60 4.43e-09 7.01 2.24e-08 7.02 9.57e-08 7.03
Approximation given by λACh (u
∗
h)
1
20 6.05e-03 – 2.39e-02 – 4.23e-02 –
30 4.77e-04 6.26 1.87e-03 6.29 2.87e-03 6.64
40 6.85e-05 6.74 2.83e-04 6.56 4.58e-04 6.38
50 1.49e-05 6.85 6.31e-05 6.72 1.05e-04 6.59
2
20 2.53e-03 – 1.05e-02 – 3.66e-02 –
30 6.17e-05 9.16 3.10e-04 8.68 5.57e-04 10.32
40 3.12e-06 10.37 1.69e-05 10.11 3.07e-05 10.08
50 2.71e-07 10.96 1.52e-06 10.79 2.77e-06 10.77
3
30 2.60e-05 – 1.59e-04 – 8.01e-04 –
40 6.30e-07 12.93 4.29e-06 12.55 2.49e-05 12.07
50 2.76e-08 14.01 1.99e-07 13,76 1.24e-06 13.45
60 1.92e-09 14.62 1.43e-08 14.45 9.25e-08 14.23
Table 3.5: History of convergence for the harmonic oscillator.
70
Figure 3.6: The HDG approximation (left) and the accuracy-enhanced approximation
(right) of the fourth eigenfunction of the harmonic oscillator equation.
3.4.5 The relevance of the adjoint-correction term
In this subsection, we first show that including the adjoint-correction term ACh(u
∗
h)
in the definition of the approximate functional is important as it enhances the quality of
the approximation. We then show that the adjoint-corrected HDG eigenvalue converges
with a rate greater than O(h2k+1) and explore the possibility of using the adjoint-
correction term ACh(uh) as an a posteriori error estimate. In the following experiments,
we choose the same local spaces for the HDG method as considered in Section 3.4.2.
The relevance of ACh(u
∗
h)
Here we show that the adjoint-correction term ACh(u
∗
h) plays an vital role in re-
ducing the error. Let us first consider the eigenvalue problem on a unit square domain
studied in Section 3.4.2. We have already listed the results of λACh (u
∗
h) = λh(u
∗
h) +
ACh(u
∗)/∥u∗h∥2L2(Ω) in Table 3.2. To see the importance of ACh(u∗h), we report the errors
and rates of convergence of λh(u
∗
h) in Table 3.6. Comparing the two tables, we observe
that the error of the adjoint-corrected eigenvalue approximation, |λ(u) − λACh (u∗h)|, is
significantly smaller than |λ(u)−λh(u∗h)|, especially for k = 2, 3. Moreover, it converges
with a faster rate of O(h4k+2).
To end, we display |ACh(u∗h)|/∥u∗h∥2L2(Ω) and ∥u∗h∥2L2(Ω)|λ(u) − λh(u∗h)|/|ACh(u∗h)|
on Table 3.7. We observe that the latter ratio is very close to 1. This shows that
71
ACh(u
∗
h)/∥u∗h∥2L2(Ω) is a good candidate for an a posteriori error estimate for λ(u) −
λh(u
∗
h).
First Second Third Fourth Fifth
k 1/h error order error order error order error order error order
1
20 2.97e-06 – 7.27e-05 – 7.27e-05 – 7.60e-04 – 5.78e-04 –
30 2.66e-07 5.95 6.85e-06 5.83 6.85e-06 5.83 6.99e-05 5.88 6.28e-05 5.47
40 4.75e-08 5.99 1.23e-06 5.96 1.23e-06 5.96 1.24e-05 6.01 1.15e-05 5.91
50 1.24e-08 6.01 3.23e-07 6.00 3.23e-07 6.00 3.22e-06 6.05 3.04e-06 5.96
2
20 8.22e-07 – 2.06e-05 – 2.06e-05 – 3.12e-05 – 2.04e-04 –
30 1.88e-07 3.64 5.61e-06 3.21 5.61e-06 3.21 1.02e-05 2.77 5.70e-05 3.14
40 6.15e-08 3.88 1.93e-06 3.71 1.93e-06 3.71 3.66e-06 3.55 2.05e-05 3.55
50 2.55e-08 3.94 8.14e-07 3.86 8.14e-07 3.86 1.57e-06 3.79 8.88e-06 3.76
3
20 2.19e-14 – 8.36e-12 – 8.36e-12 – 2.49e-11 – 1.12e-09 –
25 2.49e-15 9.74 1.33e-13 18.56 1.33e-13 18.56 3.65e-13 18.93 3.66e-11 15.33
30 4.17e-16 9.81 8.95e-14 2.17 8.95e-14 2.17 2.40e-13 2.30 1.40e-11 5.28
35 9.11e-17 9.86 2.72e-14 7.73 2.72e-14 7.73 7.06e-14 7.93 2.94e-12 10.12
Table 3.6: History of convergence of λh(u
∗
h) on the unit square domain.
First Second Third Fourth Fifth
k 1/h
ACh(u
∗
h)
∥u∗
h
∥2L2(Ω)
ratio
ACh(u
∗
h)
∥u∗
h
∥2L2(Ω)
ratio
ACh(u
∗
h)
∥u∗
h
∥2L2(Ω)
ratio
ACh(u
∗
h)
∥u∗
h
∥2L2(Ω)
ratio
ACh(u
∗
h)
∥u∗
h
∥2L2(Ω)
ratio
1
20 2.98e-06 1.00 7.40e-05 0.98 7.40e-05 0.98 7.63e-04 1.00 6.09e-04 0.95
30 2.67e-07 1.00 6.96e-06 0.98 6.96e-06 0.98 7.01e-05 1.00 6.57e-05 0.96
40 4.76e-08 1.00 1.25e-06 0.98 1.25e-06 0.98 1.24e-05 1.00 1.20e-05 0.96
50 1.25e-08 1.00 3.29e-07 0.98 3.29e-07 0.98 3.23e-06 1.00 3.17e-06 0.96
2
20 8.22e-07 1.00 2.06e-05 1.00 2.06e-05 1.00 3.12e-05 1.00 2.04e-04 1.00
30 1.88e-07 1.00 5.61e-06 1.00 5.61e-06 1.00 1.02e-05 1.00 5.70e-05 1.00
40 6.15e-08 1.00 1.93e-06 1.00 1.93e-06 1.00 3.66e-06 1.00 2.05e-05 1.00
50 2.55e-08 1.00 8.14e-07 1.00 8.14e-07 1.00 1.57e-06 1.00 8.88e-06 1.00
3
20 2.19e-14 1.00 8.36e-12 1.00 8.36e-12 1.00 2.49e-11 1.00 1.12e-09 1.00
25 2.49e-15 1.00 1.33e-13 1.00 1.33e-13 1.00 3.65e-13 1.00 3.66e-11 1.00
30 4.17e-16 1.00 8.95e-14 1.00 8.95e-14 1.00 2.40e-13 1.00 1.40e-11 1.00
35 9.11e-17 1.00 2.72e-14 1.00 2.72e-14 1.00 7.06e-14 1.00 2.94e-12 1.00
Table 3.7: The adjoint-correction term |ACh(u∗h)|/∥u∗h∥2L2(Ω) and the ratio ∥u∗h∥2L2(Ω)|λ−
λh(u
∗
h)|/ACh(u∗h)
The relevance of ACh(uh)
As discussed in Section 3.2.1, for the HDG method, we have that
ACh(uh) = −⟨uh − ûh, τ(uh − ûh)⟩∂Th ,
is a non-zero term. To uncover its importance, let us first consider the eigenvalue
problem on a unit square domain studied in Section 3.4.2. In Table 3.8, we report the
values of |ACh(uh)|/∥uh∥2L2(Ω) and the ∥uh∥2L2(Ω)|λ(u)−λh(uh)|/|ACh(uh)|. We observe
72
that the ratio is very close to 1. This indicates that ACh(uh)/∥uh∥2L2(Ω) provides a
good error estimate of λ− λh(uh), and, since it is computable, that it could be used to
enhance the approximation to the eigenvalue.
To see this, we report the history of convergence of the adjoint-corrected eigenvalue
λACh (uh) = λh(uh) +ACh(uh)/∥uh∥2L2(Ω) in Table 3.9. Comparing it with Table 3.2, we
see that the error |λ(u)−λACh (uh)| is several orders of magnitude smaller than the error
|λ(u)− λh|. We also observe that the rate of convergence of λACh (uh) is clearly of order
O(h2k+2) for k = 2, 3, which is an additional order higher than that of λh.
To end, let us consider the L-shaped domain studied in Section 3.4.2. We report
the values of |ACh(uh)|/∥uh∥2L2(Ω) and the ratio ∥uh∥2L2(Ω)|λ(u)− λh(uh)|/|ACh(uh)| in
Table 3.10. As in our previous experiments, we see that the ratio is very close to one for
the third eigenvalue. However, it increases as h decreases for the first eigenvalue. This
seems to be a reflection of the lack of smoothness that the first eignefunction displays,
as its behavior is singular at the origin, unlike the third eigenfunction.
This implies that the adjoint-correction term ACh(uh) cannot always be taken as
an a posteriori error estimate. However, if we plot the adjoint-correction term element-
by-element as shown in Figure 3.7, we observe that, it captures well the location of the
biggest approximation error, that is, the corner at which the exact eigenfunction has a
singularity.
First Second Third Fourth Fifth
k 1/h
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
1
20 3.40e-04 1.01 3.08e-03 1.03 3.08e-03 1.03 9.14e-03 0.42 1.78e-02 0.95
30 1.00e-04 1.02 8.70e-04 1.06 8.70e-04 1.06 1.96e-03 0.77 4.42e-03 1.13
40 4.21e-05 1.02 3.62e-04 1.06 3.62e-04 1.06 7.48e-04 0.90 1.79e-03 1.14
50 2.16e-05 1.02 1.85e-04 1.06 1.85e-04 1.06 3.67e-04 0.96 9.02e-04 1.13
2
20 8.25e-08 1.04 2.67e-06 1.14 2.67e-06 1.14 5.27e-06 1.14 2.99e-05 1.32
30 1.09e-08 1.03 3.54e-07 1.10 3.54e-07 1.10 6.97e-07 1.10 3.97e-06 1.22
40 2.59e-09 1.02 8.42e-08 1.07 8.42e-08 1.07 1.66e-07 1.07 9.45e-07 1.16
50 8.51e-10 1.02 2.76e-08 1.06 2.76e-08 1.06 5.44e-08 1.06 3.10e-07 1.13
3
20 1.04e-11 1.03 1.34e-09 1.11 1.34e-09 1.11 2.66e-09 1.11 3.40e-08 1.25
25 2.19e-12 1.02 2.81e-10 1.09 2.81e-10 1.09 5.60e-10 1.09 7.15e-09 1.20
30 6.12e-13 1.02 7.85e-11 1.08 7.85e-11 1.08 1.56e-10 1.07 2.00e-09 1.17
35 2.08e-13 1.02 2.67e-11 1.06 2.67e-11 1.06 5.32e-11 1.06 6.81e-10 1.14
Table 3.8: The adjoint-correction term |ACh(uh)|/∥uh∥2L2(Ω) and the ratio ∥uh∥2L2(Ω)|λ−
λh(uh)|/|ACh(uh)|
73
First Second Third Fourth Fifth
k 1/h error order error order error order error order error order
1
20 3.89e-06 – 7.94e-05 – 7.94e-05 – 5.32e-03 – 8.33e-04 –
30 1.88e-06 1.79 5.53e-05 0.89 5.53e-05 0.89 4.52e-04 6.08 5.58e-04 0.99
40 7.61e-07 3.15 2.27e-05 3.09 2.27e-05 3.09 7.43e-05 6.28 2.50e-04 2.79
50 3.51e-07 3.47 1.05e-05 3.47 1.05e-05 3.47 1.60e-05 6.89 1.17e-04 3.40
2
20 3.32e-09 – 3.86e-07 – 3.86e-07 – 7.34e-07 – 9.61e-06 –
30 2.93e-10 5.98 3.41e-08 5.98 3.41e-08 5.98 6.68e-08 5.91 8.53e-07 5.97
40 5.22e-11 6.00 6.09e-09 5.99 6.09e-09 5.99 1.20e-08 5.96 1.52e-07 5.99
50 1.37e-11 6.01 1.60e-09 5.99 1.60e-09 5.99 3.17e-09 5.98 4.00e-08 5.99
3
20 3.21e-13 – 1.50e-10 – 1.50e-10 – 2.86e-10 – 8.47e-09 –
25 5.41e-14 7.98 2.53e-11 7.98 2.53e-11 7.98 4.92e-11 7.90 1.43e-09 7.97
30 1.26e-14 7.98 5.91e-12 7.98 5.91e-12 7.98 1.16e-11 7.93 3.34e-10 7.98
35 3.68e-15 7.99 1.72e-12 7.99 1.72e-12 7.99 3.40e-12 7.95 9.75e-11 7.99
Table 3.9: History of convergence of λACh (uh) on the unit square domain.
First Third
k 1/h
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
|ACh(uh)|
∥uh∥2L2(Ω)
ratio
1
10 1.65e-03 8.37 2.91e-03 0.91
20 2.68e-04 22.03 3.40e-04 1.01
30 9.44e-05 37.06 1.00e-04 1.02
40 4.55e-05 52.85 4.21e-05 1.02
2
10 2.16e-04 36.43 2.61e-06 1.08
20 4.45e-05 71.06 8.25e-08 1.04
30 1.75e-05 105.60 1.09e-08 1.03
40 9.02e-06 140.12 2.59e-09 1.02
Table 3.10: The adjoint-correction term |ACh(uh)|/∥uh∥2L2(Ω) and the ratio
∥uh∥2L2(Ω)|λ− λh(uh)|/ACh(uh) in an L-shaped domain.
First Third
|λ1 − λAC1,h (uh)| |λ3 − λAC3,h (uh)|
k 1/h error order error order
1
10 1.55e-02 – 2.59e-04 –
20 6.16e-03 1.33 3.89e-06 6.06
30 3.59e-03 1.33 1.88e-06 1.79
40 2.45e-03 1.33 7.61e-07 3.15
2
10 8.07e-03 – 2.04e-07 –
20 3.21e-03 1.33 3.32e-09 5.94
30 1.87e-03 1.33 2.93e-10 5.98
40 1.27e-03 1.33 5.23e-11 5.99
Table 3.11: History of convergence of λACh (uh) on an L-shaped domain.
74
Figure 3.7: The element-wise plot of the adjoint-correction term |ACh(uh)|K = ⟨uh −
ûh, (q̂h − qh) · n⟩∂K on an L-shaped domain when k = 1 and h = 0.1. It captures well
the location of the singularity of the eigenfucntion.
3.5 Conclusion
In this chapter, we present a new method for computing high-order accurate approx-
imations of eigenvalues in terms of Galerkin approximations. We detail the definition
of the method and provide the approximation error expression in Section 3.2. Our ex-
tensive numerical results in Section 3.4 show that the approximate eigenvalues obtained
by our method converge with a rate of O(h4k+2) when the eigenfunction is smooth
enough. This is significantly faster than the classic finite element methods in the lit-
erature, which converge at most O(h2k+2). Although we consider a simple self-adjoint
eigenvalue problem, the framework of our method is general enough to be applied to
a wide variety of eigenvalue problems and other non-linear functional problems. Let
us sketch how we can extend our method to non self-adjoint eigenproblems such as
the convection-diffusion problem studied in [48]. Suppose the primal non self-adjoint
75
eigenproblem we consider is to seek eigenpair (λp, up) ∈ C× V such that
a(up, v) = λp b(up, v) ∀v ∈ V,
in the complex plane C and some function space V , where a(·, ·) and b(·, ·) are sesquilin-
ear forms. The adjoint (dual) eigenproblem is defined as to seek eigenpair (λd, ud) ∈
C× V such that
a(v, ud) = λdb(v, ud) ∀v ∈ V,
where (·) denotes complex conjugate. The primal and adjoint eigenvalues are connected
via λp = λd. Now let uph, u
d
h ∈ V are numerical approximations of up and ud, respectively.
We can define
λh =
a(uph, u
d
h)
b(uph, u
d
h)
as an approximation to λp.
Then simple calculation shows that
λp − λh =
λpb(uph − up, udh − ud)− a(up − uph, ud − udh)
b(uph, u
d
h)
.
So roughly speaking, we have
|λp − λh| ≤ C
∣∣∣∣∣∣uph − up∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣udh − ud∣∣∣∣∣∣∣∣∣
for some norm |||·|||. Therefore, following Algorithm 1, we can use the accuracy-enhancing
technique for the primal and dual approximate eigenfunctions and use the formula
defined above to obtain a high order accurate eigenvalue approximation.
Our preliminary numerical experiments show that, for the particular HDG method
on rectangular meshes we considered in Section 3.4.5 , the adjoint-correction term is
asymptotically exact a posteriori error estimate when the eigenfunction is very smooth.
Since, to compute the adjoint-correction term, no refined mesh or higher-order approx-
imation spaces are needed for solving the adjoint problem, this represents an advantage
over the standard approach. However, further investigation is needed to see if this
behavior also holds for other HDG methods.
We also observe that the adjoint-correction term can identify the crucial elements
76
that need to be refined even when the eigenfunction is not smooth. This suggests that
the adjoint-correction term could be used to drive adaptive algorithms in an efficient
way. The study of this approach to adaptivity is presented in the next chapter.
Chapter 4
An adjoint-based adaptive error
approximation of functionals by
the HDG method for
second-order elliptic equations
4.1 Introduction
In this chapter, we present a novel, fully computable error approximation and mesh
adaptation for functionals defined by second-order elliptic equations. Unlike most
adjoint-based error estimations in the literature, the novelty of our method is that
the error approximation is obtained without requiring an auxiliary finer mesh or higher
order approximation spaces for solving the adjoint problem. This reduces the compu-
tational cost and eases the implementation of the problem.
In Chapter 2, under suitable regularity assumptions, a super-convergent method of
order O(h4k) is studied for approximating linear functionals with Galerkin method of
polynomial degree k. Its extension to eigenvalue problems is explored in Chapter 3.
However, no such high order convergence is observed when the solutions are not smooth
enough. The error approximation method proposed in this chapter is a stepping-stone
to incorporate adaptive mesh refinements with the super-convergent method in [27] to
77
78
handle non-smooth solutions.
Output functionals, such as the normal force on a surface, the mean value of a solu-
tion, and the flux across certain boundary, are of great interest in scientific and engineer-
ing applications, see [40, 53, 67, 61, 59]. The output functional J(u) is often determined
by one or several field variables u, which is governed by some model partial differential
equations (PDEs). Therefore, one can get an approximation J(uh) by approximating u
numerically. To compute functionals as accurately as possible within given computing
resources, a posteriori error estimate of J(u)−J(uh) and adaptive mesh refinement are
usually employed. Based on the pioneering work of Eriksson, Estep, Hansbo and John-
son in [35], Becker and Rannacher developed the dual-weighted-residual (DWR) method
for error control and mesh optimization for computing functionals within the context
of finite element methods. See the thorough review [5] and the reference therein. At
around the same time, Pierce and Giles propose the adjoint-based recovery method in
[62] that extends the super-convergence results automatically inherent in many finite el-
ement methods to any numeric method, including the finite difference and finite volume
methods. Later, Giles and Su¨li explore a posteriori error analysis of the adjoint-based
method for linear and non-linear partial differential equations in [41]. Venditti and Dar-
mofal apply the adjoint-based error correction technique to viscous flows and inviscid
flows in [73, 72]. In recent decades, considerable number of studies on the develop-
ment of reliable a posteriori error estimate and mesh adaption are carried out, see, for
example, [44, 45, 43, 52, 68, 34]. In particular, see the review papers of Hartmann and
Houston [46] for fluid dynamic applications focusing on discontinuous Galerkin methods.
For general discretizations, Fidkowski and Darmofal provide a comprehensive literature
review in [37].
More recently, the adjoint-based error estimation for the HDG methods are stud-
ied. In [74], Woopen, May and Schu¨tz apply the HDG method and the adjoint-based
approach to nonlinear convection–diffusion problems. In [15], a priori estimates are de-
veloped for the convergence of output functionals, error estimates and the localizations
of the estimate for the HDG method. In [36], a mesh optimization approach is coupled
with the HDG method.
Roughly speaking, the key ingredient to estimate the functional error is to approx-
imate the weight terms u − uh and v − vh, where u (uh) and v (vh) are the exact
79
(approximate) solutions of the primal and adjoint problem, respectively. Since the ex-
act solutions u and v are usually unavailable, approximations u˜ and v˜ obtained with a
higher-order method or on a finer mesh are often used. However, this results in addi-
tional computational cost. As remarked by Becker and Rannacher in [5], this approach
may not be feasible for complex high dimensional problems. An alternative approach
is to reconstruct higher-order approximations u∗h and v
∗
h evaluated from uh and vh. For
example, a patch-wise, higher-order interpolation method for structured quadrilateral
meshes is studied in [4]. This approach is further extended to unstructured meshes of
simplices in [6, 14]. In [32], a least-squares reconstruction of the solution is considered.
The error approximation method we propose here takes advantage of a local post-
processing technique of the HDG method for second-order elliptic equations. Thanks
to this technique, the primal problem and the adjoint problem are solved in the same
manner without requiring an auxiliary finer mesh or higher order approximation spaces.
The HDG method is first introduced in the framework of steady-state diffusion in [21]
and the local post-processing technique is studied in [24, 26] and later in [19]. It has
been proved that with suitable chosen local spaces, the local post-processed solution u∗h
converges one order faster than the HDG solution uh. What’s more, the post-processing
procedure is completely local so that it can be implemented in parallel.
The remainder of this chapter is organized as follows. In Section 4.2, we recall
the adjoint-based method and present our error approximation of functionals when
the adjoint-based method is applied to the HDG method. Next, in Section 4.3, we
consider the error approximation for three types of functionals. Then in Section 4.4, we
carry out numerical experiments to test our error approximation for these three type
of functionals. Finally, in Section 4.5, we conclude with a short discussion of possible
extensions and a description of our forthcoming work.
4.2 The adjoint-based error approximation of functionals
In this section, we first recall the adjoint-based method and the HDG method in
an abstract form. Next, we present our error approximation method for approximating
functionals.
To illustrate the method, let us consider a linear functional J(·) that is determined
80
by the solutions of a second-order elliptic problem:
q +∇u = 0 in Ω,
∇ · q = f in Ω,
u = uD on ∂Ω,
(4.2.1)
where f ∈ L2(Ω) and uD ∈ L2(∂Ω). For non-linear functionals or non-linear PDEs, our
method can be applied with a linearization of the functional or PDEs. We refer to the
discussion in [5, 43, 37].
4.2.1 The adjoint-based method in an abstract form
Let us start by recalling the adjoint-based method within the context of Galerkin
methods, see also [62, 41]. We use the following notation because it will unify our
analysis in the next section.
Let J(u, q, u|Fh) be the linear functional we are interested in, where Th is a con-
forming triangulation of Ω and Fh is the set of all the faces F ∈ ∂K for all the elements
K ∈ Th. Let us assume we discretize the model problem (4.2.1) by some Galerkin
method as follows: find (uh, qh, ûh) ∈Wh × V h ×Mh such that
A(uh, qh, ûh;w, r, µ) = F (w, r, µ) ∀ (w, r, µ) ∈Wh × V h ×Mh. (4.2.2)
Here uh and qh are approximations of u and q, and ûh approximates the traces of u on
all the faces. Furthermore, let us assume the discretization (4.2.2) is consistent, i.e. the
exact solution (u, q, u|Fh) ∈W × V ×M satisfies the following equation:
A(u, q, u|Fh ;w, r, µ) = F (w, r, µ) ∀ (w, r, µ) ∈W × V ×M. (4.2.3)
Here W × V × M are suitably chosen function spaces containing the exact solution
(u, q, u|Fh) and the finite dimensional discrete spaces Wh × V h ×Mh are subspaces of
W × V ×M .
We define J(uh, qh, ûh) as the standard approximation of the functional J(u, q, u|Fh).
Next, we consider the adjoint problem as follows: find (v,p, v|Fh) ∈ W × V ×M such
81
that
A(w, r, µ; v,p, v|Fh) = J(w, r, µ) ∀ (w, r, µ) ∈W × V ×M. (4.2.4)
Therefore, we have
J(u, q, u|Fh)− J(uh, qh, ûh)
=J(u− uh, q − qh, u− ûh)
=A(u− uh, q − qh, u− ûh; v,p, v|Fh) by (4.2.4)
In the last equality above, we can add and subtract any choices of (vh,ph, v̂h) ∈ W ×
V ×M and we have
J(u, q, u|Fh)− J(uh, qh, ûh) = ACh + Eh, (4.2.5)
where
ACh = A(u− uh, q − qh, u− ûh; vh,ph, v̂h), (4.2.6)
and
Eh = A(u− uh, q − qh, u− ûh; v − vh,p− ph, v − v̂h). (4.2.7)
We refer (4.2.6) as the adjoint-correction term and (4.2.7) as the error term. Note that
the adjoint-correction term is computable since
ACh =A(u− uh, q − qh, u− ûh; vh,ph, v̂h)
=A(u, q, u|Fh ; vh,ph, v̂h)−A(uh, qh, ûh; vh,ph, v̂h)
=F (vh,ph, v̂h)−A(uh, qh, ûh; vh,ph, v̂h) by (4.2.3).
As we can see, if we correct the approximation J(uh, qh, ûh) by adding the adjoint-
correction term ACh, then the error |J(u, q, u|Fh)− J(uh, qh, ûh)−ACh| becomes
|Eh| ≤ |||(u, q, u|Fh)− (uh, qh, ûh)||||||(v,p, v|Fh)− (vh,ph, v̂h)|||,
for some norm |||·|||. This shows that if we choose (vh,ph, v̂h) to be the approximation
of the adjoint problem, the convergence rate of the corrected functional approximation
J(uh, qh, ûh) + ACh is the sum of the convergence rates of the primal and adjoint
82
solutions. This is the so-called adjoint-based recovery method.
Let us end this section by noting that the adjoint-based method shows above can
be applied to any numerical method, not limited to Galerkin methods. However, ACh
is zero if (uh, qh, ûh) are the solutions of (4.2.2) and (vh,ph, v̂h) ∈ Wh × V h × Mh.
Therefore, for many finite element methods, the approximation J(uh, qh, ûh) super-
converges automatically without correction.
4.2.2 The HDG method
Next, let us introduce a class of Galerkin methods, called the hybridizable discon-
tinuous Galerkin method (HDG). We use the HDG method here because an element-
by-element local post-processing technique can be applied to the HDG solution uh and
we can obtain a super-convergent approximation u∗h. This is the key for the functional
error approximation we will discuss in the next section.
Let us use notation
(u, v)Th :=
∑
K∈Th
∫
K
uv dx and ⟨v , w⟩∂Th :=
∑
K∈Th
∫
∂K
vu ds.
The HDGmethod seeks an approximation, (uh, qh, ûh), of the exact solution, (u|Ω, q|Ω, , u|Fh),
of the model problem (4.2.1) such that the following equations are satisfied:
(qh, r)Th − (uh,∇ · r)Th + ⟨ûh , r · n⟩∂Th = 0, (4.2.8a)
−(qh,∇w)Th + ⟨q̂h · n , w⟩∂Th = (f, w)Th , (4.2.8b)
⟨ûh , µ⟩∂Ω =< uD, µ >∂Ω, (4.2.8c)
⟨q̂h · n , µ⟩∂Th\∂Ω = 0, (4.2.8d)
for all (r, w, µ) ∈ V h ×Wh ×Mh, where
q̂h · n := qh · n+ τ(uh − ûh) on ∂Th.
Here τ is called the stabilization function and the finite dimensional spacesWh×V h×Mh
83
are defined as
Wh := {w ∈ L2(Th) :w|K ∈W(K) ∀K ∈ Th},
V h := {v ∈ L2(Th) :v|K ∈ V (K) ∀K ∈ Th},
Mh := {µ ∈ L2(Fh) :µ|F ∈M(F ) ∀F ∈ Fh}.
We can obtain different methods by taking particular choices of the local spaces V (K),
W (K) and M(F ), and the stabilization function τ .
It is known that with suitable chosen local spaces, the HDG method has the super-
convergence property. In particular, we can apply an element-by-element post-processing
procedure to the HDG solution uh and obtain a new approximation u
∗
h that converges
one order faster than uh in the L
2 norm. Let us define the post-processing technique.
Following [69, 24, 26, 19], we define a new approximation to u, u∗h in the space
W ∗h := {v ∈ L2(Th) : v|K ∈W ∗(K) ∀K ∈ Th},
such that on each element K ∈ Th, u∗h satisfies
(∇u∗h,∇v)K = −(qh,∇v)K ∀v ∈W ∗(K), (4.2.9a)
(u∗k, 1)K = (uh, 1)K . (4.2.9b)
For the mesh Th consisting of simplexes in n-dimension, we can have the super-
convergence property by simply choosing V (K) = P k(K), W (K) = P k(K), M(F ) =
P k(F ) and W ∗(K) = P k+1(K), where P k(D) denotes the space of polynomials of
total degree k defined on domain D and P k(D) = [P k(D)]n is the vector space in n-
dimensional space. For other types of meshes, such as squares, cubes and prisms, we
refer the reader to [26, 19] for the choices of local spaces.
Let us end this section by noting that as shown in the Appendix (A.4), we can
rewrite the HDG equations (4.2.8) as the formulation (4.2.2) presented in the previous
84
section, where
A(uh, qh, ûh;w, r, µ)
=− (qh, r)Th + (∇ · qh, w)Th + (uh,∇ · r)Th
+ ⟨τuh , w⟩∂Th − ⟨τ ûh , µ⟩∂Th
− ⟨ûh , r · n+ τ(w − µ)⟩∂Th\∂Ω − ⟨qh · n+ τ(uh − ûh) , µ⟩∂Th\∂Ω
(4.2.10)
and
F (w, r, µ) = (f, w)Th + ⟨uD , r · n+ τ(w − µ)⟩∂Ω. (4.2.11)
Note that for our model problem, the bilinear form A is symmetric i.e.
A(uh, qh, ûh;w, r, µ) = A(w, r, µ;uh, qh, ûh).
4.2.3 Functional error approximation
Now we are ready to present our method to approximate the functional error term
Eh. Let (uh, qh, ûh) and (vh,ph, v̂h) be the HDG solutions of the model problem (4.2.3)
and the adjoint problem (4.2.4), respectively. In this case, we know that the adjoint-
correction term ACh = 0. Based on the defintion (4.2.10), a simple calculation (see
lemma A.5 in the Appendix) shows that the error term Eh in (4.2.7) becomes:
Eh =− (q − qh,p− ph)Th
+ (f −∇ · qh, v − vh)Th + (u− uh,∇ · p−∇ · ph)Th
+ ⟨τ(ûh − uh) , v − vh⟩∂Th + ⟨u− uh , τ(v̂h − vh)⟩∂Th
− ⟨τ(uh − ûh) , vh − v̂h⟩∂Th .
(4.2.12)
Here we use the fact that ∇·q = f . Let us assume that ∇·p = g is also known since it is
the source term of the adjoint problem, which is usually given data from the functional.
See the examples in Section 4.3. We can see that the two terms (f −∇ · qh, v− vh)Th +
(u− uh, g −∇ · ph)Th measure the oscillation from the data.
Let us also note that the error term Eh obtained here is equivalent to the definition
in our previous work [27, 28] except that the last term −⟨τ(uh− ûh) , vh− v̂h⟩∂Th , which
is computable, is put in the adjoint-correction term ACh in our previous work. For
85
more details, see lemma A.6 in the Appendix.
As we can see, to approximate Eh and obtain a computable quantity, we need to
replace the unknowns (u, q) and (v,p) by suitable approximations. To do that, let u∗h
and v∗h be the HDG local post-processed solutions defined by (4.2.9). Then we can
replace (u, q) and (v,p) by (u∗h,−∇u∗h) and (v∗h,−∇vh∗), respectively. Therefore, we
have an error approximation Eh as
Eh =− (∇u∗h + qh,∇v∗h + ph)Th
+ (f −∇ · qh, v∗h − vh)Th + (u∗h − uh, g −∇ · ph)Th
+ ⟨τ(ûh − uh) , v∗h − v̂h⟩∂Th + ⟨u∗h − ûh , τ(v̂h − vh)⟩∂Th
− ⟨τ(uh − ûh) , vh − v̂h⟩∂Th .
(4.2.13)
So we have
J(u, q, u|Fh)− J(uh, qh, ûh) ≈ Eh.
With this error approximation, we can correct the error and define a new approximation
J(uh, qh, ûh) + Eh. But more importantly, we can use the error approximation to drive
adaptive algorithms. To be more specific, we can localize the error approximation Eh
and define a local error indicator
ηK =− (∇u∗h + qh,∇v∗h + ph)K
+ (f −∇ · qh, v∗h − vh)K + (u∗h − uh, g −∇ · ph)K
+ ⟨τ(ûh − uh) , v∗h − v̂h⟩∂K + ⟨u∗h − ûh , τ(v̂h − vh)⟩∂K
− ⟨τ(uh − ûh) , vh − v̂h⟩∂K .
Next, we define a localized error bound Eh,|·| =
∑
K∈Th |ηK | and we have
|J(u, q, u|Fh)− J(uh, qh, ûh)| ≈ |
∑
K∈Th
ηK | ≤ Eh,|·|.
Following the classic adaptive mesh refinement paradigm
Solve → Estimate → Mark → Refine,
86
we can design the following adaptive algorithms to approximate the functional J(u, q, u|Fh)
in an accurate and efficient way.
Algorithm 2: The adaptive algorithm to approximate the functional
J(u, q, u|Fh)
(1) Construct an initial mesh Th.
(2) Compute the HDG solutions (uh, qh, ûh) and (vh,ph, v̂h) for model
problem (4.2.3) and the adjoint problem (4.2.4), respectively.
(3) Compute the local post-processed solutions u∗h and v
∗
h based on (4.2.9).
(4) Compute the approximation J(uh, qh, ûh) and the error approximation
Eh =
∑
K ηK .
(5)Check if Eh,|·| =
∑
K |ηK | ≤ TOL, where TOL is a given error tolerance,
then STOP.
(6) Otherwise, refine (and coarsen) mesh elements based on the size of |ηK |
and generate a new mesh Th. Go to Step 2.
Let us end this section by noting that so far the theoretical analysis of the efficiency
and reliability of our adaptive method is not clear. Unlike the standard error estimate,
which is carried out in terms of some error norms, the main difficulty here is to provide
estimates of the inner-product type of terms in Eh with respect to Eh. It’s not clear
how we can connect these inner-product type of terms with inequalities. Further study
needs to be carried out.
4.3 Examples of approximating functionals
In this section, we present examples of approximating 3 types of functionals. The
first two are linear functionals, discussed under the framework in the Section 4.2. The
last one is a non-linear functional that is not covered above, but we show that a similar
idea can be applied.
87
4.3.1 Volume integral
The first type of functional we consider is a volume integral
J1(u) =
∫
Ω
u g dx,
where g ∈ L2(Ω). For example, one may be interested in the mean value of the solution
u over some domain, which is represented by the first type of functional. We define
J(w, r, µ) = (w, g)Th for the adjoint problem (4.2.4) and we can check that it corresponds
to the following equations:
p+∇v = 0 in Ω,
∇ · p = g in Ω,
v = 0 on ∂Ω.
(4.3.1)
Therefore, we define J(uh, qh, ûh) as the approximation of J1(u) obtained by the
HDG method and we have the error approximation Eh in (4.2.13) with ∇ · p = g.
4.3.2 Normal flux
The second type of functional is a boundary integral, the weighted normal flux of u
over ∂Ω,
J2(u) =
∫
∂Ω
∂u
∂n
ψds,
where n is the unit outward normal vector and the weight function ψ ∈ H1/2(Ω).
For example, the heat flux across the boundary in thermodynamics and the lift and
drag forces in fluid dynamics can be considered as this type of functional. Now we
define J(w, r, µ) = −⟨ψ , r · n + τ(w − µ)⟩∂Ω in (4.2.4). It’s easy to check that
J(u, q, u|Fh) = −⟨ψ , q ·n⟩∂Ω = ⟨ψ , ∇u ·n⟩∂Ω = J2(u). The adjoint-problem defined in
(4.2.4) corresponds to the following equations:
p+∇v = 0 in Ω,
∇ · p = 0 in Ω,
v = −ψ on ∂Ω.
(4.3.2)
Therefore, we define J(uh, qh, ûh) = −⟨ψ , qh·n+τ(uh−ûh)⟩∂Ω as the approximation
88
and we have the error approximation Eh as (4.2.13) with ∇ · p = 0.
4.3.3 Eigenvalues
The third type of functional we consider is the eigenvalue of the following second-
order elliptic eigenvalue problem:
q +∇u = 0 in Ω,
∇ · q = λu in Ω,
u = 0 on ∂Ω.
(4.3.3)
Based on the Rayleigh quotient formula, we consider the non-linear functional
J3(u) =
(∇u,∇u)Ω − ⟨u,∇u · n⟩∂Ω
∥u∥2L2(Ω)
.
We see that the if u is the corresponding eigenfunction of the eigenvalue λ, then
J3(u) = λ.
Following [42], the HDG discretization for (4.3.3) reads
A(uh, qh, ûh;w, r, µ) = λh(uh, w)Th ∀(w, r, µ) ∈Wh × V h ×Mh, (4.3.4)
where A is the same as defined in (4.2.10). We can then define an approximation to
J3(u) as
J(uh, qh, ûh) =
A(uh, qh, ûh;uh, qh, ûh)
(uh, uh)Th
. (4.3.5)
We can also see that if (uh, qh, ûh) is the solution of (4.3.4), then J(uh, qh, ûh) = λh.
Now let u be the eigenfunction corresponding to the eigenvalue λ. Let uh and λh be
89
the approximation of u and λ, respectively, obtained by the HDG method. We have
(λ− λh)(uh, uh)Th
= λ(uh, uh)Th −A(uh, qh, ûh;uh, qh, ûh)
= λ(uh − u, uh − u) + λ(uh, u)Th + λ(u, uh)Th − λ(u, u)Th −A(uh, qh, ûh;uh, qh, ûh)
= λ(uh − u, uh − u)−A(u− uh, q − qh, u− ûh;u− uh, q − qh, u− ûh)
In the last step, we use the consistency of the discretization λ(u,w)Th = A(u, q, u;w, r, µ)
and the fact that HDG discretization for (4.3.3) is symmetric.
Notice that
A(u− uh, q − qh, u− ûh;u− uh, q − qh, u− ûh)
=− (q − qh, q − qh)Th
+ (λu−∇ · qh, u− uh)Th + (u− uh, λu−∇ · qh)Th
+ ⟨τ(ûh − uh) , u− uh⟩∂Th + ⟨u− uh , τ(ûh − uh)⟩∂Th
− ⟨τ(uh − ûh) , uh − ûh⟩∂Th .
Therefore, following the idea in Section 4.2.3, let u∗h be the local post-processed approx-
imation of the HDG solution uh. We can define Θh as an approximation:
Θh :=− (∇u∗h + qh,∇u∗h + qh)Th
+ (λhu
∗
h −∇ · qh, u∗h − uh)Th + (u∗h − uh, λhu∗h −∇ · qh)Th
+ ⟨τ(ûh − uh) , u∗h − uh⟩∂Th + ⟨u∗h − uh , τ(ûh − uh)⟩∂Th
− ⟨τ(uh − ûh) , uh − ûh⟩∂Th .
Then the error λ− λh is approximated by
Eh = λh(uh − u
∗
h, uh − u∗h)−Θh
(uh, uh)Th
.
4.4 Numerical experiments
In this section, we present the numerical experiments in 2 dimensional space for
the three types of functionals {Ji(u)}3i=1 discussed in Section 4.3. We consider solutions
90
with different smoothness and show that our error approximation is efficient and reliable
to drive adaptive algorithms. To alleviate the notation, let us denote the approximation
J(uh, qh, ûh) simply by Jh. For all the experiments, we use triangular meshes and we
choose the suitable local spaces for the HDG method as discussed in Section 4.2.2. The
order of convergence is computed as
−2 log(e1/e2)
log(N1/N2)
,
where ei is the error |J(u)−Jh| and Ni is the total number of degrees of freedom (DOF)
for i = 1, 2. Let us also define the effectivity index Ieff, the ratio between the error
approximation and the true error:
Ieff = | Eh
J(u)− Jh. |
and the localized effectivity index:
I loceff =
Eh,|·|
|J(u)− Jh|.
If values of Ieff ≈ 1, it indicates that the error approximation Eh approximates the true
error very well. If values of I loceff ≥ 1 , it shows that the localized error bound Eh,|·| is
sharp in the sense that it overestimates the true error.
4.4.1 Uniform Refinement
Let us first test our method with smooth solutions and uniform refinements. We
consider the unit square domain Ω = [0, 1]2 with an initial mesh Th of size h0 = 0.25.
We choose u(x, y) = sin(πx) sin(πy) and define f and uD accordingly for the model
problem (4.2.1).
The first functional
For the first functional J1(u) = (u, g)Ω, we choose g(x, y) = 2π
2 sin(πx) sin(πy) so
that the solution of the adjoint problem (4.3.1) is also smooth. We report the history
of convergence in Table 4.1 and Fig. 4.1. As we can see, the error approximation Eh is
91
fairly close to the true error and the effectivity index Ieff is around 1 for all the meshes
and k = 1, 2, 3. We also see that the index I loceff is slightly larger than 1, which shows
Eh,|·| is a good upper bound of the true error. If we correct the approximation Jh by
adding the error approximation Eh, we observe that the convergence rate of the error
|J(u)− (Jh + Eh)| is order 2k + 2, one order faster than that of |J(u)− Jh|.
k Mesh DOF |J(u)− Jh| order |Eh| Ieff Eh,|·| I loceff |J(u)− Jh − Eh| order
1
1 3.68e+02 8.33e-02 - 9.13e-02 1.09 1.05e-01 1.27 7.92e-03 -
2 1.50e+03 1.19e-02 2.77 1.23e-02 1.04 1.33e-02 1.12 4.84e-04 3.97
3 6.08e+03 1.56e-03 2.90 1.59e-03 1.02 1.68e-03 1.07 2.97e-05 4.00
4 2.44e+04 1.99e-04 2.96 2.01e-04 1.01 2.10e-04 1.06 1.83e-06 4.00
2
1 6.96e+02 1.30e-03 - 1.43e-03 1.10 1.75e-03 1.35 1.35e-04 -
2 2.83e+03 4.54e-05 4.78 4.75e-05 1.05 5.24e-05 1.15 2.13e-06 5.91
3 1.14e+04 1.48e-06 4.91 1.52e-06 1.02 1.66e-06 1.12 3.31e-08 5.97
4 4.59e+04 4.73e-08 4.96 4.78e-08 1.01 5.19e-08 1.10 5.15e-10 5.99
3
1 1.12e+03 1.22e-05 - 1.33e-05 1.09 1.68e-05 1.38 1.13e-06 -
2 4.54e+03 1.05e-07 6.79 1.10e-07 1.04 1.23e-07 1.18 4.48e-09 7.89
3 1.83e+04 8.54e-10 6.91 8.71e-10 1.02 9.84e-10 1.15 1.75e-11 7.96
4 7.35e+04 7.29e-12 6.85 7.36e-12 1.01 7.70e-12 1.06 7.44e-14 7.86
Table 4.1: History of convergence of approximating functional J(u) = (u, g)Ω uniformly
on the unit square Ω by the HDG method. Here Eh is the approximation of the error
and Eh,|·| is the localized error bound.
92
1.2 1.4 1.6 1.8 2 2.2
-5
-4
-3
-2
-1
0
4
1
3
1
1 1.5 2 2.5 3 3.5 4
Mesh
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
1.4 1.6 1.8 2 2.2 2.4
-10
-9
-8
-7
-6
-5
-4
-3
-2
5
6
1
1
1 1.5 2 2.5 3 3.5 4
Mesh
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
1.6
1.5 2 2.5
-14
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
1
7
1
8
1 1.5 2 2.5 3 3.5 4
Mesh
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 4.1: Plots of the error and the error approximation against the DOF in loga-
rithmic scale (left column) and plots of the effectivity index Ieff and I
loc
eff (right column)
when the functional J(u) = (u, g)Ω is approximated uniformly with k = 1 (top row),
k = 2 (middle row), k = 3 (bottom row). Here the solutions are smooth.
93
The second functional
For the second functional J2(u) =< ∇u ·n, ψ >∂Ω, we choose ψ = sin(πx)eπy on ∂Ω
so that the solution of the adjoint problem (4.3.2) is v = −ψ, smooth. The history of
convergence for this test case is reported in Table 4.2 and Fig. 4.2. We observe that the
effectivity index Ieff is approaching to 1 as we refine the mesh for k = 1, 2, 3. This shows
that the error approximation Eh is asymptotically exact in this case. We also see that
Eh,|·| is an upper bound of the true error. Furthermore, we observe that the corrected
approximation Jh + Eh converges one order faster than Jh.
k Mesh DOF |J(u)− Jh| order |Eh| Ieff Eh,|·| I loceff |J(u)− Jh − Eh| order
1
1 3.68e+02 5.77e-02 - 3.94e-02 0.68 1.56e-01 2.71 1.83e-02 -
2 1.50e+03 8.17e-03 2.78 7.03e-03 0.86 1.80e-02 2.21 1.14e-03 3.94
3 6.08e+03 1.06e-03 2.93 9.90e-04 0.94 2.19e-03 2.06 6.87e-05 4.02
4 2.44e+04 1.34e-04 2.97 1.30e-04 0.97 2.69e-04 2.01 4.18e-06 4.02
2
1 6.96e+02 1.26e-03 - 9.03e-04 0.72 2.18e-03 1.73 3.56e-04 -
2 2.83e+03 3.75e-05 5.01 3.22e-05 0.86 6.44e-05 1.72 5.35e-06 5.98
3 1.14e+04 1.14e-06 5.02 1.05e-06 0.93 1.96e-06 1.73 8.17e-08 6.00
4 4.59e+04 3.48e-08 5.01 3.36e-08 0.96 6.03e-08 1.73 1.26e-09 6.00
3
1 1.12e+03 6.64e-06 - 6.37e-06 0.96 1.36e-05 2.05 2.66e-07 -
2 4.54e+03 4.61e-08 7.10 4.55e-08 0.99 9.84e-08 2.13 6.37e-10 8.62
3 1.83e+04 3.38e-10 7.06 3.38e-10 1.00 7.37e-10 2.18 4.25e-13 10.52
Table 4.2: History of convergence of approximating functional J(u) =< ∇u · n, ψ >∂Ω
on the unit square Ω by the HDG method. Here Eh is the approximation of the error
and Eh,|·| is the localized error bound.
94
1.2 1.4 1.6 1.8 2 2.2
-5.5
-5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
1
3
4
1
1 1.5 2 2.5 3 3.5 4
Mesh
0
0.5
1
1.5
2
2.5
3
3.5
4
1.4 1.6 1.8 2 2.2 2.4
-9
-8
-7
-6
-5
-4
-3
-2
5
6
1
1
1 1.5 2 2.5 3 3.5 4
Mesh
0
0.5
1
1.5
2
2.5
3
1.5 1.6 1.7 1.8 1.9 2 2.1 2.2
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
1
1
8
7
1 1.5 2 2.5 3
Mesh
0
0.5
1
1.5
2
2.5
3
Figure 4.2: Plots of the error and the error approximation against the DOF in loga-
rithmic scale (left column) and plots of the effectivity index Ieff and I
loc
eff (right column)
when the functional J(u) =< ∇u ·n, ψ >∂Ω is approximated uniformly with k = 1 (top
row), k = 2 (middle row), k = 3 (bottom row). Here the solutions are smooth.
95
The third functional
For the third functional
J3(u) = ((∇u,∇u)Ω − ⟨u,∇u · n⟩Ω)/(u, u)Ω,
let us consider the eigenvalue problem (4.3.3) on the unit square Ω. We know that the
eigenvalues are λm,n = (m2+n2)π2 form,n ∈ N+ and the corresponding eigenfunctions
are {um,n = sin(mπx) sin(nπy)}. Let us consider the 1st eigenvalue, namely, J3(u1,1) =
λ1,1 = 2π2. We first solve the eigenvalue problem by the HDG method following [42] and
obtain the approximation (u1,1h , q
1,1
h , û
1,1
h ), and we define Jh as in (4.3.5). As reported
in Table 4.3 and Fig. 4.3, the results match what we observe for the previous two
functionals. We see that Eh approximates the error very well, Eh,|·| is a good upper
bound of the error, and the corrected approximation Jh + Eh converges with order
2k + 2 for k = 1, 2, 3, which is one order faster than Jh.
k Mesh DOF |λ1,1 − Jh| order |Eh| Ieff Eh,|·| I loceff |λ1,1 − Jh − Eh| order
1
1 3.68e+02 5.67e-01 - 5.63e-01 0.99 5.84e-01 1.03 3.99e-03 -
2 1.50e+03 6.19e-02 3.15 6.35e-02 1.03 6.54e-02 1.06 1.66e-03 1.24
3 6.08e+03 7.10e-03 3.10 7.22e-03 1.02 7.48e-03 1.05 1.18e-04 3.79
4 2.44e+04 8.50e-04 3.05 8.58e-04 1.01 8.89e-04 1.05 7.38e-06 3.98
2
1 6.96e+02 7.73e-03 - 8.35e-03 1.08 9.09e-03 1.18 6.21e-04 -
2 2.83e+03 2.19e-04 5.08 2.28e-04 1.04 2.41e-04 1.10 8.91e-06 6.05
3 1.14e+04 6.50e-06 5.04 6.63e-06 1.02 7.10e-06 1.09 1.35e-07 6.01
4 4.59e+04 1.98e-07 5.02 2.00e-07 1.01 2.15e-07 1.08 2.08e-09 6.00
3
1 1.12e+03 6.56e-05 - 7.05e-05 1.08 8.12e-05 1.24 4.94e-06 -
2 4.54e+03 4.84e-07 7.01 5.02e-07 1.04 5.45e-07 1.13 1.86e-08 7.97
3 1.83e+04 3.66e-09 7.01 3.73e-09 1.02 4.13e-09 1.13 7.12e-11 7.99
Table 4.3: History of convergence of approximating the 1st eigenvalue λ1,1 of (4.3.3)
uniformly by the HDG method. Here Eh is the approximation of the error and Eh,|·| is
the localized error bound.
96
1.2 1.4 1.6 1.8 2 2.2
-6
-5
-4
-3
-2
-1
0
1
1
3
4
1 1.5 2 2.5 3 3.5 4
Mesh
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
1.4 1.6 1.8 2 2.2 2.4
-9
-8
-7
-6
-5
-4
-3
-2
1
1
6
5
1 1.5 2 2.5 3 3.5 4
Mesh
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
1.5 2 2.5
-13
-12
-11
-10
-9
-8
-7
-6
-5
-4
1
7
8
1
1 1.5 2 2.5 3 3.5 4
Mesh
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 4.3: Plots of the error and the error approximation against the DOF in loga-
rithmic scale (left column) and plots of the effectivity index Ieff and I
loc
eff (right column)
when the 1st eigenvalue λ1,1 of (4.3.3) is approximated uniformly with k = 1 (top row),
k = 2 (middle row), k = 3 (bottom row). Here the eigenfunctions are smooth.
97
4.4.2 Adaptive Refinement
In this section, let us test our error estimate method with adaptive mesh refinement.
For all the experiments in this section, we follow the adaptive algorithm stated in
Algorithm 2 except that we stop the algorithm after some number of refinements. We
use the bulk marking strategy proposed by Do¨rfler [33]. We order the elements according
to the size of the local error indicator |ηK | and define the marking set M such that∑
K∈M
|ηK | > θ
∑
K∈Th
|ηK |,
where θ ∈ (0, 1). In our experiments, we choose θ = 70%. We refine the marked
elements in M with the red-green-blue refinement strategy. Triangulation by red-green-
blue refinement maintains the minimum angle condition and avoids over-refinements.
Details about the red-green-blue procedure and its implementation are discussed in
[16, 38].
The first functional
For the first functional J1(u) = (u, g)Ω, let us first consider the same problem in
Section 4.4.1, where the solutions are smooth on the unit square domain Ω = [0, 1]2.
We start with an initial mesh Th of size h0 = 0.25 and report the adaptive results with 8
refinements in Fig. 4.4. We see that the logarithm of the true error |J(u)−Jh| converges,
roughly, in a linear way with the optimal convergence rate 2k + 1, for k = 1, 2, 3. We
also see that the effectivity index Ieff gets closer to one as we adaptively refine the mesh,
which shows that the error approximation behaves quite well in this case.
98
1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8
-5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
-1
-0.5
3
1
1 2 3 4 5 6 7 8
Mesh
0.8
0.85
0.9
0.95
1
1.05
1.1
1.15
1.2
1.4 1.6 1.8 2 2.2 2.4 2.6 2.8
-9
-8
-7
-6
-5
-4
-3
-2
5
1
1 2 3 4 5 6 7 8
Mesh
0
0.5
1
1.5
2
1.4 1.6 1.8 2 2.2 2.4 2.6
-12
-11
-10
-9
-8
-7
-6
-5
-4
1
7
1 2 3 4 5 6 7 8
Mesh
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
1.5
Figure 4.4: Plots of the error and the error approximation against the DOF in logarith-
mic scale for k = 1, 2, 3 (left column) and plots of the corresponding effectivity index
Ieff (right column) when approximating J1(u) = (u, g)Ω and the solutions are smooth.
99
Next, let us consider a singular solution for the model problem (4.2.1):
u = r2/3 sin(
2
3
θ), (4.4.1)
and we define f and uD accordingly. Note that u is singular at the origin and u ∈
H5/3−ϵ(Ω) for any ϵ > 0. We simply choose g = 1 and we start with an initial mesh
Th of size h0 = 0.25. In Fig. 4.5, we display some meshes obtained with our adaptive
procedure. Note that the refinement is concentrated at the origin, where the singularity
of the exact solution is. Note also that for around the same accuracy, the use of poly-
nomial of degree 3 results in a mesh with 60% fewer triangles than the one obtained
by using polynomial of degree 2. It also reduces the number of degrees of freedom by
30%. In Fig. 4.6, we plot the true error |J(u) − Jh|, the error approximation Eh and
the error bound Eh,|·| against the total number of degrees of freedom in logarithmic
scale for k = 1, 2, 3. We observe that Eh,|·| is a upper bound of the error and the error
approximation Eh captures the trend of the error J(u)−Jh well. We can also see in Fig.
4.6 that the effectivity index Ieff is approaching to 1 as we refine. Also note that Eh,|·|
decreases almost in a linear way and the slope is around 2k + 1, for k = 1, 2, 3. This
shows that with our adaptive method, the error converges with the optimal convergence
rate.
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 4.5: Example of adaptive meshes for the first functional when u is singular. Left:
k =2, 254 elements, the number of DOF is 5649 and error is 2.27e-10; Right: k=3, 112
elements, the number of DOF is 3964 and error is 2.49e-10.
100
1 1.5 2 2.5 3
-10
-9
-8
-7
-6
-5
-4
3
1
0 2 4 6 8 10 12
Mesh
0.5
1
1.5
2
2.5
3
3.5
4
4.5
1.4 1.6 1.8 2 2.2 2.4 2.6
-12
-11
-10
-9
-8
-7
-6
-5
5
1
0 2 4 6 8 10 12
Mesh
0
0.5
1
1.5
2
2.5
1.5 1.6 1.7 1.8 1.9 2 2.1 2.2
-13
-12
-11
-10
-9
-8
-7
-6
-5
1
7
0 2 4 6 8 10
Mesh
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
I
eff
Figure 4.6: Plots of the error and the error approximation against the DOF in logarith-
mic scale for k = 1, 2, 3 (left column) and plots of the corresponding effectivity index
Ieff (right column) when approximating J1(u) = (u, g)Ω and u is singular.
101
The second functional
For the second functional J2(u) =< ∇u·n, ψ >∂Ω, let u be the same singular solution
(4.4.1) and we choose ψ to be a discontinuous function on ∂Ω:
ψ =
(x− 0.2)(x− 0.8) + 1 if 0.2 < x < 0.8 and y = 0,0 otherwise.
We start with an initial mesh Th of size h0 = 0.2. In Fig. 4.7 we report sample
adaptive meshes with polynomial of degree 2 and 3 . We observe that our algorithm is
able to properly locate the singularity point at the origin and the discontinuity points
(0.2, 0) and (0.8, 0). Note also that, for about the same accuracy, the use of polynomial
of degree 3 results in a mesh with 85% less triangles than the one obtained by using
polynomial of degree 2. The number of degrees of freedom is also reduced 77%, changing
from 48144 for k = 1 to 11076 for k = 3. Let us also plot the adaptive results and the
effectivity index for k = 1, 2, 3 in Fig. 4.8. We observe that Eh,|·| bounds the error
|J(u) − Jh| and the error approximation Eh stays close to the error. The convergence
rate of the error is around 2k + 1 for k = 1, 2, 3. We can also see the effectivity index
Ieff is bounded in all cases. As we refine, it is approaching to 1 for k= 1 and 2 .
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 4.7: Example of adaptive meshes for the second functional. Left: k =2, 2148
elements, the number of DOF is 48144, and error is 1.12e-10; Right: k=3, 310 elements,
the number of DOF is 11076, and error is 1.97e-10.
102
1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3
-9
-8
-7
-6
-5
-4
-3
-2
3
1
0 2 4 6 8 10 12
Mesh
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
1.5 2 2.5
-11
-10
-9
-8
-7
-6
-5
-4
-3
-2
5
1
0 2 4 6 8 10 12
Mesh
0.2
1
2
3
4
5
6
7
1.6 1.7 1.8 1.9 2 2.1 2.2 2.3
-14
-12
-10
-8
-6
-4
-2
1
7
1 2 3 4 5 6 7 8 9 10
DOF
0.1
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
I
eff
Figure 4.8: Plots of the error and the error approximation against the DOF in logarith-
mic scale for k = 1, 2, 3 (left column) and plots of the corresponding effectivity index
Ieff(right column) when approximating J2(u) =< ∇u · n, ψ >∂Ω.
103
The third functional
For the third functional
J3(u) = ((∇u,∇u)Ω − ⟨u,∇u · n⟩Ω)/(u, u)Ω,
let us consider the eigenvalue problem (4.3.3) on a L-shaped domain Ω = [0, 2]2\[1, 1]×
[2, 2]. Since Ω has a reentrant corner at the point (1, 1), the first eigenfunction is
singular. Let (u1h, q
1
h, û
1
h) be the approximations of the first eigenfunction obtained by
the HDG method and then we define Jh as in (4.3.5). Numerical results in [42, 30] show
that if we refine the mesh uniformly, the convergence rate of the first eigenvalue is only
of order 4/3.
We start with an initial mesh Th of size h0 = 0.25. We report a sample of adaptive
meshes in Fig. 4.9. We can see that the refinement is concentrated at the corner (1, 1),
where the singularity is. In addition, we observe that, for about the same error, the use
of higher degree polynomials results in much fewer elements and degrees of freedom.
For example, to achieve the error around 1.0e-3 , we need 2905 elements for k = 1 while
we only need 170 elements for k = 3, which is a remarkable 94% reduction of elements.
The number of degrees of freedom is also reduced 82%, changing from 34719 for k = 1
to 6044 for k = 3. This shows how advantageous is the use of higher-degree polynomials
with our adaptive method. Next, as shown in Fig. 4.10, for k = 1, 2, 3, we can see that
logarithm of the error |λ− λh| decreases in a linear way and that the convergence rate
is at least 2k + 1 for k = 1, 2, 3. It actually converges faster than 2k + 1 for k = 3. We
can also see that the error approximation Eh is almost parallel to the error, but Eh is
clearly smaller than the error and Eh,|·| does not bound the true error. The effectivity
index Ieff lies between 0.2 and 0.6 and it is roughly a constant for each k. This shows
that we can still use our error approximation to guide the mesh adaption even though
it may not be very close to the true error. Our study of the error approximation clearly
needs to be refined though.
104
0 0.5 1 1.5 2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0 0.5 1 1.5 2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0 0.5 1 1.5 2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Figure 4.9: Examples of adaptive meshes for the third functional. Top left: k =1,
2905 elements, the number of DOF is 34719, and error is 8.01e-04; Top right: k=2, 244
elements, the number of DOF is 5427, and error is 9.01e-04; Bottom: k=3, 170 elements,
the number of DOF is 6044 and error is 1.00e-3.
105
1.4 1.6 1.8 2 2.2 2.4 2.6
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
-1
3
1
1 2 3 4 5 6 7 8
Mesh
0.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0.6
0.65
0.7
1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 2.05
-5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
1
5
1 2 3 4 5 6 7 8
Mesh
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
1.75 1.8 1.85 1.9 1.95 2 2.05
-5.5
-5
-4.5
-4
-3.5
-3
-2.5
-2
-1.5
1
7
1 2 3 4 5 6 7 8
Mesh
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Figure 4.10: Plots of the error and the error approximation against the DOF in loga-
rithmic scale for k = 1, 2, 3 (left column) and plots of the corresponding effectivity index
Ieff (right column) when approximating the 1st eigenvalue on the L-shaped domain.
106
4.5 Concluding Remarks
In this chapter, we propose a novel error approximation for functionals approxi-
mated by the HDG method. We detailed the definition of the method in Section 4.2
and we present numerical tests in Section 4.4 to show our error approximation can be
used to drive adaptive algorithms in a efficient and robust way. Although our error
approximation depends on the HDG local post-processing technique for second-order
elliptic equations, the framework we described in Section 4.2 can be applied to other
PDEs and to other linear/nonlinear functionals if a super-convergent post-processing
method is obtained.
Incorporating the adaptive strategy proposed here into the the convolution filtering
technique developed in our previous work [27, 28, 30] is presented in the next chapter.
Chapter 5
An adjoint-based adaptive
method with convolution filtering
5.1 Introduction
In this chapter, we present a novel method that can achieve highly accurate approx-
imations even when the exact solution has corner singularities or when the given data
is discontinuous. This method takes advantage of the convolution filtering technique
of the adjoint-based method as well as the error approximation and mesh adaptation
technique.
As presented in chapter 2 and 3, under suitable regularity assumptions, the adjoint-
based method converges with order at least O(h4k) for approximating linear functionals
and eigenvalues with Galerkin method of polynomial degree k. However, we lose such
high order convergence rate when the exact solution or the functional are not smooth
enough. To overcome this difficulty, in Chapter 4, a fully computable error approx-
imation technique is is introduced to drive an adaptive strategy. As our numerical
experiments show, this allows us to handle the lack of smoothness as well as achieve
super-convergence.
The main idea of our method is in twofold. First, we use the convolution filtering
technique to enhance the accuracy of the Galerkin solution for the subdomains where the
solution is smooth enough. In these subdomains, the functional error super-converges.
Second, we use the error approximation technique to adaptively refine the elements
107
108
where the solution is not smooth. With several refinements, the functional error for the
non-smooth region is reduced to a comparable size as the error for the smooth region.
Hence the total error super-converges.
The remainder of this chapter is organized as follows. In Section 5.2, we introduce
our method. Next, in Section 5.3, we present numerical results to show the 4k super-
convence of our method. Last, in Section 5.4, we conclude this chapter with a discussion
of future works.
5.2 The method
To illustrate the method, let us consider a linear functional J(u) that is determined
by the solution u of the model problem (4.2.1). Let us begin by describing the assump-
tions on the meshes.
Assumption 6. Let {Th}h>0 be a family of shape-regular meshes of Ω. There exist
subdomains Ω0 ⊂⊂ Ω′0 ⊂ Ω such that
(1) Th ∩ Ω′0 is a translation-invariant mesh with element size h,
(2) The solution u is smooth in the subdomain Ω′0,
(3) Any element K ∈ Th is fully contained in either Ω0 or Ω1 := Ω\Ω0.
The assumption (1) and (2) make sure that we can use the convolution filtering
technique to enhance the numerical solution in the subdomain Ω0. Let us denote by
T0,h := Th ∩ Ω0. For the subdomain Ω1, since we don’t apply the convolution filtering
technique there, we can simply take T1,h(0) := Th ∩ Ω1 or we can define a different
(unstructured) mesh T1,h(0) with size h on Ω1 as long as assumption (3) is satisfied.
Here the superscript 0 stands for the initial outer mesh without any adaptivity.
Now we are ready to define our adaptive method. We follow the idea of the adjoint-
based method presented in Chapter 2 but for the subdomain Ω1, we use the adaptive
approach described in Chapter 4. In our numerical results, we observe that when we
take the adaptive steps M to be 6 to 8 steps, we usually obtain the best results.
109
Algorithm 3: The adjoint-based adaptive algorithm with convolution filtering
(1) Construct an initial mesh Th.
(2) Compute the HDG solutions (uh, qh, ûh) and (vh,ph, v̂h) for model
problem (4.2.3) and the adjoint problem (4.2.4), respectively, with polynomial
degree k as discussed in Section 4.2.1.
(3) Compute the accuracy-enhanced solutions u∗h and v
∗
h in Ω0 using the
convolution filtering technique in Section 2.2.
(4) Apply the adaptive algorithm 2 described in Section 4.2.3 in Ω1 with
initial mesh T1,h(0) and polynomial degree 2k for M steps. We obtain u2kh and
v2kh on mesh T1,h(M).
(5) Compute Jh(u
∗
h, v
∗
h) = J(u
∗
h) +ACh as the adjoint-based approximation
to J(u), where
u∗h :=
K2kh ∗ ukh in K ∈ T0,h,u2kh in K ∈ T1,h(M),
and
v∗h :=
K2kh ∗ vkh in K ∈ T0,h,v2kh in K ∈ T1,h(M).
Let us end this section by remarking that our method assumes that we have priori
knowledge of the solution so that we can choose suitable subdomain Ω0 to satisfy item
(2) in Assumption 6. For cases when we do not know the behavior of the solution,
further studies are needed to determine how to choose subdomain Ω0.
5.3 Numerical results
In this section, we present numerical examples to show the advantage of our adaptive
method. We consider the model problem (4.2.1) on the unit square Ω = [0, 1] × [0, 1].
The functional J(u) we consider is defined with non-smooth solution or discontinuous
data. In the adaptive method, the subdomain Ω0 is chosen to be 2k elements away from
the boundary and 2k + 1 elements away from the non-smooth region. The adaptive
mesh refinement strategy we use here is the same as the one we use in Section 4.4.2.
110
Let us denote by eHDG := |J(u)− J(uh)−ACh| and eAdap := |J(u)− J(u∗h)−AC∗h| the
functional error obtained by the HDG method and our adaptive method, respectively.
Similarly, let tHDG and tAdap be the CPU time needed to approximate the functional
using the two methods, and let NHDG and NAdap be the number of degrees of freedom
(DOFs) of the two methods.
5.3.1 Volume integral with corner singularity
The first example is to approximate a volume integral
J(u) =
∫
Ω
u g dx,
where u = r2/3 sin(23θ) in polar coordinates, g = 1. Note that u is singular at the origin
(0, 0). More specifically, u ∈ H5/3−ϵ(Ω) for any ϵ > 0. In our numerical experiment,
we consider 3 meshes with initial mesh size h0 = 0.1, 0.05, 0.025, respectively. For the
adaptive method, we choose the number of adaptive steps M = 6 in Algorithm 3. In
Figure 5.1, we show the last adaptive mesh T1,h(6) for k = 1 and k = 2 with the initial
mesh size h0 = 0.05. We can see that the algorithm captures the singularity very well
and the refinement is concentrated at the origin. In Table 5.1 and Table 5.2, we report
a comparison of the functional error, the CPU time, and the number of degrees of
freedom for the two methods. We can observe that the error from the adaptive method
is significantly smaller than the error obtained from the HDG method. Although the
CPU time needed for the adaptive method is 10 times longer than the HDG method
since the algorithm needs to adaptively solve the PDEs on Ω1 several times, the error for
the adaptive method is 2 to 5 order of magnitude smaller than the error from the HDG
method. We can also see that as the mesh gets finer (so the area of Ω1 gets smaller),
the CPU time ratio tAdap/tHDG also gets smaller. These results indicate that it is worth
using the adaptive method to obtain more accurate functional approximation.
111
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
k=1
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
k=2
Figure 5.1: Example of adaptive meshes when u is singular with the initial mesh size
h0 = 0.05. Left: k = 1; Right: k = 2.
k Mesh eHDG eAdap eHDG/eAdap
1
1 1.25E-03 1.11E-07 1.13E+04
2 1.49E-05 3.36E-08 4.43E+02
3 2.77E-05 4.58E-09 6.05E+03
2
1 6.86E-04 3.25E-09 2.11E+05
2 1.16E-04 7.61E-10 1.52E+05
3 1.89E-05 1.58E-10 1.20E+05
3
1 - - -
2 2.85E-05 8.98E-10 3.17E+04
3 4.55E-06 1.41E-10 3.23E+04
Table 5.1: Functional error comparison of the HDG method and adjoint-based adaptive
method for the first functional.
112
k Mesh tAdap/tHDG NAdap/NHDG
1
1 36.2 5.0
2 15.9 2.0
3 6.6 1.3
2
1 42.3 2.5
2 29.7 1.9
3 11.8 1.5
3
1 - -
2 55.6 2.4
3 20.4 1.8
Table 5.2: The CPU time and DOFs comparison of the two methods for the first
functional.
5.3.2 Boundary integral with discontinuity
The second example we consider here is a boundary integral
J(u) =< ∇u · n, ψ >∂Ω,
where u(x, y) = sin(πx) sin(πy) is smooth but ψ is chosen to be a discontinuous function
on ∂Ω:
ψ =
(x− 0.2)(x− 0.8) + 1 if 0.2 < x < 0.8 and y = 0,0 otherwise.
In our numerical experiment, we consider 3 meshes with initial mesh size h0 = 0.1, 0.05, 0.025,
respectively. For the adaptive method, we choose the number of adaptive steps M = 6
in Algorithm 3. In Figure 5.2, we show the last adaptive mesh T1,h(6) for k = 1 and
k = 2. We can see that the algorithm captures the two discontinuities well and the
adaptive refinement is concentrated at these two points. We report a comparison of the
functional error, the CPU time, and the number of DOFs for the two methods in Table
5.2 and Table 5.3. We observe similar results as we see in the previous example. Even
though the adaptive method takes 10s times of CPU time, the error of the adaptive
method is reduced by hundreds, even thousands of times compared the HDG method.
113
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
k=1
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
k=2
Figure 5.2: Example of adaptive meshes when ψ is discontinuous with the initial mesh
size h0 = 0.05. Left: k = 1; Right: k = 2.
k Mesh eHDG eAdap eHDG/eAdap
1
1 1.48E-02 2.72E-04 5.44E+01
2 1.69E-03 1.09E-06 1.55E+03
3 2.00E-04 9.79E-08 2.04E+03
2
1 5.89E-05 2.31E-08 2.55E+03
2 2.00E-06 5.20E-09 3.85E+02
Table 5.3: Functional error comparison of the HDG method and adjoint-based adaptive
method for the second functional.
k Mesh tAdap/tHDG NAdap/NHDG
1
1 46.5 14.3
2 29.8 6.4
3 10.3 2.8
2
1 58.8 8.0
2 34.4 4.0
Table 5.4: The CPU time and DOFs comparison of the two methods for the second
functional
114
5.4 Concluding Remarks
In this chapter, we present a novel method that combines the convolution filtering
technique and the error approximation and mesh adaptation technique studied in pre-
vious chapters. Our numerical results show that for non-smooth problems, our new
method can help enhance the accuracy by hundreds, even thousands of times compared
to the standard Galerkin methods.
Chapter 6
Conclusion and Discussion
We study adjoint-based methods for approximating functionals in a highly accurate
manner. We rigorously prove the convergence of the adjoint-based method proposed by
Cockburn and Wang in [27] under suitable regularity assumptions. We then extend this
method to a nonlinear functional problem, the eigenvalue problem. To handle the lack
of smoothness of the solution, we further develop an adjoint-based error approximation
and mesh adaptation algorithm. Let us discuss three lines of possible future work.
The first one is to devise a mechanism that identifies the smooth region Ω0 without
priori knowledge of the solution. In Algorithm 3, we assume that we can choose a
suitable subdomain Ω0 to satisfy Assumption 6 based on our priori knowledge of the
solution. But this will not be the case for more complicated problems. So we need a
practical algorithm to choose the subdomain Ω0 wisely.
The second direction is to determine the number of adaptive steps automatically
instead of manually picking a number. In our numerical tests, we observe that after
certain number of adaptive steps, the error at the subdomain Ω1 will be small enough so
that the error at the inner subdomain Ω0 will dominate. In such case, further refinement
will not reduce the whole functional error. Besides, if the solution is smooth, there is
no need for adaptive refinements. Therefore, an intelligent algorithm to determine the
sufficient number of adaptive steps will be helpful to avoid unnecessary computations.
The last direction is to bypass the translation-invariance requirement of the meshes
needed for the convolution filtering technique. This is a difficult task, but it would allow
the methods studied here to be applied to general meshes.
115
References
[1] L. Advanpix, Multiprecision computing toolbox for matlab, YoNohama, Japan,
(2019). https://www.advanpix.com.
[2] P. F. Antonietti, A. Buffa, and I. Perugia, Discontinuous Galerkin approx-
imation of the laplace eigenproblem, Comput. Methods Appl. Mech. Engrg., 195
(2006), pp. 3483–3503.
[3] I. Babusˇka and J. Osborn, Eigenvalue problems, Handbook of numerical anal-
ysis, 2 (1991), pp. 641–787.
[4] W. Bangerth and R. Rannacher, Adaptive finite element methods for differ-
ential equations, Birkha¨user, 2013.
[5] R. Becker and R. Rannacher, An optimal control approach to a posteriori
error estimation in finite element methods, Acta Numerica, 10 (2001), pp. 1–102.
[6] R. Bermejo and J. Carpio, A space-time adaptive finite element algorithm based
on dual weighted residual methodology for parabolic equations, SIAM J. Sci. Com-
put., 31 (2009), pp. 3324–3355.
[7] T. Betcke and L. N. Trefethen, Reviving the method of particular solutions,
SIAM review, 47 (2005), pp. 469–491.
[8] D. Boffi, F. Brezzi, and L. Gastaldi, On the problem of spurious eigenvalues
in the approximation of linear elliptic problems in mixed form, Math. Comp., 69
(2000), pp. 121–140.
116
117
[9] J. Bramble and A. Schatz, Higher order local accuracy by averaging in the finite
element method, Math. Comp., 31 (1977), pp. 94–111.
[10] F. Brezzi, J. Douglas, and L. D. Marini, Two families of mixed finite elements
for second order elliptic problems, Numer. Math., 47 (1985), pp. 217–235.
[11] A. Buffa and I. Perugia, Discontinuous Galerkin approximation of the maxwell
eigenproblem, SIAM J. Numer. Anal., 44 (2006), pp. 2198–2226.
[12] A. Cangiani, Z. Dong, and E. H. Georgoulis, hp-version discontinu-
ous Galerkin methods on essentially arbitrarily-shaped elements, arXiv preprint
arXiv:1906.01715, (2019).
[13] A. Cangiani, E. Georgoulis, and Y. Sabawi, Adaptive discontinuous Galerkin
methods for elliptic interface problems, Math. Comp., 87 (2018), pp. 2675–2707.
[14] J. Carpio, J. L. Prieto, and R. Bermejo, Anisotropic “goal-oriented” mesh
adaptivity for elliptic problems, SIAM J. Sci. Comput., 35 (2013), pp. A861–A885.
[15] H. A. Carson, D. L. Darmofal, M. C. Galbraith, and S. R. Allmaras,
Analysis of output-based error estimation for finite element methods, Appl. Numer.
Math., 118 (2017), pp. 182–202.
[16] C. Carstensen, An adaptive mesh-refining algorithm allowing for an h1 stable
l2 projection onto courant finite element spaces, Constructive Approximation, 20
(2004), pp. 549–564.
[17] B. Cockburn and G. Fu, Superconvergence by M-decompositions. Part II: Con-
struction of two-dimensional finite elements, ESAIM Math. Model. Numer. Anal.,
51 (2017), pp. 165–186.
[18] , Superconvergence by M-decompositions. Part III: Construction of three-
dimensional finite elements, ESAIM Math. Model. Numer. Anal., 51 (2017),
pp. 365–398.
[19] B. Cockburn, G. Fu, and F. Sayas, Superconvergence by M-decompositions.
Part I: General theory for HDG methods for diffusion, Math. Comp., 86 (2017),
pp. 1609–1641.
118
[20] B. Cockburn and J. Gopalakrishnan, A characterization of hybridized mixed
methods for second order elliptic problems, SIAM J. Numer. Anal., 42 (2004),
pp. 283–301.
[21] B. Cockburn, J. Gopalakrishnan, and R. Lazarov, Unified hybridization of
discontinuous Galerkin, mixed, and continuous Galerkin methods for second order
elliptic problems, SIAM J. Numer. Anal., 47 (2009), pp. 1319–1365.
[22] B. Cockburn, J. Gopalakrishnan, F. Li, N.-C. Nguyen, and J. Peraire,
Hybridization and postprocessing techniques for mixed eigenfunctions, SIAM J. Nu-
mer. Anal., 48 (2010), pp. 857–881.
[23] B. Cockburn, J. Gopalakrishnan, and F.-J. Sayas, A projection-based error
analysis of HDG methods, Math. Comp., 79 (2010), pp. 1351–1367.
[24] B. Cockburn, J. Guzma´n, and H. Wang, Superconvergent discontinuous
Galerkin methods for second-order elliptic problems, Math. Comp., 78 (2009), pp. 1–
24.
[25] B. Cockburn, M. Luskin, C.-W. Shu, and E. Su¨li, Enhanced accuracy by
post-processing for finite element methods for hyperbolic equations, Math. Comp.,
72 (2003), pp. 577–606.
[26] B. Cockburn, W. Qiu, and K. Shi, Conditions for superconvergence of HDG
methods for second-order elliptic problems, Math. Comp., 81 (2012), pp. 1327–1353.
[27] B. Cockburn and Z. Wang, Adjoint-based, superconvergent Galerkin approxi-
mations of linear functionals, J. Sci. Comput., 73 (2017), pp. 644–666.
[28] B. Cockburn and S. Xia, An a priori error analysis of adjoint-based super-
convergent Galerkin approximations of linear functionals, IMA J. Num. Anal.,
(2021).
[29] , An adjoint-based adaptive error approximation of functionals by the hybridiz-
able discontinuous Galerkin method for second-order elliptic equations, J. Comput.
Phys., 457 (2022), p. 111078.
119
[30] , An adjoint-based super-convergent Galerkin approximation of eigenvalues, J.
Comput. Phys., 449 (2022), p. 110816.
[31] S. Curtis, R. M. Kirby, J. K. Ryan, and C.-W. Shu, Postprocessing for the
discontinuous Galerkin method over nonuniform meshes, SIAM J. Sci. Comput.,
30 (2007), pp. 272–289.
[32] V. Dolejsi, G. May, A. Rangarajan, and F. Roskovec, A goal-oriented
high-order anisotropic mesh adaptation using discontinuous Galerkin method for
linear convection-diffusion-reaction problems, SIAM J. Sci. Comput., 41 (2019),
pp. A1899–A1922.
[33] W. Do¨rfler, A convergent adaptive algorithm for poisson’s equation, SIAM J.
Numer. Anal., 33 (1996), pp. 1106–1124.
[34] B. Endtmayer, U. Langer, and T. Wick, Two-side a posteriori error estimates
for the dual-weighted residual method, SIAM J. Sci. Comput., 42 (2020), pp. A371–
A394.
[35] K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Introduction to adaptive
methods for differential equations, Acta Numerica, 4 (1995), pp. 105–158.
[36] K. J. Fidkowski and G. Chen, Output-based mesh optimization for hybridized
and embedded discontinuous Galerkin methods, Internat. J. Numer. Methods En-
grg., 121 (2020), pp. 867–887.
[37] K. J. Fidkowski and D. L. Darmofal, Review of output-based error estimation
and mesh adaptation in computational fluid dynamics, AIAA journal, 49 (2011),
pp. 673–694.
[38] S. A. Funken and A. Schmidt, Adaptive mesh refinement in 2d – an efficient
implementation in matlab, Comput. Methods Appl. Math., 20 (2020), pp. 459–479.
[39] F. Gardini, Mixed approximation of eigenvalue problems: a superconvergence re-
sult, Mode´l. Math. Anal. Nume´r., 43 (2009), pp. 853–865.
[40] M. Giles, M. G. Larson, J. M. Levenstam, and E. Su¨li, Adaptive error
control for finite element approximations of the lift and drag coefficients in viscous
120
flow, Tech. Rep. NA-97/06, Oxford University Computing Laboratory, Wolfson
Building, Parks Road, Oxford OX1 3QD, 1997.
[41] M. B. Giles and E. Su¨li, Adjoint methods for PDEs: a posteriori error analysis
and postprocessing by duality, Acta Numerica, 11 (2002), pp. 145–236.
[42] J. Gopalakrishnan, F. Li, N.-C. Nguyen, and J. Peraire, Spectral approx-
imations by the HDG method, Math. Comp., 84 (2015), pp. 1037–1059.
[43] R. Hartmann, Multitarget error estimation and adaptivity in aerodynamic flow
simulations, SIAM J. Sci. Comput., 31 (2008), pp. 708–731.
[44] R. Hartmann and P. Houston, Adaptive discontinuous Galerkin finite ele-
ment methods for the compressible Euler equations, J. Comput. Phys., 183 (2002),
pp. 508–532.
[45] , Symmetric interior penalty DG methods for the compressible Navier–Stokes
equations ii: Goal–oriented a posteriori error estimation, Internat. J. Numer. Anal.
Model., 3 (2006), pp. 141–162.
[46] R. Hartmann and P. Houston, Error estimation and adaptive mesh refine-
ment for aerodynamic flows, in VKI LS 2010-01: 36th CFD/ADIGMA course on
hp-adaptive and hp-multigrid methods, Oct. 26-30, 2009, H. Deconinck, ed., Von
Karman Institute for Fluid Dynamics, Rhode Saint Gene`se, Belgium, 2009.
[47] J. S. Hesthaven and T. Warburton, High–order nodal discontinuous Galerkin
methods for the maxwell eigenvalue problem, Philos. Trans. Roy. Soc. Lond. Ser. A:
Math. Phys. Eng. Sci., 362 (2004), pp. 493–524.
[48] V. Heuveline and R. Rannacher, A posteriori error control for finite element
approximations of elliptic eigenvalue problems, Adv. Comput. Math., 15 (2001),
pp. 107–138.
[49] L. Ji, P. Van Slingerland, J. K. Ryan, and K. Vuik, Superconvergent
error estimates for position-dependent smoothness-increasing accuracy-conserving
(siac) post-processing of discontinuous Galerkin solutions, Math. Comp., (2014),
pp. 2239–2262.
121
[50] J. King, H. Mirzaee, J. K. Ryan, and R. M. Kirby, Smoothness-increasing
accuracy-conserving (siac) filtering for discontinuous Galerkin solutions: improved
errors versus higher-order accuracy, J. Sci. Comput., 53 (2012), pp. 129–149.
[51] X. Li, J. K. Ryan, R. M. Kirby, and K. Vuik, Smoothness-increasing accuracy-
conserving (siac) filtering for discontinuous Galerkin solutions over nonuniform
meshes: Superconvergence and optimal accuracy, J. Sci. Comput., (2019), pp. 1–
31.
[52] A. Loseille, A. Dervieux, and F. Alauzet, Fully anisotropic goal-oriented
mesh adaptation for 3d steady Euler equations, J. Comput. Phys., 229 (2010),
pp. 2866–2897.
[53] L. Machiels, J. Peraire, and A. Patera, A posteriori finite-element output
bounds for the incompressible Navier–Stokes equations: application to a natural
convection problem, J. Comput. Phys., 172 (2001), pp. 401–425.
[54] B. Mercier, J. Osborn, J. Rappaz, and P.-A. Raviart, Eigenvalue approxi-
mation by mixed and hybrid methods, Math. Comp., 36 (1981), pp. 427–453.
[55] H. Mirzaee, L. Ji, J. K. Ryan, and R. M. Kirby, Smoothness-increasing
accuracy-conserving (siac) postprocessing for discontinuous Galerkin solutions over
structured triangular meshes, SIAM J. Numer. Anal., 49 (2011), pp. 1899–1920.
[56] H. Mirzaee, J. King, J. K. Ryan, and R. M. Kirby, Smoothness-increasing
accuracy-conserving filters for discontinuous Galerkin solutions over unstructured
triangular meshes, SIAM J. Sci. Comput., 35 (2013), pp. A212–A230.
[57] H. Mirzaee, J. K. Ryan, and R. M. Kirby, Efficient implementation of
smoothness-increasing accuracy-conserving (siac) filters for discontinuous Galerkin
solutions, J. Sci. Comput., 52 (2012), pp. 85–112.
[58] , Smoothness-increasing accuracy-conserving (siac) filters for discontinuous
Galerkin solutions: Application to structured tetrahedral meshes, J. Sci. Comput.,
58 (2014), pp. 690–704.
122
[59] P. Monk and E. Su¨li, The adaptive computation of far-field patterns by a pos-
teriori error estimation of linear functionals, SIAM J. Numer. Anal., 36 (1998),
pp. 251–274.
[60] J. A. Nitsche and A. H. Schatz, Interior estimates for Ritz-Galerkin methods,
Math. Comp., 28 (1974), pp. 937–958.
[61] M. Paraschivoiu, J. Peraire, and A. T. Patera, A posteriori finite element
bounds for linear-functional outputs of elliptic partial differential equations, Com-
put. Methods Appl. Mech. Engrg., 150 (1997), pp. 289–312.
[62] N. A. Pierce and M. B. Giles, Adjoint recovery of superconvergent functionals
from PDE approximations, SIAM Review, 42 (2000), pp. 247–264.
[63] , Adjoint and defect error bounding and correction for functional estimates, J.
Comput. Phys., 200 (2004), pp. 769–794.
[64] J. Ryan and C.-W. Shu, On a one-sided post-processing technique for the dis-
continuous Galerkin methods, Methods and Applications of Analysis, 10 (2003),
pp. 295–308.
[65] J. K. Ryan and B. Cockburn, Local derivative post-processing for the discon-
tinuous Galerkin method, J. Comput. Phys., 228 (2009), pp. 8642–8664.
[66] J. K. Ryan, X. Li, R. M. Kirby, and K. Vuik, One-sided position-dependent
smoothness-increasing accuracy-conserving (siac) filtering over uniform and non-
uniform meshes, J. Sci. Comput., 64 (2015), pp. 773–817.
[67] J. Sarrate, J. Peraire, and A. Patera, A posteriori finite element error
bounds for non-linear outputs of the helmholtz equation, Internat. J. Numer. Meth-
ods Fluids, 31 (1999), pp. 17–36.
[68] M. Sharbatdar and C. Ollivier-Gooch, Adjoint-based functional correction
for unstructured mesh finite volume methods, J. Sci. Comput., 76 (2018), pp. 1–23.
[69] R. Stenberg, Postprocessing schemes for some mixed finite elements, RAIRO
Mode´l. Math. Anal. Nume´r., 25 (1991), pp. 151–167.
123
[70] V. Thome´e, High order local approximations to derivatives in the finite element
method, Math. Comp., 31 (1977), pp. 652–660.
[71] P. van Slingerland, J. K. Ryan, and C. Vuik, Position-dependent
smoothness-increasing accuracy-conserving (siac) filtering for improving discontin-
uous Galerkin solutions, SIAM J. Sci. Comput., 33 (2011), pp. 802–825.
[72] D. A. Venditti and D. L. Darmofal, Grid adaptation for functional out-
puts: Application to two-dimensional inviscid flows, J. Comput. Phys., 176 (2002),
pp. 40–69.
[73] D. A. Venditti and D. L. Darmofal, Anisotropic grid adaptation for func-
tional outputs: application to two-dimensional viscous flows, J. Comput. Phys.,
187 (2003), pp. 22–46.
[74] M. Woopen, G. May, and J. Schu¨tz, Adjoint-based error estimation and mesh
adaptation for hybridized discontinuous Galerkin methods, Internat. J. Numer.
Methods Fluids, 76 (2014), pp. 811–834.
Appendix A
Supplements
A.1 Special projections and their properties
A.1.1 The M-decompositions
The concept of M -decomposition was introduced in [19] to systematically construct
HDG and hybridized mixed methods with superconvergence properties on unstructured
meshes. In what follows, for each element K, we consider local finite dimensional spaces
V ⊂H(div,K), W ⊂ H1(K) and M ⊂ L2(∂K).
Definition 2. We say V ×W admits an M -decomposition if
1. tr(V ×W ) ⊂M ,
and there exists a subspace V˜ × W˜ ⊂ V ×W such that
2. ∇W ×∇ · V ⊂ V˜ × W˜ ,
3. tr : V˜
⊥ × W˜⊥ →M is an isomorphism.
Here V˜
⊥
and W˜⊥ are the L2-orthogonal complements of V˜ in V and W˜ in W , respec-
tively, and tr is the trace operator defined as follows
tr : V ×W → L2(∂K),
(v, w) 7→ (v · n+ w)|∂K .
124
125
Note that the choice of V˜ × W˜ is not unique. The canonical (orthogonal) M -
decomposition is
V˜ = ∇W ⊕ V sbb, and W˜ = ∇ · V ,
where V sbb := {v ∈ V : ∇ · v = 0, v · n = 0}.
A.1.2 The HDG Projection
In this section we define the HDG projection and state its approximation properties.
Since we deal with a general domain Ω, we may have two types of elements. The first
type is a regular polygonal or polyhedral element and the second type is a boundary
polygonal or polyhedral element with one curved face. Let us start with a regular
polygonal or polyhedral element K and define the standard HDG projection.
Definition 3. Let (q, u) be smooth enough so that their boundary traces are in L2(∂K).
Assume that local spaces V ×W admits an M -decomposition on K. Then the HDG
projection Πh(q, u) := (ΠV q, ΠWu) is defined by the following equations:
(u−ΠWu,w)K = 0 ∀ w ∈ W˜ , (A.1.1a)
(q −ΠV q,v)K = 0 ∀ v ∈ V˜ , (A.1.1b)
⟨(q −ΠV q) · n+ τ(u−ΠWu) , µ⟩∂K = 0 ∀ µ ∈M. (A.1.1c)
Next, we define the HDG projection for elements that have one curved face.
Definition 4. Let K be a polygonal or polyhedral element with one curved face E.
With the same assumption in Definition 3, we define the HDG projection Πh(q, u) :=
(ΠV q, ΠWu) by the following equations:
(u−ΠWu,w)K = 0 ∀ w ∈ W˜ , (A.1.2a)
(q −ΠV q,v)K = 0 ∀ v ∈ S, (A.1.2b)
⟨(q −ΠV q) · n+ τ(u−ΠWu) , µ⟩F = 0 ∀ µ ∈M(F ) F ̸= E (A.1.2c)
(∇ · (q −ΠV q), w) + ⟨τ(u−ΠWu) , w⟩∂K = 0 ∀ w ∈W, (A.1.2d)
126
where S := S1 ⊕ S2, S1 = ∇W˜⊥ and
S2 = {v ∈ V |∇ · v = 0 and v · nF = 0 on F ̸= E and (v,p) = 0 ∀p ∈ S1}.
For both cases, the HDG projection depends on the stabilization parameter τ . Let
us introduce the quantity
τ
W˜⊥ :=
infµ∈W˜⊥\{0}
⟨τµ,µ⟩∂K
∥µ∥2∂K
if W˜⊥ ̸= 0,
∞ if W˜⊥ = 0.
It can be shown that, if τ
W˜⊥ > 0, then the HDG projections (A.1.1) and (A.1.2) are
well-defined. Let us show the approximation properties for the HDG projections. T he
proof for the regular element case is given in [19] and the proof for the curved element
case is displayed in Appendix A.1.3.
Theorem A.1. Suppose V ×W admits an M -decomposition and let the stabilization
function τ satisfy the condition
τ
W˜⊥ > 0.
Let the element K satisfy Assumption 4. Then the HDG projection given by (A.1.1) or
(A.1.2) is well-defined and has the following approximation properties
∥q −ΠV q∥K ≤ ∥(Id− PV )q∥K + Ch
1
2
K∥(Id− PV )q · n∥∂K +Θ,
∥u−ΠWu∥K ≤ ∥(Id− PW )u∥K +Θ,
where Θ = Ch
1
2
K∥(Id− PW )u∥∂K +ChK∥(Id− PW˜ )∇ · q∥K and P is the L2-projection
operator.
A.1.3 The HDG projection for curved elements
In this section, we prove that the HDG projection for curved elements is well-defined
and has the optimal approximation properties as stated in Theorem A.1.
For anyK ∈ Th, let the local space V (K)×W (K) admit anM(∂K)−decomposition.
If K is a polygonal or polyhedral element with one curved face E, we define the HDG
projection by equations (A.1.2).
127
We assume that the stabilization function τ is nonnegative on each face F ∈ ∂K
such that
w ∈ W˜⊥ : ⟨τ(w) , w⟩∂K = 0 =⇒ w = 0. (A.1.3)
We can see that τ
W˜⊥ > 0 implies (A.1.3).
Remark 1. If W˜ =W (K) and τ = 0, the projection becomes:
(ΠWu,w)K = (u,w)K ∀ w ∈W (K),
(∇ ·ΠV q, w) = (∇ · q, w) ∀ w ∈W (K),
⟨ΠV q · n , µ⟩F = ⟨q · n , µ⟩F ∀ µ ∈M(F ) F ̸= E
(ΠV q,v)K = (q,v)K ∀ v ∈ S.
Here S := {v ∈ V (K)|∇ · v = 0, v · nF = 0 on F ̸= E}. This is the same projection
proposed in [10] for the BDM mixed method.
Remark 2. In general, even if E is flat, the projection is different from the standard
HDG projection. However, if the local spaces satisfy the assumption that W˜⊥|E =
M(E), then the projection is equivalent to the standard HDG projection. For example,
it is easy to show that for simplicial elements, the projection is the same as the standard
HDG projection. In this case, we choose V (K) := [Pk(K)]d, W (K) := Pk(K) and
M(F ) := PF (K) for all faces F ∈ ∂K and V˜ := Pk−1(K) and W˜ := Pk−1(K). Then,
the assumption W˜⊥|E = P⊥k (K)|E = Pk(E) =M(E) is satisfied.
Let us write the projection (A.1.2) in a more detailed manner:
(u−ΠWu,w)K = 0 ∀ w ∈ W˜ , (A.1.4a)
(q −ΠV q,∇w)K = 0 ∀ w ∈ W˜⊥, (A.1.4b)
(q −ΠV q,v)K = 0 ∀ v ∈ S2, (A.1.4c)
⟨(q −ΠV q) · n+ τ(u−ΠWu) , µ⟩F = 0 ∀ µ ∈M(F ) F ̸= E, (A.1.4d)
(∇ · (q −ΠV q), w) + ⟨τ(u−ΠWu) , w⟩∂K = 0 ∀ w ∈W (K). (A.1.4e)
128
Well-posedness
Theorem A.2. Suppose V ×W admits an M−decomposition and let the stabilization
function τ satisfy the condition (A.1.3). Then the HDG projection (A.1.4) is well-
defined.
Proof. First notice that we can decouple the projection component ΠWu. If we take
w ∈ W˜⊥ in (A.1.4e), then the equations (A.1.4a) and (A.1.4e) read
(u−ΠWu,w)K = 0 ∀ w ∈ W˜ , (A.1.5a)
⟨τ(u−ΠWu) , w⟩∂K = −(∇ · q, w) ∀ w ∈ W˜⊥. (A.1.5b)
Clearly, this is a square system for ΠWu. So to show the existence and uniqueness, we
only need to show that, if q = 0 and u = 0, then ΠWu = 0. From the first equation,
we know that ΠWu ∈ W˜⊥ and from the second, that ΠWu = 0 on ∂K. Since V ×W
admits an M-decomposition, we also know that ∇ · (∇ΠWu) ∈ W˜ . This implies that
(∇ΠWu,∇ΠWu)K = ⟨ΠWu,∇ΠWu · n⟩∂K − (ΠWu,∇ · (∇ΠWu)))K = 0.
Therefore, ΠWu is a constant and vanishes on ∂K. So ΠWu = 0.
Next, let us consider the component ΠV q. The equations (A.1.4e), (A.1.4d),
(A.1.4b) and (A.1.4c), read
(∇ · (q −ΠV q), w) + ⟨τ(u−ΠWu) , w⟩∂K = 0 ∀ w ∈ W˜ , (A.1.6a)
⟨(q −ΠV q) · n+ τ(u−ΠWu) , µ⟩F = 0 ∀ µ ∈M(F ) F ̸= F0 (A.1.6b)
(q −ΠV q,∇w)K = 0 ∀ w ∈ W˜⊥. (A.1.6c)
(q −ΠV q, v)K = 0 ∀ v ∈ S2. (A.1.6d)
Using the fact that ∇ · V = W˜ , we can easily see this is a square system for ΠV q.
129
Now set q = 0 and u = 0. Since we already proved that ΠWu = 0, we get that
(∇ ·ΠV q, w) = 0 ∀ w ∈ W˜ ,
⟨ΠV q · n , µ⟩F = 0 ∀ µ ∈M(F ) F ̸= F0
(ΠV q,∇w)K = 0 ∀ w ∈ W˜⊥.
(ΠV q, v)K = 0 ∀ v ∈ S2.
Since, by definition, S2 := {v ∈ V (K)|∇ · v = 0, v · nF = 0 on F ̸= F0 and (v,∇w) =
0 ∀w ∈ W˜⊥}, we can see that ΠV q = 0. Therefore, the projection is well-posed.
Approximation properties
Let us first define the following constants
∥τ∥ = sup
λ,µ∈L2(∂K)
⟨τ(λ) , µ⟩∂K
∥λ∥∂K∥µ∥∂K ,
C
W˜⊥ = sup
0̸=w∈W˜⊥
h
− 1
2
K ∥w∥K
∥w∥∂K .
We are now ready to show the proof of Theorem A.1 for the HDG projection for
curved elements.
Proof. Let us first prove the estimate for ΠWu. Let PWu be the standard L
2-orthogonal
projection of u on W . Since
∥u−ΠWu∥K ≤ ∥(Id− PW )u∥K + ∥PWu−ΠWu∥K ,
we only need to estimate δk := PWu−ΠWu.
By the first of the equations (A.1.5), (δk, w)K = (PWu−ΠWu,w)K = (u−ΠWu)k =
0 ∀w ∈ W˜ , and so δk ∈ W˜⊥. This means that we have
∥δk∥K ≤ CW˜⊥h
1
2
K∥δk∥∂K ,
130
and that we can take w := δk in the second of the equations (A.1.5) to get
⟨τδk , δk⟩∂K = (∇ · q, δk)K + ⟨τ(u− PWu) , δk⟩∂K .
Then,
∥δk∥2∂K ≤
1
τ
W˜⊥
⟨τδk , δk⟩∂K
≤ 1
τ
W˜⊥
(∥∇ · (q)− P
W˜
∇ · q∥K∥δk∥K + ∥τ∥∥u− PWu∥∂K∥δk∥∂K)
≤ CW˜⊥
τ
W˜⊥
h
1
2
K∥(Id− PW˜ )∇ · q∥K∥δk∥∂K +
∥τ∥
τ
W˜⊥
∥u− PWu∥∂K∥δk∥∂K ,
where P
W˜
is the L2 projection onto W˜ . Thus
∥δk∥∂K ≤
C
W˜⊥
τ
W˜⊥
h
1
2
K∥(Id− PW˜ )∇ · q∥K +
∥τ∥
τ
W˜⊥
∥u− PWu∥∂K , (A.1.7)
and
∥δk∥K ≤
C2
W˜⊥
τ
W˜⊥
hK∥(Id− PW˜ )∇ · q∥K +
C
W˜⊥∥τ∥
τ
W˜⊥
h
1
2
K∥u− PWu∥∂K
This completes the proof of the second inequality.
Now let us prove the first inequality. Let PV q be the standard L
2-orthogonal pro-
jection of q on V . Then
∥q −ΠV q∥K ≤ ∥(Id− PV )q∥K + ∥PV q −ΠV q∥K .
To estimate the second term, let us define norm
|||v|||K := ∥PS1(v)∥K + ∥PS2(v)∥K + hK∥∇ · v∥K +
∑
F ̸=E
h
1
2
K∥v · n∥F .
By the definition of space S1 and S2, it is easy to see that |||·||| is a norm on Pk(K). Since
all norms are equivalent in finite dimensional spaces, we have that ∥v∥K ≤ C|||v|||K
131
for any v ∈ V . By a simple scaling argument, we know C > 0 is independent of hK .
Therefore, we only need to estimate |||ΠV q − PV q|||K .
If we set ξk := ΠV q − PV q, the equations (A.1.6) can be written as follows:
(∇ · ξk, w)K = (∇ · (q − PV q), w)K + ⟨τ(u−ΠWu) , w⟩∂K
= ⟨(q − PV q) · n , w⟩∂K + ⟨τ(u−ΠWu) , w⟩∂K ∀ w ∈W
⟨ξk · n , µ⟩F = ⟨(q − PV q) · n+ τ(u−ΠWu) , µ⟩F ∀ µ ∈ Pk(F ) F ̸= E,
(ξk,v)K = 0 ∀ v ∈ S = S1 ∪ S2.
Taking w := ∇ · ξk in the first equation, µ := ξk ·n in the second, and using an inverse
estimate, we get
∥∇ · ξk∥K ≤ Ch−
1
2
K ∥(q − PV q) · n∥∂K + Ch
− 1
2
K ∥τ∥∥u−ΠWu∥∂K ,
∥ξk · n∥F ≤ ∥(q − PV q) · n∥F + ∥τ∥∥u−ΠWu∥F ,
∥PS1(ξk)∥K = 0
∥PS2(ξk)∥K = 0
Therefore, we have
|||ξk||| = ∥PS1(ξk)∥K + ∥PS2(ξk)∥K + hK∥∇ · ξk∥K +
∑
F ̸=E
h
1
2
K∥ξk · n∥F
≤ Ch
1
2
K∥(Id− PV )q · n∥∂K + Ch
1
2
K∥τ∥∥u−ΠWu∥∂K .
By (A.1.7), we have
∥u−ΠWu∥∂K ≤ ∥u− PWu∥∂K + ∥δk∥∂K ≤ C∥(Id− PW )u∥∂K + Ch
1
2
K∥(Id− PW˜ )∇ · q∥∂K ,
and we get
∥ξk∥K ≤ C(h
1
2
K∥(Id− PV )q · n∥∂K + hK∥(Id− PW˜ )∇ · q∥∂K) + h
1
2
K∥(Id− PW )u∥∂K)
This completes the proof.
132
A.1.4 Estimate of the jump εuh − εûh
Here, we give the following estimate of the jump εuh − εûh.
Lemma A.3. Let K be an element of the mesh Th. Assume that the local space V (K)×
W (K) admits anM(∂K)-decomposition. Then for the approximate solution of the HDG
method (2.2.2), we have the following local stability estimate
∥εuh − εûh∥∂K ≤ Ch
1
2
K∥q − qh∥K .
In the proof, we are going to use the adjoint HDG Projection which is defined as
follows.
Definition 5 ([19]). Let K be a regular polygonal or polyhedral element. Assume that
V (K) × W (K) admits an M(∂K)-decomposition. Let d := (dV , dw, dµ) ∈ V (K) ×
W (K)×M(∂K). Then the adjoint-HDG projection Π∗hd := (Π∗V d,Π∗Wd) ∈ V ×W is
defined by equations
(Π∗V d, v)K = (dV , v)K ∀ v ∈ V˜ (K),
(Π∗Wd,w)K = (dw, w)K ∀ w ∈ W˜ (K),
⟨Π∗V d · n− τΠ∗Wd , µ⟩F = ⟨dµ , µ⟩F ∀ µ ∈M(F ).
Definition 6. With the same assumption in Definition 5, let K be an element with one
curved face E. Then the adjoint-HDG projection Π∗hd := (Π
∗
V d,Π
∗
Wd) ∈ V (K)×W (K)
is defined by equations
(Π∗V d, v)K = (dV , v)K ∀ v ∈ S(K),
(Π∗Wd,w)K = (dw, w)K ∀ w ∈ W˜ (K),
⟨Π∗V d · n− τΠ∗Wd , µ⟩F = ⟨dµ , µ⟩F ∀ µ ∈M(F ), ∀F ̸= E,
⟨Π∗V d · n− τΠ∗Wd , w⟩∂K − (Π∗V d,∇w)K = ⟨dµ , w⟩∂K − (dV ,∇w)K ∀ w ∈W (K).
Proof. We rewrite the first two error equations (2.4.2a) and (2.4.2b) on the element K
133
as
(∇εuh,v)K − ⟨εuh − εûh , v · n⟩∂K = (qh − q,v),
(∇ · εqh, w)K + ⟨τ(εuh − εûh) , w⟩∂K = 0,
for all (v, w) ∈ V (K) ×W (K). Adding these equalities and using the fact that the
stabilization parameter τ is self-adjoint, we get
(∇ · εqh, w)K + (∇εuh,v)K − ⟨εuh − εûh , v · n− τw⟩∂K = (qh − q,v)K .
Now we take (v, w) := (Π∗V d,Π∗Wd), where d = −(0, 0, εuh − εûh) and Π∗ is the adjoint-
HDG projection. Note that for any interior face F , we have εuh − εûh ∈ M(F ) and
so
⟨εuh − εûh , Π∗V d·n− τΠ∗Wd⟩F = ⟨εuh − εûh , εuh − εûh⟩F ,
by definition. If K is a regular element, we have
−⟨εuh − εûh , εuh − εûh⟩∂K = (qh − q,Π∗V d)
On the other hand, notice that εuh − εûh = εuh on ∂Ω. Let K be a boundary element
with one curved face E and then we have
⟨εuh , Π∗V d·n− τΠ∗Wd⟩E = ⟨εuh , Π∗V d·n− τΠ∗Wd⟩∂K − ⟨εuh , Π∗V d·n− τΠ∗Wd⟩∂K\E
= ⟨dµ , εuh⟩∂K − (dV −Π∗V d,∇εuh)K − ⟨dµ , εuh⟩∂K\E
= ⟨dµ , εuh⟩E − (dV −Π∗V d,∇εuh)K
= ⟨εuh − εûh , εuh − εûh⟩E + (Π∗V d,∇εuh)K .
Therefore,
⟨εuh − εûh , Π∗V d·n− τΠ∗Wd⟩∂K = ⟨εuh − εûh , εuh − εûh⟩∂K + (Π∗V d,∇εuh)K .
So if K is an element with one curved face, we have
(∇εuh,Π∗V d)K − ⟨εuh − εûh , εuh − εûh⟩∂K − (Π∗V d,∇εuh)K = (qh − q,Π∗V d)K .
134
Hence, in both cases, we get
⟨εuh − εûh , εuh − εûh⟩∂K = (q − qh,Π∗V d)K
By the approximation properties of the adjoint HDG projection, we have that for d =
−(0, 0, εuh − εûh),
∥Π∗V d∥K ≤ Ch
1
2
K∥dµ∥∂K .
This estimate was shown in [19] for the case of regular elements K. For the case of
elements with a curved face, the proof is similar.
Therefore, after a simple application of the Cauchy-Schwarz inequality, we have
∥εuh − εûh∥∂K ≤ Ch
1
2
K∥q − qh∥K .
This completes the proof.
A.2 Auxiliary Lemmas
Lemma A.4. The HDG method defined in (2.2.2) can be rewritten as
A(uh, qh, ûh;w, r, µ) = F (w, r, µ) ∀ (w, r, µ) ∈Wh × V h ×Mh,
where the bilinear form A and linear form F are defined in (4.2.10) and (4.2.11), re-
spectively.
Proof. Let us rewrite the HDG equations (2.2.2) as follows:
(qh, r)Th − (uh,∇ · r)Th + ⟨ûh , r · n⟩∂Th\∂Ω = −⟨uD , r · n⟩∂Ω, (A.2.1a)
(∇ · qh, w)Th + ⟨τuh , w⟩∂Th − ⟨τ ûh , w⟩∂Th\∂Ω = (f, w)Th + ⟨τuD , w⟩∂Ω, (A.2.1b)
⟨ûh , µ⟩∂Ω =< uD, µ >∂Ω, (A.2.1c)
⟨qh · n+ τ(uh − ûh) , µ⟩∂Th\∂Ω = 0, (A.2.1d)
Now let (A.2.1b) - (A.2.1a) - (A.2.1d) - τ(A.2.1c) and use the fact that ⟨ûh , µ⟩∂Ω =
135
⟨ûh , µ⟩∂Th − ⟨ûh , µ⟩∂Th\∂Ω , we get
(∇ · qh, w)Th + (uh,∇ · r)Th − (qh, r)Th + ⟨τuh , w⟩∂Th − ⟨τ ûh , µ⟩∂Th
− ⟨ûh , r · n+ τ(w − µ)⟩∂Th\∂Ω − ⟨qh · n+ τ(uh − ûh) , µ⟩∂Th\∂Ω
=(f, w)Th + ⟨uD , r · n+ τ(w − µ)⟩∂Ω
Following the notation in Section 4.2.1, we have
A(uh, qh, ûh;w, r, µ)
=(∇ · qh, w)Th + (uh,∇ · r)Th − (qh, r)Th
+ ⟨τuh , w⟩∂Th − ⟨τ ûh , µ⟩∂Th
− ⟨ûh , r · n+ τ(w − µ)⟩∂Th\∂Ω − ⟨qh · n+ τ(uh − ûh) , µ⟩∂Th\∂Ω
and
F (w, r, µ) = (f, w)Th + ⟨uD , r · n+ τ(w − µ)⟩∂Ω
Lemma A.5. The error term (4.2.7) for the HDG method is the same as (4.2.12).
Proof. By the definition of the bilinear form A for the HDG method (4.2.10), we have
Eh =A(u− uh, q − qh, u− ûh; v − vh,p− ph, v − v̂h)
=(∇ · (q − qh), v − vh)Th + (u− uh,∇ · (p− ph))Th − (q − qh,p− ph)Th
+ ⟨τ(u− uh) , v − vh⟩∂Th − ⟨τ(u− ûh) , v − v̂h⟩∂Th
− ⟨u− ûh , (p− ph) · n+ τ(v̂h − vh)⟩∂Th\∂Ω
− ⟨(q − qh) · n+ τ(ûh − uh) , v − v̂h⟩∂Th\∂Ω
=(∇ · (q − qh), v − vh)Th + (u− uh,∇ · (p− ph))Th − (q − qh,p− ph)Th
+ ⟨τ(u− uh) , v − vh⟩∂Th − ⟨τ(u− ûh) , v − v̂h⟩∂Th
In the last step above , we used the definiton of q̂h · n and p̂h · n and the fact that
u, ûh, v, v̂h are single values on ∂Th\∂Ω. Note that ∇ · q = f . Compared with (4.2.12),
136
we only need to show
⟨τ(u− uh) , v − vh⟩∂Th − ⟨τ(u− ûh) , v − v̂h⟩∂Th
=⟨u− uh , τ(v̂h − vh)⟩∂Th + ⟨τ(ûh − uh) , v − vh⟩∂Th − ⟨τ(uh − ûh) , vh − v̂h⟩∂Th .
Indeed, with simple algebraic manipulation, we have
⟨τ(u− uh) , v − vh⟩∂Th − ⟨τ(u− ûh) , v − v̂h⟩∂Th
=⟨τ(u− ûh + ûh − uh) , v − vh⟩∂Th − ⟨τ(u− ûh) , v − v̂h⟩∂Th
=⟨τ(u− ûh) , v̂h − vh⟩∂Th + ⟨τ(ûh − uh) , v − vh⟩∂Th
=⟨τ(u− uh) , v̂h − vh⟩∂Th + ⟨τ(uh − ûh) , v̂h − vh⟩∂Th + ⟨τ(ûh − uh) , v − vh⟩∂Th
=⟨u− uh , τ(v̂h − vh)⟩∂Th + ⟨τ(ûh − uh) , v − vh⟩∂Th
−⟨τ(uh − ûh) , vh − v̂h⟩∂Th
This completes the proof.
Lemma A.6. The error term Eh in (4.2.12) is equivalent to A˜Ch + E˜h as defined in
[27, 28] for the HDG method.
Proof. Let us recall the definition of the adjoint-correction term A˜Ch and E˜h in our
previous work [27, 28]:
A˜Ch :=(f, vh)Th + (qh,∇vh)Th − ⟨q̂h · n , vh⟩∂Th
+ (qh +∇uh,ph)Th − ⟨uh − ûh , ph · n⟩∂Th
+ ⟨q̂h · n , v̂h⟩∂Th\∂Ω
+ ⟨uh − ûh , (ph − p̂h) · n⟩∂Th .
and
E˜h :=(q − qh,p− ph)Th
+ (q − qh,ph +∇vh)Th + (qh +∇uh,p− ph)Th
+ ⟨(q̂h − q) · n , vh − v̂h⟩∂Th + ⟨uh − ûh , (p̂h − p) · n⟩∂Th .
For the HDGmethod, we see that the first three lines of A˜Ch vanish by (2.2.2a), (2.2.2b),
137
(2.2.2d). Hence we have
A˜Ch = −⟨uh − ûh , τ(vh − v̂h) · n⟩∂Th
by the definition of p̂h · n.
Next, note that for the error term E˜h, we have
E˜h =− (q − qh,p− ph)Th
+ (q − qh,p− ph)Th + (q − qh,ph +∇vh)Th + ⟨(q̂h − q) · n , vh − v̂h⟩∂Th
+ (q − qh,p− ph)Th + (qh +∇uh,p− ph)Th + ⟨uh − ûh , (p̂h − p) · n⟩∂Th .
Notice that
(q − qh,p− ph)Th + (q − qh,ph +∇vh)Th
=(q − qh,−∇v +∇vh)Th since p = −∇v
=⟨(q − qh) · n , −v + vh⟩∂Th + (∇ · (q − qh), v − vh)Th
and
⟨(q̂h − q) · n , vh − v̂h⟩∂Th = ⟨(q̂h − q) · n , vh − v⟩∂Th
since q̂h · n and v̂h are single values on ∂Th and v̂h = v on ∂Ω. Therefore, we have
(q − qh,p− ph)Th + (q − qh,ph +∇vh)Th + ⟨(q̂h − q) · n , vh − v̂h⟩∂Th
=⟨(q̂h − qh) · n , vh − v⟩∂Th + (∇ · (q − qh), v − vh)Th
=⟨τ(ûh − uh) , v − vh⟩∂Th + (f −∇ · qh), v − vh)Th
Similarly, we have
(q − qh,p− ph)Th + (qh +∇uh,p− ph)Th + ⟨uh − ûh , (p̂h − p) · n⟩∂Th
=⟨u− uh , τ(v̂h − vh)⟩∂Th + (u− uh,∇ · p−∇ · ph)Th .
This completes the proof.