Computational lower bounds of the Maxwell eigenvalues

by   Dietmar Gallistl, et al.

A method to compute guaranteed lower bounds to the eigenvalues of the Maxwell system in two or three space dimensions is proposed as a generalization of the method of Liu and Oishi [SIAM J. Numer. Anal., 51, 2013] for the Laplace operator. The main tool is the computation of an explicit upper bound to the error of the Galerkin projection. The error is split in two parts: one part is controlled by a hypercircle principle and an auxiliary eigenvalue problem. The second part requires a perturbation argument for the right-hand side replaced by a suitable piecewise polynomial. The latter error is controlled through the use of the commuting quasi-interpolation by Falk–Winther and computational bounds on its stability constant. This situation is different from the Laplace operator where such a perturbation is easily controlled through local Poincaré inequalities. The practical viability of the approach is demonstrated in two-dimensional test cases.



There are no comments yet.


page 1

page 2

page 3

page 4


Information inequalities for the estimation of principal components

We provide lower bounds for the estimation of the eigenspaces of a covar...

Guaranteed Lower Eigenvalue Bound of Steklov Operator with Conforming Finite Element Methods

For the eigenvalue problem of the Steklov differential operator, by foll...

Guaranteed a posteriori bounds for eigenvalues and eigenvectors: multiplicities and clusters

This paper presents a posteriori error estimates for conforming numerica...

Asymptotic analysis and numerical computation of the Laplacian eigenvalues using the conformal mapping

We consider the Dirichlet and Neumann eigenvalues of the Laplacian for a...

Direct guaranteed lower eigenvalue bounds with optimal a priori convergence rates for the bi-Laplacian

An extra-stabilised Morley finite element method (FEM) directly computes...

Electromagnetic Stekloff eigenvalues: approximation analysis

We continue the work of [Camano, Lackner, Monk, SIAM J. Math. Anal., Vol...

Uniform Error Estimates for the Lanczos Method

The Lanczos method is one of the most powerful and fundamental technique...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

This paper is devoted to the computation of guaranteed lower bounds of the Maxwell eigenvalues. The Maxwell eigenvalue problem over a suitable bounded domain in dimension or seeks eigenpairs with nontrivial such that


Here, is the outer unit normal to and is the tangential trace of . The usual rotation (or curl) operator is denoted by , while its formal adjoint is denoted by ; precise definitions are given below. The operator has an infinite-dimensional kernel containing all admissible gradient fields, leading to an eigenvalue of infinite multiplicity. Sorting out this eigenvalue in numerical computations requires the incorporation of a divergence-free constraint. In the setting of -conforming finite elements, the divergence constraint is necessarily imposed in a discrete weak form because simultaneous and -conformity may lead to non-dense and thus wrong approximations [Cos91, Mon03]. In conclusion, a variational form of (1.1) in an energy space is in general approximated with a nonconforming discrete space and no monotonicity principles are applicable for a comparison of discrete and true eigenvalues. Upper eigenvalue bounds can be expected from discontinuous Galerkin (dG) schemes [BP06], but the practically more interesting question of guaranteed lower bounds has remained open until the contributions [BBB14, BBB17]. For a detailed exposition of the eigenvalue problem and its numerical approximation, the reader is referred to [Hip02, Mon03, Bof10] and the references therein.

In the finite element framework for coercive operators (e.g., the Laplacian), guaranteed lower eigenvalue bounds were successfully derived by the independent contributions [LO13] and [CG14b], which basically follow the same reasoning, illustrated here for the first eigenpair of a variational eigenvalue problem

with inner products and and corresponding norms and . For a (possibly nonconforming) discretization with the first discrete eigenpair , the discrete Rayleigh–Ritz principle [WS72] implies

for any . Given an -orthogonal projection operator (assuming is defined on the sum ), this and some elementary algebraic manipulations show

Assuming the normalization so that , it turns out that explicit control of by yields a computational lower bound. In [LO13] is the standard Galerkin projection while in [CG14b] is the interpolation operator in a Crouzeix–Raviart method.

In this paper we aim at extending the idea of [LO13] to the Maxwell eigenvalue problem (1.1) discretized with lowest-order Nédélec (edge) elements [Mon03]. The main novelty in contrast to [LO13] is the guaranteed computational control of the Galerkin error in a linear Maxwell system with right-hand side . In [LO13] such bound is achieved for the Laplacian by splitting in a piecewise polynomial part and some remainder. The first part of the error is quantified through a hypercircle principle [Bra07] and an auxiliary global eigenvalue problem. This idea goes back to [LO13]. The second part of the error is —in the case of the Laplacian, where is simply the piecewise mean of — easily controlled because it reduces to element-wise Poincaré inequalities whose constants can be explicitly bounded [PW60]. In the present case of the Maxwell system, the situation is more involved. The operator maps the Nédélec space to the divergence-free Raviart–Thomas elements [BBF13], and the -orthogonal projection to the latter is nonlocal and explicit bounds on that projection are unknown. In order to obtain a computable bound, we make use of recent developments of Finite Element Exterior Calculus [AFW06, Arn18], namely the Falk–Winther projection [FW14]. This family of operators commutes with the exterior derivative and is locally defined, so that it is actually computable. Practical implementations of the operator have been used in the context of numerical homogenization [GHV18, HP19, HP20]

. In this work, the advantage of the local construction is that the involved stability constant can be computationally bounded from above. In a perturbation argument for the Maxwell system, this tool replaces the Poincaré inequality from the Laplacian case. The bounds are achieved by solving local discrete eigenvalue problems combined with standard estimates.

The main result is the guaranteed lower bound

from Theorem 4.1 for the th Maxwell eigenvalue . The quantities on the left-hand side are the th discrete eigenvalue and the computable mesh-dependent number . In particular, the computation of involves guaranteed control over the bound for the Galerkin error in a linear Maxwell problem and the local stability constants of the Falk–Winther interpolation. Guaranteed computational bounds for the quantities entering are carefully described in this paper.

The remaining parts of this article are organized as follows. §2 lists preliminaries on the Maxwell problem, discrete spaces, and the Falk–Winther interpolation. Guaranteed computational bounds on the Galerkin error are presented in §3. The lower eigenvalue bounds are shown in §4 and practically computed in the numerical experiments of §5. The remarks of §6 conclude this paper. An appendix describes the concrete computation of the involved constants.

2 Preliminaries

2.1 Notation

Let for be a bounded and open polytopal Lipschitz domain, which we assume to be simply connected. The involved differential operators read

For the formal adjoint operators we write

so that the integration-by-parts formula

holds for sufficiently regular vector fields

, with vanishing tangential trace over .

Standard notation on Lebesgue and Sobolev spaces is employed throughout this paper. Given any open set , the inner product is denoted by with the norm . The usual -based first-order Sobolev space is denoted by and is the subspace with vanishing trace over . The space of vector fields over with weak divergence in is denoted by ; and the subspace of divergence-free vector fields reads . The space of vector fields with weak rotation in is denoted by while its subspace with vanishing tangential trace is denoted by .

In the context of eigenvalue problems, the inner product is also denoted by and the norm is denoted by .

On , we define the bilinear form

Let . On , the form is an inner product [Mon03, Corollary 4.8] and the seminorm is a norm on . Given a divergence-free right-hand side , the linear Maxwell problem seeks such that


It is well known [Mon03] and needed in some arguments of this article that (2.1) is even satisfied for all test functions from the larger space .

Let be a regular simplicial triangulation of . The diameter of any is denoted by and is the maximum mesh size. Given any , the space of first-order polynomial functions over is denoted by . The lowest-order standard finite element space (with or without homogeneous Dirichlet boundary conditions) is denoted by

The space of lowest-order edge elements [Mon03, BBF13] reads

and we denote

The approximation of (2.1) with edge elements uses the space

The elements of are weakly divergence-free and need in general not be elements of , i.e., . It is known [Mon03] that is an inner product on . The finite element system seeks such that


Given , its approximation is called the Galerkin projection. This terminology is justified by the fact that the approach is conforming when viewed in a saddle-point setting. In particluar, since (2.1) is satisfied for all , the following “Galerkin orthogonality” is valid


The Raviart–Thomas finite element space [BBF13] is defined as

2.2 Falk–Winther interpolation

Given a regular triangulation and any element , the element patch built by the simplices having nontrivial intersection with is defined as

There is a projection with local stability in the sense that there exist constants , such that for any and any we have


Furthermore, there is a projection where

with constants , such that for any and any we have


The crucial property is that these two operators commute with the exterior derivative in the sense that . The corresponding commuting diagram is displayed in Figure 1.

Figure 1: Commuting diagram of the Falk–Winther operator.

3 Bounds on the Galerkin projection

The goal of this section is a fully computable bound on the Galerkin error.

3.1 error control

From elliptic regularity theory, it is known that there exists a mesh-dependent (but -independent) number such that


On convex domains, is proportional to the mesh size while in general, reduced regularity [CD00] implies that is proportional to with some .

Given some , the Galerkin approximation is usually not divergence-free and therefore possesses a nontrivial orthogonal decomposition


with and . The inclusion furthermore shows that . The following lemma states an error estimate. The proof uses the classical Aubin–Nitsche duality technique.

Lemma 3.1.

The divergence-free part of the error satisfies the following error estimate


The error is divergence-free and thus . There exists a unique solution satisfying

Since , we infer with the symmetry of and that

The Galerkin property (2.3) shows that is -orthogonal to any element of . Thus

The application of (3.1) to with right-hand side reveals that the Galerkin error is bounded by , which implies the asserted bound. ∎

3.2 Perturbation of the right-hand side

In this section, we quantify the error that arises in the solution of the linear Maxwell system when the right-hand side is replaced by a piecewise polynomial approximation .

Let . The regular decomposition [CM10] implies that there exists such that


We denote the overlap constant of element patches by

is an overlap constant.

Given and its element patch , the Poincaré inequality states for any function with vanishing average over that with a constant proportional to . By we denote the smallest constant such that

holds for all such functions uniformly in .

Recall that, due to its commutation property, the Falk–Winther interpolation maps to .

Lemma 3.2.

Let and let be its Falk–Winther interpolation. Let and denote the solution to (2.1) with right-hand side and , respectively. Let furthermore and denote the solution to (2.2) with right-hand side and , respectively. Then

for the constant


The proof is only shown for and ; the case of and is completely analogous. Abbreviate . The solution properties imply

From the regular decomposition (3.3) and the commuting property of the operators , we obtain

Thus, integration by parts and the homogeneous boundary conditions of imply

The combination with the above chain of identities implies


Let be arbitrary. Since locally preserves constants, we obtain for the patch average

that the difference can be split with the triangle inequality and the inclusion as follows

Estimate (2.5) with followed by the Poincaré inequality on with constant thus reveal

This local result generalizes to the whole domain as follows

We estimate with and obtain

The combination with (3.4) concludes the proof. ∎

3.3 Error bound on the Galerkin projection

Let . Given , let



The next lemma states that, for discrete data, the Galerkin error can be quantified through .

Lemma 3.3.

Let , and let and be solutions to (2.1) and (2.2), respectively, with right-hand side . Then the following computable error estimate holds


Let be arbitrary and let . Integration by parts implies the orthogonality

Thus the following hypercircle identity holds



The left-hand side is minimal for amongst all elements , whence

where the last estimate follows from the definition (3.5). ∎

For possibly non-discrete , the Galerkin error is bounded by the following perturbation argument.

Lemma 3.4.

Given , let solve (2.1) and let solve (2.2). Then, the following error bound is satisfied


Let denote the Falk–Winther interpolation of and denote as in Lemma 3.3 by and the solutions with respect to . The triangle inequality leads to


The second term on the right-hand side is bounded through Lemma 3.3 as follows

The local stability (2.4) of (note that ) and the overlap of element patches imply



The sum of the first and the third term on the right-hand side of (3.7) is bounded through Lemma 3.2 by . The combination of this bound with (3.7) and (3.8) concludes the proof. ∎

4 Eigenvalue problem and lower bound

The Maxwell eigenvalue problem seeks pairs with such that


The condition implies that and, thus, the non-compact part of the spectrum of the

operator (corresponding to gradient fields as eigenfunctions) is sorted out in this formulation. It is well known

[Mon03] that the eigenvalues to (4.1) form an infinite discrete set

The discrete counterpart seeks with such that


We now focus on the th eigenpair and its approximation . Recall from (3.1), which is practically bounded by Lemma 3.4.

Theorem 4.1.

Let with be the th eigenpair to (4.1) and with be the th discrete eigenpair to (4.2). The following lower bound

is valid for . For it is valid provided the separation condition holds.


We focus on the case . Recall the Helmholtz decomposition (3.2) with the divergence-free part of . The discrete Rayleigh–Ritz principle [WS72] states

and therefore


We expand the square on the left-hand side, use that , and use the Young inequality with an arbitrary to infer

(note that ). The combination with Lemma 3.1 results in


The Galerkin orthogonality (2.3) in the -product and the identity for the right-hand side of (4.3) result in


The choice and the combination of (4.3)–(4.5) results in

The proof for proceeds analogously to [CG14b] where the separation condition guarantees that the space spanned by the first eigenfunctinos is mapped to a -dimensional space by the Galerkin projection. The details are omitted. The reader is instead referred to [CG14b, Lemma 5.2] and [Gal14, Lemma 3.13]. ∎

5 Numerical results

In this section, we present numerical results for the two-dimensional case. The Python implementation uses the FEnics library [ABH15] and the realization of the Falk–Winther operator from [HP19].

The two domains we consider are the unit square and the L-shaped domain . The initial triangulations are displayed in Figure 2. We consider uniform mesh refinement (red-refinement).

The computational bound on (see §5.3 for actual values) is achieved with the techniques described in Appendix A.

Figure 2: Initial triangulations.

5.1 Results on the square domain

For the square domain, it is known that the first two eigenvalues are given by and . The numerical results are displayed in Table 1. The lower bounds are guaranteed, but their values become practically relevant after a moderate number of refinements only when is sufficiently small.

lower bd. lower bd.
0.1443 12.4297 9.6000 0.0065 20.2871 0.0065
0.0721 6.2134 9.8305 0.0258 20.0235 0.0259
0.0356 3.0934 9.8612 0.1034 19.8205 0.1040
0.0180 1.5520 9.8676 0.3984 19.7601 0.4066
0.0089 0.7734 9.8691 1.4298 19.7445 1.5415
0.0044 0.3852 9.8695 4.0048 19.7405 5.0242
0.0023 0.1941 9.8696 7.1949 19.7395 11.3218
0.0011 0.0972 9.8696 9.0280 19.7393 16.6375
0.0006 0.0484 9.8696 9.6462 19.7392 18.8653
Table 1: Numerical results on the square domain.

5.2 Results on the L-shaped domain

On the L-shaped domain, we use the reference value from [CDMV03] for the first eigenvalue for comparison. The numerical results are displayed in Table 2. As in the previous example, the lower bounds take values of practical significance after a couple of refinement steps.

lower bd.
0.1355 12.1702 1.3180 0.0067
0.0709 6.1780 1.4157 0.0257
0.0356 3.0934 1.4526 0.0975
0.0180 1.5526 1.4667 0.3234
0.0090 0.7762 1.4721 0.7801
0.0045 0.3876 1.4743 1.2070
0.0022 0.1941 1.4751 1.3975
0.0011 0.0972 1.4754 1.4550
Table 2: Numerical results on the L-shaped domain.

5.3 Values of the individual constants

Table 3 displays the values of the constants for the unit square and the L-shaped domain entering the computation of the upper bound

to from Lemma 3.4. All these constants coincide for the two domains, which have a similar local mesh geometry.

constant value
Table 3: Values of some of the relevant constants for the square and the L-shaped domain.

6 Conclusive remarks

The theory on computational lower bounds to the Maxwell eigenvalues applies to the case of two or three space dimensions. The sharpness of the bounds critically depends on the actual value of on coarse meshes. We remark that in two dimensions the eigenvalues coincide with the Laplace–Neumann eigenvalues, so our computations should be rather seen as a proof of concept. We succeeded in bounding such that meaningful lower bounds could be achieved on moderately fine meshes. The novel ingredient is the explicit computational stability bound on the Falk–Winther operator that makes a full quantification of the Galerkin error possible. Computations in the three-dimensional case are beyond the scope of this work. Furthermore, the methodology does not immediately generalize to adaptive meshes in the sense that is expected to scale like (some power of) the maximum mesh size.

Appendix A Appendix: Bounds on the involved constants

The appendix describes how the critical constants are computed or computationally bounded.

a.1 Computation of

In this paragraph we briefly sketch the numerical computation of . The reasoning is similar to [LO13], and we illustrate the extension of their approach to the Maxwell operator. Recall from §3.3 the spaces and for given . Let furthermore denote the solution to the linear problem (2.1) with right-hand side . Expanding squares and integration by parts then shows for any and any the hypercircle identity (3.6) because . Thus, the optimization problem

is equivalent to

The minimizer to the first term is given by the solution to the primal problem (2.2) while the minimizer to the second term is given by as part of the solution pair to the following (dual) saddle-point problem

This can be shown with arguments analogous to [Bra07, III§9, Lemma 9.1].

Let denote the solution operator to the primal problem and let , denote the components of the solution operator to the dual problem. The operators and are symmetric, and straightforward calculations involving the definitions of , , prove that the error can be represented as

The divergence-free constraint in the primal problem is practically implemented in a saddle-point fashion. In what follows, we identity the piecewise constant function with its vector representation. The matrix structure of the discrete problem is as follows

The solution can be expressed by

where denotes the relevant rows and columns of for the computation of the first component of the solution.

The dual system can be implemented by introducing Lagrange multipliers related to the interior hyper-faces to enforce normal-continuity. In terms of matrices, the dual system reads

The solution can be expressed by