# A discontinuous least squares finite element method for time-harmonic Maxwell equations

We propose and analyze a discontinuous least squares finite element method for solving the indefinite time-harmonic Maxwell equations. The scheme is based on the L^2 norm least squares functional with the weak imposition of the continuity across the interior faces. We minimize the functional over the piecewise polynomial spaces to seek numerical solutions. The method is shown to be stable without any constraint on the mesh size. We prove the convergence orders under both the energy norm and the L^2 norm. Numerical results in two and three dimensions are presented to verify the error estimates.

## Authors

• 15 publications
• 2 publications
• 5 publications
05/05/2021

### A discontinuous least squares finite element method for Helmholtz equations

We propose a discontinuous least squares finite element method for solvi...
11/13/2019

### A Least-Squares Finite Element Method Based on the Helmholtz Decomposition for Hyperbolic Balance Laws

In this paper, a least-squares finite element method for scalar nonlinea...
03/04/2020

### Reconstructed Discontinuous Approximation to Stokes Equation in A Sequential Least Squares Formulation

We propose a new least squares finite element method to solve the Stokes...
01/14/2021

### Lebesgue integration. Detailed proofs to be formalized in Coq

To obtain the highest confidence on the correction of numerical simulati...
09/18/2019

### First-order system least squares finite-elements for singularly perturbed reaction-diffusion equations

We propose a new first-order-system least squares (FOSLS) finite-element...
09/30/2021

### Least-Squares Finite Element Method for Ordinary Differential Equations

We consider the least-squares finite element method (lsfem) for systems ...
05/22/2020

### A short note on plain convergence of adaptive least-squares finite element methods

We show that adaptive least-squares finite element methods driven by the...
##### 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

The time-harmonic Maxwell equations are often encountered in many engineering applications such as antenna design, microwaves, and satellites [13, 29]. These applications require us to carry out numerical studies of the time-harmonic Maxwell equations. The Maxwell’s operator is strongly indefinite especially for the case of the high wave number, which brings many difficulties in the numerical simulation and the error analysis [25]. Despite these difficulties, there are plenty of studies on numerical methods for this problem, such as finite difference methods, spectral methods, and finite element methods.

There are a variety of finite element methods for solving the time-harmonic Maxwell equations and related problems. A common choice is to employ the -conforming elements, known as edge elements. We refer to [28, 27, 26, 13, 11, 19] for more details of these conforming methods. We note that the implementation of edge elements is still challenging especially in the higher-order case [7].

Using discontinuous functions to approximate the solution, the discontinuous Galerkin (DG) methods have been applied in the numerical simulation of the Maxwell equations. The DG methods can offer great flexibility in the mesh structure which allows the elements of different shapes and can easily handle the irregular non-conforming meshes. Additionally, the implementation of discontinuous elements is very straightforward and there is no need to use the curl-conforming elemental mappings [25]. We refer to [25, 12, 15, 32, 29, 9, 30, 31, 10] and the references therein for those DG methods.

The least squares finite element method (LSFEM) is a general technique in numerical PDEs which is based on the minimization of the norm of the residual over an approximation space, and we refer to the paper [6] for an overview. The LSFEM has also been applied to solve the Maxwell system and we refer to [7, 8, 18, 16, 17] for such least squares methods. One of the advantages of LSFEM is that the resulting linear system is always symmetric and positive definite. Recently, the LSFEM has been extended for the numerical approximation to some classical PDEs using discontinuous elements, and this method is referred to as the discontinuous LSFEM. We refer the readers to [2, 3, 5, 4, 23, 24] for more details.

Noticing that the Maxwell’s operator is indefinite while the LSFEM always provides a positive definite linear system, we are motivated to develop a stable numerical method, based on the least squares functional and discontinuous elements, to solve the time-harmonic Maxwell equations. For this purpose, we define an norm least squares functional involving the proper penalty terms which weakly enforce the tangential continuity across the interior faces as well as the boundary conditions. The functional is then minimized over the discontinuous piecewise polynomial spaces to seek the numerical solutions. The method is easy to be implemented and the resulting linear system is shown to be symmetric and positive definite. The method combines the attractive features of DG methods and LSFEMs, giving us a new discontinuous least squares finite element method for the time-harmonic Maxwell equations.

We estimate the error of the new method to derive the convergence rates for all solution variables under both the norms and the energy norms. It is interesting that the proposed method is shown to be unconditionally stable without making any assumptions about the mesh size. We carry out a series of numerical tests in two and three dimensions to verify the theoretical predictions. Noticing that the least squares functional can serve as an a posteriori estimator, we particularly present a low-regularity example, which is solved by the -adaptive refinement strategy using the least squares functional as an adaptive indicator.

The rest of this paper is organized as follows. In Section 2, we introduce the notations and define the first-order system to the time-harmonic Maxwell equations. In Section 3, we present our least squares method and define the least squares functional. The error estimates are also proven in this section. In Section 4, we present a series of numerical examples to illustrate the accuracy of the proposed method.

## 2. Preliminaries

Let be an open, bounded polygonal (polyhedral) domain with the boundary in , . We denote by a regular and shape-regular partition over the domain into triangles (tetrahedrons). We let be the collection of all dimensional interior faces with respect to the partition , and we let be the collection of all dimensional faces that are on the boundary , and then we set . In particular, for the three-dimensional case we denote by the set of all dimensional interior edges of all elements in , and by the set of all dimensional boundary edges, and we let . For the element and the face , we let and be their diameters, respectively, and we denote as the mesh size of . Then, the shape-regularity of is in the sense of that there exists a constant such that for any element ,

 hKρK≤C,

where denotes the diameter of the largest disk (ball) inscribed in .

Next, we introduce the following trace operators that are commonly used in the DG framework. Let be an interior face shared by two adjacent elements and

with the unit outward normal vectors

and on , respectively. For the piecewise smooth scalar-valued function and the piecewise smooth vector-valued function , we define the jump operator on as

 [[n×v]]:=n+×v|K++n−×v|K−,[[n×q]]:=n+×q|K++n−×q|K−.

For a boundary face , we let be the element that contains and the jump operator on is defined as

 [[v]]:=n×v|K,[[n×q]]:=n×q|K,

where is the unit outward normal on .

We note that the capital with or without subscripts are generic positive constants, which are possibly different from line to line, are independent of the mesh size but may depend on the wave number and the domain . Given a bounded domain , we will follow the standard notations for the Sobolev spaces , , and with the regular exponent . We also use the standard definitions of their corresponding inner products, semi-norms and norms. Throughout the paper, we mainly use the notation for the three-dimensional case. For the case , it is natural to identity the space with the plane in . Specifically, in two dimensions for the vector-valued function , the curl of reads

 ∇×u=∂u2∂x−∂u1∂y,

and for the scalar-valued function , we let be the formal adjoint, which reads

 ∇×q=(∂q∂y,−∂q∂x)T.

For the problem domain , we define the space

 H(curl) :={v∈L2(Ω) | ∇×v∈L2(Ω)2}, for scalar-valued functions% ,{v∈L2(Ω)2 | ∇×v∈L2(Ω)}, for vector-valued functions,d=2, H(curl) :={v∈L2(Ω)3 | ∇×v∈L2(Ω)3},d=3.

Further, we denote by the space of functions in with vanishing tangential trace,

 H0(curl):={v∈H(curl) | n×v=0, on ∂Ω}.

For the partition , we will use the notations and the definitions for the broken Sobolev space , , and with the exponent and their associated inner products and norms [1].

The problem under our consideration is to find a numerical approximation to the time-harmonic Maxwell equations in a lossless medium with an inhomogeneous boundary condition, which seeks the electric field such that

 (1) ∇×μ−1r∇×u−k2εru =f,in Ω, n×u =g,on ∂Ω.

Here is an external source file and is the wave number with the assumption

is not an eigenvalue of the Maxwell system

[21, 26, 14]. and are the relative magnetic permeability and the relative electric permittivity of the medium. For simplicity, we set and .

Below let us propose a least squares method for the equations (LABEL:eq_Maxwell) based on the discontinuous approximation. We first introduce an auxiliary variable to rewrite the Maxwell equations LABEL:eq_Maxwell into an equivalent first-order system:

 (2) ∇×p−ku=˜f, in Ω, ∇×u−kp=0, in Ω, n×u=g on ∂Ω,

where . We note that to rewrite the problem into a first-order system is the fundamental idea in the modern least squares finite element method [6], and our discontinuous least squares method is then based on the system (2).

## 3. Discontinuous Least Squares Method for Time-Harmonic Equations

Aiming to construct a discontinuous least squares finite element method for (2), we first define a least squares functional based on (2), which reads:

 (3) Jh(u, p):=∑K∈Th(∥∇×p−ku−˜f∥2L2(K)+∥∇×u−kp∥2L2(K)) +

where is a positive parameter and will be specified later on. We note that in two dimensions the variable in (3) is scalar-valued and in three dimensions the variable is a vector in . Then, we define two approximation spaces and for the variables and , respectively,

 Vmh :=(Vmh)d,Σmh:=(Vmh)2d−3,

where is the standard piecewise polynomial space,

 Vmh :={vh∈L2(Ω) | vh|K∈Pm(K), ∀K∈Th}.

Clearly, the functions in and can be discontinuous across inter-element faces. We seek the numerical solution by minimizing the functional (3) over the space , which reads:

 (4) (uh,ph)=argmin(vh,qh)∈Vmh×ΣmhJh(vh,qh).

To solve the minimization problem (4), we can write the corresponding Euler-Lagrange equation, which takes the form: find such that

 (5) ah(uh,ph;vh,qh)=lh(vh,qh),∀(vh,qh)∈Vmh×Σmh,

where the bilinear form and the linear form are defined as

 (6) ah(uh,ph;vh,qh) :=∑K∈Th∫K(∇×ph−kuh)⋅(∇×qh−kvh)dx +∑K∈Th∫K(∇×uh−kph)⋅(∇×vh−kqh)dx +∑f∈Fbhμhf∫f(n×uh)⋅(n×vh)ds,

and

 lh(vh,qh) :=−∑K∈Th∫Kf⋅vhdx+∑K∈Th∫K∇×qh⋅˜fdx +∑f∈Fbhμhf∫fn×vh⋅gds.

Then we focus on the error estimate to the problem (5). To do so, we first define the spaces and as

 Vh=Vmh+H0(curl),Σh=Σmh+H(curl).

We introduce the following energy norms for both the spaces:

 ∥u∥2u =∑K∈Th(∥u∥2L2(K)+∥∇×u∥2L2(K))+∑f∈Fh1hf∥[[n×u]]∥2L2(f),∀u∈Vh,

and

 ∥p∥2p =∑K∈Th(∥p∥2L2(K)+∥∇×p∥2L2(K))+∑f∈Fih1hf∥[[n×p]]∥2L2(f),∀p∈Σh,

and we define the energy norm as

 |||(u,p)|||2=∥u∥2u+∥p∥2p,∀(u,p)∈Vh×Σh.

It is easy to see that indeed defines a norm on and indeed defines a norm on . As a result, defines a norm on . Then we state the continuity result of the bilinear form with respect to the norm .

###### Lemma 1.

Let the bilinear form be defined as (6) with any positive , there exists a constant such that

 (7) |ah(u,p;v,q)|≤C|||(u,p)||||||(v,q)|||,

for any .

###### Proof.

Using the Cauchy-Schwartz inequality, we have that

 ∑K∈Th∫Kk2u⋅vdx ≤k2⎛⎝∑K∈Th∥u∥2L2(K)⎞⎠12⎛⎝∑K∈Th∥v∥2L2(K)⎞⎠12 ≤k2|||(u,p)||||||(v,q)|||.

Other terms that appear in the bilinear form (6) can be bounded similarly, which gives us the inequality (7) and completes the proof. ∎

In order to prove the coercivity of the bilinear form , we require the following lemmas.

###### Lemma 2.

For any , there exists a piecewise polynomial such that

 (8) ∥uh−vh∥2L2(Ω)≤C∑f∈Fhhf∥[[n×uh]]∥2L2(f), ∥uh−vh∥2u≤C∑f∈Fh1hf∥[[n×uh]]∥2L2(f),

and for any , there exists a piecewise polynomial such that

 (9) ∥ph−wh∥2L2(Ω)≤C∑f∈Fihhf∥[[n×ph]]∥2L2(f), ∥ph−wh∥2p≤C∑f∈Fih1hf∥[[n×ph]]∥2L2(f).
###### Proof.

We prove this lemma by using the similar techniques as those in [15, 20, 22]. In two dimensions, is a scalar-valued piecewise polynomial function. For scalar-valued functions, the space is equal to the space . By [20, Theorem 2.1], there exists a piecewise polynomial satisfying the estimate (9). For the vector-valued function , we will use Nédélec ’s elements of the second type in two dimensions to prove the result and the estimate (8). In three dimensions, and are both vector-valued piecewise polynomials and we will apply 3D Nédélec ’s elements of the second type to prove the two results. We primarily deduce that for the three-dimensional case and it is almost trivial to extend the proof to 2D for vector-valued functions.

We first give some properties about Nédélec ’s edge elements of the second type, which was introduced by Nédélec [28]. For a bounded domain , we define the polynomial space as , where is the space of homogeneous polynomials of degree . For a tetrahedral element and a polynomial

, three types of moments (degrees of freedom) of the Nédélec ’s elements associated with edges of

, faces of and itself are defined as follows,

 MeK(t) ={∫e(t⋅τe)qds,∀q∈Pm(e)},%foranyedgee∈∂K, MfK(t) ={1|f|∫ft⋅qds,∀q∈Dm−1(f)},for any face f∈∂K, MbK(t) ={∫Kt⋅qdx,∀q∈Dm−2(K)},

where denotes the unit vector along the edge. Further, for any polynomial , we define , and as the corresponding moments of . We let , , are the Lagrange bases of the space with respect to the moments, respectively. Then, any polynomial can be uniquely expressed as

 v=∑e∈E(K)Ne∑i=1viK,eϕiK,e+∑f∈F(K)Nf∑i=1viK,fϕiK,f+Nb∑i=1viK,bϕiK,b,

where and are the sets of edges and faces of the element , respectively.

Then we state the following estimates. For an element and any polynomial , there exists a constant such that

 (10) h−2K∥v∥2L2(K) +∥∇×v∥2L2(K)≤ Ch−1K ⎛⎝∑e∈E(K)Ne∑i=1(viK,e)2+∑f∈F(K)Nf∑i=1(viK,f)2+Nb∑i=1(viK,b)2⎞⎠.

For an interior face shared by and . For any polynomial and , there exists a constant such that

 (11) Nf∑i=1(viK1,f−viK2,f)2+∑e∈E(f)Ne∑i=1(viK1,e−viK2,e)2≤C∫f|nf×(v1−v2)|2ds,

where is the set of edges belonging to the face and is the unit outward normal on . For a boundary face , we let be the element such that . Then, there exists a constant such that

 (12)

The estimates (10), (11) and (12) are obtained from the equivalence of norms over finite dimensional spaces and the scaling argument. We refer to [15, 26] for the details of their proofs.

Then we construct two new piecewise polynomials and satisfying the estimates (9) and (8), respectively. To construct , for the element , we define its moments , and as follows,

 (13) wiK,e=∑˜K∈N(e)1|N(e)|pi˜K,e,∀e∈Eh,1≤i≤Ne,

and

 (14) wiK,f=∑˜K∈N(f)1|N(f)|pi˜K,f,∀f∈Fh,1≤i≤Nf,

and

 wiK,b=piK,b,1≤i≤Nb.

Here we let , be the sets of elements that contain the edge and the face , respectively, and we denote their cardinalities as and . Clearly, the above moments will yield a piecewise polynomial . Then we focus on the error between and . From (10), we deduce that

 (15) hα−1K∥ph−wh∥2L2(K)+hα+1K∥∇×(ph−wh)∥2L2(K)≤ ChαK⎛⎝∑e∈E(K)Ne∑i=1(piK,e−wiK,e)2+∑f∈F(K)Nf∑i=1(piK,f−wiK,f)2⎞⎠,α=±1,

on the element . By (13), we only need to consider the error for the edge which satisfies . For such an edge , we let be shared by a sequence of elements with and are two neighbouring elements. Then the edge will be shared by faces, which are where are the adjacent faces of the element , and we further have that . Then from (10) and the mesh regularity, we deduce that

 (16) Ne∑i=1hαK(piK,e−wiK,e)2 ≤C|N(e)|∑j=2Ne∑i=1hαK(piKe,1,e−wiKe,j,e)2 ≤C|N(e)|−1∑j=1Ne∑i=1hαK(piKe,j,e−wiKe,j+1,e)2 ≤C|N(e)|∑j=2∫fe,jhαfe,j∥[[n×ph]]∥2ds

where are the set of faces that contain the edge .

By (14), we consider the errors on interior faces, and from (11) we obtain that

 (17) Nf∑i=1hαK(piK,f−wiK,f)2 ≤C∑K′∈N(f)Nf∑i=1hαK(piK,f−piK′,f)2 ≤Chαf∥[[n×ph]]∥2L2(f),

for the face . Combining the estimates (15), (16) and (17), and summing over all elements arrives at the estimate (9).

The construction of the piecewise polynomial is similar. We define its corresponding moments , and as

 viK,e={∑˜K∈N(e)1|N(e)|ui˜K,e,∀e∈Eih,0,∀e∈Ebh,,1≤i≤Ne,

and

 viK,f={∑˜K∈N(f)1|N(f)|ui˜K,f,∀f∈Fih,0,∀f∈Fbh,,1≤i≤Nf,

and

 viK,b=uiK,b,1≤i≤Nb.

The above moments obviously imply that . We note that the moments of and are only different on boundary edges and faces. Hence the estimates (17) and (16) also hold for on interior edges and faces. On the element , we have that

 Ne∑i=1hαK(uiK,e−viK,e)≤C∑f∈F(e)hαf∥[[n×uh]]∥2L2(f),

for an interior edge , and

 Nf∑i=1hαK(uiK,f−viK,f)2≤Chαf∥[[n×uh]]∥2L2(f),

for an interior face . For boundary edges and faces, we directly apply the estimate (12) to obtain the analogous inequalities, which bring us that

 hα−1K∥uh −vh∥2L2(K)+hα+1K∥∇×(uh−vh)∥2L2(K) ≤C⎛⎝∑e∈E(K)∑f∈F(e)hαf∥[[n×uh]]∥2L2(f)+∑f∈F(K)hαf∥[[n×uh]]∥2L2(f)⎞⎠,α=±1.

We finally arrive at the estimate (8) by summing over all elements. This completes the proof. ∎

###### Lemma 3.

There exists a constant such that

 (18) ∥u∥u+∥p∥p≤C(∥∇×u−kp∥L2(Ω)+∥∇×p−ku∥L2(Ω)),

for any and any .

###### Proof.

We let

 ∇×u−kp=f1,∇×p−ku=f2.

For any , we have

 (∇×u,∇×ψ)L2(Ω)−k(p,∇×ψ)L2(Ω)=(f1,∇×ψ)L2(Ω),

and

 (∇×p,ψ)L2(Ω)−k(u,ψ)=(f2,ψ)L2(Ω).

Using the integration by parts to obtain that

 (∇×u,∇×ψ)L2(Ω)−k2(u,ψ)L2(Ω) =(f1,∇×ψ)L2(Ω)+k(f2,ψ)L2(Ω).

By [13, Theorem 5.2], we derive that

 ∥u∥H(curl) ≤Csupψ∈H0(curl)(∇×u,∇×ψ)L2(Ω)−k2(u,ψ)L2(Ω)∥ψ∥H(curl) ≤Csupψ∈H0(curl)(f1,∇×ψ)L2(Ω)+k(f2,ψ)L2(Ω)∥ψ∥H(curl) ≤C(∥f1∥L2(Ω)+∥f2∥L2(Ω)).

Since , we have that . Further, we deduce that

 ∥p∥p =∥p∥H(curl)≤C(∥p∥L2(Ω)+∥∇×p∥L2(Ω)) ≤C(∥f1∥L2(Ω)+∥∇×u∥L2(Ω)+∥f2∥L2(Ω)+∥u∥L2(Ω)) ≤C(∥f1∥L2(Ω)+∥f2∥L2(Ω)),

which gives us the estimate (18) and completes the proof. ∎

Now we are ready to state that the bilinear form is coercive with respect to the energy norms.

###### Lemma 4.

Let the bilinear form be defined as (6) with any positive , there exists a constant such that

 (19) ah(uh,ph;uh,ph)≥C|||(uh,ph)|||2,

for any .

###### Proof.

Clearly, we have that

 ah(uh,ph;uh, ph)=∑K∈Th(∥∇×ph−kuh∥2L2(K)+∥∇×uh−kph∥2L2(K)) +

By Lemma 2, for and