Received March 31, 2019, accepted May 6, 2019, date of publication May 13, 2019, date of current version

May 29, 2019.

Digital Object Identifier 10.1109/ACCESS.2019.2916394

Time Domain Finite Element Method

for Maxwell’s Equations
ASAD ANEES AND LUTZ ANGERMANN
Institute of Mathematics, Clausthal University of Technology, 38678 Clausthal-Zellerfeld, Germany
Corresponding author: Asad Anees (asad.anees@tu-clausthal.de)
This work was partially supported by the German Academic Exchange Service (DAAD) through a Funding Program under
Grant 57129429, and by the Open Access Publishing Fund of the Clausthal University of Technology.

ABSTRACT In this paper, we discuss a time domain finite element method for the approximate solution of
Maxwell’s equations. A weak formulation is derived for the electric and magnetic fields with appropriate
initial and boundary conditions, and the problem is discretized both in space and time. In space, Nédeléc
curl-conforming and Raviart–Thomas div-conforming finite elements are used to discretize the electric and
magnetic fields, respectively. The backward Euler and symplectic schemes are applied to discretize the
problem in time. For this system, we prove an error estimate. In addition, computational experiments are
presented to validate the method, the electric and magnetic fields are visualized. The method also allows
treating complex geometries of various physical systems coupled to electromagnetic fields in 3D.

INDEX TERMS Backward Euler method, error estimates, Maxwell’s equations, time domain finite element
methods, simulation, symplectic method, visualization.

I. INTRODUCTION magnetic field formulations, respectively. In the literature,

The solution of Maxwell’s equations in the time domain a lot of papers have been devoted to the solution of sec-
formulation is involved in many engineering and industrial ond order wave equations [1]–[15]. In [1], Monk presented
problems, e.g. RF, radar, mixed signal integrated circuits semi-discrete error estimates in the energy and L2 norms
(ICs), diffraction of electromagnetic waves, plasma physics, by employing Nédeléc curl-conforming elements for the
acoustic or seismic wave propagation, radiation, scattering, vector wave equations. The paper [2] presents TDFEM for
environmental and medical imaging, and microwave devices. the second order vector wave equations and hyperbolic sys-
In the presence of complex media or geometries, finite ele- tems using node-based and edge-based elements. Various
ment methods in either continuous or discontinuous vari- numerical experiments were performed to investigate the
ants are the main numerical approaches. In the literature advantages or disadvantages of the mass lumping scheme.
an abundance of work about convergence analysis, semi- Comparisons of various TDFEM and semi-discrete error
discrete, fully-discrete error estimates, and numerical simula- estimates are presented in [3] for the electric field hyper-
tions for the time dependent Maxwell’s equations exists. The bolic equation in anisotropic and inhomogeneous media with
Galerkin time domain finite element methods (TDFEM) can respect to test and trial spaces, explicit and implicit formu-
be grouped into two classes, one class of schemes directly lations. Moreover, the convergence of a fully discrete finite
deals with the system of Maxwell’s equations, whereas the element scheme is analyzed for the second order electric
other class solves the second order wave equations (classical field equation (vector wave equation) in [4], and optimal
approach). error estimates are obtained in the L2 norm. Furthermore,
In the classical approach, second order wave equa- the only simulations of the second order vector wave equa-
tions or hyperbolic system are obtained either by eliminating tion are described in [5], [12], [14], [15]. In the papers [4]
the electric or magnetic field from the system of Maxwell’s and [5], the vector wave equation is discretized by Nédeléc
equations. The resulting problems are called electric or curl-conforming finite elements in space, and second order
backward and central finite differences, respectively are used
The associate editor coordinating the review of this manuscript and to discretize in time. The scheme described by White and
approving it for publication was Gang Mei. Stowell [5] is second order accurate in space and time.

2169-3536
2019 IEEE. Translations and content mining are permitted for academic research only.
63852 Personal use is also permitted, but republication/redistribution requires IEEE permission. VOLUME 7, 2019
See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

Edge finite elements (for the magnetic vector potential) in discretization. The L2 error is estimated for a semi-discrete
the time domain are presented in [12] to address the prob- scheme in [26]. We proposed a splitting approach for the
lem of inductive and capacitive effects. A TDFEM (curl-curl Maxwell’s equation that splits the system of Maxwell’s equa-
electric field equation) forward solver transmitting loop is tions into two uncoupled system, to deal with ε and µ as
presented in [13] to address for a complex shaped domain. matrix function of space (complex material) [36]. The split-
A local time stepping method (LTS) based on explicit Runge- ting method allows to solve the uncoupled systems inde-
Kutta schemes having arbitrary accuracy in time for wave pendently, and is proved to be stable and convergent at the
propagation is demonstrated in [14]. semi-discrete level. The operator form of uncoupled systems
The vector wave equation and magnetic vector potential and semi-discrete error estimates are presented in [32], but
approaches allow to address the Maxwell’s equations in time fully-discrete error estimation and computational experiment
domain, an easy implementation for analysis, error estimation are not given in [26], [32], [36]. In [18], semi-discrete and
and simulations. However, the simulation, analysis and error fully-discrete error estimates are obtained from the operator
estimation of the wave equation (electric field formulation form of the system of Maxwell’s equations, and the time
and magnetic field formulation) cause spurious and non- discretizations by rational approximations of the exponential
physical solutions that are linearly raising corresponding to function are investigated, but no computational experiments
time [6], [9]. were given. Moreover, the only simulations for the system of
In [16]–[32], a number of mixed time domain finite ele- Maxwell’s equations are given in [21], [35]
ment methods are explained to deal with the system of In the article [19], the Maxwell’s equations are discretized
Maxwell’s equations. An abundance of mixed time domain in space by a node-based and edge-based finite element
finite element methods for the direct application to Maxwell’s method, and an efficient solver is described both for the
system is available, where the electric and magnetic fields frequency and time domain formulations. The paper [24]
are discretized in space by discontinuous and Nédeéc curl- describes a general way to investigate the stability of temporal
conforming spaces, respectively, see e.g. [16], [20], [22], discretization schemes such as backward difference, forward
[23], [29] and [31]. In the work [16], error estimates difference and central difference methods in electromagnet-
are demonstrated for a semi-discrete problem, but com- ics. The stability is determined by analyzing the root locus
putational experiments and fully-discrete error estimates map of a characteristic equation and evaluating the spectral
are not provided. Both semi-discrete, and fully-discrete radius of finite element system matrix. Stability properties are
(using a Crank-Nicolson discretization) point-wise super- given in [3] for simulations of transient electromagnetic phe-
convergence results are obtained for Maxwell’s equations in nomena for the Maxwell’s equations. In the article [25], time
metamaterials for nonuniform cublic and rectangular edge domain finite element methods based on Whitney elements
elements in [33]. Lee and Madsen [20] also demonstrated a are proposed for solving transient response problems on tetra-
mixed time domain finite element simulation for Maxwell’s hedral meshes. One of the proposed schemes is uncondition-
equations and employed a explicit leapfrog time integration ally stable, another scheme is explicit but does not require
scheme. A mixed semi-discrete and a fully-discrete error matrix inversions. An explicit time domain finite element
analysis for Maxwell’s equations in double negative mate- algorithm is presented for the Maxwell’s equations in [30],
rial are given in [22], but computational experiments were for complex media in [28] (only numerical result). An energy
not performed. In [34], a variable time step method for conserving method for 3D Maxwell’s equations is obtained
time domain Maxwell’s equations is presented. Fully discrete in [37], based on an exponential operator splitting approach.
error estimates and computational experiments for Maxwell’s Many papers have been written about time domain discon-
equations simultaneously for dispersive and Lorentz meta- tinuous Galerkin (TDDG) methods in computational electro-
materials are presented in [23], [31]. There the temporal magnetics [38]–[49]. Other time domain methods to solve
discretization is done by means of first order backward finite either the Maxwell’s equations or the vector wave equation
differences, and computational experiments are performed can be found in [50]–[62]. Several books for electromagnet-
for 2D situations. Error estimates are also presented in [27], ics [63]–[68] are available for analysis and simulation.
where the problem is discretized in space and time by means The TDFEM proposed in [18], [27], [32], [34], [36] also
of Nédeléc curl-conforming finite elements and backward cause spurious and non-physical solutions because these
finite differences, respectively. The scheme is called decou- methods do not figure out quantities from the system of
pled or explicit magnetic field scheme and causes spurious Maxwell’s equations directly. It is well known that H1 con-
solutions. Furthermore, other semi-discrete theoretical and forming finite elements for electromagnetics may result in
numerical results for the time dependent Maxwell’s equation spurious and non-physical solutions. The degrees of free-
in composite material are presented in [29], where the mate- dom for Nédeléc curl-conforming and Raviart-Thomas div-
rial parameters ε, µ and σ are 3 × 3 positive definite matrices conforming finite elements are related to the edges and faces
depending on the spatial variables. of the meshes, respectively, and not to the mesh nodes. These
In [18], [21], [26], [32], [35], time domain finite ele- finite elements also avoid the appearance of non-physical,
ment methods for Maxwell’s equations are discussed using spurious and divergent solutions [69], [70]. These are good
curl-conforming and div-conforming elements for spatial reasons to use Nédeléc curl-conforming and Raviart-Thomas

VOLUME 7, 2019 63853

 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

div-conforming finite elements. To date most of contributions A perfect conducting boundary condition on  is assumed,
have been based on Nédeléc curl-conforming (edge elements) that is,
for the time domain solution of Maxwell’s equation and a few
are using (Nédeléc and Raviart-Thomas) edge and face ele- n×E=0 on 0 × (0, T ).
ments with spurious solutions. The technique of error estima- Finally, initial conditions have to be specified so that
tion and simulation we present is motivated by the last three
decades works. We demonstrate error estimates for fully dis- E(x, 0) = E0 (x) and H(x, 0) = H0 (x) for all x ∈ , (3)
cretized Maxwell’s equations based on a time domain finite
where E0 and H0 are given functions on 0, and H0 satisfies
element approach, and simulations using various solvers and
visualizations of the computed electric and magnetic fields. ∇ · (µH0 ) = 0 in , H0 · n = 0 on 0.
In our approach, we deal with the system of Maxwell’s
equations rather than the second order vector wave equation. We assume that the solution (E, H) of the system (1) − (3)
Additionally, the electric and magnetic fields figure out exists and is unique, for details see [74]. The Maxwell’s
directly both in the error estimates and numerical experi- equations with piecewise constant coefficient have been
ments. For simplicity of presentation, the material parame- investigated in [16], [75].
ters ε and µ are considered as time independent, piecewise The paper is structured as follows. The next section
constant scalar functions, but the results can be general- gives an overview over the basic function spaces and nota-
ized to more complicated material parameters, e.g. positively tion. Section III describes the weak formulation and an
definite tensors. The electric and magnetic fields are dis- error estimate for the backward Euler semi-discrete method
cretized by Nédeléc curl-conforming and Raviart-Thomas (Rothe method). The spatial discretization is discussed in
div-conforming finite elements in space, respectively. The Section IV. In Section V we describe and investigate the
properties of these finite elements have been investigated in full discretization. Finally Section VI presents a collection of
many articles [71], [72] and [73]. In addition, the problem is numerical examples.
discretized in time by backward Euler and sympletic meth-
ods. The error analysis of the mixed finite element method II. SPACES AND NOTATION
for the fully discrete problem is given in the case of only the For a real number p ≥ 1, the space L p () consists of equiva-
backward Euler method, and computational results are given lence classes of Lebesgue-measurable functions u :  → R
in both cases. such that
Z
Our proposed schemes deal with the system of Maxwell’s |u|p dx < ∞.
equations directly in 3D, which cause no spurious solution. 
Fully-discrete error estimates and simulation results with If u ∈ L p (), we define its L p ()-norm as follows :
visualizations of the electromagnetic quantities are given. Z 1/p
Similar results for fully-discrete error estimation could be
kukL p () = p
|u(x)| dx .
obtained also for our previous results about the decou- 
pled Maxwell’s equations [32], [36]. These are intermediate ∞
Furthermore, the space L () consists of the equiva-
results that provide a starting point for the development and
lence classes of essentially bounded measurable functions
theoretical-numerical investigation of TDFEM for nonlinear
u :  → R equipped with the norm
problems in optics and photonics. These include energy-
conserving methods in 3D. kukL ∞ () = ess sup |u(x)|.
To begin with the problem formulation, let  ⊂ R3 be a x∈
simply connected domain with a smooth boundary 0 and unit The analogous spaces of vector fields u :  → R3 are
outward normal n. The symbols E = E(x, t) and H = H(x, t) denoted by Lp () := [L p ()]3 .
denote the electric and magnetic field intensities respectively, In what follows we have to deal with weighted function
where the time variable t belongs to some finite interval spaces. Given a weight ω :  → R, where the values of
(0, T ), 0 < T < ∞. Given a current density function ω are positive a.e. on , we define a weighted inner product
J = J(x, t), specifying the applied current, Maxwell’s and a weighted norm by
equations state that Z
εEt − ∇ × H = J in  × (0, T ), (1) (u, v)ω := ω u · v dx and,

µHt + ∇ × E = 0 in  × (0, T ). (2)
p
kukω := kukL2ω () := (u, u)ω ,
The material parameters ε and µ do not depend on time
and the space L2ω () consists of vector fields u :  → R3
and are piecewise constant. In addition, there exist positive
with Lebesgue-measurable components and such that
constants εmin , εmax , µmin , µmax such that, for all x ∈ ,
0 < εmin ≤ ε(x) ≤ εmax < ∞, kukω < ∞.
0 < µmin ≤ µ(x) ≤ µmax < ∞. In the case ω = 1, the subscript is omitted.

63854 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

As transient problems are addressed, we will work with Furthermore, we need the following Hilbert spaces that are
functions that depend on time and have values in certain related to the (weak) rotation and divergence operators:
Banach spaces. If u = u(x, t) is a vector field of the space
variable x and the time variable t, it is suitable to separate H(curl; ) := {u ∈ L2 () : ∇ × u ∈ L2 ()},
these variables in such a way that u(t) = u(·, t) is considered H0 (curl; ) := {u ∈ H(curl; ) : u × n|0 = 0},
as a function of t with values in a Banach space, say X , with H(div; ) := {u ∈ L2 () : ∇ · u ∈ L2 ()},
the norm k · kX . That is, for any t ∈ (0, T ), the mapping x 7 →
H0 (div; ) := {u ∈ H(div; ) : u · n|0 = 0}.
u(x, t) is interpreted as a parameter-dependent element u(t)
of X . In this sense we will write E(t) = E(·, t), H(t) = H(·, t) These Hilbert spaces are equipped with the norms (resp.
and so on. induced norms)
The space Cm (0, T , X ), m ∈ N ∪ {0}, consists of all 1/2
kukH(curl;) := kuk20 + k∇ × uk20 ,

continuous functions u : (0, T ) → X that have continuous
derivatives up to order m on (0, T ). It is equipped with the 2 1/2
2
.

kukH(div;) := kuk0 + k∇ · uk0
norm
m
X We refer to [73], [75], [78] and [79] for details about these
sup ku(j) (t)kX . spaces.
j=0 t∈(0,T )

For the sake of consistency in the notation we will write III. WEAK FORMULATION
C(0, T , X ) := C0 (0, T , X ). Given J ∈ C(0, T , L2ε−1 ()), and the weak solution
The space Lp (0, T , X ) with 1 ≤ p < ∞ contains C(0, T , H0 (curl; )) ∩ C1 (0,
T , L2ε ()) ×


(E, H) ∈
(equivalent classes of) strongly measurable functions u : C(0, T , H(div; )) ∩ C1 (0, T , L2µ ()) of the
(0, T ) → X such that system (1)–(2) satisfies
Z T
p
ku(t)kX dt < ∞ (εEt , 9) − (H, ∇ × 9) = (J, 9) ∀ 9 ∈ H0 (curl; ),
0 (4)
(for the definition of strongly measurable functions we refer
(µHt , 8) + (∇ × E, 8) = 0 ∀ 8 ∈ H(div; ), (5)
to [76]). The norm on Lp (0, T , X ) is defined by
Z T 1/p for t ∈ (0, T ) with the initial conditions (3).
p
kukLp (0,T ,X ) := ku(t)kX dt . Theorem 1: Let J ∈ C(0, T , L2ε−1 ()), and the solution
0
C(0, T , H0 (curl; )) ∩
C1 (0, T , L2ε ()) ×


(E, H) ∈
These spaces can be equipped with a weight, too. In particu- C(0, T , H(div; ))∩C1 (0, T , L2µ ()) of the system (1)–(2)
lar, we will write satisfies
Z T Z 1/2
√
 
kukL2 (0,T ,L2ω ()) := |u(t)|2 ω dxdt . kEkε +kHkµ ≤ 2 kE0 kε + kH0 kµ +kJkL1 0,T ,L2 ()
.
0  ε−1
Next we introduce the Sobolev spaces of functions with weak Proof: Take test functions 9 := E and 8 := H
spatial derivatives of maximal order r ∈ N in L p (), where in (4)–(5) respectively, then
α is a multi-index:
W r,p () := u ∈ L p () : ∂ α u ∈ L p () ∀ |α| ≤ r .
 (εEt , E) − (H, ∇ × E) = (J, E), (6)
(µHt , H) + (∇ × E, H) = 0. (7)
For p ≥ 1, the norms and semi-norms are defined by
XZ Adding the equations (6)–(7), we have
k∂ α ukp dx,
p
kukW r,p () :=
|α|≤r  (εEt , E) + (µHt , H) = (J, E),
XZ
k∂ α ukp dx.
p
|u|W r,p () := therefore
|α|=r 
dh i
The modifications for p = ∞ are obvious. If p = 2, we write kEk2ε + kHk2µ = 2(J, E) ≤ 2kJkε−1 kEkε .
dt
H r () := W r,2 () and k · kH r () := k · kW r,2 () .
Integrating both sides from 0 to T , we obtain,
The space H01 () is defined as the closure of C0∞ () with
respect to the norm k·kH 1 () , where C0∞ () denotes the space kEk2ε + kHk2µ
of all arbitrarily often differentiable functions with compact Z T
support on . It is well known that H01 () is a closed subspace ≤ kE0 k2ε + kH0 k2µ + 2 kJkε−1 kEkε ds
of H 1 () and consists of elements u such that u = 0 on 0 0
Z T
in the sense of traces [77]. As in the case of the L p -spaces,
q
≤ kE0 k2ε + kH0 k2µ + 2 kJkε−1 kEk2ε + kHk2µ ds.
we shall write Wr,p () := [W r,p ()]3 and so on. 0

VOLUME 7, 2019 63855

 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

Then it follows from the Gronwall–Ou-Iang’s inequality Using this estimate in (12), we obtain
(see, e.g., [80]) that
kEn k2ε − kEn−1 k2ε + kHn k2µ − kHn−1 k2µ
q q Z T
kEk2ε + kHk2µ ≤ kE0 k2ε + kH0 k2µ + kJkε−1 ds. ≤ 1tkEn k2ε + 1tkJn k2ε−1 .
0
Summing up from n = 1 to N , we arrive at
Since
q √ q kEN k2ε + kHN k2µ
kEk2ε + kHk2µ ≤ kEkε + kHkµ ≤ 2 kEk2ε + kHk2µ , XN N
X
≤ 1t kEn k2ε + 1t kJn k2ε−1 + kE0 k2ε + kH0 k2µ
the statement follows. 
n=1 n=1
Let us now turn to time discretizations for the Maxwell N
X N
X
system (1)–(2). The time interval (0, T ) is divided into N ∈ N ≤ 1t kEn k2ε + 1t kJn k2ε−1 + kE0 k2ε + kH0 k2µ .
equally spaced subintervals by using nodal points n=0 n=0

0 =: t 0 < t 1 < t 2 < · · · < t N := T , Now we are ready to apply Gronwall’s inequality (Lemma 2)
with δ := 1t ≥ 0, g0 := kE0 k2ε + kH0 k2µ ≥ 0, an :=
with t n = n1t, 1t > 0, n = 0, 1, 2, . . . , N . kEn k2ε + kHn k2µ ≥ 0, bn := 0, cn := kJn k2ε−1 ≥ 0, and
Replacing the time derivatives in (4)–(5) at t n by the back- γn := 1 ≥ 0. Then the condition γn δ < 1 corresponds to
ward finite difference quotient, that is 1t < 1 and the final estimate follows from the observation
that (n + 1)1t ≤ T + 1t ≤ T + 1/2:
E(t n ) − E(t n−1 )
Et (t n ) ≈ etc.,
1t kEN k2ε + kHN k2µ
and get a sequence of problems of the type N
X
≤ 1t kJn k2ε−1 + kE0 k2ε + kH0 k2µ


 En − En−1 
ε , 9 − (Hn , ∇ × 9) = (Jn , 9) n=0
1t N
∀ 9 ∈ H0 (curl; ),
 X 
(8) × exp 1t (1 − 1t)−1
 Hn − Hn−1 
µ , 8 + (∇ × En , 8) = 0 n=0
1t N
!
∀ 8 ∈ H(div; ),
X
(9) ≤ kJn k2ε−1 1t + kE0 k2ε + kH0 k2µ exp (2T + 1) .
n=0
where (En , Hn ) ∈ H0 (curl; ) × H(div; ) are to be
determined for n ∈ {1, . . . , N } (as approximations to N
X
(E(t n ), H(t n ))), (E0 , H0 ) ∈ H0 (curl; )×H(div; ) are given It remains to note that the term kJn k2ε−1 1t is an approx-
(as approximations to (E0 , H0 )), and Jn := J(t n ) ∈ L2ε−1 (). Z T
n=0
Theorem 2: For 0 < 1t < 1/2, there exists a constant imation to kJ(t)k2ε−1 dt = kJk2L2 (0,T ,L2 and thus
())
C > 0 independent of 1t (but dependent on T ) such that 0 ε−1
bounded. 
kEN k2ε + kHN k2µ ≤ C. Next we want to prove an estimate of the error in time. To
Proof: In (8)–(9), we choose the test functions 9 := do so, we introduce the errors
21tEn and 8 := 21tHn . Then
ζcn := E(t n ) − En , ξcn := H(t n ) − Hn . (13)
n
2(ε(E −E n−1
), E )−21t(H , ∇ × E ) = 21t(J , E ),
n n n n n

(10) Theorem 3: If
(E, H) ∈ C(0, T , H0 (curl; )) ∩

n
2(µ(H −H n−1
), H )+21t(∇ × E , H ) = 0.
n n n
(11) C2 (0, T , L2ε ()) × C(0, T , H(div; )) ∩ C2 (0, T , L2µ ())
and if the time step 1t is sufficiently small, then there exists
Adding the equations (10) and (11), we get a constant C > 0 independent of 1t (but dependent on T )
2(ε(En −En−1 ), En )+2(µ(Hn −Hn−1 ), Hn ) = 21t(Jn , En ). such that
kζcN k2ε + kξcN k2µ ≤ C 1t + kζc0 kε + kξc0 kµ .


The identity (1) from Lemma 1 implies that
Proof: From Taylor’s formula with integral remainder
kEn k2ε + kEn − En−1 k2ε − kEn−1 k2ε + kHn k2µ it follows that
Z t
+ kHn −Hn−1 k2µ − kHn−1 k2µ = 21t(Jn , En ). (12) n n n
E(t) = E(t ) + Et (t )(t − t ) + (t − s)Ett (s)ds,
tn
The right-hand side is estimated similarly to the proof of the
inequality (2) from Lemma 1: hence
 E(t n ) − E(t n−1 ) 
2(Jn , En ) = 2(ε −1/2 Jn , ε1/2 En ) ≤ kε −1/2 Jn k2 ε , 9 = (εEt (t n ), 9) + (εRnE , 9),
1t
+ kε 1/2 En k2 = kJn k2ε−1 +kEn k2ε . (14)

63856 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

where It remains to estimate the sum terms. Since

tn
Z Z n 2
1 1 t
RnE := (t n−1 − s)Ett (s)ds. kRnE k2ε =

(t n−1
− s)E (s)ds.

1t t n−1 (1t)2 t n−1
tt
ε

An analogous relation can be obtained with respect to H. Z tn Z tn
1
Making use of (4)–(5), we get ≤ (s − t n−1 )2 ds kEtt (s)k2ε ds
(1t)2 t n−1 t n−1
Z n
 E(t n ) − E(t n−1 )  1t t
ε , 9 − (H(t n ), ∇ × 9) = kEtt (s)k2ε ds ,
1t 3 t n−1
= (Jn + εRnE , 9) ∀ 9 ∈ H0 (curl; ), (15) it follows that
 H(t n ) − H(t n−1 )  N
(1t)2 T
Z
µ , 8 + (∇ × E(t n ), 8)
X
kRnE k2ε 1t ≤ kEtt (s)k2ε ds
1t 3 0
n=0
= (µRnH , 8) ∀ 8 ∈ H(div; ). (16)
(1t)2
= kEtt k2L2 (0,T ,L2 ()) .
Now the equations (8)-(9) are subtracted from the 3 ε

equations (15)-(16). Together with (13) the result is: In summary, we get
 ζ n − ζ n−1
kζcN k2ε + kξcN k2µ

ε c c
, 9 − ξcn , ∇ × 9 = (εRnE , 9)


(17)
1t (1t)2 
kEtt k2L2 (0,T ,L2 ()) + kHtt k2L2 (0,T ,L2 ())

 ξ n − ξ n−1  ≤
µ c c
, 8 + ∇ × ζcn , 8 = (µRnH , 8).


(18) 3 ε µ
1t + kζc0 k2ε + kξc0 k2µ exp (2T + 1) .

Taking 9 := 21tζcn and 8 := 21tξcn in equations (17) and

(18), by the same reasoning as in the proof of Theorem 2, 
from this we obtain the estimate
IV. SPATIAL DISCRETIZATION
kζcn k2ε − kζcn−1 k2ε + kξcn k2µ − kξcn−1 k2µ In this section we introduce families of finite dimensional
h i subspaces Uh ⊂ H(curl; ) and Vh ⊂ H(div; ) to discretize
≤ 1t kζcn k2ε + kξcn k2µ + 1t kRnE k2ε + kRnH k2µ .
 
the problem (8)–(9) in space.
Let Th be an arbitrary member of a family of triangula-
The summation leads to
tions of  consisting of geometric elements K . Each element
kζcN k2ε + kξcN k2µ K ∈ Th is assumed to be an open tetrahedron if K has no
N face or edge on 0. The elements that have an edge or face
on 0 are allowed to have one curved edge or one curved
h X
≤ 1t
 n 2
kζc kε + kξcn k2µ


n=1
face, respectively. These elements are called boundary ele-
N h i i ments [16]. If K is a boundary element we assign a standard
tetrahedron K̃ by connecting the four vertices of K by straight
X
+ 1t kRnE k2ε + kRnH k2µ + kζc0 k2ε + kξc0 k2µ .
n=1
edges.
In addition, all triangulations should be compatible with
Next we apply Gronwall’s inequality (Lemma 2) with δ := the discontinuities of the coefficients ε and µ, that is their
1t ≥ 0, g0 := kζc0 k2ε + kξc0 k2µ ≥ 0, an := kζcn k2ε + kξcn k2µ ≥ discontinuities lie on the boundaries of the elements of the
0, bn := 0, cn := kRnE k2ε + kRnH k2µ ≥ 0 (n ≥ 1), c0 := 0, triangulations only. Moreover we assume that the family of
and γn := 1 ≥ 0. Then the condition γn δ < 1 corresponds to triangulations is quasi-uniform. That is, there exist constants
1t < 1 and we get cF > 0, cF > 0 independent of K and Th such that
kζcN k2ε + kξcN k2µ cF h ≤ hK ≤ cF ρK ∀K ∈ Th ∀Th ,
N h
where ρK is the maximum diameter of the largest ball con-
h X i
≤ 1t kRnE k2ε + kRnH k2µ + kζc0 k2ε + kξc0 k2µ


tained in K or K̃ , hK is the diameter of K and h :=
n=0
N
maxK ∈Th hK [78].

i Now let h be the interior of the set
X
× exp 1t (1 − 1t)−1
[
n=0
K̃ .
N h
K ∈ Th
h X i
kRnE k2ε + kRnH k2µ 1t + kζc0 k2ε + kξc0 k2µ


≤
n=0 Let Pk be the space of scalar real-valued polynomials in three
variables of maximal degree k, and P̃k be the space of scalar
i
× exp (2T + 1) .
real-valued homogeneous polynomials of exact degree k.
VOLUME 7, 2019 63857
 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

The following estimate holds for rh ( [71], Thm. 2): If v ∈

Hk+1 () then
kv − rh vkH(curl;) ≤ Chk kvkHk+1 () .
Definition 2 (Dk -Unisolvent and Div-Conforming Dofs):
Let K be a tetrahedron in R3 with faces denoted by f with
normal n. Let v ∈ W1,p (K ) for some p > 2. We define the
following two set of moments of v on K :
For the four facesf of K
Z 
FIGURE 1. A tetrahedron K : t is an edge tangent vector, n is a face normal.
Mf (v) := v · n q, ds ∀q ∈ Pk−1 (f ) ,
f
Z 
MK (v) := v · q dx ∀q ∈ Pk−2 (K ) .
K
For any k ∈ N, we define the following subspaces of Pk :=
[Pk ]3 (for details see [65]), [71] or [72]: This set of moments is Dk -unisolvent and div-conforming as
proved in [71], Thm. 3. For any v ∈ W1,p (K ), we define a
Dk := Pk−1 ⊕ x P̃k−1 , local interpolant wK v ∈ Dk as follows:
Sk := {p ∈ [P̃k ]3 : x · p(x) = 0}, Mf (v − wK v) = MK (v − wK v) = {0}.
Rk := Pk−1 ⊕ Sk . The global interpolant wh v ∈ Vh := {w ∈ H(div; ) :
w|K ∈ Dk ∀K ∈ Th } is defined element-wise:
Obviously, Sk ⊂ Pk and Rk ⊂ Pk .
Here we describe the so-called first family of Nédélec edge wh v|K := wK (v|K ) ∀K ∈ Th .
elements. We mention that in the two definitions below it is An error estimate also holds for wh (see [18], equation (19),
assumed that a standard reference tetrahedron K̂ is used and and [71], Thm. 4): If v ∈ Hk (), then
that an affine transformation between K̂ and K is applied. kv − wh vkε ≤ Chk kvkHk () .
Definition 1 (Rk -Unisolvent and Curl-Conforming Dofs):
According to [65], Lemma 5.40, the spaces Uh and Vh are
Let K be a tetrahedron in R3 with faces denoted by f and
related via
edges denoted by e. t and n in Fig (1) represent the unit
vectors along the edge e and perpendicular to the face f , ∇ × Uh ⊂ Vh , (19)
respectively. Let v ∈ W1,p (K ) for some p > 2. We define the and the interpolation operators rh and wh are linked together
following three sets of moments of v on K : as follows: ∇ × rh v = wh (∇ × v) for all v such that both the
interpolants rh v and wh (∇ × v) are defined.
For the six edges of e of K In the subsequent error analysis, we will make use of
Z 
Me (v) := v · t q de ∀q ∈ Pk−1 (e) , two projection operators into U0h := Uh ∩ H0 (curl; )
e and Vh , respectively. The first projection operator 5h :
for the four facesf of K
Z  H0 (curl; ) → U0h is defined by
v × n · q ds ∀q ∈ [Pk−2 (f )]2 , ∇ × 5h u, 9h = ∇ × u, 9h ∀9h ∈ Vh ,



Mf (v) :=
f
Z  (5h u, ∇ph ) = (u, ∇ph ) ∀ph ∈ Shk , (20)
MK (v) := v · q dx ∀q ∈ Pk−3 (K ) . where Shk is the standard space of continuous finite elements
K
on Th :
Remark 1: The dofs Mf are given here in the original Shk := {v ∈ H01 () : v|K ∈ Pk ∀K ∈ Th },
version of [71], Def. 4. However, from the point of view of see [17].
affine equivalence it may be useful to use a different, but An error bound for the projection 5h can be derived as
equivalent representation [65], Remark 5.31. in [18]: If u ∈ Hk+1 (), then
This set of moments is Rk -unisolvent and curl-conforming as
proved in [71], Thm. 1. For any v ∈ W1,p (K ), we define a ku − 5h ukε ≤ Chk kukHk+1 () . (21)
local interpolant rK v ∈ Rk such that The second projection operator is the L2 -orthogonal projec-
tion of H(div; ) onto Vh :
Me (v − rK v) = Mf (v − rK v) = MK (v − rK v) = {0}
(Ph H, 8h ) = (H, 8h ) ∀8h ∈ Vh . (22)
The global interpolant rh v ∈ Uh := {w ∈ H(curl; ) : By the help of similar arguments as in [18], the following
w|K ∈ Rk ∀K ∈ Th } is defined element-wise: estimate can be shown for v ∈ Hk ():

rh v|K := rK (v|K ) ∀K ∈ Th . kPh v − vkµ ≤ Chk kukHk () . (23)

63858 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

V. FULL DISCRETIZATION USING THE BACKWARD Now we set 9h := −21tηhn and 8h := −21tθhn :
EULER METHOD
2 ε(ηhn −ηhn−1 ), ηhn −21t(θhn , ∇ ×ηhn ) = 2 ε(ηn −ηn−1 ), ηhn



The fully discrete electric and magnetic fields (Enh , Hnh ) ∈
U0h × Vh satisfy −21t(θ n , ∇ ×ηhn )−21t(εRnE , ηhn )
2 µ(θhn −θhn−1 ), θh +21t ∇ ×ηhn , θhn = 2 µ(θ n −θ n−1 ), θhn
n




 En − En−1 
ε h h
, 9h − (Hnh , ∇ × 9h ) + 21t ∇ ×ηn , θhn −21t(µRnH , θhn ).


1t
= (Jn , 9h ) ∀ 9h ∈ U0h , (24) Thanks to (20), (22) and (19), the middle terms on the right-
 Hn − Hn−1 
µ h h
, 8h + (∇ × Enh , 8h ) hand sides vanish. Adding the resulting equations, we get
1t
= 0 ∀ 8h ∈ Vh , 2 ε(ηhn − ηhn−1 ), ηhn + 2 µ(θhn − θhn−1 ), θhn



(25)
= 2 ε(ηn − ηn−1 ), ηhn + 2 µ(θ n − θ n−1 ), θhn



where (Enh , Hnh ) ∈ U0h × Vh are to be determined for n ∈
{1, . . . , N } and (E0h , H0h ) ∈ U0h × Vh are given (and will − 21t(εRnE , ηhn ) − 21t(µRnH , θhn ).
be specified later). Before formulating the theorem, we still
Now the identity (1) from Lemma 1 allows to rewrite the left-
introduce a few terms. The full error of the electric field will
hand side as,
be denoted by
kηhn k2ε + kηhn − ηhn−1 k2ε − kηhn−1 k2ε + kθhn k2µ
ζ n := E(t n ) − Enh = ηn − ηhn , (26)
+ kθhn − θhn−1 k2µ − kθhn−1 k2µ
where
= 2 ε(ηn − ηn−1 ), ηhn + 2 µ(θ n − θ n−1 ), θhn



ηn := E(t n ) − 5h E(t n ), ηhn := Enh − 5h E(t n ). (27) − 21t(εRnE , ηhn ) − 21t(µRnH , θhn ). (30)
Similarly for the magnetic field: The first two terms on the right-hand side are treated by
ξ := H(t
n n
) − Hnh =θ n
− θhn (28) means of formula (14). Namely, replacing there E(t n ) by
(I − 5h )E(t n ), we obtain
with
ε(ηn − ηn−1 ), ηhn = 1t(ε(I − 5h )Et (t n ), ηhn )


θ n := H(t n ) − Ph H(t n ), θhn := Hnh − Ph H(t n ). (29) + 1t(ε(I − 5h )RnE , ηhn ).
Theorem 4: Let (E, H) be the solution of (4)–(5) such that, In the same way we get
for some k ∈ N,
µ(θ n − θ n−1 ), θhn = 1t(µ(I − Ph )Ht (t n ), θhn )


E ∈ C(0, T , H0 (curl; )) ∩ C 2 (0, T , L2ε () ∩ Hk+1 ()),
+ 1t(µ(I − Ph )RnH , θhn ).
H ∈ C 2 (0, T , L2µ () ∩ Hk ()),
Using this in (30) and applying the estimate (2) from
and (Enh , Hnh ) be the fully discrete solution of (24)-(25). Then, Lemma 1 to each of the resulting terms on the right-hand side,
for sufficiently small 1t and h, the following error estimate we obtain
holds:
kηhn k2ε + kηhn − ηhn−1 k2ε − kηhn−1 k2ε + kθhn k2µ
kζ N kε + kξ N kµ ≤ C 1t + hk + hk 1t ,
 
+ kθhn − θhn−1 k2µ − kθhn−1 k2µ
where the constant C > 0 does not depend on 1t and h (the ≤ 1tk(I−5h )Et (t n )k2ε +1tk(I−5h )RnE k2ε +1tkRnE k2ε
structure of C will be seen from the proof).
+ 31tkηhn k2ε +1tk(I−Ph )Ht (t n )k2µ +1tk(I−Ph )RnH k2µ
Proof: Taking 9 = 9h and 8 = 8h in (15)-(16),
subtracting the system (24)-(25) from the system (4)–(5) and + 1tkRnH k2µ + 31tkθhn k2µ
using the definitions (26), (28), we obtain: Summing up from n = 1 to N and ignoring the second and
 ζ n − ζ n−1  fifth terms on the left-hand side, we arrive at
ε , 9h − (ξ n , ∇ × 9h ) = (εRnE , 9h ),
 ξ n −1tξ n−1
N
X
kηhN k2ε + kθhN k2µ ≤ 31t
 n 2
kηh kε + kθhn k2µ
 
µ , 8h + (∇ × ζ n , 8h ) = (µRnH , 8h ).
1t n=1
Furthermore, using the decompositions (27), (29), after a N h
X
simple rearrangement we arrive at + 1t k(I−5h )Et (t n )k2ε +k(I−5h )RnE k2ε +kRnE k2ε
 (ηn − ηn−1 ) − (ηn − ηn−1 )  n=1
i
ε h h
, 9h − (θ n − θhn , ∇ × 9h ) + k(I − Ph )Ht (t n )k2µ + k(I − Ph )RnH k2µ + kRnH k2µ
1t
= (εRnE , 9h ),
+ kηh0 k2ε + kθh0 k2µ .
 (θ n − θ n−1 ) − (θ n − θ n−1 ) 
µ h h
, 8h + (∇ × (ηn − ηhn ), 8h ) Now we are ready to apply Gronwall’s inequality (Lemma 2)
1t
= (µRnH , 8h ). with δ := 1t ≥ 0, g0 := kηh0 k2ε + kθh0 k2µ ≥ 0,

VOLUME 7, 2019 63859

 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

an := kηhn k2ε +kθhn k2µ ≥ 0, bn := 0, cn := k(I−5h )Et (t n )k2ε + we conclude that,

k(I − 5h )RnE k2ε + kRnE k2ε + k(I − Ph )Ht (t n )k2µ + k(I −
kηhN kε + kθhN kµ ≤ C 1t + hk + hk 1t .
 
Ph )RnH k2µ +kRnH k2µ ≥ 0, and γn := 3 ≥ 0. Then the condition
γn δ < 1 corresponds to 1t < 1/3 and we obtain, say for The terms kηN kε and kθ N kµ are estimated by (21) and (23),
1t < 1/6, respectively:
kηhN k2ε + kθhN k2µ kηN kε = k(I − 5h )E(t N )kε ≤ Chk kE(t N )kHk+1 () ,
N h
kθ N kµ = k(I − Ph )H(t N )kµ ≤ Chk kH(t N )kHk () .
 X
≤ 1t k(I − 5h )Et (t n )k2ε + k(I − 5h )RnE k2ε
n=0 The triangle inequality yields the stated result:

+ kRnE k2ε + k(I − Ph )Ht (t n )k2µ + k(I − Ph )RnH k2µ kE(t N ) − EN N N

h kε + kH(t ) − Hh kµ

i  ≤ kηN kε + kηhN kε + kθ N kµ + kθhN kµ

+ kRnH k2µ + kηh0 k2ε + kθh0 k2µ exp (6T + 1) .
≤ C 1t + hk + hk 1t .
 

From the end of the proof of Thm. 3 it is known that 

N
VI. NUMERICAL RESULTS
X
1t kRnE k2ε + kRnH k2µ
 

n=0
Denoting by en , hn and jn the representation vectors of Enh , Hnh
and the L2 -orthogonal projection of Jn onto U0h , the
(1t)2  method (24)–(25) can be written as follows:
kEtt k2L2 (0,T ,L2 ()) + kHtt k2L2 (0,T ,L2 ()) .

≤
3 ε µ
en − en−1
Furthermore, from (21) we see that, Mε = G> hn + jn , (32)
1t
k(I − 5h )Et (t n )kε ≤ Chk kEt (t n )kHk+1 () , hn − hn−1
Mµ = −Gen , (33)
1t
k(I − 5h )RnE kε ≤ Chk kRnE kHk+1 () ,
where Mε is the positively definite mass matrix of size
and (23) yields the estimates, dim U0h × dim U0h for the material parameter ε, Mµ is the
positively definite mass matrix of size dim Vh × dim Vh for
k(I − Ph )Ht (t n )kµ ≤ Chk kHt (t n )kHk ()
the material parameter µ, and G is a discrete representation
k(I − Ph )RnH kµ ≤ Chk kRnH kHk () . of −curl with size dim Vh × dim U0h .
Therefore,
A. IMPLEMENTATION OF THE BACKWARD EULER
N
X METHOD
1t k(I − 5h )Et (t n )k2ε + k(I − Ph )Ht (t n )k2µ
 
The formal algorithm for the system (32)-(33) reads as
n=0 follows:
N
≤ Ch2k
X 
kEt (t n )k2Hk+1 () + kHt (t n )k2Hk () 1t
 Compute the number of time steps:
n=0
T − t0
nstep :=
Z T 1t
≤ Ch2k kEt (t)k2Hk+1 () + kHt (t)k2Hk () dt Compute the initial values for the electric and magnetic
 
0
fields:
= Ch kEt k2L2 (0,T ,Hk+1 ()) + kHt k2L2 (0,T ,Hk ()) ,
2k
 
e0 ← E0
and h0 ← H0
N
X Loop over time steps:
1t k(I − 5h )RnE k2ε + k(I − Ph )RnH k2µ
 
for n = 1 to nstep do:
n=0
N
X Begin integration method update:
kRE kHk+1 () + kRnH k2Hk () 1t
2k
 n 2 
≤ Ch ein ← en−1
n=0
hin ← hn−1
≤ Ch2k (1t)2 kEtt k2L2 (0,T ,Hk+1 ()) + kHtt k2L2 (0,T ,Hk ()) .
 
Update the electric and magnetic fields values:
Finally, if we take eout ← ein + 1tM−1 > n


ε G hout + j
E0h := 5h E0 , H0h := Ph H0 , (31) hout ← hin + 1tM−1
µ Geout

63860 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

Update the electric and magnetic fields values for this

time step:
en ← eout
hn ← hout
end for
Completion:
eN ← enstep
hN ← hnstep

B. IMPLEMENTATION OF THE SYMPLECTIC METHOD

We have also tested the application of symplectic time inte-
gration methods. A 4th order symplectic integration algo-
rithm for the time dependent Maxwell’s equations with
parameters
FIGURE 2. The energy of the system remains constant if the symplectic
1
− 13 time integration method is applied.
2+2 +2 3
β1 = , α1 = 0,
6
1 1
1 − 2 3 − 2− 3 1 Completion:
β2 = , α2 = 1
,
6 2−2 3 eN ← enstep
1
− 13
1−2 −2 3 1 hN ← hnstep
β3 = , α3 = 2
,
6 1−2 3
1 Krylov solvers were used to invert the mass matrices Mε
− 13
2+2 +2 3 1 and Mµ . A preconditioned conjugate gradient solver was also
β4 = , α4 = 1
6 2 − 23 implemented.
(see [81], [82]) is given as follows:
C. ENERGY CONSERVATION
Compute the number of time steps: In order to verify a correct physical behaviour of the numer-
T − t0
nstep := ical methods, we considered the energy evolution, too. The
1t discrete instantaneous energy is the total energy that is stored
Compute the initial values for the electric and magnetic
in the discrete electric and magnetics fields. It is computed as
fields:
e0 ← E0 1 T
e Mε e + h T Mµ h .


Energy =
h0 ← H0 2
Loop over time steps: It is an substantial advantage of the symplectic method that it
conserves the energy, see Fig. 2.
for n = 1 to nstep do:
Begin integration method update: D. SIMULATION RESULTS, VALIDATIONS AND
ein ← en−1 DISCUSSION
hin ← hn−1 A number of numerical experiments were performed to
Update the electric and magnetic fields values: approximate solutions of time dependent Maxwell’s prob-
lems. We visualized the electromagnetic fields for cases
for j = 1 to 4 do
where the exact solution is known, but also for cases with
eout ← ein + αj 1tM−1 > n


ε G hin + j unknown analytical solution, and checked the stability and
hout ← hin + βj 1tM−1
µ Geout
convergence properties in problems with complicated geome-
ein ← eout tries. The main object of these simulations was to validate the
code. The simulations are conditionally stable in the case of
hin ← hout the symplectic time integration method.
end for EXAMPLE 1: This test example is characterized by the
Update the electric and magnetic fields values for this following parameters, where a symplectic time integration
time step: method
√ is applied to a Fichera mesh. The frequency is f =
3
2 0 c Hz, where c0 denotes the speed of light in vacuum, and
en ← eout
the wave length is λ = 1.1547 m. The angular frequency is
hn ← hout ω = 2πf (rad·s−1 ). The permittivity and the permeability are
end for equal to the constant vacuum values ε = ε0 and µ = µ0 as

VOLUME 7, 2019 63861

 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

FIGURE 3. The snapshot is taken at the final time step (N = 100) for the electric and magnetic fields, by employing the backward
Euler method for the beam tetrahedron. The time step size is 1t = 0.03125. (a) Electric field. (b) Magnetic field.

FIGURE 4. The snapshot is taken at the first time step (n = 1) for the electric and magnetic fields by taking the
projections (31), and employing the backward Euler method for the Fichera mesh. The time step size is 1t = 0.0005.
(a) Electric field. (b) Magnetic field.

in [83]. The exact electric and magnetic fields are given as obtained by taking the projections (31) of the exact electric
and magnetic fields, where the exact fields given by
E1 (t) = − cos(π x) sin(π y) sin(π z) cos(ωt),
 T
E2 (t) = 0, E= sin(2t − 3z), sin(2t − 3x), sin(2t − 3y) ,
E3 (t) = sin(π x) sin(π y) cos(π z) cos(ωt),  T
π H= sin(2t − 3y), sin(2t − 3z), sin(2t − 3x) ,
H1 (t) = − sin(π x) cos(πy) cos(πz) sin(ωt),
ω  T
2π J= sin(2t − 3z), sin(2t − 3x), sin(2t − 3y) .
H2 (t) = cos(π x) sin(π y) cos(π z) sin(ωt),
ω
π
H3 (t) = − cos(π x) cos(π y) sin(πz) sin(ωt). Theorem 1 states that the problem is well-posed. An anal-
ω ogous result for the Rothe method is given in Theorem 2.
For the case of the symplectic time integration method, The Theorems 3 and 4 present a priori estimates of the abso-
an upper bound of the largest stable time step (CFL) is given lute error. These results show that we get optimal solutions
by [2], [5], [81], as within the selected finite element spaces, whereas many exist-
ing methods exhibit spurious solutions e.g. [1]–[15], [18],
2
1t ≤ q
,
[27], [32], [34], [36] also because they solve other than the
ρ M−1 > −1
ε G Mµ G direct Maxwell’s problem (1)-(3). We measured the L2 -norm
of the error for a sequence of successively refined meshes
where ρ is the spectral radius function. The largest eigenvalue starting from a uniform coarse mesh. The refinement level at
can be efficiently determined by [84] or the power method. l = 1 shows the initial geometry of the mesh, and the levels
EXAMPLE 2: Here the backward Euler method is con- at l = 2, l = 3, and l = 4 show the uniform refinement
sidered. The permittivity and the permeability are constant: simultaneously at the subsequent 2nd, 3rd and 4th steps.
ε = 2, µ = 1.5. The initial electric and magnetic fields are We summarize the obtained absolute errors for the fourth

63862 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

FIGURE 5. The scale shows the values of electric fields at the final time step (N = 50), by employing the backward
Euler method. The time step size is 1t = 0.0005. (a) Snapshot of top. (b) Snapshot of bottom. (c) Snapshot of top and
front. (d) Snapshot of back and bottom. (e) Snapshot of back, left and bottom. (f) Snapshot of right and front.

TABLE 1. Absolute error.

order symplectic integration method in Table 1. The Table 1 time step by taking the projections (31), and by employing the
shows that the symplectic method is conditionally stable and backward Euler method for the Fichera mesh (3D L-shaped
its order of the convergence is approximately equal to 4. domain). The time step size is 1t = 0.0005 in Fig. 4. Differ-
Fig 2 illustrates that the symplectic method conserves the ent orientations of the Fichera mesh are illustrated in Figs. 5.
energy. This is an additional aspect to underline the good Furthermore, the electric and magnetic fields at time step
accuracy of our numerical results. N = 50 are shown in Fig. 5, by employing the backward
The snapshot of the electric and magnetic fields depicted Euler for 1t = 0.0005. Both the absolute errors and the con-
in Fig. 3 is taken at the final time step N = 100, using the servation property of energy are determined and are strictly
time step size 1t = 0.03125, by employing the backward fulfilled as the electric and magnetic fields visualization in 3D
Euler method for the beam tetrahedron meshes. Fig. 4 shows underline for our cases, in contrast to [45]. This is clear from
the initial values of the electric and magnetic fields at the first Table 1 and the conservation property of energy, see Fig. 2.

VOLUME 7, 2019 63863

 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

FIGURE 6. The scale shows the values of magnetic fields at the final time step (N = 50), by employing the
backward Euler method. The time step size is 1t = 0.0005. (a) Snapshot of top. (b) Snapshot of bottom.
(c) Snapshot of top and front. (d) Snapshot of back and bottom. (e) Snapshot of back, left and bottom.
(f) Snapshot of right and front.

In addition, the electric and magnetic fields are visualized at VII. CONCLUSION
initial and final time steps for beam and Fichera meshes. The paper summarizes some time domain finite element
The backward Euler method is unconditionally stable methods for the system of Maxwell’s equations in three
and computationally expensive. In contrast to the backward dimensions, where the electric and magnetic fields are dis-
Euler method the symplectic method is conditionally sta- cretized by means of different finite element spaces. The
ble. We conclude that our proposed time domain finite ele- time domain mixed finite element methods have the advan-
ment method methods possess good accuracy and no spuri- tage of being substantially more powerful and reliable than
ous solutions in 3D complex geometries, and allow to treat FDTD or other existing methods with respect to error
the systems of Maxwell’s equations directly and more effi- estimates and numerical experiment, because they directly
ciently than many existing methods such as A-formulation, solve the system for electric and magnetic field intensities.
A − 8 method, operator-form, electric field formulation, Moreover the Rothe method and fully-discrete error esti-
magnetic field formation and decoupled scheme (explicit mation yield optimal solutions. To the best of our knowl-
magnetic field), for details see [1]–[15], [18], [27], [32], [34], edge, fully discrete error estimates, simulations and visu-
[36]. Moreover the proposed methods are a good basis for alizations of the type presented here, using the Nédeléc
the development of energy conserving methods for nonlin- curl-conforming and Raviart Thomas div-conforming ele-
ear problems in Optics and Photonics. The A-formulation ments with backward Euler temporal discretization for
scheme [9] also caused spurious solution in the nonlinear the system of Maxwell’s equations, were not yet avail-
case. Moreover, our proposed method could replace the able. Numerical examples are given for the fully discrete
A-formulation [9] to solve the fully non-linear system of problems.
High-Power microwave air breakdown without spurious solu- The presented symplectic time integration method is
tions. We have also obtained some first parallel results, accurate up to fourth order in time, conserves the energy
see [85]. and is conditionally stable. The backward Euler method is

63864 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

unconditionally stable. The computed electric and magnetic [4] P. Ciarlet, Jr., and L. Zou, ‘‘Fully discrete finite element approaches for
fields are visualized at intermediate and final time steps. time-dependent Maxwell’s equations,’’ Numerische Mathematik, vol. 82,
no. 2, pp. 193–219, Apr. 1999.
[5] D. A. White and M. Stowell, ‘‘Full-wave simulation of electromagnetic
APPENDIX coupling effects in RF and mixed-signal ICs using a time-domain finite-
SOME INEQUALITIES element method,’’ IEEE Trans. Microw. Theory Techn., vol. 52, no. 5,
pp. 1404–1413, May 2004.
Lemma 1: Let X be a real Hilbert space with inner product [6] W. A. Artuzi, ‘‘Improving the Newmark time integration scheme in finite
(·, ·). Then the following relations are valid for all u, v ∈ X : element time domain methods,’’ IEEE Microw. Wireless Compon. Lett.,
vol. 15, no. 12, pp. 898–900, Dec. 2005.
1) 2(u − v, u) = kuk2 + ku − vk2 − kvk2 ,
α 1 [7] M. Hochbruck and C. Stohrer, ‘‘Finite element heterogeneous multi-
2) |(u, v)| ≤ kuk2 + kvk2 for all α > 0. scale method for time-dependent Maxwell’s equations,’’ in Proc. Spectral
2 2α High Order Methods Partial Differ. Equ. (ICOSAHOM). Rio de Janeiro,
Proof: (1) follows from Brazil: Springer, Aug. 2017, pp. 269–281.
[8] A. Akbarzadeh-Sharbaf and D. D. Giannacopoulos, ‘‘A stable and effi-
kvk2 = (v − u + u, v − u + u) cient direct time integration of the vector wave equation in the finite-
element time-domain method for dispersive media,’’ IEEE Trans. Antennas
= kv − uk2 + 2(v − u, u) + kuk2 Propag., vol. 63, no. 1, pp. 314–321, Jan. 2015.
= ku − vk2 − 2(u − v, u) + kuk2 . [9] S. Yan and J.-M. Jin, ‘‘A fully coupled nonlinear scheme for
time-domain modeling of high-power microwave air breakdown,’’
(2) Obviously, 0 ≤ ku ± vk2 ≤ kuk2 ± 2(u, v) + kvk2 , hence IEEE Trans. Microw. Theory Techn., vol. 64, no. 9, pp. 2718–2729,
Sep. 2016.
2|(u, v)| ≤ kuk2 + kvk2 . √ √ [10] N. Liu, G. Cai, C. Zhu, Y. Huang, and Q. H. Liu, ‘‘The mixed finite-element
Replacing in this inequality u by αu and v by v/ α, method with mass lumping for computing optical waveguide modes,’’
the statement follows.  IEEE J. Sel. Topics Quantum Electron., vol. 22, no. 2, Mar./Apr. 2016,
Art. no. 4400709.
Lemma 2: Let δ ≥ 0, g0 ≥ 0 and (an ), (bn ), (cn ) and (γn ) [11] J. Diaz and M. J. Grote, ‘‘Multi-level explicit local time-stepping methods
be sequences of nonnegative numbers such that for second-order wave equations,’’ Comput. Methods Appl. Mech. Eng.,
vol. 291, pp. 240–265, Jul. 2015.
n n n
X X X [12] S. Ho, Y. Zhao, W. N. Fu, and P. Zhou, ‘‘Application of edge ele-
an + δ bj ≤ δ γj aj + δ cj + g0 ments to 3-D electromagnetic field analysis accounting for both inductive
j=0 j=0 j=0 and capacitive effects,’’ IEEE Trans. Magn., vol. 52, no. 3, Mar. 2016,
Art. no. 7400504.
for all n = 0, 1, 2 . . . (34) [13] J. Li, X. Lu, C. G. Farquharson, and X. Hu, ‘‘A finite-element time-
domain forward solver for electromagnetic methods with complex-
Assume that γj δ < 1 for all j, and set σj := (1 − γj δ)−1 . Then shaped loop sources,’’ Geophysics, vol. 83, no. 3, pp. E117–E132,
it holds, for all n ≥ 0: May 2018.
    [14] M. J. Grote, M. Mehlin, and T. Mitkova, ‘‘Runge–Kutta-based explicit
n n n local time-stepping methods for wave propagation,’’ SIAM J. Sci. Comput.,
X X X vol. 37, no. 2, pp. A747–A775, 2015.
an + δ bj ≤ δ cj + g0  exp δ σj γj  . [15] J. Bai, Y. Cao, X. He, H. Liu, and X. Yang, ‘‘Modeling and an immersed
j=0 j=0 j=0 finite element method for an interface wave equation,’’ Comput. Math.
Proof: Since the assumption (34) can be rewritten as Appl., vol. 76, no. 7, pp. 1625–1638, Oct. 2018.
[16] P. B. Monk, ‘‘A mixed method for approximating Maxwell’s equations,’’
n
X n−1
X n
X SIAM J. Numer. Anal., vol. 28, no. 6, pp. 1610–1634, 1991.
(1 − γn δ)an + δ bj ≤ δ γj aj + δ cj + g0 , [17] P. Monk, ‘‘An analysis of Nédélec’s method for the spatial discretization of
Maxwell’s equations,’’ J. Comput. Appl. Math., vol. 47, no. 1, pp. 101–121,
j=0 j=0 j=0
Jun. 1993.
it takes – after some further simple manipulations – a form [18] C. Makridakis and P. Monk, ‘‘Time-discrete finite element schemes for
Maxwell’s equations,’’ ESAIM: Math. Model. Numer. Anal., vol. 29, no. 2,
as assumed in [86, Lemma 2.1], and the results follows from pp. 171–197, 1995.
that Lemma.  [19] J.-M. Jin, M. Zunoubi, K. C. Donepudi, and W. C. Chew, ‘‘Frequency-
domain and time-domain finite-element solution of Maxwell’s equations
using spectral Lanczos decomposition method,’’ Comput. Methods Appl.
ACKNOWLEDGMENT Mech. Eng., vol. 169, nos. 3–4, pp. 279–296, Feb. 1999.
The authors would like to thank the anonymous review- [20] R. L. Lee and N. K. Madsen, ‘‘A mixed finite element formulation for
ers for their helpful and constructive comments that greatly Maxwell’s equations in the time domain,’’ J. Comput. Phys., vol. 88, no. 2,
pp. 284–304, Jun. 1990.
contributed to improving the final version of the paper. [21] G. Rodrigue and D. White, ‘‘A vector finite element time-domain method
We would also like to thank the Associate Editor for the for solving Maxwell’s equations on unstructured hexahedral grids,’’ SIAM
generous comments and support during the review process. J. Sci. Comput., vol. 23, no. 3, pp. 683–706, 2001.
[22] J. Li, ‘‘Error analysis of fully discrete mixed finite element schemes for
3-D Maxwell’s equations in dispersive media,’’ Comput. Methods Appl.
REFERENCES Mech. Eng., vol. 196, nos. 33–34, pp. 3081–3094, Jul. 2007.
[1] P. Monk, ‘‘Analysis of a finite element method for Maxwell’s equations,’’ [23] J. Li, ‘‘Error analysis of mixed finite element methods for wave propaga-
SIAM J. Numer. Anal., vol. 29, no. 3, pp. 714–729, Jun. 1992. tion in double negative metamaterials,’’ J. Comput. Appl. Math., vol. 209,
[2] J.-F. Lee, R. Lee, and A. Cangellaris, ‘‘Time-domain finite-element no. 1, pp. 81–96, Dec. 2007.
methods,’’ IEEE Trans. Antennas Propag., vol. 45, no. 3, pp. 430–442, [24] D. Jiao and J.-M. Jin, ‘‘A general approach for the stability analysis
May 1997. of the time-domain finite-element method for electromagnetic simula-
[3] P. Neittaanmäki and J. Saranen, ‘‘Semi-discrete Galerkin approximation tions,’’ IEEE Trans. Antennas Propag., vol. 50, no. 11, pp. 1624–1632,
method applied to initial boundary value problems for Maxwell’s equations Nov. 2002.
in anisotropic, inhomogeneous media,’’ Proc. Roy. Soc. Edinburgh Sect. A, [25] J.-F. Lee and Z. Sacks, ‘‘Whitney elements time domain (WETD) meth-
Math., vol. 89, nos. 1–2, pp. 125–133, 1981. ods,’’ IEEE Trans. Magn., vol. 31, no. 3, pp. 1325–1329, May 1995.

VOLUME 7, 2019 63865

 A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

[26] J. Zhao, ‘‘Analysis of finite element approximation for time-dependent [48] B. Wang, Z. Xie, and Z. Zhang, ‘‘Space-time discontinuous Galerkin
Maxwell problems,’’ Math. Comput., vol. 73, no. 247, pp. 1089–1105, method for Maxwell equations in dispersive media,’’ Acta Math. Scientia,
2004. vol. 34, no. 5, pp. 1357–1376, Sep. 2014.
[27] C. Ma, ‘‘Finite-element method for time-dependent Maxwell’s equations [49] M. Hochbruck and T. Pažur, ‘‘Implicit Runge–Kutta methods and discon-
based on an explicit-magnetic-field scheme,’’ J. Comput. Appl. Math., tinuous Galerkin discretizations for linear Maxwell’s equations,’’ SIAM J.
vol. 194, no. 2, pp. 409–424, Oct. 2006. Numer. Anal., vol. 53, no. 1, pp. 485–507, 2015.
[28] F. L. Teixeira, ‘‘Time-domain finite-difference and finite-element methods [50] M. Dehghan and R. Salehi, ‘‘A meshless local Petrov–Galerkin method for
for Maxwell equations in complex media,’’ IEEE Trans. Antennas Propag., the time-dependent Maxwell equations,’’ J. Comput. Appl. Math., vol. 268,
vol. 56, no. 8, pp. 2150–2166, Aug. 2008. pp. 93–110, Oct. 2014.
[29] Y. Zhang, L.-Q. Cao, and Y.-S. Wong, ‘‘Multiscale computations for 3D [51] P. Monk and E. Süli, ‘‘A convergence analysis of Yee’s scheme on
time-dependent Maxwell’s equations in composite materials,’’ SIAM J. Sci. nonuniform grids,’’ SIAM J. Numer. Anal., vol. 31, no. 2, pp. 393–412,
Comput., vol. 32, no. 5, pp. 2560–2583, 2010. 1994.
[30] J. Kim and F. L. Teixeira, ‘‘Parallel and explicit finite-element time-domain [52] T. Kang, K. Van Bockstal, and R. Wang, ‘‘The reconstruction of a time-
method for Maxwell’s equations,’’ IEEE Trans. Antennas Propag., vol. 59, dependent source from a surface measurement for full Maxwell’s equations
no. 6, pp. 2350–2356, Jun. 2011. by means of the potential field method,’’ Comput. Math. Appl., vol. 75,
[31] W. Yang, Y. Huang, and J. Li, ‘‘Developing a time-domain finite element no. 3, pp. 764–786, Feb. 2018.
method for the Lorentz metamaterial model and applications,’’ J. Sci. [53] P. Monk, ‘‘A comparison of three mixed methods for the time-dependent
Comput., vol. 68, no. 2, pp. 438–463, Aug. 2016. Maxwell’s equations,’’ SIAM J. Sci. Stat. Comput., vol. 13, no. 5,
[32] A. Anees and L. Angermann, ‘‘Mixed finite element methods pp. 1097–1122, 1992.
for the Maxwell’s equations with matrix parameters,’’ in Proc. [54] R. Hiptmair, ‘‘Finite elements in computational electromagnetism,’’ Acta
Int. Appl. Comput. Electromagn. Soc. Symp. (ACES), Mar. 2018, Numer., vol. 11, pp. 237–339, Jan. 2002.
pp. 1–2. [55] S. Labbé, ‘‘Fast computation for large magnetostatic systems adapted for
[33] Y. Huang, J. Li, and Q. Lin, ‘‘Superconvergence analysis for time- micromagnetism,’’ SIAM J. Sci. Comput., vol. 26, no. 6, pp. 2160–2175,
dependent Maxwell’s equations in metamaterials,’’ Numer. Methods Par- 2005.
tial Differ. Equ., vol. 28, no. 6, pp. 1794–1816, Nov. 2012. [56] H. Fu et al., ‘‘A parallel finite-element time-domain method for transient
[34] J. Li and Y. Lin, ‘‘A priori and posteriori error analysis for time-dependent electromagnetic simulation,’’ Geophysics, vol. 80, no. 4, pp. E213–E224,
Maxwell’s equations,’’ Comput. Methods Appl. Mech. Eng., vol. 292, Jul. 2015.
pp. 54–68, Aug. 2015. [57] L. D. Angulo, J. Alvarez, F. L. Teixeira, M. F. Pantoja,
[35] A. Anees and L. Angermann, ‘‘Time domain finite element methods for and S. G. Garcia, ‘‘A nodal continuous-discontinuous Galerkin time-
Maxwell’s equations in three dimensions,’’ in Proc. Int. Appl. Comput. domain method for Maxwell’s equations,’’ IEEE Trans. Microw. Theory
Electromagn. Soc. Symp. (ACES), Mar. 2018, pp. 1–2. Techn., vol. 63, no. 10, pp. 3081–3093, Oct. 2015.
[36] A. Anees and L. Angermann, ‘‘A mixed finite element methods approx- [58] K. S. Atia, A. Heikal, and S. S. A. Obayya, ‘‘Efficient smoothed finite
imation for the Maxwell’s equations in electromagnetics,’’ in Proc. element time domain analysis for photonic devices,’’ Opt. Express, vol. 23,
IEEE/ACES Int. Conf. Wireless Inf. Technol. Syst. (ICWITS) Appl. Comput. no. 17, pp. 22199–22213, Aug. 2015.
Electromagn. (ACES), Mar. 2016, pp. 1–2.
[59] M. C. Pinto, M. Mounier, and E. Sonnendrücker, ‘‘Handling the divergence
[37] J. Cai, J. Hong, Y. Wang, and Y. Gong, ‘‘Two energy-conserved split- constraints in Maxwell and Vlasov–Maxwell simulations,’’ Appl. Math.
ting methods for three-dimensional time-domain Maxwell’s equations Comput., vol. 272, pp. 403–419, Jan. 2016.
and the convergence analysis,’’ SIAM J. Numer. Anal., vol. 53, no. 4,
[60] Y. Shao, J. J. Yang, and M. Huang, ‘‘A review of computational electromag-
pp. 1918–1940, 2015.
netic methods for graphene modeling,’’ Int. J. Antennas Propag., vol. 2016,
[38] R. Abedi and S. Mudaliar, ‘‘An asynchronous spacetime discontinu- Apr. 2016, Art. no. 7478621.
ous Galerkin finite element method for time domain electromagnetics,’’
[61] H. Cai, X. Hu, J. Li, M. Endo, and B. Xiong, ‘‘Parallelized 3D CSEM
J. Comput. Phys., vol. 351, pp. 121–144, Dec. 2017. [Online]. Available:
modeling using edge-based finite element with total field formula-
http://www.sciencedirect.com/science/article/pii/S0021999117306514
tion and unstructured mesh,’’ Comput. Geosci., vol. 99, pp. 125–134,
[39] M. Dawson, R. Seville, and K. Morgan, ‘‘The application of a high-
Feb. 2017.
order discontinuous Galerkin time-domain method for the computation of
[62] J. Zheng, H. Huang, S. Zhang, and Z. Deng, ‘‘A general method to
electromagnetic resonant modes,’’ Appl. Math. Model., vol. 55, pp. 94–108,
simulate the electromagnetic characteristics of HTS Maglev systems by
Mar. 2018.
finite element software,’’ IEEE Trans. Appl. Supercond., vol. 28, no. 5,
[40] E. T. Chung and B. Engquist, ‘‘Optimal discontinuous Galerkin meth-
Aug. 2018, Art. no. 3600808.
ods for wave propagation,’’ SIAM J. Numer. Anal., vol. 44, no. 5,
[63] J.-M. Jin, The Finite Element Method in Electromagnetics. Hoboken, NJ,
pp. 2131–2158, 2006.
USA: Wiley, 2015.
[41] V. Dolean, S. Lanteri, and R. Perrussel, ‘‘A domain decomposition method
for solving the three-dimensional time-harmonic Maxwell equations dis- [64] J. Li and Y. Huang, Time-Domain Finite Element Methods for Maxwell’s
cretized by discontinuous Galerkin methods,’’ J. Comput. Phys., vol. 227, Equations in Metamaterials, vol. 43. New York, NY, USA: Springer, 2012.
no. 3, pp. 2044–2072, Jan. 2008. [65] P. Monk, Finite Element Methods for Maxwell’s Equations. New York, NY,
[42] H. Fahs, ‘‘High-order discontinuous Galerkin method for time-domain USA: Oxford Univ. Press, 2003.
electromagnetics on non-conforming hybrid meshes,’’ Math. Comput. [66] D. B. Davidson, Computational Electromagnetics for RF and Microwave
Simul., vol. 107, pp. 134–156, Jan. 2015. Engineering. Cambridge, U.K.: Cambridge Univ. Press, 2010.
[43] M. J. Grote and T. Mitkova, ‘‘Explicit local time-stepping methods for [67] S. R. H. Hoole, Finite Elements, Electromagnetics and Design. Amster-
MaxwellâĂŹs equations,’’ J. Comput. Appl. Math., vol. 234, no. 12, dam, The Netherlands: Elsevier, 1995.
pp. 3283–3302, Oct. 2010. [68] L. Angermann and V. V. Yatsyk, Resonant Scattering and Generation
[44] J. Li, C. Shi, and C.-W. Shu, ‘‘Optimal non-dissipative discontinuous of Waves: Cubically Polarizable Layers. Cham, Switzerland: Springer,
Galerkin methods for Maxwell’s equations in Drude metamaterials,’’ Com- 2018.
put. Math. Appl., vol. 73, no. 8, pp. 1760–1780, Apr. 2017. [69] A. Bossavit, ‘‘Whitney forms: A class of finite elements for three-
[45] C. Shi, J. Li, and C.-W. Shu, ‘‘Discontinuous Galerkin methods for dimensional computations in electromagnetism,’’ IEE A, Phys. Sci.,
Maxwell’s equations in Drude metamaterials on unstructured meshes,’’ Meas. Instrum., Manage. Educ., Rev., vol. 135, no. 8, pp. 493–500,
J. Comput. Appl. Math., vol. 342, pp. 147–163, Nov. 2018. Nov. 1988.
[46] S. Shields, J. Li, and E. A. Machorro, ‘‘Weak Galerkin methods for time- [70] A. Bossavit, ‘‘Solving Maxwell equations in a closed cavity, and the
dependent Maxwell’s equations,’’ Comput. Math. Appl., vol. 74, no. 9, question of ‘spurious modes’,’’ IEEE Trans. Magn., vol. 26, no. 2,
pp. 2106–2124, Nov. 2017. pp. 702–705, Mar. 1990.
[47] D. Sármány, M. A. Botchev, and J. J. van der Vegt, ‘‘Time- [71] J. C. Nédélec, ‘‘Mixed finite elements in R3 ,’’ Numerische Mathematik,
integration methods for finite element discretisations of the second-order vol. 35, no. 3, pp. 315–341, Sep. 1980.
Maxwell equation,’’ Comput. Math. Appl., vol. 65, no. 3, pp. 528–543, [72] J. C. Nédélec, ‘‘A new family of mixed finite elements in R3 ,’’ Numerische
Feb. 2013. Mathematik, vol. 50, no. 1, pp. 57–81, Jan. 1986.

63866 VOLUME 7, 2019

A. Anees, L. Angermann: TDFEM for Maxwell’s Equations

[73] P.-A. Raviart and J. M. Thomas, ‘‘A mixed finite element method for ASAD ANEES received the M.Sc. degree in
2-nd order elliptic problems,’’ in Mathematical Aspects of Finite Element applied mathematics from the University of
Methods. Berlin, Germany: Springer, 1977, pp. 292–315. Engineering and Technology at Lahore, Lahore,
[74] R. Leis, Initial Boundary Value Problems in Mathematical Physics. Pakistan, in 2009. He is currently pursuing the
Hoboken, NJ, USA: Wiley, 1988. Ph.D. degree in applied mathematics with the
[75] J. Lions and G. Duvaut, Inequalities in Mechanics and Physics. Berlin, Institute of Mathematics, Clausthal University of
Germany: Springer, 1976. Technology, Germany, where he is also a Member
[76] A. Kufner, O. John, and S. Fucik, Function spaces. Prague, Czech Repub-
of a research project.
lic: Academic, 1977.
Since 2011, he has been a tenured Lecturer in
[77] R. Adams and J. Fournier, ‘‘Sobolev spaces,’’ in Pure and Applied Math-
ematics (Amsterdam), vol. 140, 2nd ed. Amsterdam, The Netherlands: mathematics with the Department of Mathematics
Elsevier, 2003. and Statistics, University of Agriculture Faisalabad, Punjab, Pakistan. He has
[78] V. Girault and P. Raviart, Finite Element Methods for Navier-Stokes Equa- been having two positions, since 2014. His current research interests include
tions: Theory and Algorithms. New York, NY, USA: Springer-Verlag, high performance computing, space-time parallel method, error estimation,
1986. and scientific visualization for electromagnetic and nonlinear optics. In par-
[79] P. G. Ciarlet, The Finite Element Method for Elliptic Problems. Philadel- ticular, his expertise includes the development of theoretical aspects of higher
phia, PA, USA: SIAM, 2002. order spatial, temporal discretization, and energy conserving time domain
[80] B. G. Pachpatte, ‘‘On a certain inequality arising in the theory of dif- finite element methods (TDFEM) with nonlinear material in computational
ferential equations,’’ J. Math. Anal. Appl., vol. 182, no. 1, pp. 143–157, electromagnetic and optics.
Feb. 1994.
[81] A. Fisher, D. White, and G. Rodrigue, ‘‘An efficient vector finite ele-
ment method for nonlinear electromagnetic modeling,’’ J. Comput. Phys.,
vol. 225, no. 2, pp. 1331–1346, 2007. LUTZ ANGERMANN received the Ph.D. degree
[82] R. Rieben, D. White, and G. Rodrigue, ‘‘High-order symplectic integration from the University of Technology at Dresden,
methods for finite element solutions to time dependent Maxwell equa- in 1987, and the Habilitation degree from the Uni-
tions,’’ IEEE Trans. Antennas Propag., vol. 52, no. 8, pp. 2190–2195, versity of Erlangen-Nürnberg, in 1995. He studied
Aug. 2004. Mathematics with the State University of Kharkiv
[83] A. Christophe, S. Descombes, and S. Lanteri, ‘‘An implicit hybridized dis- (now V. N. Karazin Kharkiv National University,
continuous Galerkin method for the 3D time-domain Maxwell equations,’’ Ukraine).
Appl. Math. Comput., vol. 319, pp. 395–408, Feb. 2015.
From 1998 to 2001, he was an Associate Profes-
[84] W. E. Arnoldi, ‘‘The principle of minimized iterations in the solution of the
matrix eigenvalue problem,’’ Quart. Appl. Math., vol. 9, no. 1, pp. 17–29,
sor of numerical mathematics with the University
1951. of Magdeburg. Since 2001, he has been a Professor
[85] A. Anees and L. Angermann, ‘‘Time domain finite element meth- of numerical mathematics with the Department of Mathematics, Clausthal
ods for Maxwell’s equations,’’ in Proc. 7th Workshop Parallel-in- University of Technology. He has authored or coauthored about 120 scientific
Time Methods, Roscoff Marine station. Roscoff, France, May 2018, papers and four books among them. His research interest includes the math-
pp. 1–32. [Online]. Available: http://www.math.univ-paris13.fr/~japhet/ ematical analysis of numerical algorithms for partial differential equations
PINT/Slides/Asad_ANEES_PinT2018.pdf with special interests in finite volume and finite element methods and their
[86] A. Quarteroni and A. Valli, Numerical Methods for Partial Differential application to problems in physics and engineering.
Equations. New York, NY, USA: Springer, 1994.

VOLUME 7, 2019 63867