capbtabboxtable[][]

# A Yee-like finite-element scheme for
Maxwell’s equations on unstructured grids

B. Radu† and H. Egger†,∗ †Johann Radon Institute for Computational and Applied
Mathematics, Linz, Austria
∗Institute for Computational Mathematics, Johannes-Kepler University Linz,
Austria<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

A novel finite element scheme is studied for solving the time-dependent
Maxwell’s equations on unstructured grids efficiently. Similar to the
traditional Yee scheme, the method has one degree of freedom for most edges
and a sparse inverse mass matrix. This allows for an efficient realization by
explicit time-stepping without solving linear systems. The method is
constructed by algebraic reduction of another underlying finite element scheme
which involves two degrees of freedom for every edge. Mass-lumping and
additional modifications are used in the construction of this method to allow
for the mentioned algebraic reduction in the presence of source terms and
lossy media later on. A full error analysis of the underlying method is
developed which by construction also carries over to the reduced scheme and
allows to prove convergence rates for the latter. The efficiency and accuracy
of both methods are illustrated by numerical tests. The proposed schemes and
their analysis can be extended to structured grids and in special cases the
reduced method turns out to be algebraically equivalent to the Yee scheme. The
analysis of this paper highlights possible difficulties in extensions of the
Yee scheme to non-orthogonal or unstructured grids, discontinuous material
parameters, and non-smooth source terms, and also offers potential remedies.

> Keywords: Maxwell’s equations, explicit time stepping, finite element
> methods, mass-lumping, Yee-like schemes
>
> AMS subject classifications: 35Q61, 65M60

##  1\. Introduction

We consider the efficient numerical simulation of electromagnetic wave
propagation through linear non-dispersive lossy media. As the basic
mathematical model, we use Maxwell’s equations in second-order form, i.e.

$\displaystyle\varepsilon\partial_{tt}{E}+\sigma\partial_{t}{E}+\operatorname{curl}(\nu\operatorname{curl}{E})$
$\displaystyle=f,\qquad\text{in }\Omega,\,t>0,$ (1) $\displaystyle
n\times(\nu\operatorname{curl}E)$ $\displaystyle=g,\qquad\text{on
}\partial\Omega,\,t>0,$ (2)

Here $\varepsilon$ and $\sigma$ denote the electric permittivity and
conductivity, and $\nu$ the magnetic reluctivity of the medium; further $E$ is
the electric field intensity, $f$, $g$ describe the excitation by volume or
surface currents, and $\Omega$ is the computational domain.

One of the standard approaches for solving (1)–(2) numerically is the Yee
scheme and its relatives, e.g., the finite difference time domain (FDTD)
method and the finite-integration technique (FIT); see [38], [34, 35], and
[36, 37] for details. The common algebraic form of all these methods reads

$\displaystyle\mathsf{M}_{\varepsilon}\frac{\mathsf{E}^{n+1}-2\mathsf{E}^{n}+\mathsf{E}^{n}}{\tau^{2}}+\mathsf{M}_{\sigma}\frac{\mathsf{E}^{n+1}-\mathsf{E}^{n-1}}{2\tau}+\mathsf{K}_{\nu}\mathsf{E}^{n}$
$\displaystyle=\mathsf{f}^{n}+\mathsf{g}^{n}.$ (3)

This allows for an efficient realization by explicit time stepping, whenever
the mass matrices $\mathsf{M}_{\varepsilon}$, $\mathsf{M}_{\sigma}$ are
diagonal. On structured orthogonal grids and for homogeneous material
distributions, such constructions are possible and lead to second-order
accurate approximations in space and time; see e.g. [6, 35, 38]. These results
remain valid for non-uniform orthogonal grids [28]; further see [6, 7] for
related higher-order methods. The extension of the Yee scheme to discontinuous
material parameters or source terms and to non-orthogonal or even unstructured
grids is, however, not straightforward; see [17, 20, 32, 33] and [2, 4, 5] for
some attempts in the latter case. To the best of our knowledge, a rigorous
convergence analysis of Yee-like schemes in such situations is not available
up to date.

An alternative approach that allows the construction of stable schemes for
structured as well as unstructured grids, and also for discontinuous material
parameters, is offered by finite element discretizations. A rigorous error
analysis for various methods can be developed; see e.g. [23, 24, 25]; we refer
to [6, 18] for an overview of different approaches, their analysis, and
further references. Together with the leap-frog time-stepping scheme, the
finite element approximation of (1)–(2) again leads to algebraic systems of
the form (3). In contrast to the finite difference schemes mentioned above,
the mass matrices $\mathsf{M}_{\varepsilon}$ and $\mathsf{M}_{\sigma}$ are,
however, usually not diagonal, and the realization of (3) thus requires the
solution of linear systems in every time step. Therefore, the use of implicit
time-stepping schemes has been advocated in the literature [22, 26].

This inherent disadvantage of finite element approximations can be overcome by
mass-lumping, which aims at replacing $\mathsf{M}_{\varepsilon}$,
$\mathsf{M}_{\sigma}$ in (3) by diagonal or block-diagonal matrices. For
orthogonal grids, corresponding schemes have been proposed in [7] and
modifications for unstructured grids have been considered in [11, 12]; we also
refer to [6] for a detailed discussion. A full convergence analysis for a
second order method has been given recently in [10, 31], and in [19] higher
order approximations of Maxwell’s equations were obtained based on staggered
grid complexes. Another approach for generating non-conforming finite element
approximations with block-diagonal mass matrices is provided by discontinuous-
Galerkin methods [8, 16]. For low-order approximations however, these methods
suffer from a substantial increase in the number of degrees of freedom; we
refer to [15] for comparison with mass-lumping schemes in the context of
elastodynamics.

### Main contributions

In this paper, we propose and analyze a fully discrete finite element
approximation for (1)–(2) which can be considered as a natural extension of
the Yee scheme to unstructured grids. This method involves only one degree of
freedom for most edges and it can be realized efficiently as an explicit time-
stepping scheme; moreover, a full convergence analysis is possible. The
approach and its analysis are based on the following key ideas, already
presented in [9]: First, a finite element method is considered involving two
degrees of freedom on each edge. Mass lumping is achieved by appropriate
numerical quadrature and a corresponding choice of basis functions, and a full
error analysis can be developed based on standard arguments. In a second step,
the method is then reduced on the algebraic level to a scheme that involves
only one degree of freedom for most edges. While a direct variational
characterization of this scheme is no longer possible, its close relation to
the first method nevertheless allows to develop a rigorous error analysis.
Compared to [9], some non-trivial modifications are required in the
construction and the analysis of the proposed methods to deal with non-trivial
conductivities $\sigma$ and inhomogeneous problem data $f$, $g$. On orthogonal
grids and for special situations, the reduced method becomes equivalent to the
traditional Yee scheme. The analysis of this paper therefore also offers some
recipes for extensions of the Yee scheme to non-orthogonal grids and
discontinuous parameters and problem data. The algebraic form of the numerical
scheme obtained after reduction also shares similarities with the approaches
of [4, 5], which will be briefly discussed at the end of the manuscript.

### Outline

In Section 2, we introduce our notation and basic assumptions. The methods and
main results are presented in Section 3. The error analysis is developed in
Section 4 and the implementation of the method is discussed in some detail in
Section 5. For an illustration of our theoretical results, some computational
tests are presented in Section 6, and we close with a short discussion.

##  2\. Notation and basic assumptions

We consider a three-dimensional setting and assume that
$\Omega\subset\mathbb{R}^{3}$ is a bounded polyhedral Lipschitz domain.
Throughout the paper, we assume that

  1. (A1)

$\mathcal{T}_{h}$ is a geometrically-conforming non-overlapping partition of
$\Omega$ into tetrahedral elements $K\in\mathcal{T}_{h}$; the mesh
$\mathcal{T}_{h}$ is shape-regular and quasi-uniform, i.e., $\gamma
h_{K}^{3}\leq|K|\leq h_{K}^{3}$ and $\delta h\leq h_{K}\leq h$ for all
$K\in\mathcal{T}_{h}$ with $\gamma,\delta>0$.

As usual, $h_{K}$ and $|K|$ denote the diameter and the volume of the element
$K\in\mathcal{T}_{h}$, and $h=\max_{K}h_{K}$ is the global mesh size; see [13]
for further details. All results are presented in detail for the three-
dimensional setting and unstructured grids, but the translation to two
dimensions and structured grids is possible and more or less straightforward;
see [31, Chapter 5] and Section 6.

The material parameters and problem data are required to be sufficiently
regular and to satisfy the usual physical bounds. Moreover, material
discontinuities shall be resolved by the mesh. For ease of presentation, we
thus assume that

  1. (A2)

$\varepsilon$, $\sigma$, $\nu\in P_{0}(\mathcal{T}_{h})$ and
$0\leq\varepsilon,\sigma,\nu\leq\overline{c}$ and
$\varepsilon,\nu\geq\underline{c}$ for some constants
$\underline{c},\overline{c}>0$.

  2. (A3)

$f:[0,T]\to L^{2}(\Omega)^{3}$ and $g:[0,T]\to L^{2}(\partial\Omega)^{3}$ are
smooth functions of time, with $f(0)=0$ and $g(0)=0$, and the initial
conditions are $E(0)=\partial_{t}E(0)=0$.

We write $P_{k}(\mathcal{T}_{h})=\\{v:v|_{K}\in P_{k}(K),\ \forall
K\in\mathcal{T}_{h}\\}$ for the space of piecewise polynomials of degree $k$
on $\mathcal{T}_{h}$. Under these assumptions, the existence of a unique
solution to (1)–(2) can be established by standard arguments; see e.g. [21,
23]. Also, more general conditions could be treated with minor modifications
to our analysis.

### Function spaces

We use standard symbols $L^{2}(\Omega)$, $H^{1}(\Omega)$, and
$H(\operatorname{curl};\Omega)$ for the spaces of square-integrable function
with square integrable weak gradients and curls, respectively; see [27] for
details. The norm of a space $X$ is denoted by $\|\cdot\|_{X}$ and we will
often write $\langle a,b\rangle=\int_{\Omega}a\cdot b\,dx$ and $\langle
a,b\rangle_{\partial}=\int_{\partial\Omega}a\cdot b\,ds(x)$ for the scalar
products of two functions in $L^{2}(\Omega)^{3}$ and
$L^{2}(\partial\Omega)^{3}$, respectively. For the error analysis, we use
$H^{k}(\mathcal{T}_{h})=\\{v\in L^{2}(\Omega):v|_{K}\in H^{k}(K),\ \forall
K\in\mathcal{T}_{h}\\}$, to denote spaces of piecewise smooth functions, and
we write $\|v\|_{H^{k}(\mathcal{T}_{h})}=(\sum_{K}\|v\|_{H^{k}(K)}^{2})^{1/2}$
for the corresponding norms. We further denote by $L^{p}(0,T;X)$ the Bochner
spaces of functions with values in $X$ whose $p$-th power is integrable in
time. The functions in $W^{k,p}(0,T;X)$ further have weak derivatives in
$L^{p}(0,T;X)$. For brevity, we will sometimes write $L^{p}(X)$ and
$W^{k,p}(X)$, and omit explicit reference to the time interval.

### Space discretization

For the spatial approximation of the electric field $E$, we first consider
Nédélec finite elements of type II, i.e.

$\displaystyle V_{h}=\\{v_{h}\in H(\operatorname{curl};\Omega):v_{h}|_{K}\in
P_{1}(K)^{3}\ \forall K\in\mathcal{T}_{h}\\}.$ (4)

This amounts to the space of piecewise linear vector-valued functions with
tangential continuity across element boundaries; see [1, 27, 29]. Let us
recall that functions in $V_{h}$ have two degrees of freedom for every edge,
and the canonical interpolation operator for the space is given by
$\Pi_{h}:H(\operatorname{curl};\Omega)\cap H^{1}(\mathcal{T}_{h})^{3}\to
V_{h}$, with

$\displaystyle\int_{e}\Pi_{h}v\cdot\tau_{e}\,p_{e}\,ds$
$\displaystyle=\int_{e}v\cdot\tau_{e}\,p_{e}\,ds\qquad\forall p_{e}\in
P_{1}(e),\ e\in\mathcal{E}_{h}.$ (5)

Here $\mathcal{E}_{h}=\\{e_{ij}:i<j\\}$ is the set of edges
$e_{ij}=(v_{i},v_{j})$ and $\tau_{e}$ is the unit tangential vector on
$e=e_{ij}$ pointing from vertex $v_{i}$ to $v_{j}$ with $i<j$.

For a given subset $\widetilde{\mathcal{E}}_{h}\subset\mathcal{E}_{h}$ of
edges, we define the corresponding subspace

$\displaystyle\widetilde{V}_{h}=\\{v_{h}\in V_{h}:v_{h}\cdot\tau_{e}\in
P_{0}(e)\ \forall e\in\widetilde{\mathcal{E}}_{h}\\}\subset V_{h},$ (6)

which consists of functions in $V_{h}$ having only constant tangential trace,
and therefore only one degree of freedom for edges
$e\in\widetilde{\mathcal{E}}_{h}$. We thus call $\widetilde{V}_{h}$ the
_reduced space_ in the following. The canonical interpolation operator for
$\widetilde{V}_{h}$ is given by
$\widetilde{\Pi}_{h}:H^{1}(\mathcal{T}_{h})^{3}\cap
H(\operatorname{curl};\Omega)\to\widetilde{V}_{h}$ with

$\displaystyle\int_{e}\widetilde{\Pi}_{h}v\cdot\tau_{e}\,p_{e}\,ds$
$\displaystyle=\int_{e}v\cdot\tau_{e}\,p_{e}\,ds\qquad\forall p_{e}\in
P_{k_{e}}(e),\ e\in\mathcal{E}_{h},$ (7)

and with degree $k_{e}=0$ for $e\in\widetilde{\mathcal{E}}_{h}$ and $k_{e}=1$
for $e\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}$. For any choice
of $\widetilde{\mathcal{E}}_{h}\subset\mathcal{E}_{h}$, we have the inclusions

$\displaystyle\mathcal{N}_{0}(\mathcal{T}_{h})\cap
H(\operatorname{curl};\Omega)\subset\widetilde{V}_{h}\subset\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}(\mathcal{T}_{h})\cap
H(\operatorname{curl};\Omega).$

Here $\mathcal{N}_{0}$ and $\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$ denote
the lowest order Nédélec elements of type I and II. This ensures good
approximation properties for both spaces; see [1, 30, 29] for details. For the
choice $\widetilde{\mathcal{E}}_{h}=\mathcal{E}_{h}$ or
$\widetilde{\mathcal{E}}_{h}=\emptyset$, one of the two inclusions becomes an
identity.

### Mass lumping

For the approximation of some of the integrals arising in the finite element
approximation of (1)–(2), we use numerical integration by the vertex rule;
this will allow for mass-lumping later on. For ease of notation, we introduce

$\displaystyle\langle
a,b\rangle_{h}=\sum\nolimits_{T}\tfrac{|T|}{4}\sum\nolimits_{v_{i}\in
T}a(v_{i})\cdot b(v_{i}).$ (8)

Here $a,b$ are assumed to be piecewise smooth vector-valued functions over the
mesh $\mathcal{T}_{h}$. Let us note that the quadrature is exact if $a\cdot
b\in P_{1}(\mathcal{T}_{h})$.

### Time discretization

Let $\tau=T/N$ and $t^{n}=n\tau$ be a sequence of uniformly spaced time steps.
Further let $(a^{n})_{n\geq 0}\subset X$ be a sequence in some vector space
$X$. Then

$\displaystyle\partial_{\tau\tau}a^{n}:=\frac{a^{n+1}-2a^{n}+a^{n-1}}{\tau^{2}}\qquad\text{and}\qquad\partial_{\tau}a^{n-1/2}:=\frac{a^{n}-a^{n-1}}{\tau},$
(9)

are used to denote the standard central difference quotients approximating the
second and first derivative at time $t=t^{n}$ and
$t=t^{n-1/2}=t^{n}-\frac{1}{2}\tau$, respectively.

##  3\. Main results

For the numerical approximation of (1)–(2) with homogeneous initial
conditions, see assumption (A3), we now consider the following fully discrete
scheme.

######  Method 3.1.

Let $\widetilde{\mathcal{E}}_{h}\subset\mathcal{E}_{h}$ and
$\widetilde{\Pi}_{h}:V_{h}\to\widetilde{V}_{h}$ denote the appropriate
projection. Find $E_{h}^{n}\in V_{h}$, $0\leq n\leq N$, with
$E_{h}^{0}=E_{h}^{1}=0$ and such that

$\displaystyle\langle(\varepsilon+\tfrac{\tau}{2}\sigma)\partial_{\tau\tau}E_{h}^{n},v_{h}\rangle_{h}+\langle\sigma\widetilde{\Pi}_{h}$
$\displaystyle\partial_{\tau}E_{h}^{n-1/2},\widetilde{\Pi}_{h}v_{h}\rangle_{h}+\langle\nu\operatorname{curl}E_{h}^{n},\operatorname{curl}v_{h}\rangle$
(10) $\displaystyle=\langle f(t^{n}),v_{h}\rangle+\langle
g(t^{n}),v_{h}\rangle_{\partial\Omega},\qquad\forall\ 1\leq n<N.$

The solution (sequence) will be abbreviated by the symbol
$E_{h}=(E_{h}^{n})_{0\leq n\leq N}$.

Independently of the choice of the set of edges
$\widetilde{\mathcal{E}}_{h}\subset\mathcal{E}_{h}$, on which the number of
degrees of freedom is reduced, the implementation of the method leads to a
finite-dimensional recursion of the form

$\displaystyle\mathsf{M}_{\varepsilon+\tau\sigma/2}\frac{\mathsf{E}^{n+1}-2\mathsf{E}^{n}+\mathsf{E}^{n-1}}{\tau^{2}}+\widehat{\mathsf{M}}_{\sigma}\frac{\mathsf{E}^{n}-\mathsf{E}^{n-1}}{\tau}+\mathsf{K}_{\nu}\mathsf{E}^{n}$
$\displaystyle=\mathsf{f}^{n}+\mathsf{g}^{n}.$ (11)

The well-posedness of the discretization scheme then follows immediately from
the regularity of the matrix $\mathsf{M}_{\varepsilon+\tau\sigma/2}$, which is
a direct consequence of Lemma 4.3 below. In Section 5, we further show that an
appropriate choice of basis functions for the space $V_{h}$, adopted to the
numerical quadrature, leads to a block-diagonal mass matrix
$\mathsf{M}_{\varepsilon+\tau\sigma/2}$, such that time stepping in (11) can
be realized efficiently.

The algebraic form (11) reveals that Method 3.1 is based on an explicit time-
stepping scheme and a restriction on the time step size $\tau$ is, therefore,
required to ensure discrete stability and convergence with $h,\tau\to 0$. We
thus assume that

  1. (A4)

the time step $\tau>0$ is chosen to satisfy for all ${v}_{h}\in{V}_{h}$ the
inequality

$\displaystyle\tfrac{{\tau}^{2}}{4}\langle\nu\operatorname{curl}{v}_{h},\operatorname{curl}{v}_{h}\rangle+\tfrac{\tau}{2}\left|\langle\sigma\widetilde{\Pi}_{h}{v}_{h},\widetilde{\Pi}_{h}{v}_{h}\rangle-\langle\sigma{v}_{h},{v}_{h}\rangle\right|\leq\tfrac{1}{2}\langle\varepsilon{v}_{h},{v}_{h}\rangle_{h}.$

For conductivity $\sigma=0$ and
$\langle\cdot,\cdot\rangle_{h}=\langle\cdot,\cdot\rangle$, this assumption
reduces to the usual CFL condition as used, e.g., in [6, 18]. Under our
assumptions on the mesh and the model parameters, one can verify that
$\tau\leq Ch$ for some appropriate constant $C>0$ is sufficient to guarantee
condition (A4); see Section 6. In practice, an appropriate time step $\tau$
satisfying (A4) can be found by performing a few vector iterations.

To guarantee good approximation properties, we further need some restriction
on the set $\widetilde{\mathcal{E}}_{h}$ of edges, on which the polynomial
order is reduced. We thus require that

  1. (A5)

$\sigma$ is continuous across edges $e\in\widetilde{\mathcal{E}}_{h}$ inside
$\Omega$, and $\sigma=0$ for all $e\in\widetilde{\mathcal{E}}_{h}$ on
$\partial\Omega$.

This condition simply means that we stay with two degrees of freedom on edges
where the conductivity $\sigma$ is either discontinuous or non-trivial at the
boundary. The reason for this restriction will become clear from the error
analysis given in the next section and its necessity will be illustrated by
numerical tests.

For ease of notation, we write $\|u_{h}\|_{\ell_{\infty}(X)}=\max_{0\leq k\leq
N-1}\|u_{h}^{n+1/2}\|_{X}$ in the following statements, and we write
$u^{n+1/2}=u(t^{n+1/2})$ for functions $u$ that are continuous in time. We
further denote by $\widehat{u}^{\,n+1/2}=\frac{1}{2}(u^{n+1}+u^{n})$ the
average at intermediate time steps. This allows us to present our first main
result as follows.

######  Theorem 3.2.

Let $E$ be a sufficiently smooth solution of (1)–(2) and let (A1)–(A5) hold.
Then Method 3.1 is well-defined and the discrete solution
$E_{h}=(E_{h}^{n})_{n}$ satisfies

$\displaystyle\|\partial_{t}E-\partial_{\tau}E_{h}\|_{\ell_{\infty}(L^{2}(\Omega))}$
$\displaystyle+\|\operatorname{curl}(E-\widehat{E}_{h})\|_{\ell_{\infty}(L^{2}(\Omega))}\leq
C(E)h+C^{\prime}(E)\tau^{2},$

with constants

$\displaystyle C(E)$
$\displaystyle=\|\partial_{t}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}+\|\operatorname{curl}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}+\|\partial_{tt}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}$
$\displaystyle\qquad\qquad\qquad+\|\operatorname{curl}\partial_{t}E\|_{L^{1}(H^{1}(\mathcal{T}_{h}))}+\|\partial_{t}E\|_{L^{1}(H^{1}(\Omega))},$
$\displaystyle C^{\prime}(E)$
$\displaystyle=\|\partial_{tttt}E\|_{L^{1}(L^{2}(\Omega))}+\|\partial_{ttt}E\|_{L^{1}(L^{2}(\Omega))}+\|\operatorname{curl}\partial_{tt}E\|_{L^{1}(H^{1}(\mathcal{T}_{h}))}.\qquad$

The implementation leads to a time-stepping scheme (11), and for an
appropriate choice of a basis, the matrix
$\mathsf{M}_{\varepsilon+\tau\sigma/2}$ is block-diagonal, while
$\widehat{\mathsf{M}}_{\sigma}$ and $\mathsf{K}_{\nu}$ are sparse.

The convergence results is proven in Section 4, while the algebraic structure
of the scheme is derived in Section 5. Let us note that the assertions hold,
in particular, for the choice $\widetilde{\mathcal{E}}_{h}=\emptyset$, for
which the projection $\widetilde{\Pi}_{h}$ drops out and the method as well as
its implementation become somewhat simpler; see [31] for details.

### The Yee-like scheme

The reason for introducing the projection $\widetilde{\Pi}_{h}$ in Method 3.1
lies in the following important observation, which is summarized as our second
main result and leads to the reduced Yee-like scheme announced in the
introduction.

######  Theorem 3.3.

Let the assumptions of Theorem 3.2 be valid and $E_{h}=(E_{h}^{n})_{n}$ denote
the solution of Method 3.1. Further define
$\widetilde{E}_{h}^{n}=\widetilde{\Pi}_{h}E_{h}^{n}$ for all $n\geq 0$. Then

$\displaystyle\|\partial_{t}E-\partial_{\tau}\widetilde{E}_{h}\|_{\ell_{\infty}(L^{2}(\Omega))}$
$\displaystyle+\|\operatorname{curl}(E-\widehat{\widetilde{E}}_{h})\|_{\ell_{\infty}(L^{2}(\Omega))}\leq
C(E)h+C^{\prime}(E)\tau^{2},$

with constants $C(E)$ and $C^{\prime}(E)$ of the same form as in Theorem 3.2.
Moreover, the coefficients of the solution $\widetilde{E}_{h}$ can be computed
by the time-stepping scheme

$\displaystyle\partial_{\tau\tau}\widetilde{\mathsf{E}}^{n}$
$\displaystyle=\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}(-\widetilde{\mathsf{M}}_{\sigma}\partial_{\tau}\widetilde{\mathsf{E}}^{n-1/2}-\widetilde{\mathsf{K}}_{\nu}\widetilde{\mathsf{E}}^{n})+\widetilde{\mathsf{F}}^{n}+\widetilde{\mathsf{G}}^{n},$
(12)

and for an appropriate choice of basis for $\widetilde{V}_{h}$, the matrices
$\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}$,
$\widetilde{\mathsf{M}}_{\sigma}$ and $\widetilde{\mathsf{K}}_{\nu}$ are
sparse and the vectors $\widetilde{\mathsf{F}}^{n}$ and
$\widetilde{\mathsf{G}}^{n}$ can be cheaply assembled from $\mathsf{f}^{n}$
and $\mathsf{g}^{n}$ in (11).

######  Remark 3.4.

If $\sigma=0$, we may choose $\widetilde{\mathcal{E}}_{h}=\mathcal{E}_{h}$ and
obtain an explicit time-stepping method for Maxwell’s equations with exactly
one degree of freedom per edge. For orthogonal grids and homogeneous data
$f,g\equiv 0$, the presented approach becomes equivalent to the Yee scheme;
see [9] for details.

The two assertions of Theorem 3.3 are again proven in the following two
sections.

######  Remark 3.5.

Let us emphasize that the matrix
$\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}$ in (12) has a sparse
inverse, but it is not a sparse matrix by itself. In contrast to (11), which
corresponds to Method 3.1, we can not give a variational characterization of
the scheme (12) in closed form. This poses a severe challenge for the analysis
of this method which can be overcome only by a somewhat non-standard analysis.

##  4\. Proof of convergence rates

In this section, we establish the convergence rates stated in Theorem 3.2 and
3.3. To be able to do so, we require a couple of auxiliary results, which are
stated first.

###  4.1. Projection operators

For later reference, we collect some well-known properties of projection
operators arising in our error analysis below.

######  Lemma 4.1.

Let (A1) hold, $\widetilde{\mathcal{E}}_{h}\subset\mathcal{E}_{h}$, and
$\Pi_{h}$, $\widetilde{\Pi}_{h}$ be defined as in Section 2. Then

$\displaystyle\begin{split}\|{E}-\widetilde{\Pi}_{h}{E}\|_{L^{2}(K)}&\leq
Ch\|{E}\|_{H^{1}(K)},\\\
\|\operatorname{curl}({E}-\widetilde{\Pi}_{h}{E})\|_{L^{2}(K)}&\leq
Ch\|\operatorname{curl}{E}\|_{H^{1}(K)},\end{split}$ (13)

for all $E\in H^{1}(\mathcal{T}_{h})$ and $K\in\mathcal{T}_{h}$ with a
constant $C$ depending only on $\gamma,\delta$ in assumption (A1). The same
estimates also hold for the projection operator $\Pi_{h}$.

The proof of the assertions follows from the arguments given in [1, Sec. 2.5].
In our analysis, we will also make use of the $L^{2}$-orthogonal projection
$\pi^{0}_{\omega}:L^{2}(\omega)\to P_{0}(\omega)$ to constants for certain
subsets $\omega\subset\Omega$, which is defined by

$\displaystyle\int_{\omega}\pi^{0}_{\omega}v\,dx=\int_{\omega}v\,dx.$ (14)

The same symbol will also be used for the projection of vector-valued
functions.

######  Lemma 4.2.

Let (A1) hold and $\omega=K$ for $K\in\mathcal{T}_{h}$ or
$\omega=\bigcup_{K\cap e=e}K$ for $e\in\mathcal{E}_{h}$. Then the projection
error can be estimated by

$\displaystyle\|v-\pi^{0}_{\omega}v\|_{L^{2}(\omega)}\leq
Ch\|v\|_{H^{1}(\omega)}$ (15)

for all $v\in H^{1}(\omega)$ with constant $C$ depending only on
$\gamma,\delta$ in assumption (A1).

The proof of this assertion is based on the Poicaré–Friedrichs inequality [14,
Thm. 1.1] and standard scaling arguments; also see [3, Ch. 4]. Here we use
that the sets $\omega$ appearing in the lemma are uniformly star-shaped with
respect to balls of size $h$, which follows from assumption (A1) on the mesh,
and hence $C$ is universal.

We will further write $\pi^{0}_{h}:L^{2}(\Omega)\to P_{0}(\mathcal{T}_{h})$
for the projection to piecewise constants over the mesh $\mathcal{T}_{h}$,
defined by $(\pi^{0}_{h}v)|_{K}=\pi^{0}_{K}(v|_{K})$, and note that

$\displaystyle\|v-\pi^{0}_{h}v\|_{L^{2}(K)}\leq Ch\|v\|_{H^{1}(K)},$ (16)

which follows immediately from the assertion of the previous lemma.

###  4.2. Properties of the quadrature rule

As a second ingredient, we now state some elementary facts about the
quadrature rule introduced in Section 2.

######  Lemma 4.3.

Let (A1) hold and $\alpha\in P_{0}(\mathcal{T}_{h})$ with $\alpha\geq 0$. Then

$\displaystyle c\,\langle\alpha v_{h},v_{h}\rangle\leq\langle\alpha
v_{h},v_{h}\rangle_{h}\leq C\,\langle\alpha v_{h},v_{h}\rangle\quad\forall
v_{h}\in V_{h},$ (17)

with uniform constants $c,C>0$ depending only on the bounds in assumption
(A1).

###### Proof.

We consider a single element $K\in\mathcal{T}_{h}$ and abbreviate $\langle
a,b\rangle_{K}=\int_{K}a\cdot b\,dx$ and $\langle
a,b\rangle_{h,K}=\frac{|K|}{4}\sum_{i}a(v_{i})\cdot b(v_{i})$. By mapping to
the reference element, using the finite dimensionality of $P_{1}(K)^{3}$, and
noting that $\alpha\geq 0$ is piecewise constant, we get

$\displaystyle c\,\langle\alpha v_{h},v_{h}\rangle_{K}\leq\langle\alpha
v_{h},v_{h}\rangle_{h,K}\leq C\,\langle\alpha v_{h},v_{h}\rangle_{K},$

with constants $c,C$ that only depend on the shape regularity of the element
$K$. The assertion of the lemma then follows by summation over all elements
$K\in\mathcal{T}_{h}$. ∎

As a direct consequence of the previous result, we obtain the following
assertions:

######  Corollary 4.4.

Let (A1) hold. Then $\|v_{h}\|_{h}^{2}=\langle v_{h},v_{h}\rangle_{h}$ defines
a norm on $V_{h}$ and

$\displaystyle c_{1}\|v_{h}\|_{L^{2}(\Omega)}\leq\|v_{h}\|_{h}^{2}\leq
c_{2}\|v_{h}\|_{L^{2}(\Omega)}^{2}\quad\forall v_{h}\in V_{h}.$ (18)

Further let $\alpha\in P_{0}(\mathcal{T}_{h})$ with
$0<\underline{\alpha}\leq\alpha\leq\overline{\alpha}$. Then $\langle\alpha
u_{h},v_{h}\rangle_{h}$ defines a continuous and elliptic symmetric bilinear
form on $V_{h}$, more precisely

$\displaystyle
c_{1}\underline{\alpha}\|u_{h}\|_{L^{2}(\Omega)}^{2}\leq\langle\alpha
u_{h},u_{h}\rangle_{h}\quad\text{and}\quad\langle\alpha
u_{h},v_{h}\rangle_{h}\leq
c_{2}\overline{\alpha}\|u_{h}\|_{L^{2}(\Omega)}\quad\forall u_{h},v_{h}\in
V_{h}.$

These properties immediately imply the well-posedness of Method 3.1. As a next
ingredient for our analysis, we analyze the quadrature error.

######  Lemma 4.5.

Let assumption (A1) hold and let

$\displaystyle\delta_{h}(\alpha u_{h},v_{h})\coloneqq\langle\alpha
u_{h},v_{h}\rangle_{h}-\langle\alpha u_{h},v_{h}\rangle$

denote the quadrature error for some $\alpha\in P_{0}(\mathcal{T}_{h})$. Then

$\displaystyle|\delta_{h}(\alpha\widetilde{\Pi}_{h}{u},{v}_{h})|\leq
Ch\|u\|_{H^{1}(\mathcal{T}_{h})}\|v_{h}\|_{L^{2}(\Omega)}\qquad\forall u\in
H^{1}(\mathcal{T}_{h})^{3},v_{h}\in V_{h},$

with constant $C$ depending only on $\|\alpha\|_{L^{\infty}(\Omega)}$ and the
constants in assumption (A1).

###### Proof.

We define the local error $\delta_{h,K}(\alpha
u_{h},v_{h})\coloneqq\langle\alpha u_{h},v_{h}\rangle_{h,K}-\langle\alpha
u_{h},v_{h}\rangle_{K}$, and split

$\displaystyle|\delta_{h,K}(\alpha\widetilde{\Pi}_{h}{u},{v}_{h})|$
$\displaystyle\leq|\alpha\delta_{h,K}(\pi^{0}_{h}{u},{v}_{h})|+|\alpha\delta_{K}(\widetilde{\Pi}_{h}{u}-\pi^{0}_{h}{u},v_{h})|=(i)+(ii),$

where we used that $\alpha$ is piecewise constant. Since the quadrature rule
integrates linear polynomials exactly, we obtain $(i)=0$. The second term can
again be bounded elementwise. We may therefore omit $\alpha$ and obtain

$\displaystyle|\delta_{K}(\widetilde{\Pi}_{h}{u}-\pi^{0}_{h}{u},{v}_{h})|$
$\displaystyle\leq(1+c_{1})\|\widetilde{\Pi}_{h}{u}-\pi^{0}_{K}{u}\|_{L^{2}(K)}\|{v}_{h}\|_{L^{2}(K)}$
$\displaystyle\leq
C\big{(}\|{u}-\pi^{0}_{h}{u}\|_{L^{2}(K)}+\|{u}-\widetilde{\Pi}_{h}{u}\|_{L^{2}(K)}\big{)}\|{v}_{h}\|_{L^{2}(K)}$
$\displaystyle\leq C^{\prime}h\,\|u\|_{H^{1}(K)}\|{v}_{h}\|_{L^{2}(K)}.$

Here we used the assertions of Corollary 4.4 and the Cauchy-Schwarz inequality
in the first step, and the projection error estimates of Lemma 4.1 and 4.2 in
the last. Scaling by the constant $\alpha$ and summation over all elements
leads to the assertion. ∎

###  4.3. Estimates for the loss term

We now present a particular approximation property, which explains why the
degree of approximation can be reduced on the edges in the set
$\widetilde{\mathcal{E}}_{h}$ satisfying assumption (A5) without decreasing
the accuracy.

######  Lemma 4.6.

Let $\lambda_{i}\in P_{1}(\mathcal{T}_{h})\cap H^{1}(\Omega)$ denote the
barycentric coordinates defined by $\lambda_{i}(v_{j})=\delta_{ij}$ for all
vertices $v_{j}$ of the mesh. To every edge $e_{ij}\in\mathcal{E}_{h}$, we
define

$\displaystyle\Phi_{ij}=\lambda_{i}\nabla\lambda_{j}\qquad\text{and}\qquad\Phi_{ji}=-\lambda_{j}\nabla\lambda_{i}.$
(19)

These basis functions are linearly independent and
$V_{h}=\operatorname{span}\\{\Phi_{ij},\Phi_{ji}:e_{ij}\in\mathcal{E}_{h}\\}$.
Hence the functions defined in (19) comprise a basis for $V_{h}$.

The assertion follows immediately from the considerations in [1, 29]. In the
subsequent analysis, we will also make use of the following norm equivalence.

######  Lemma 4.7.

Let (A1) hold and let $v_{h}\in V_{h}$, i.e.,
$v_{h}=\sum_{e_{ij}\in\mathcal{E}_{h}}\mathsf{v}_{ij}\Phi_{ij}+\mathsf{v}_{ji}\Phi_{ji}$
for appropriate coefficients $\mathsf{v}_{ij},\mathsf{v}_{ji}\in\mathbb{R}$.
Further let $\omega(e_{ij})=\bigcup_{K\cap e_{ij}=e_{ij}}K$ denote the patch
of elements containing the edge $e_{ij}$. Then

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{h}^{2}\coloneqq\sum\nolimits_{e_{ij}\in\mathcal{E}_{h}}\mathsf{v}_{ij}^{2}\|\Phi_{ij}\|^{2}_{L^{2}(\omega(e_{ij}))}+\mathsf{v}_{ji}^{2}\|\Phi_{ji}\|^{2}_{L^{2}(\omega(e_{ij}))}$

defines a norm on $V_{h}$ which is equivalent to $\|\cdot\|_{L^{2}(\Omega)}$.
More precisely, one has

$\displaystyle
c_{1}^{\prime}\|v_{h}\|_{L^{2}(\Omega)}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{h}\leq
c_{2}^{\prime}\|v_{h}\|_{L^{2}(\Omega)}\quad\forall v_{h}\in V_{h}$

with uniform constants $c_{1}^{\prime},c_{2}^{\prime}>0$ depending only on the
bounds in assumption (A1).

The result again follows by scaling arguments and the equivalence of norms on
finite-dimensional spaces. Using these observations, we can now prove the
following.

######  Lemma 4.8.

Let (A1)–(A2) and (A5) hold. Then for any $u\in H^{1}(\Omega)^{3}$, we have

$\displaystyle\langle\sigma\widetilde{\Pi}_{h}u,\widetilde{\Pi}_{h}{v}_{h}-{v}_{h}\rangle\leq
C\,h\,\|u\|_{H^{1}(\Omega)}\|v_{h}\|_{L^{2}(\Omega)}\qquad\forall v_{h}\in
V_{h},$ (20)

with a uniform constant $C$ depending only on the bounds in the assumptions.

###### Proof.

We start by considering a special test function $v_{h}=\Phi_{ij}$, where
$\Phi_{ij}$ is one of the basis functions introduced in Lemma 4.6. We then
split

$\displaystyle\langle\sigma\widetilde{\Pi}_{h}u,\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij}\rangle$
$\displaystyle=\langle\sigma\pi_{\omega}^{0}u,\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij}\rangle+\langle\sigma(\widetilde{\Pi}_{h}u-\pi^{0}_{\omega}u),\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij}\rangle$
$\displaystyle=(i)+(ii),$ (21)

where $\pi^{0}_{\omega}:L^{2}(\omega)\to P_{0}(\omega)$ is the
$L^{2}$-projection onto constants on the support $\omega=\omega(e_{ij})$ of
the basis function $\Phi_{ij}$. If
$e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}$, we have
$\Phi_{ij}=\widetilde{\Pi}_{h}\Phi_{ij}$, which means that $(i)=0$ in this
case. If $e_{ij}\in\widetilde{\mathcal{E}}_{h}$, on the other hand, then we
deduce from assumptions (A2) and (A5) that $\sigma$ is constant on the patch
$\omega=\omega(e_{ij})$. We can thus find a vector $b_{1}\in
P_{1}(\omega)^{3}$ such that
$\operatorname{curl}b_{1}=\sigma\pi^{0}_{\omega}u$, and evaluate

$\displaystyle(i)$
$\displaystyle=\langle\operatorname{curl}b_{1},\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij}\rangle_{\omega}$
$\displaystyle=\langle
b_{1},\operatorname{curl}(\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij})\rangle_{\omega}+\langle
b_{1},n\times(\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij})\rangle_{\partial\omega}=(iii)+(iv).$

By elementary computations, see Lemma 5.1 below, one can verify that

$\displaystyle\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij}=\nabla(\lambda_{i}\lambda_{j}),\text{
and hence }\operatorname{curl}(\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij})=0,$
(22)

which in turn implies $(iii)=0$. By definition of the basis functions, the
tangential components of $\Phi_{ij}$ and also that of
$\widetilde{\Pi}_{h}\Phi_{ij}$ vanish on $\partial\omega(e_{ij})$, unless
$e_{ij}$ is a boundary edge, which is excluded by assumption (A5). Hence
$(iv)=0$, and as a consequence, we see that $(i)=0$. Since we assumed $u\in
H^{1}(\Omega)$, we may further bound

$\displaystyle(ii)\leq
C\,h\,\|u\|_{H^{1}(\omega)}\|\Phi_{ij}\|_{L^{2}(\omega)}$

by employing the estimates of Lemma 4.1 and 4.2. In summary, we thus have

$\displaystyle\langle\sigma\widetilde{\Pi}_{h}u,\widetilde{\Pi}_{h}\Phi_{ij}-\Phi_{ij}\rangle\leq
C\,h\,\|u\|_{H^{1}(\omega)}\|\Phi_{ij}\|_{L^{2}(\omega)}.$

The same estimate is obtained for the basis functions $\Phi_{ji}$. By Lemma
4.6, any test function in $V_{h}$ can be expanded as
$v_{h}=\sum_{e_{ij}}\mathsf{v}_{ij}\Phi_{ij}+\mathsf{v}_{ji}\Phi_{ji}$, and by
splitting and summing over all elements, we immediately obtain

$\displaystyle\langle\sigma\widetilde{\Pi}_{h}u,\widetilde{\Pi}_{h}{v}_{h}-{v}_{h}\rangle$
$\displaystyle\leq
Ch\sum\nolimits_{e_{ij}}\|u\|_{H^{1}(\omega(e_{ij}))}\left(\mathsf{v}_{ij}\|\Phi_{ij}\|_{L^{2}(\omega(e_{ij}))}+\mathsf{v}_{ji}\|\Phi_{ji}\|_{L^{2}(\omega(e_{ij}))}\right)$
$\displaystyle\leq
Ch\Big{(}\sum\nolimits_{e_{ij}}\|u\|^{2}_{H^{1}(\omega(e_{ij}))}\Big{)}^{1/2}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{h}\leq
c^{\prime\prime}h\|u\|_{H^{1}(\Omega)}\|v_{h}\|_{L^{2}(\Omega)}.$

Here we used the Cauchy-Schwarz inequality in the second step, the finite
overlap of the patches, and the norm equivalence of Corollary 4.4 in the last.
∎

###  4.4. Discrete stability

We now derive discrete stability estimates for solutions of Method 3.1. To
simplify the presentation, we introduce the short-hand notation

$\displaystyle\|{v}_{h}\|^{2}_{\alpha}\coloneqq\langle\alpha{v}_{h},{v}_{h}\rangle,\qquad\|{v}_{h}\|^{2}_{h,\alpha}\coloneqq\langle\alpha{v}_{h},{v}_{h}\rangle_{h},$

for non-negative piecewise constant parameters $\alpha\in
P_{0}(\mathcal{T}_{h})$. Furthermore, we use

$\displaystyle\widehat{u}_{h}^{\,n+1/2}=\tfrac{1}{2}(u_{h}^{n}+u_{h}^{n+1})\qquad\text{and}\qquad\partial_{\tau}\widehat{u}_{h}^{\,n}=\tfrac{1}{2\tau}(u_{h}^{n+1}-u_{h}^{n-1}).$
(23)

As a final ingredient, we introduce a discrete energy functional, defined by

$\displaystyle\mathscr{E}_{h}(u_{h}^{n},u_{h}^{n+1})$
$\displaystyle\coloneqq\|\partial_{\tau}u_{h}^{n+1/2}\|_{h,\varepsilon}^{2}+\|\operatorname{curl}\widehat{u}_{h}^{\,n+1/2}\|^{2}_{\nu}$
(24)
$\displaystyle\qquad-\tfrac{\tau^{2}}{4}\|\operatorname{curl}\partial_{\tau}u_{h}^{n+1/2}\|^{2}_{\nu}-\tfrac{\tau}{2}\left(\|\widetilde{\Pi}_{h}\partial_{\tau}u_{h}^{n+1/2}\|_{h,\sigma}^{2}-\|\partial_{\tau}u_{h}^{n+1/2}\|_{h,\sigma}^{2}\right).$

We begin with some elementary auxiliary observations.

######  Lemma 4.9.

Let (A1)–(A2) and (A4) hold and $u_{h}^{n},u_{h}^{n+1}\in V_{h}$ be given.
Then

$\displaystyle\tfrac{2}{3}\mathscr{E}_{h}(u_{h}^{n},u_{h}^{n+1})\leq\|\partial_{\tau}u_{h}^{n+1/2}\|_{h,\varepsilon}^{2}+\|\operatorname{curl}\widehat{u}_{h}^{\,n+1/2}\|^{2}_{\nu}\leq
2\mathscr{E}_{h}(u_{h}^{n},u_{h}^{n+1}).$ (25)

###### Proof.

The result follows directly from the definition of the energy functional and
the CFL condition (A4), which was specifically tailored to obtain this result.
∎

We can now establish the required stability estimate for solutions of Method
3.1.

######  Lemma 4.10 (Discrete stability).

Let (A1)–(A2) and (A4) hold. Further let $\xi_{h}^{n}$, ${r}_{h}^{n}\in
V_{h}$, $n\geq 0$ be given such that

$\displaystyle\langle(\varepsilon+\tfrac{\tau}{2}\sigma)\partial_{\tau\tau}\xi_{h}^{n},v_{h}\rangle_{h}+\langle\sigma\widetilde{\Pi}_{h}\partial_{\tau}\xi_{h}^{n-1/2},\widetilde{\Pi}_{h}v_{h}\rangle_{h}+\langle\nu\operatorname{curl}\xi_{h}^{n},\operatorname{curl}v_{h}\rangle=\langle
r_{h}^{n},v_{h}\rangle,\qquad$ (26)

for all $v_{h}\in V_{h}$ and $n\geq 0$. Then for all time steps $0\leq n<N$,
there holds

$\displaystyle\mathscr{E}_{h}(\xi_{h}^{n},\xi_{h}^{n+1})\leq\mathscr{E}_{h}(\xi_{h}^{0},\xi_{h}^{1})+2\sum\nolimits_{k=1}^{n}\tau({r}_{h}^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}).$

###### Proof.

We begin by setting
$v_{h}=\partial_{\tau}\widehat{\xi}_{h}^{\,n}:=\frac{1}{2\tau}(\xi_{h}^{n+1}-\xi_{h}^{n-1})$
in the identity (26). For the first term on the left-hand side, this results
in

$\displaystyle(\varepsilon\partial_{\tau\tau}\xi_{h}^{n},\partial_{\tau}\widehat{\xi}_{h}^{\,n})_{h}$
$\displaystyle=\frac{1}{2\tau}\big{(}\|\partial_{\tau}\xi_{h}^{n+1/2}\|_{h,\varepsilon}^{2}-\|\partial_{\tau}\xi_{h}^{n-1/2}\|_{h,\varepsilon}^{2}\big{)}.$
(27)

For the third term in (26), we see in a similar manner that

$\displaystyle(\nu\operatorname{curl}\xi_{h}^{n},\operatorname{curl}\partial_{\tau}\widehat{\xi}_{h}^{\,n})$
$\displaystyle=\tfrac{1}{2\tau}\Big{(}\|\operatorname{curl}\widehat{\xi}_{h}^{\,n+1/2}\|^{2}_{\nu}-\|\operatorname{curl}\widehat{\xi}_{h}^{\,n-1/2}\|^{2}_{\nu}$
(28)
$\displaystyle\qquad-\tfrac{\tau^{2}}{4}\|\operatorname{curl}\partial_{\tau}\xi_{h}^{n+1/2}\|^{2}_{\nu}+\tfrac{\tau^{2}}{4}\|\operatorname{curl}\partial_{\tau}\xi_{h}^{n-1/2}\|^{2}_{\nu}\Big{)}.$

The loss terms involving $\sigma$ have to be treated more carefully. Here we
use that

$\displaystyle\langle\tfrac{\tau}{2}\sigma$
$\displaystyle\partial_{\tau\tau}\xi_{h}^{n},\partial_{\tau}\widehat{\xi}_{h}^{\,n}\rangle_{h}+\langle\sigma\widetilde{\Pi}_{h}\partial_{\tau}\xi_{h}^{n-1/2},\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,n}\rangle_{h}$
(29)
$\displaystyle=\langle\tfrac{\tau}{2}\sigma\partial_{\tau\tau}\xi_{h}^{n},\partial_{\tau}\widehat{\xi}_{h}^{\,n}\rangle_{h}-\langle\tfrac{\tau}{2}\sigma\partial_{\tau\tau}\widetilde{\Pi}_{h}\xi_{h}^{n},\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,n}\rangle_{h}+\langle\sigma\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,n},\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,n}\rangle_{h}$
$\displaystyle=\tfrac{\tau}{2}\tfrac{1}{2\tau}\Big{(}\|\partial_{\tau}\xi_{h}^{n+1/2}\|_{h,\sigma}^{2}-\|\partial_{\tau}\xi_{h}^{n-1/2}\|_{h,\sigma}^{2}$
$\displaystyle\qquad\qquad\qquad\qquad-\|\widetilde{\Pi}_{h}\partial_{\tau}\xi_{h}^{n+1/2}\|_{h,\sigma}^{2}+\|\widetilde{\Pi}_{h}\partial_{\tau}\xi_{h}^{n-1/2}\|_{h,\sigma}^{2}\Big{)}+\|\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,n}\|_{h,\sigma}^{2}.$

By summing up the identities (27)–(29), we can deduce that

$\displaystyle\mathscr{E}_{h}(\xi_{h}^{n},\xi_{h}^{n+1})$
$\displaystyle\leq\mathscr{E}_{h}(\xi_{h}^{n},\xi_{h}^{n+1})+2\tau\|\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,n}\|_{\sigma}^{2}=\mathscr{E}_{h}(\xi_{h}^{n-1},\xi_{h}^{n})+2\tau\langle{r}_{h}^{n},\partial_{\tau}\widehat{\xi}_{h}^{\,n}\rangle.$

The assertion of the theorem then follows by induction over $n$. ∎

###  4.5. Error estimates

Let us start with proving the estimate of Theorem 3.2. Following standard
practice, we split

$\displaystyle
E(t^{n})-E_{h}^{n}=-(\widetilde{\Pi}_{h}E(t^{n})-E(t^{n}))+(\widetilde{\Pi}_{h}E(t^{n})-E_{h}^{n})\eqqcolon-\eta^{n}+\xi_{h}^{n},$
(30)

into an interpolation error $\eta^{n}$ and a discrete error component
$\xi_{h}^{n}$. The use of the particular projection $\widetilde{\Pi}_{h}$ in
this splitting will become important below. With the interpolation error
estimates of Lemma 4.1, we immediately obtain

$\displaystyle\max\limits_{0\leq n<N}\Big{(}\|\partial_{\tau}\eta^{n+1/2}$
$\displaystyle\|_{L^{2}}^{2}+\|\operatorname{curl}\widehat{\eta}^{\,\,n+1/2}\|_{L^{2}}^{2}\Big{)}$
$\displaystyle\leq
Ch\left(\|\partial_{t}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}+\|\operatorname{curl}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}\right).$

For estimating the discrete error $\xi_{h}^{n}$, we first note that
$\xi_{h}^{0}=0$. Moreover, $E_{h}^{1}=0$ by assumption, and by Taylor
expansion, we further see that

$\displaystyle E(\tau)$
$\displaystyle=E(0)+\tau\partial_{t}E(0)+\tfrac{\tau^{2}}{2}\partial_{tt}E(0)+\tfrac{\tau^{3}}{6}\partial_{ttt}E(s_{3}),$

for some $0<s_{3}<\tau$. From assumption (A3) and using equation (1), one can
see that $E(0)=\partial_{t}E(0)=\partial_{tt}E(0)=0$, and hence
$\xi_{h}^{1}=\frac{\tau^{3}}{6}\widetilde{\Pi}_{h}\partial_{ttt}E(s_{3})$.
Alternatively, we could also get
$\xi_{h}^{1}=\frac{\tau^{2}}{2}\widetilde{\Pi}_{h}\partial_{tt}E(s_{2})$ for
some $0<s_{2}<\tau$ by truncating the Taylor series earlier. In summary, this
yields

$\displaystyle\begin{split}\tfrac{2}{3}\mathscr{E}_{h}(\xi_{h}^{0},\xi_{h}^{1})&\leq\|\partial_{\tau}\xi_{h}^{1/2}\|_{L^{2}}^{2}+\|\operatorname{curl}\widehat{\xi}_{h}^{\;1/2}\|_{L^{2}}^{2}=\tfrac{1}{\tau^{2}}\|\xi_{h}^{1}\|_{L^{2}}^{2}+\tfrac{1}{4}\|\operatorname{curl}\xi_{h}^{1}\|_{L^{2}}^{2}\\\
&\leq
C\tau^{4}\Big{(}\|\partial_{ttt}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}^{2}+\|\partial_{tt}E\|_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}^{2}\Big{)}.\end{split}$
(31)

For the last step, we used the formulas for $\xi_{h}^{0}$ and $\xi_{h}^{1}$
derived above, and the stability estimates
$\|\widetilde{\Pi}_{h}v\|_{L^{2}}\leq C\|v\|_{H^{1}}$ and
$\|\operatorname{curl}\widetilde{\Pi}_{h}v\|_{L^{2}}\leq C\|v\|_{H^{1}}$ for
the projection.

As a next step, let us observe that the continuous solution $E(t)$ of (1)–(2)
satisfies

$\displaystyle\langle\varepsilon\partial_{tt}{E}(t^{n}),v_{h}\rangle+\langle\sigma\partial_{t}{E}(t^{n}),v_{h}\rangle$
$\displaystyle+\langle\nu\operatorname{curl}{E}(t^{n}),\operatorname{curl}v_{h}\rangle=\langle
f(t^{n}),v_{h}\rangle+\langle g(t^{n}),v_{h}\rangle_{\partial\Omega}$

for all $v_{h}\in V_{h}$ and $n>0$. By combination with (10), one can then see
that the discrete error $\xi_{h}^{n}=\widetilde{\Pi}_{h}E(t^{n})-E_{h}^{n}$
thus satisfies the error equation (26) with

$\displaystyle\langle{r}_{h}^{n},v_{h}\rangle=\langle{r}_{h,t}^{n},v_{h}\rangle+\langle{r}_{h,s}^{n},v_{h}\rangle+\langle{r}_{h,q}^{n},v_{h}\rangle+\langle{r}_{h,p}^{n},v_{h}\rangle,$

where the four partial residuals are defined by

$\displaystyle\langle{r}_{h,t}^{n},v_{h}\rangle$
$\displaystyle=\langle\varepsilon(\partial_{\tau\tau}E^{n}-\partial_{tt}E^{n}),v_{h}\rangle+\langle\sigma(\partial_{\tau}\widehat{E}^{n}-\partial_{t}E^{n}),v_{h}\rangle,$
$\displaystyle\langle{r}_{h,s}^{n},v_{h}\rangle$
$\displaystyle=\langle\varepsilon\partial_{tt}\eta^{n},v_{h}\rangle+\langle\sigma\partial_{t}\eta^{n},v_{h}\rangle+\langle\nu\operatorname{curl}\eta^{n},\operatorname{curl}v_{h}\rangle,$
$\displaystyle\langle{r}_{h,q}^{n},v_{h}\rangle$
$\displaystyle=\delta_{h}((\varepsilon+\tfrac{\tau}{2}\sigma)\widetilde{\Pi}_{h}\partial_{\tau\tau}E^{n},v_{h})+\delta_{h}(\sigma\widetilde{\Pi}_{h}\partial_{\tau}E^{n-1/2},v_{h}),$
$\displaystyle\langle{r}_{h,p}^{n},v_{h}\rangle$
$\displaystyle=\langle\widetilde{\Pi}_{h}\sigma\partial_{\tau}E^{\,n-1/2},\widetilde{\Pi}_{h}v_{h}-v_{h}\rangle,$

which represent the temporal, spatial, quadrature, and projection errors,
respectively. In the third term, we again used $\delta_{h}(\alpha
u_{h},v_{h})=\langle\alpha u_{h},v_{h}\rangle_{h}-\langle\alpha
u_{h},v_{h}\rangle$ to abbreviate the quadrature error. We now estimate the
four residuals independently.

First residual. By summation over the time steps, we get

$\displaystyle\sum\nolimits_{k}\tau\langle{r}_{h,t}^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle=\sum\nolimits_{k}\tau$
$\displaystyle\langle\varepsilon(\partial_{\tau\tau}E^{k}-\partial_{tt}E^{k}),\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle$
$\displaystyle+\sum\nolimits_{k}\tau\langle\sigma(\partial_{\tau}\widehat{E}^{k}-\partial_{t}E^{k}),\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle=(i)+(ii).$

The first term in this expansion can be estimated by Taylor expansion, giving

$\displaystyle|(i)|$
$\displaystyle\leq\sum\nolimits_{k}c\tau^{2}\|\partial_{tttt}E\|_{L^{1}(t^{k-1},t^{k};L^{2})}\|\partial_{\tau}\widehat{\xi}_{h}^{\,k}\|_{L^{2}}$
$\displaystyle\leq
C\tau^{4}\|\partial_{tttt}E\|_{L^{1}(L^{2})}^{2}+\tfrac{1}{56}\|\partial_{\tau}\xi_{h}\|_{\ell_{\infty}(L^{2})}^{2}$
$\displaystyle\leq
C\tau^{4}\|\partial_{tttt}E\|_{L^{1}(L^{2})}^{2}+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1}).$

With similar arguments, the second term can be bounded by

$\displaystyle|(ii)|$ $\displaystyle\leq
C\tau^{4}\|\partial_{ttt}E\|_{L^{1}(L^{2})}^{2}+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1})$

Second residual. For the spatial errors, we use

$\displaystyle\sum\nolimits_{k}\tau\langle{r}_{h,s}^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle$
$\displaystyle=\sum\nolimits_{k}\tau\langle\varepsilon\partial_{tt}\eta^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle+\sum\nolimits_{k}\tau\langle\sigma\partial_{t}\eta^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle$
$\displaystyle\qquad+\sum\nolimits_{k}\tau\langle\nu\operatorname{curl}\eta^{k},\operatorname{curl}\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle=(iii)+(iv)+(v).$

By the interpolation error estimates of Lemma 4.1 and Young’s inequality, we
obtain

$\displaystyle|(iii)+(iv)|$ $\displaystyle\leq
Ch^{2}(\|\partial_{t}E\|^{2}_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}+\|\partial_{tt}E\|^{2}_{L^{\infty}(H^{1}(\mathcal{T}_{h}))})+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1}).$

The factor $\tfrac{1}{28}$ in the last term was simply chosen small enough for
later on. For the remaining term in the second residual, we use summation by
parts to arrive at

$\displaystyle|(v)|\leq\langle\operatorname{curl}$
$\displaystyle\eta^{n},\operatorname{curl}\widehat{\xi}_{h}^{\;n+1/2}\rangle-\langle\operatorname{curl}\eta^{1},\operatorname{curl}\widehat{\xi}_{h}^{\;1/2}\rangle$
$\displaystyle-\sum\nolimits_{k}\langle\nu\operatorname{curl}\partial_{\tau}\eta^{k-1/2},\operatorname{curl}\widehat{\xi}_{h}^{\;k-1/2}\rangle=(vi)+(vii)+(viii).$

By the interpolation error estimates and Young’s inequality, we obtain

$\displaystyle|(vi)+(vii)|$ $\displaystyle\leq
Ch^{2}\|\operatorname{curl}{E}\|^{2}_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}+\tfrac{1}{56}\|\operatorname{curl}\widehat{\xi}_{h}\|^{2}_{\ell_{\infty}(L^{2})}$
$\displaystyle\leq
Ch^{2}\|\operatorname{curl}{E}\|^{2}_{L^{\infty}(H^{1}(\mathcal{T}_{h}))}+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1}).$

By Taylor expansion and arguments similar to before, we further obtain

$\displaystyle|(viii)|$ $\displaystyle\leq
Ch^{2}\|\operatorname{curl}\partial_{t}E\|_{L^{1}(H^{1}(\mathcal{T}_{h}))}^{2}+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1}).$

The prefactors in these estimates were again simply chosen sufficiently small.

Third residual. Using Lemma 4.5 and Taylor expansion in time, we can estimate
the residuals caused by the quadrature errors via

$\displaystyle\sum\nolimits_{k}\tau\langle{r}_{h,q}^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle$
$\displaystyle\leq
Ch^{2}(\|\partial_{tt}E\|^{2}_{L^{1}(H^{1}(\mathcal{T}_{h}))}+\|\partial_{t}E\|^{2}_{L^{1}(H^{1}(\mathcal{T}_{h}))})+\tfrac{1}{56}\|\partial_{\tau}\xi_{h}\|^{2}_{\ell_{\infty}(L^{2})}$
$\displaystyle\leq
Ch^{2}(\|\partial_{tt}E\|^{2}_{L^{1}(H^{1}(\mathcal{T}_{h}))}+\|\partial_{t}E\|^{2}_{L^{1}(H^{1}(\mathcal{T}_{h}))})+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1}).$

Fourth residual. For the projection errors, we again use Lemma 4.8 and obtain

$\displaystyle\sum\nolimits_{k}\tau\langle{r}_{h,p}^{k},\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle$
$\displaystyle=\sum\nolimits_{k}\tau\langle\sigma\widetilde{\Pi}_{h}\partial_{\tau}E^{\,k-1/2},\widetilde{\Pi}_{h}\partial_{\tau}\widehat{\xi}_{h}^{\,k}-\partial_{\tau}\widehat{\xi}_{h}^{\,k}\rangle$
$\displaystyle\leq
Ch^{2}\|\partial_{t}E\|_{L^{1}(H^{1})}^{2}+\tfrac{1}{56}\|\partial_{\tau}\xi_{h}\|^{2}_{\ell_{\infty}(L^{2})}$
$\displaystyle\leq
Ch^{2}\|\partial_{t}E\|_{L^{1}(H^{1})}^{2}+\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}_{h}(\xi_{h}^{k},\xi_{h}^{k+1}).$

Let us note that the appearance of the projection $\widetilde{\Pi}_{h}$ in the
residual $r_{h,p}^{k}$ was essential here, in order to be able to apply Lemma
4.8.

Conclusion. By summation of all the individual estimates for the respective
residuals and application of Lemma 4.10, we finally obtain

$\displaystyle\mathscr{E}(\xi_{h}^{n},\xi_{h}^{n+1})$
$\displaystyle\leq\mathscr{E}(\xi_{h}^{0},\xi_{h}^{1})+2\left(7\cdot\tfrac{1}{28}\max_{0\leq
k<N}\mathscr{E}(\xi_{h}^{k},\xi_{h}^{k+1})+C(E)h^{2}+C^{\prime}(E)\tau^{2}\right).$

Taking the maximum over all $n$, using (25), and the bound (31) for
$\mathscr{E}(\xi_{h}^{0},\xi_{h}^{1})$, we thus obtain the required estimates
for the discrete error component. By combination with the interpolation error
estimate, we arrive at the estimate of Theorem 3.2. ∎

### Estimate of Theorem 3.3

We now split the corresponding error by

$\displaystyle
E(t^{n})-\widetilde{\Pi}_{h}E_{h}^{n}=-(\widetilde{\Pi}_{h}E(t^{n})-E(t^{n}))+(\widetilde{\Pi}_{h}E(t^{n})-\widetilde{\Pi}_{h}E_{h}^{n})=:-\eta^{n}+\widetilde{\xi}_{h}^{n}$

The first term is the same as before, and for the discrete error, we can use
the following arguments: By the commuting diagram property of the projectors,
we see that
$\operatorname{curl}\widetilde{\xi}_{h}^{n}=\operatorname{curl}\xi_{h}^{n}$,
which allows using the bounds of the previous proof to handle the
$\operatorname{curl}$-terms in the estimate of the discrete error
$\widetilde{\xi}_{h}^{n}$. For the $L^{2}$-terms, we use

$\displaystyle\|\partial_{\tau}\widetilde{\xi}_{h}^{n+1/2}\|_{L^{2}}$
$\displaystyle=\|\widetilde{\Pi}_{h}(\Pi_{h}\partial_{\tau}E^{n+1/2}-\partial_{\tau}E_{h}^{n+1/2})\|_{L^{2}}\leq
C\|\Pi_{h}\partial_{\tau}E^{n+1/2}-\partial_{\tau}E_{h}^{n+1/2}\|_{L^{2}}$
$\displaystyle\leq
C(\|\widetilde{\Pi}\partial_{\tau}E^{n+1/2}-\partial_{\tau}E_{h}^{n+1/2}\|_{L^{2}}+\|\Pi_{h}\partial_{\tau}E^{n+1/2}-\widetilde{\Pi}_{h}\partial_{\tau}E^{n+1/2}\|_{L^{2}}).$

In the first inequality, we made use of the bound
$\|\widetilde{\Pi}_{h}v_{h}\|_{L^{2}}\leq C\|v_{h}\|_{L^{2}}$, which follows
by the usual scaling arguments. The first term in the above estimate only
involves the discrete error
$\partial_{\tau}\xi_{h}^{n+1/2}=\widetilde{\Pi}\partial_{\tau}E^{n+1/2}-\partial_{\tau}E_{h}^{n+1/2}$,
which was already analyzed before, and the second term can be bounded by the
interpolation error estimates for the two projectors. This then already yields
the estimate of Theorem 3.3. ∎

##  5\. Algebraic properties

We now discuss the implementation of the proposed methods and the algebraic
properties stated in Theorems 3.2 and 3.3. Let us recall the basis functions

$\displaystyle\Phi_{ij}=\lambda_{i}\nabla\lambda_{j}\qquad\text{and}\qquad\Phi_{ji}=-\lambda_{j}\nabla\lambda_{i}$

associated with the edges $e_{ij}\in\mathcal{E}_{h}$. By basic computations,
one can see the following.

######  Lemma 5.1.

Let $v_{h}=\sum_{e_{ij}}\mathsf{v}_{ij}\Phi_{ij}+\mathsf{v}_{ji}\Phi_{ji}\in
V_{h}$ be given. Then

$\displaystyle\widetilde{\Pi}_{h}v_{h}=\sum_{e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}}\mathsf{v}_{ij}\Phi_{ij}+\mathsf{v}_{ji}\Phi_{ji}+\sum_{e_{ij}\in\widetilde{\mathcal{E}}_{h}}\tfrac{1}{2}(\mathsf{v}_{ij}+\mathsf{v}_{ji})(\Phi_{ij}+\Phi_{ji}).$
(32)

Hence
$\widetilde{\Pi}_{h}v_{h}=\sum_{e_{ij}\in\mathcal{E}_{h}}\widehat{\mathsf{v}}_{ij}\Phi_{ij}+\widehat{\mathsf{v}}_{ji}\Phi_{ji}$
with coefficients given by
$\widehat{\mathsf{v}}_{ij}\coloneqq\mathsf{v}_{ij}$,
$\widehat{\mathsf{v}}_{ji}\coloneqq\mathsf{v}_{ji}$ for the edges
$e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}$ and
$\widehat{\mathsf{v}}_{ij}=\widehat{\mathsf{v}}_{ji}\coloneqq\frac{1}{2}(\mathsf{v}_{ij}+\mathsf{v}_{ji})$
for the edges $e_{ij}\in\widetilde{\mathcal{E}}_{h}$.

For the following considerations, we always assume that the degrees of freedom
are sorted edge-wise. This allows us to make the following statements.

######  Lemma 5.2.

The equation (10) is equivalent to the algebraic system (11) with

$\displaystyle[\mathsf{M}_{\varepsilon+\tau\sigma/2}]_{ij,kl}$
$\displaystyle=\langle(\varepsilon+\tau\sigma/2)\Phi_{kl},\Phi_{ij}\rangle_{h},\qquad$
$\displaystyle[\mathsf{f}^{n}]_{ij}$ $\displaystyle=\langle
f(t^{n}),\Phi_{ij}\rangle,$ $\displaystyle[\mathsf{K}_{\nu}]_{ij,kl}$
$\displaystyle=\langle\nu\operatorname{curl}\Phi_{kl},\operatorname{curl}\Phi_{ij}\rangle,\qquad$
$\displaystyle[\mathsf{g}^{n}]_{ij}$ $\displaystyle=\langle
g(t^{n}),\Phi_{ij}\rangle_{\partial\Omega}.$

Both matrices are sparse, and $\mathsf{M}_{\varepsilon+\tau\sigma/2}$ is block
diagonal, with one block per vertex of the mesh coupling all degrees of
freedom not vanishing at this vertex. Furthermore

$\displaystyle\widehat{\mathsf{M}}_{\sigma}=\mathsf{Q}^{\top}\mathsf{M}_{\sigma}\mathsf{Q}\qquad\text{with}\qquad[\mathsf{M}_{\sigma}]_{ij,kl}=(\sigma\Phi_{kl},\Phi_{ij})_{h},$

and the projection matrix $\mathsf{Q}$ is block diagonal with $2\times 2$
blocks of the form

$\displaystyle\frac{1}{2}\begin{pmatrix}1&1\\\
1&1\end{pmatrix}\qquad\text{and}\qquad\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$

for the edges $e_{ij}\in\widetilde{\mathcal{E}}_{h}$ and
$e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}$, respectively.

The assertions follow immediately from the properties of the basis functions;
for details, see [9, 31]. This already proves the second claim of Theorem 3.2.
∎

### Algebraic reduction

We next derive the reduced scheme (12). To do so, we start with the following
observation, which again follows from elementary arguments.

######  Lemma 5.3.

The function $\widetilde{\Pi}_{h}v_{h}$ can be expressed equivalently as

$\displaystyle\widetilde{\Pi}_{h}v_{h}=\sum_{e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}}\widetilde{\mathsf{v}}_{ij}\Phi_{ij}+\widetilde{\mathsf{v}}_{ji}\Phi_{ji}+\sum_{e_{ij}\in\widetilde{\mathcal{E}}_{h}}\widetilde{\mathsf{v}}_{ij}\widetilde{\Phi}_{ij},\qquad\widetilde{\Phi}_{ij}=\Phi_{ij}+\Phi_{ji},$

with $\widetilde{\mathsf{v}}_{ij}=\mathsf{v}_{ij}$,
$\widetilde{\mathsf{v}}_{ji}=\mathsf{v}_{ji}$ for
$e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}$ and
$\widetilde{\mathsf{v}}_{ij}=\frac{1}{2}(\mathsf{v}_{ij}+\mathsf{v}_{ji})$ for
$e_{ij}\in\widetilde{\mathcal{E}}_{h}$, and the collection
$\\{\Phi_{ij},\Phi_{ji}:e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}\\}\cup\\{\widetilde{\Phi}_{ij}:=\Phi_{ij}+\Phi_{ji}):e_{ij}\in\widetilde{\mathcal{E}}_{h}\\}$
is a basis for $\widetilde{V}_{h}$.

Also for the reduced space, the degrees of freedom are sorted edgewise. Then
the relation between the coefficients $\mathsf{v}_{ij}$,
$\widehat{\mathsf{v}}_{ij}$ and $\widetilde{\mathsf{v}}_{ij}$ can be expressed
as follows.

######  Lemma 5.4.

Let $\mathsf{P}$ be the block-diagonal (prolongation) matrix with blocks

$\displaystyle\begin{pmatrix}1\\\
1\end{pmatrix}\qquad\text{and}\qquad\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$
(33)

for edges $e_{ij}\in\widetilde{\mathcal{E}}_{h}$ and
$e_{ij}\in\mathcal{E}_{h}\setminus\widetilde{\mathcal{E}}_{h}$, respectively.
Then
$\mathsf{Q}=\mathsf{P}(\mathsf{P}^{\top}\mathsf{P})^{-1}\mathsf{P}^{\top}$ and
$\mathsf{P}^{\top}\mathsf{P}$ is diagonal with entries $1/2$ and $1$,
respectively. Furthermore,

$\displaystyle\widehat{\mathsf{v}}=\mathsf{P}\,\widetilde{\mathsf{v}}\qquad\text{and}\qquad\widetilde{\mathsf{v}}=\mathsf{R}\,\mathsf{v}\quad\text{with}\quad\mathsf{R}=(\mathsf{P}^{\top}\mathsf{P})^{-1}\mathsf{P}^{\top}.$
(34)

This lemma allows us to express the coefficients of $\widetilde{\Pi}_{h}v_{h}$
in the basis of $\widetilde{V}_{h}$ by the coefficients of the expansion in
the basis of $V_{h}$. With the help of this result, we can now derive the
algebraic form (12) of the reduced scheme of Theorem 3.3.

######  Lemma 5.5.

Let $(\mathsf{E}^{n})_{n}$ be a solution of (11). Then
$\widetilde{\mathsf{E}}^{n}=\mathsf{R}\mathsf{E}^{n}$ satisfies (12) with

$\displaystyle\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}$
$\displaystyle=\mathsf{R}\mathsf{M}_{\varepsilon+\tau\sigma/2}^{-1}\mathsf{R}^{\top},\qquad\widetilde{\mathsf{M}}_{\sigma}=\mathsf{P}^{\top}\mathsf{M}_{\sigma}\mathsf{P},\qquad\widetilde{\mathsf{K}}_{\nu}=\mathsf{P}^{\top}\mathsf{K}_{\nu}\mathsf{P}$

and right-hand sides

$\displaystyle\widetilde{\mathsf{F}}^{n}=\mathsf{R}\mathsf{M}_{\varepsilon+\tau\sigma/2}^{-1}\mathsf{f}^{n}\qquad\widetilde{\mathsf{G}}^{n}=\mathsf{R}\mathsf{M}_{\varepsilon+\tau\sigma/2}^{-1}\mathsf{g}^{n},$
(35)

where $\mathsf{f}^{n},\mathsf{g}^{n}$ are defined as in Lemma 5.2. Hence, the
projection $\widetilde{E}_{h}^{n}=\widetilde{\Pi}_{h}E_{h}^{n}$ of the
solution of Method 3.1 is given by (12).

###### Proof.

Since $\mathsf{M}_{\varepsilon+\tau\sigma/2}$ is regular, the scheme (11) is
equivalent to

$\displaystyle\partial_{\tau\tau}\mathsf{E}^{n}$
$\displaystyle=\mathsf{M}_{\varepsilon+\tau\sigma/2}^{-1}\left(-\widehat{\mathsf{M}}_{\sigma}\partial_{\tau}\mathsf{E}^{n-1/2}-\mathsf{K}_{\nu}\mathsf{E}^{n}+\mathsf{f}^{n}+\mathsf{g}^{n}\right).$

We further multiply this equation from the left by
$\mathsf{R}^{\top}=\mathsf{P}(\mathsf{P}^{\top}\mathsf{P})^{-1}$ and note that
$\mathsf{P}\mathsf{R}=\mathsf{P}(\mathsf{P}^{\top}\mathsf{P})^{-1}\mathsf{P}^{\top}=\mathsf{Q}$,
which follows from the definition of the matrices. We see

$\displaystyle\widehat{\mathsf{M}}_{\sigma}=\mathsf{R}^{\top}\mathsf{P}^{\top}\mathsf{M}_{\sigma}\mathsf{P}\mathsf{R}\qquad\text{and}\qquad\mathsf{K}_{\nu}=\mathsf{Q}^{\top}\mathsf{K}_{\nu}\mathsf{Q}=\mathsf{R}^{\top}\widetilde{\mathsf{K}}_{\nu}\mathsf{R},$

where we used the algebraic form of the commuting diagram property of the
projection $\widetilde{\Pi}_{h}$, see (22). Using the definitions of
$\widetilde{\mathsf{E}}^{n}$, $\widetilde{\mathsf{F}}^{n}$, and
$\widetilde{\mathsf{G}}^{n}$, we thus conclude that

$\displaystyle\partial_{\tau\tau}\widetilde{\mathsf{E}}^{n}=\partial_{\tau\tau}\mathsf{R}\mathsf{E}^{n}$
$\displaystyle=\mathsf{R}\mathsf{M}_{\varepsilon+\tau\sigma/2}^{-1}(-\mathsf{R}^{\top}\widetilde{\mathsf{M}}_{\sigma}\mathsf{R}\,\partial_{\tau}\mathsf{E}^{n-1/2}-\mathsf{R}^{\top}\widetilde{\mathsf{K}}_{\nu}\mathsf{R}\mathsf{E}^{n}+\mathsf{f}^{n}+\mathsf{g}^{n})$
$\displaystyle=\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}(-\widetilde{\mathsf{M}}_{\sigma}\partial_{\tau}\widetilde{\mathsf{E}}^{n-1/2}-\widetilde{\mathsf{K}}_{\nu}\widetilde{\mathsf{E}}^{n})+\widetilde{\mathsf{F}}^{n}+\widetilde{\mathsf{G}}^{n}$

This already yields the algebraic form given in (12) and concludes the proof.
∎

######  Remark 5.6.

The choice (35) for the right-hand sides $\widetilde{\mathsf{F}}^{n}$,
$\widetilde{\mathsf{G}}^{n}$ in the reduced method makes use of the vectors
$\mathsf{f}^{n}$, $\mathsf{g}^{n}$ with two degrees of freedom for each edge.
An alternative and more direct choice would be

$\displaystyle\widetilde{\mathsf{F}}^{n}:=\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}\widetilde{\mathsf{f}}^{n}\qquad\text{and}\qquad\widetilde{\mathsf{G}}^{n}:=\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}\widetilde{\mathsf{g}}^{n},$
(36)

which now only depends on the vectors $\widetilde{\mathsf{f}}^{n}$ and
$\widetilde{\mathsf{g}}^{n}$ assembled on the reduced finite element space
$\widetilde{V}_{h}$. This modification can be included in our analysis, if one
replaces assumption (A5) by choosing
$\widetilde{\mathcal{E}}_{h}\subset\mathcal{E}_{h}$ such that

  1. (A5∗)

$\sigma$ and $f$ are continuous across edges $e\in\widetilde{\mathcal{E}}_{h}$
lying in the interior $\Omega$, and $\sigma=0$, $g=0$ for all edges
$e\in\widetilde{\mathcal{E}}_{h}$ on the boundary $\partial\Omega$.

This allows us to apply a variant of Lemma 4.8 which covers the additional
consistency error introduced by the representation of $f$; the boundary term
$g$ is already fully covered by the original assumption (A5) since it vanishes
along all edges where no reduction is applied. The time-stepping scheme
resulting from this modified choice of the right-hand sides can be stated as

$\displaystyle\partial_{\tau\tau}\widetilde{\mathsf{E}}^{n}=\widetilde{\mathsf{M}}_{\varepsilon+\tau\sigma/2}^{-1}(-\widetilde{\mathsf{M}}_{\sigma}\partial_{\tau}\widetilde{\mathsf{E}}^{n-1/2}-\widetilde{\mathsf{K}}_{\nu}\widetilde{\mathsf{E}}^{n}+\widetilde{\mathsf{f}}^{n}+\widetilde{\mathsf{g}}^{n}),$
(37)

and used as an alternative to (12). In our numerical tests, we will compare
the two versions (12) and (37), in particular, also highlighting the
importance of the additional conditions in assumption (A5∗) for the latter.

##  6\. Numerical validation

As a model problem for our numerical tests, we consider the scattering of a
plane electromagnetic wave from a cylinder. Under the usual symmetry
assumptions, this can be modeled by Maxwell’s equations in two space
dimensions. All results from the three-dimensional setting translate almost
verbatim. The main differences are that $\mathcal{T}_{h}$ now corresponds to a
triangular mesh and that the space $\widetilde{V}_{h}$ is a subspace of
$P_{1}(K)^{2}\cap H(\operatorname{curl};\Omega)$. Moreover, the triangular
vertex rule

$\displaystyle(a,b)_{h}=\sum\nolimits_{T}\tfrac{|T|}{3}\sum\nolimits_{v_{i}\in
T}a(v_{i})\cdot b(v_{i})$ (38)

is used for the numerical quadrature.

### Model problem.

The geometric setup of our test problem is illustrated in Figure 1. The
computational domain $\Omega=(-1,1)^{2}$ is split into two parts: The circular
region $\Omega_{S}$ of radius $r_{S}=0.3$ contains the scatterer with material
parameters $\varepsilon=\nu=1$ and $\sigma=100$. For the surrounding medium
$\Omega\setminus\Omega_{S}$ we take $\varepsilon=\nu=1$ and $\sigma=0$.

$\Omega_{S}$ $\Omega_{R\setminus S}$ $(-1,-1)$$(1,-1)$$(-1,1)$$(1,1)$ Figure
1. Computational domain for the wave scattering problem.

### Test scenario

For excitation of the fields, the initial and boundary conditions are chosen
to match the corresponding traces of the plane wave

$\displaystyle\widetilde{E}(x,y,t)=\begin{pmatrix}-k_{2}\\\
k_{1}\end{pmatrix}a(k_{1}x+k_{2}x-t)$

with envelope $a(x)=2e^{-10(x+3)^{2}}$ and
$k=(k_{1},k_{2})=\frac{1}{\sqrt{2}}(1,1)$. The function $\widetilde{E}$
satisfies (1)–(2) with $\varepsilon=\nu=1$, $\sigma=0$, and data $f=0$ and
$g=n\times(\nu\operatorname{curl}\widetilde{E})$. As a consequence, the plane
wave will first propagate freely through $\Omega\setminus\Omega_{S}$, but then
get scattered at the inclusion $\Omega_{S}$ with high conductivity. The
scattered wave will then travel back through the free region
$\Omega\setminus\Omega_{S}$ and finally be artificially reflected at the outer
boundary $\partial\Omega$. In Figure 2, some corresponding snapshots of the
numerical solution $E_{h}(t)$ are depicted.

Figure 2. Absolute values of the simulated electric field $E_{h}(t)$ at time
$t=1.5$ (left), $t=2$ (middle), and $t=2.5$ (right).

### Discretization setup

To examine the convergence of the methods, we consider a sequence of meshes
$\mathcal{T}_{h}$ obtained by uniform refinement of an initial unstructured
coarse mesh. The time step $\tau$ is chosen following the CFL condition (A4).
Let us note that $\sigma$ is continuous everywhere but on the edges along
$\partial\Omega_{S}$ and vanishes on the outer boundary $\partial\Omega$. To
illustrate our theoretical results, we consider Method 3.1 for three different
choices of $\widetilde{\mathcal{E}}_{h}$, namely

  * •

$\widetilde{\mathcal{E}}_{h}=\emptyset$: This leads to
$\widetilde{V}_{h}=\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$ with two degrees
of freedom for every edge;

  * •

$\widetilde{\mathcal{E}}_{h}\subseteq\mathcal{E}_{h}$ such that (A5) is
satisfied, i.e.
$\mathcal{N}_{0}\subset\widetilde{V}_{h}\subset\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$:
in particular, only one degree of freedom is used for all edges not lying on
the boundary $\partial\Omega_{S}$;

  * •

$\widetilde{\mathcal{E}}_{h}=\mathcal{E}_{h}$: This amounts to
$\widetilde{V}_{h}=\mathcal{N}_{0}$ and violates assumption (A5).

The corresponding finite element solutions are denoted by
$E_{h}^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}$, $E_{h}^{\mathcal{N}^{+}}$ and
$E_{h}^{\mathcal{N}}$. For measuring the error between two functions $E_{h}$
and $E_{h}^{*}$, we use

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|}E_{h}-E_{h}^{*}{|\kern-1.07639pt|\kern-1.07639pt|}\coloneqq\frac{\|\partial_{\tau}(E_{h}-E_{h}^{*})\|_{\ell_{\infty}(L^{2}(\Omega))}}{\|\partial_{\tau}E_{h}^{*}\|_{\ell_{\infty}(L^{2}(\Omega))}}$
$\displaystyle+\frac{\|\operatorname{curl}(\widehat{E}_{h}-\widehat{E}_{h}^{*})\|_{\ell_{\infty}(L^{2}(\Omega))}}{\|\operatorname{curl}\widehat{E}_{h}^{*}\|_{\ell_{\infty}(L^{2}(\Omega))}}$

### Computational results

We first study the error and convergence rate of the
$\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$ method. Since the exact solution to
our model problem is somewhat cumbersome to compute, we estimate the
discretization errors by comparing the numerical solutions obtained on two
different nested meshes.

$h\approx$ | ${|\kern-1.07639pt|\kern-1.07639pt|}E^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}_{h}-E^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}_{2h}{|\kern-1.07639pt|\kern-1.07639pt|}$ | eoc | dofs
---|---|---|---
$2^{-3}$ | $0.779603$ | — | 2.296
$2^{-4}$ | $0.373586$ | $1.06$ | 9.056
$2^{-5}$ | $0.185312$ | $1.01$ | 35.968
$2^{-6}$ | $0.093906$ | $0.98$ | 143.360
$2^{-7}$ | $0.046257$ | $1.02$ | 572.416
$2^{-8}$ | $0.023158$ | $1.00$ | 2.287.616
$2^{-9}$ | $0.011612$ | $1.00$ | 9.146.368
Table 1. Errors, the estimated order of convergence (eoc) and the number of
degrees of freedom (dofs) for a multitude of mesh sizes $h$ with fixed time
step sizes $\tau=0.28\,h$.

As one can infer from Table 1, first-order convergence of the
$\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$-method is observed, which is in
perfect agreement with our previous studies [9, 31].

In the following tests, we evaluate the convergence rates for the
$\mathcal{N}_{0}^{+}$ and $\mathcal{N}_{0}$ methods by comparing them to the
simulations obtained with the $\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$
method.

$h\approx$ | ${|\kern-1.07639pt|\kern-1.07639pt|}E_{h}^{\mathcal{N}^{+}}-E_{h}^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}{|\kern-1.07639pt|\kern-1.07639pt|}$ | eoc | dofs | ${|\kern-1.07639pt|\kern-1.07639pt|}E_{h}^{\mathcal{N}}-E_{h}^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}{|\kern-1.07639pt|\kern-1.07639pt|}$ | eoc | dofs
---|---|---|---|---|---|---
$2^{-3}$ | $0.098385$ | — | 1.164 | $0.105887$ | — | 1.148
$2^{-4}$ | $0.049076$ | $1.00$ | 4.560 | $0.053302$ | $1.00$ | 4.528
$2^{-5}$ | $0.024133$ | $1.02$ | 18.048 | $0.027561$ | $0.95$ | 17.984
$2^{-6}$ | $0.011923$ | $1.01$ | 71.808 | $0.015084$ | $0.87$ | 71.680
$2^{-7}$ | $0.005916$ | $1.01$ | 286.464 | $0.008739$ | $0.79$ | 286.208
$2^{-8}$ | $0.002945$ | $1.00$ | 1.144.320 | $0.005371$ | $0.70$ | 1.143.808
$2^{-9}$ | $0.001469$ | $1.00$ | 4.574.208 | $0.003474$ | $0.63$ | 4.573.184
Table 2. Errors, the estimated order of convergence (eoc) and the number of
degrees of freedom (dofs) for a multitude of mesh sizes $h$ with fixed time
step sizes $\tau=0.28\,h$.

For the $\mathcal{N}_{0}^{+}$-method, assumption (A5) is valid, and we observe
first-order convergence as predicted by our theoretical results. Recall that
this method has two degrees of freedom (only) for all edges at the interface
$\partial\Omega_{S}$ where $\sigma$ has a jump. The $\mathcal{N}_{0}$-method,
on the other hand, uses algebraic reduction also at the interface, thus
violating assumption (A5). This indeed results in a reduction of the
convergence rate, as could be expected from our analysis. Condition (A5)
therefore seems necessary to obtain the optimal convergence rates in the
presence of discontinuous conductivities. In further tests, we also verified
that assumption (A5∗) is necessary if the right-hand sides are chosen
according to (36).

### Additional consistency errors

To better understand the convergence breakdown when violating the consistency
conditions (A5) and (A5∗), we illustrate in Figure 3 the local contributions
of the errors
${|\kern-1.07639pt|\kern-1.07639pt|}E_{h}^{\mathcal{N}}-E_{h}^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}{|\kern-1.07639pt|\kern-1.07639pt|}$
for $h=2^{-4}$ for the two choices (35) and (36) of the right-hand sides.

Figure 3. Local contributions of the error
${|\kern-0.96873pt|\kern-0.96873pt|}E_{h}^{\mathcal{N}}-E_{h}^{\mathcal{N\mkern-2.0muC\mkern-2.0mu}}{|\kern-0.96873pt|\kern-0.96873pt|}$
for $h=2^{-4}$ with different implementations of the right-hand sides. The red
color indicates which part of the mesh the error dominates.

As can be expected, the dominating error contributions come from the interface
and the boundary, where the consistency conditions (A5) respectively (A5∗) are
violated.

### Inspection of the CFL condition

As we have observed in assumption (A4), the CFL condition is somewhat non-
standard, introducing an additional term depending on the projection
$\widetilde{\Pi}_{h}$. To illustrate the effect of this additional term, we
compute the largest time step $\tau=\tau_{\max}$ for which

$\displaystyle\frac{\tau_{\max}^{2}}{4}\mathsf{v}^{\top}\mathsf{K}_{\nu}\mathsf{v}-\frac{1}{2}\mathsf{v}^{\top}\mathsf{M}_{\varepsilon}\mathsf{v}+\gamma\frac{\tau_{\max}}{2}\Big{|}\mathsf{v}^{\top}\Big{(}\widehat{\mathsf{M}}_{\sigma}-\mathsf{M}_{\sigma}\Big{)}\mathsf{v}\Big{|}\geq
0$

for all vectors $\mathsf{v}$, where $\gamma=0$ for the
$\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$ method and $\gamma=1$ for the
$\mathcal{N}_{0}^{+}$ method. By Taylor expansion, we note that $\tau_{\max}$
behaves approximately like

$\displaystyle\tau_{\max}\sim\sqrt{\frac{\varepsilon}{\nu}}h+\gamma\frac{\sigma}{2\nu}h^{2}.$
(39)

For $\gamma=0$ or $\sigma=0$, we obtain the classic linear dependence
$\tau\leq Ch$. In Table 3, we display the values of $C=\tau_{\max}/h$ for each
of the methods for the model parameters used in our simulations.

$h$ | $\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$ | $\mathcal{N}_{0}^{+}$
---|---|---
$2^{-3}$ | $0.391076$ | $0.210000$
$2^{-4}$ | $0.381624$ | $0.317041$
$2^{-5}$ | $0.377977$ | $0.377977$
$2^{-6}$ | $0.376846$ | $0.376846$
$2^{-7}$ | $0.376528$ | $0.376528$
Table 3. CFL constants for both the $\mathcal{N\mkern-2.0muC\mkern-2.0mu}_{1}$
and the $\mathcal{N}_{0}^{+}$ method on a sequence of uniformly refined
meshes.

First note that the CFL constant $C$ is uniformly bounded from above for both
methods. While for larger values of $h$, the CFL constant is somewhat more
stringent for the $\mathcal{N}_{0}^{+}$ method, the constant $C$ behaves
almost the same for small $h$, which is in perfect agreement with (39). The
algebraic reduction hence does not have a severe effect on the maximal
admissible time step.

##  7\. Discussion

In this paper, we proposed and analyzed finite-element schemes which can be
seen as a natural extension of the Yee scheme to unstructured grids and
inhomogeneous lossy media. While the number of degrees of freedom can be
reduced to one for almost all edges, the necessity incorporation of two
degrees of freedom at interfaces or boundaries with material jumps has been
illustrated. A full convergence analysis of the schemes could be provided and
optimal convergence could be proven under reasonable assumptions. Alternative
extensions of the Yee scheme to triangular and tetrahedral elements based on
dual cell complexes have been proposed previously in [4, 5, 19]. While in some
settings, the algebraic form of these methods is similar or even identical to
ours, the constructions in these papers follow a rather geometric approach
which complicates the error analysis. The results obtained in our paper may be
useful to gain further insight also into these related methods.

## References

  * [1] D. Boffi, F. Brezzi, and M. Fortin.  Mixed finite element methods and applications, volume 44 of Springer Series in Computational Mathematics.  Springer, Heidelberg, 2013.
  * [2] A. Bossavit and L. Kettunen.  Yee-like schemes on a tetrahedral mesh, with diagonal lumping.  Int. J. Numer. Model., 12:129–142, 1999.
  * [3] S. C. Brenner and L. R. Scott.  The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics.  Springer, New York, third edition, 2008.
  * [4] L. Codecasa, B. Kapidani, R. Specogna, and F. Trevisan.  Novel FDTD technique over tetrahedral grids for conductive media.  IEEE Transactions on Antennas and Propagation, 66:5387–5396, 2018\.
  * [5] L. Codecasa and M. Politi.  Explicit, consistent, and conditionally stable extension of FDTD to tetrahedral grids by FIT.  IEEE Transactions on Magnetics, 44:1258–1261, 2008.
  * [6] G. Cohen.  Higher-Order Numerical Methods for Transient Wave Equations.  Springer, Heidelberg, 2002.
  * [7] G. Cohen and P. Monk.  Efficient edge finite element schemes in computational electromagnetism.  In Mathematical and numerical aspects of wave propagation, pages 250–259. SIAM, Philadelphia, PA, 1995.
  * [8] G. Cohen and S. Pernet.  Finite Element and Discontinuous Galerkin Methods for Transient Wave Equations.  Scientific Computation. Springer Netherlands, 2016.
  * [9] H. Egger and B. Radu.  A mass-lumped mixed finite element method for Maxwell’s equations.  In Scientific computing in electrical engineering, volume 32 of Math. Ind., pages 15–24. Springer, Cham, 2020.
  * [10] H. Egger and B. Radu.  A second-order finite element method with mass lumping for Maxwell’s equations on tetrahedra.  SIAM J. Numer. Anal., 59:864–885, 2021.
  * [11] A. Elmkies and P. Joly.  Éléments finis d’arête et condensation de masse pour les équations de Maxwell: le cas 2D.  Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, 324:1287–1293, 1997.
  * [12] A. Elmkies and P. Joly.  Éléments finis d’arête et condensation de masse pour les équations de Maxwell: le cas de dimension 3.  Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, 325:1217–1222, 1997.
  * [13] A. Ern and J.-L. Guermond.  Theory and practice of finite elements, volume 159 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2004.
  * [14] R. Farwig and V. Rosteck.  Note on Friedrichs’ inequality in $N$-star-shaped domains.  J. Math. Anal. Appl., 435:1514–1524, 2016.
  * [15] S. Geevers, W. Mulder, and J. van der Vegt.  New higher-order mass-lumped tetrahedral elements for wave propagation modelling.  SIAM Journal on Scientific Computing, 40:A2830–A2857, 2018.
  * [16] J. S. Hesthaven and T. Warburton.  Nodal discontinuous Galerkin methods, volume 54 of Texts in Applied Mathematics.  Springer, New York, 2008.  Algorithms, analysis, and applications.
  * [17] R. Holland.  Finite-difference solution of Maxwell’s equations in generalized nonorthogonal coordinates.  IEEE Trans. on Nuclear Science, NS-30:4589–4591, 1983.
  * [18] P. Joly.  Variational methods for time-dependent wave propagation problems.  In Topics in Computational Wave Propagation, volume 31 of LNCSE, pages 201–264. Springer, 2003.
  * [19] B. Kapidani, L. Codecasa, and J. Schöberl.  An arbitrary-order Cell Method with block-diagonal mass-matrices for the time-dependent 2D Maxwell equations.  J. Comput. Phys., 433:Paper No. 110184, 20, 2021.
  * [20] J. F. Lee, R. Pandalech, and R. Mittra.  Modeling three-dimensional discontinuities in waveguides using nonorthogonal FDTD algorithm.  IEEE Trans. on Microwave Theory and Techniques, 40:346–352, 1992\.
  * [21] R. Leis.  Initial Boundary Value Problems in Mathematical Physics.  John Wiley, New York, 1988.
  * [22] C. G. Makridakis and P. Monk.  Time-discrete finite element schemes for Maxwell’s equations.  RAIRO Modél. Math. Anal. Numér., 29:171–197, 1995.
  * [23] P. Monk.  Analysis of a finite element methods for Maxwell’s equations.  SIAM J. Numer. Anal., 29:714–729, 1992.
  * [24] P. Monk.  A comparison of three mixed methods for the time-dependent Maxwell’s equations.  SIAM J. Sci. Statist. Comput., 13:1097–1122, 1992.
  * [25] P. Monk.  An analysis of Nédélec’s method for the spatial discretization of Maxwell’s equations.  J. Comput. Appl. Math., 47:101–121, 1993.
  * [26] P. Monk.  Finite element time domain methods for Maxwell’s equations.  In Second International Conference on Mathematical and Numerical Aspects of Wave Propagation (Newark, DE, 1993), pages 380–389. SIAM, Philadelphia, PA, 1993.
  * [27] P. Monk.  Finite element methods for Maxwell’s equations.  Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2003.
  * [28] P. B. Monk and E. Süli.  A convergence analysis of Yee’s scheme on nonuniform grids.  SIAM J. Numer. Anal., 31:393–412, 1994.
  * [29] J. C. Nédélec.  A new family of mixed finite elements in $\mathbb{R}^{3}$.  Numer. Math., 50:57–81, Jan 1986.
  * [30] J. C. Nédélec.  Mixed finite elements in $\mathbb{R}^{3}$.  Numer. Math., 35:315–341, 1980.
  * [31] B. Radu.  Finite element mass lumping for H(div) and H(curl).  PhD thesis, Technische Universität Darmstadt, Darmstadt, 2022.
  * [32] R. Schuhmann and T. Weiland.  FDTD on nonorthogonal grids with triangular fillings.  IEEE Trans. Magn., pages 1470–1473, 1998.
  * [33] R. Schuhmann and T. Weiland.  A stable interpolation technique for FDTD on nonorthogonal grids.  Int. J. Numer. Model., 11:299–306, 1998.
  * [34] A. Taflove.  Application of the finite-difference time-domain method to sinusoidal steady-state electromagnetic-penetration problems.  IEEE Transactions on Electromagnetic Compatibility, EMC-22:191–202, Aug 1980.
  * [35] A. Taflove and S. C. Hagness.  Computational electrodynamics: the finite-difference time-domain method.  Artech House, Norwood, 3rd edition, 2005.
  * [36] T. Weiland.  A discretization model for the solution of Maxwell’s equations for six-component fields.  Archiv Elektronik und Uebertragungstechnik, 31:116–120, Mar. 1977\.
  * [37] T. Weiland.  Finite integration method and discrete electromagnetism.  In P. Monk, C. Carstensen, S. Funken, W. Hackbusch, and R. H. W. Hoppe, editors, Computational Electromagnetics, pages 183–198, Berlin, Heidelberg, 2003. Springer Berlin Heidelberg.
  * [38] K. Yee.  Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media.  IEEE Transactions on Antennas and Propagation, 14:302–307, May 1966\.