# Computational lower bounds of the Maxwell eigenvalues

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.

## Authors

• 4 publications
• 1 publication
05/14/2020

### Information inequalities for the estimation of principal components

We provide lower bounds for the estimation of the eigenspaces of a covar...
01/24/2020

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

For the eigenvalue problem of the Steklov differential operator, by foll...
08/10/2020

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

This paper presents a posteriori error estimates for conforming numerica...
12/21/2021

### 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...
05/04/2021

### 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...
09/02/2019

### Electromagnetic Stekloff eigenvalues: approximation analysis

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

### 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

 (−1)d−1Curlrotu=λu in Ωandu∧ν=0 on ∂Ω. (1.1)

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

 a(u,v)=λb(u,v)for all v∈V

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

 λh∥vh∥2b≤∥vh∥2a

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

 λh∥Ghu∥2b≤∥Ghu∥2a≤∥u∥2a−∥u−Ghu∥2a.

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

 λh1+M2hλh≤λ

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

 rotv=∂1v2−∂2v1 for d=2androtv=⎛⎜⎝∂2v3−∂3v2∂3v1−∂1v3∂1v2−∂2v1⎞⎟⎠ for d=3.

For the formal adjoint operators we write

 Curlϕ=(−∂2ϕ∂1ϕ) for d=2andCurl=rot for d=3

so that the integration-by-parts formula

 ∫Ωϕrotvdx=(−1)d−1∫ΩCurlϕ⋅vdx

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

 a(v,w):=(rotv,rotw)L2(Ω)% for any v,w∈H0(rot,Ω).

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

 a(u,v)=b(f,v)for all v∈V. (2.1)

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

 S1(T):={v∈H1(Ω):∀T∈T,v|T∈P1(T)} and S10(T):=H10(Ω)∩S1(T).

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

 N0(T)={v∈H(rot,Ω):∀T∈T∃αT∈Rd∃βT∈P1(T)d∀x∈T: v(x)=αT+βT(x) and βT(x)⋅x=0}

and we denote

 N0,D(T):=N0(T)∩H0(rot,Ω).

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

 Vh={ψh∈N0,D(T):∀vh∈S10(T)(∇vh,ψh)L2(Ω)=0}.

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

 a(uh,vh)=b(f,vh)for all vh∈Vh. (2.2)

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

 a(u−uh,vh)=0for all vh∈Vh. (2.3)

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

 RT0(T):={v∈H(div,Ω): ∀T∈T∃(αT,βT)∈Rd×R ∀x∈T,v|T(x)=αT+βTx}.

### 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

 ωT:=int(∪{K∈T:K∩T≠∅}).

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

 ∥πdivv∥L2(T)≤C1,div∥v∥L2(ωT)+hTC2,div∥divv∥L2(ωT). (2.4)

Furthermore, there is a projection where

 H(Curl,Ω)={H1(Ω;R2)% if d=2H(rot,Ω)if d=3andSh={S1(T)if d=2N0(T)if d=3

with constants , such that for any and any we have

 ∥πCurlv∥L2(T)≤C1,Curl∥v∥L2(ωT)+hTC2,Curl∥Curlv∥L2(ωT). (2.5)

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.

## 3 Bounds on the Galerkin projection

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

### 3.1 L2 error control

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

 ∥u−Ghu∥a≤Mh∥f∥b. (3.1)

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

 Ghu=∇ϕ+R (3.2)

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

 ∥u−R∥b≤Mh∥u−Ghu∥a.
###### Proof.

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

 a(z,v)=b(e,v)for all v∈V.

Since , we infer with the symmetry of and that

 ∥e∥2b=b(e,e)=a(z,e)=a(e,z)=a(u−Ghu,z).

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

 ∥e∥2b=a(u−Ghu,z−Ghz)≤∥u−Ghu∥a∥z−Ghz∥a.

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

 Curlβ=fand∥Dβ∥L2(Ω)≤CRD∥f∥L2(Ω). (3.3)

We denote the overlap constant of element patches by

 COL=maxT∈Tcard{K∈T:T⊆¯¯¯¯¯¯¯ωK}.

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

 ∥v∥L2(ωT)≤hT~c∥Dv∥L2(ωT)

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

 max{∥u−~u∥a,∥uh−~uh∥a}≤√2COLhmax^C∥f∥L2(Ω)

for the constant

 ^C:=√((1+C1,Curl)~cCRD)2+C22,Curl.
###### Proof.

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

 ∥u−~u∥2a=a(u−~u,e)=b(f−~f,e)=b(f−πdivf,e).

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

 f−πdivf=Curlβ−πdivCurlβ=Curl(β−πCurlβ).

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

 b(f−πdivf,e) =b(Curl(β−πCurlβ),e) =(−1)d−1b(β−πCurlβ,rote)≤∥β−πCurlβ∥b∥e∥a.

The combination with the above chain of identities implies

 ∥u−~u∥a≤∥β−πCurlβ∥b. (3.4)

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

 ¯β:=\fintωTβdx

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

 ∥β−πCurlβ∥L2(T)≤∥β−¯β∥L2(ωT)+∥πCurl(β−¯β)∥L2(T).

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

 ∥β−πCurlβ∥L2(T)≤hT~c(1+C1,Curl)∥Dβ∥L2(ωT)+C2,CurlhT∥f∥L2(ωT).

This local result generalizes to the whole domain as follows

 ∥β−πCurlβ∥2b =∑T∈T∥β−πCurlβ∥2L2(T) ≤∑T∈T(hT~c(1+C1,Curl)∥Dβ∥L2(ωT)+C2,CurlhT∥f∥L2(ωT))2 ≤2∑T∈Th2T((~c(1+C1,Curl)∥Dβ∥L2(ωT))2+(C2,Curl∥f∥L2(ωT))2) ≤2COLh2max((~c(1+C1,Curl)∥Dβ∥L2(Ω))2+(C2,Curl∥f∥L2(Ω))2).

We estimate with and obtain

The combination with (3.4) concludes the proof. ∎

### 3.3 Error bound on the Galerkin projection

Let . Given , let

 Sfh:={vh∈Sh:Curlvh=(−1)d−1fh}.

Define

 κh:=maxfh∈Xh∖{0}minvh∈Vhminτh∈Sfh∥τh−rotvh∥L2(Ω)∥fh∥L2(Ω). (3.5)

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

 ∥~u−~uh∥a≤κh∥fh∥b.
###### Proof.

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

 (τh−rot~u,rot(~u−vh))L2(Ω)=0.

Thus the following hypercircle identity holds

 ∥τh−rot~u∥2L2(Ω)+∥rot(~u−vh)∥2L2(Ω)=∥τh−rotvh∥2L2(Ω). (3.6)

Thus,

 ∥~u−vh∥a=∥rot(~u−vh)∥L2(Ω)≤minτ∈Sfh∥τh−rotvh∥L2(Ω).

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

 ∥~u−~uh∥a≤minvh∈Vhminτ∈Sfh∥τh−rotvh∥L2(Ω)≤κh∥fh∥L2(Ω)

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

 ∥u−uh∥a≤(2√2hmax^C+κhC1,div)√COL∥f∥L2(Ω).
###### Proof.

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

 ∥u−uh∥a≤∥u−~u∥a+∥~u−~uh∥a+∥~uh−uh∥a. (3.7)

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

 ∥~u−~uh∥a≤κh∥fh∥L2(Ω)=κh∥πdivf∥L2(Ω).

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

 ∥πdivf∥2L2(Ω)=∑T∈T∥πdivf∥2L2(T)≤∑T∈TC21,div∥f∥2L2(ωT)≤COLC21,div∥f∥2L2(Ω).

Thus,

 ∥~u−~uh∥a≤κh√COLC1,div∥f∥L2(Ω). (3.8)

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

 a(u,v)=λb(u,v)for all v∈V. (4.1)

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

 0<λ1≤λ2≤…with limj→∞λj=∞.

The discrete counterpart seeks with such that

 a(uh,vh)=λhb(uh,vh)for all vh∈Vh. (4.2)

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

 λh1+M2hλh≤λ

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

###### Proof.

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

 λh=minvh∈Vh∖{0}∥vh∥2a∥vh∥2b

and therefore

 λh∥R∥2b≤λh∥Ghu∥2b≤∥Ghu∥2a. (4.3)

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

 ∥R∥2b=∥R−u∥2b+1+2b(R−u,u)≥(1−δ−1)∥R−u∥2b+(1−δ)

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

 λh∥R∥2b≥λh((1−δ−1)M2h∥Ghu−u∥2a+(1−δ)). (4.4)

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

 ∥Ghu∥2a≤λ−∥u−Ghu∥2a. (4.5)

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

 λh1+M2hλh≤λ.

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.

### 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.

### 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.

### 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

 Mh≤hmax2√(1+C1,Curl)2~c2C2RD+C22,Curl√2COL+κh√COLC1,div

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

## 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 κh

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

 minvh∈N0,D(T)minτh∈Sfh∥τh−rotvh∥2L2(Ω)

is equivalent to

 minvh∈N0,D(T)∥rot(~u−vh)∥2L2(Ω)+minτh∈Sfh∥τh−rot~u∥2L2(Ω).

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

 (σh,τh)L2(Ω)+(−1)d(Curlτh,ρh)L2(Ω) =0 for all τh∈Sh (−1)d(Curlσh,ϕh)L2(Ω) =−(fh,ϕh)L2(Ω) for all ϕh∈Xh.

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

 ∥rotuh−σh∥2L2(Ω)=−(fh,Thfh)L2(Ω)+(fh,Rhfh)L2(Ω).

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

 L(w1w2)=(Bfh0)where L:=(DFTF−0).

The solution can be expressed by

 w1=L−1(∗,∗)Bfh=HBfh,

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

 K⎛⎜⎝z1z2z3⎞⎟⎠=⎛⎜⎝0Mfh0⎞⎟⎠where K:=⎛⎜⎝A−GT0G−0C0−CT0⎞⎟⎠.

The solution can be expressed by

 z2