# Deep Fourier Residual method for solving time-harmonic Maxwell’s equations

Jamie M. Taylor  CUNEF University, Madrid, Spain.
<EMAIL_ADDRESS>Manuela Bastidas  University of the Basque Country
(UPV/EHU), Leioa, Spain.  David Pardo  University of the Basque Country
(UPV/EHU), Leioa, Spain.  Basque Centre for Applied Mathematics (BCAM),
Bilbao, Spain.  Ikerbasque (Basque Foundation For Sciences), Bilbao, Spain.
Ignacio Muga  Pontificia Universidad Católica de Valparaíso, Chile.

###### Abstract

Solving PDEs with machine learning techniques has become a popular alternative
to conventional methods. In this context, Neural networks (NNs) are among the
most commonly used machine learning tools, and in those models, the choice of
an appropriate loss function is critical. In general, the main goal is to
guarantee that minimizing the loss during training translates to minimizing
the error in the solution at the same rate. In this work, we focus on the
time-harmonic Maxwell’s equations, whose weak formulation takes
${H}_{0}(\mathrm{curl},\Omega)$ as the space of test functions. We propose a
NN in which the loss function is a computable approximation of the dual norm
of the weak-form PDE residual. To that end, we employ the Helmholtz
decomposition of the space ${H}_{0}(\mathrm{curl},\Omega)$ and construct an
orthonormal basis for this space in two and three spatial dimensions. Here, we
use the Discrete Sine/Cosine Transform to accurately and efficiently compute
the discrete version of our proposed loss function. Moreover, in the numerical
examples we show a high correlation between the proposed loss function and the
$H$(curl)-norm of the error, even in problems with low-regularity solutions.

##  1 Introduction

Partial Differential Equations (PDEs) are essential tools for modeling and
simulating various scientific and industrial problems, in particular, they
form the backbone of modern physics. Herein, we concentrate our attention on
Maxwell’s equations. This set of equations describes the propagation of
electromagnetic waves through different media. Commonly used methods to
calculate the solution of Maxwell’s problems range from exact methods [19, 21]
to numerical approximations, such as Finite Differences, Finite Elements
(FEM), and Discontinuous Galerkin [8, 9, 26, 29, 36]. Nevertheless, exact
solutions can only be obtained in rare scenarios, and proposing a proper
numerical technique that is conformal, accurate, and efficient is often
challenging.

Using machine learning techniques to approximate the solution of Maxwell’s
equations is an attractive alternative to classical approaches. Some of the
most used machine learning tools in the context of PDEs are Neural Networks
(NNs) [6, 13, 27, 28, 30, 33]. These architectures have shown promising
results when tackling complex nonlinear systems of equations that underlay
physical phenomena. For instance, we highlight their potential to solve
parametric PDEs [20, 25], enhance classical numerical methods [2] and solve
problems in presence of singularities or sharp gradients [7, 35].

Some popular methods like PINNs (Physics-Informed Neural Networks) [22, 30]
enforce a NN to satisfy the strong formulation of a PDE by implementing a
numerically tractable norm of the strong-form residual as a loss function.
Even though PINNs can approximate the solutions of many physical problems,
they are inaccurate when the weak solution does not satisfy the strong-form
equation. As a result, PINNs are severely constrained in many applications
that naturally produce low-regularity solutions. Maxwell’s equations, for
example, may admit a solution in the functional space
${H}(\mathrm{curl},\Omega)\setminus H^{2}(\Omega)$ (or even in
${H}(\mathrm{curl},\Omega)\setminus H^{1}(\Omega)$ as shown in, e.g., [32])
since the smoothness of the solution depends on the regularity of the domain,
the sources, and the boundary conditions. In these situations, PINNs may not
be applicable.

In contrast, VPINNs (Variational Physics-Informed Neural Networks) [23, 24]
use the residual of the weak formulation of a PDE in the loss function. In
this approach, it is essential to select an appropriate set of test functions
and define a computable loss that controls the error of the solution, which is
generally non-trivial.

An ideal loss function would be the energy norm of the error. Since the error
function is unavailable in practice (one would need the exact solution), one
typically resorts to minimizing the dual norm of the weak-form PDE residual,
that is $\mathcal{R}:H\to H^{\prime}$, where $H$ is a Hilbert space of test
functions and $H^{\prime}$ its dual. In many practical cases of linear PDEs,
if $u^{*}$ is the exact solution of a PDE, then, there exist constants
$0<\gamma<M$, such that

$\frac{1}{M}\|\mathcal{R}(u)\|_{H^{\prime}}\leq\|u-u^{*}\|_{H}\leq\frac{1}{\gamma}\|\mathcal{R}(u)\|_{H^{\prime}},$
(1.1)

where $u$ is an approximation to the solution $u^{*}$ [37]. It is clear from
(1.1) that the dual norm of the residual $\|\mathcal{R}(u)\|_{H^{\prime}}$ is
equivalent to the energy norm of the error, and so it is an appropriate choice
of loss function, i.e., $\mathcal{L}(u):=\|\mathcal{R}(u)\|_{H^{\prime}}$.
Unfortunately, it is in general very challenging to evaluate
$\|\cdot\|_{H^{\prime}}$.

For similar problems based on $H^{1}$ test function spaces, the authors of
[37] proposed the Deep Fourier Residual (DFR) method. They employed a
numerical method for approximating the dual norm of residuals corresponding to
PDEs with $H^{1}$ test function spaces via a spectral representation of the
dual norm, which may be implemented using the Fast Fourier Transform (FFT).

This paper follows the same spirit; we extend the ideas of [37] and
numerically implement the dual norm as a loss function to solve PDEs with
${H}_{0}(\mathrm{curl},\Omega)$ test function spaces, motivated by the time-
harmonic Maxwell’s equations. Similarly to the DFR method for PDEs with
$H^{1}$ test function spaces, here the key challenge is to construct an
appropriate orthonormal basis for ${H}_{0}(\mathrm{curl},\Omega)$. To find
such a basis, we use the classical Helmholtz decomposition of a vector field
[29] and construct a complete set of basis functions for the space
${H}_{0}(\mathrm{curl},\Omega)$ on product domains in two and three spatial
dimensions. The strength of the DFR method lies in solving problems with low
regularity, where PINNs-like methods based on strong formulations fail.
Moreover, our choice of loss function bounds the $H$(curl)-norm of the error
of the solution. In this way, minimizing the loss during training implies the
reduction of the error solution at the same rate. We provide numerical
examples that demonstrate the correlation between the proposed loss function
and the $H$(curl)-norm of the error. In the specific case of Maxwell’s
equations, the DFR method produces accurate results on heterogeneous and
discontinuous media, as well as strong correlations between the loss and
$H$(curl)-norm of the error during training.

Notice that, besides our work, other studies employ the combination of Fourier
basis functions and machine learning techniques. This is usually called
Physics-informed Spectral Learning (PiSL) [15, 16]. Nevertheless, these
techniques are associated with data analysis and in contrast to our study,
they are situated within the framework of physics-informed statistical
learning.

Our proposed DFR method faces significant challenges that are similar to those
encountered in [37]. Here, we only consider product domains, to profit from
the results in [10], and with Dirichlet-type boundary conditions. Different
strategies must be used to tackle problems with non-trivial geometries and
general boundary conditions. Specifically, using the DFR method when the
problem involves general geometries is challenging since constructing an
orthonormal basis for the space ${H}_{0}(\mathrm{curl},\Omega)$ is not
straightforward. Another important limitation of our technique is that the
constants bounding the norm of the error, in (1.1), diverge in certain cases,
e.g., when one considers materials with a large variation in parameters or
frequencies close to a resonant frequency of the system. These technical
complications are seen as stability issues intrinsic to the PDE similar to
those encountered in traditional approaches such as FEM [12], where the
accuracy of the error estimates deteriorate under the same conditions.

The remainder of this work is organized as follows. In Section 2, we state
Maxwell’s equations and the weak formulation of the problem. There, we also
motivate the choice of the dual norm of the residual operator as a natural
loss function for training NNs in the case of Maxwell’s equations. Later, we
discuss the relevance of constructing an appropriate orthonormal basis for the
space ${H}_{0}(\mathrm{curl},\Omega)$. The details of generating such an
orthonormal basis, using Helmholtz decomposition in two and three spatial
dimensions, are described in Section 3, with the technical details deferred to
Appendix A. In Section 3, we also present the test functions in 2D and 3D for
the simple case of squared and cubic domains and the corresponding
calculations are presented in the Appendix B. In Section 4, we describe the
architecture of a NN and explain the core of the DFR method, which consists of
constructing a computable discretized loss function. Finally, in Sections 5
and 6, we present numerical experiments, conclusions, and directions for
future research.

##  2 Problem statement

Consider a domain $\Omega\subset{\mathbb{R}}^{n}$, with $n=2$ or $3$, whose
boundary $\Gamma:=\partial\Omega$ is polyhedral and connected. Given an
impressed field $\bm{E}^{I}$, and an electric density current source $\bm{J}$,
we look for an electric field $\bm{E}$, and a magnetic field $\bm{H}$ solving
the so-called macroscopic linear Maxwell’s equations in a time-harmonic form

$\displaystyle\mathrm{curl}(\bm{E})-i\omega\mu\bm{H}$ $\displaystyle=0$
$\displaystyle\text{in }\Omega,\qquad\text{(Faraday's Law)},$ (2.1)
$\displaystyle\mathrm{curl}(\bm{H})+i\omega\epsilon\bm{E}$
$\displaystyle=\bm{J}$ $\displaystyle\text{in }\Omega,\qquad\text{(Ampere's
Law)},$ $\displaystyle\bm{E}\times\bm{\mathrm{n}}$ $\displaystyle=\bm{E}^{I}$
$\displaystyle\text{on }\Gamma,$

where $i$ is the imaginary unit, $\bm{\mathrm{n}}$ denotes the outward unit
normal vector of $\Omega$, $\omega\in{\mathbb{R}}$ is the angular frequency,
and $\mu$ and $\epsilon$ are space-dependent functions standing for the
magnetic permeability and electrical permittivity, respectively. These
functions generally may be tensor-valued, but we consider here only the scalar
case.

The formulation (2.1) arises by considering the time-dependent Maxwell’s
equations under the following Ansatz on the electric and magnetic fields

$\bm{E}(x,t)=\mathfrak{R}(e^{i\omega t}\bm{E}(x)),\text{ and
}\bm{H}(x,t)=\mathfrak{R}(e^{i\omega t}\bm{H}(x)).$

In (2.1), the dependency of the time $t>0$ is implicit. The model (2.1) is
completed with the following Gauss’ Laws for magnetic and electric fields

$\displaystyle\mathrm{div}(\mu\bm{H})$ $\displaystyle=0$
$\displaystyle\text{in }\Omega,$ (2.2)
$\displaystyle\mathrm{div}(\epsilon\bm{E})$ $\displaystyle=\rho$
$\displaystyle\text{in }\Omega,$

where $\rho$ is the density of free charge.

For the sake of simplicity, in this work we only consider the case when
$\omega$ is non-zero. Notice that when $\omega\neq 0$ the equations in (2.2)
are consequences of (2.1) and the continuity equation

$i\omega\rho+\mathrm{div}(\bm{J})=0,$

which relates the rate of change of the charge density to the divergence of
the current density.

We point out that we have only considered Dirichlet-type boundary conditions.
Other strategies must be implemented for different boundary conditions. Whilst
the following discussion will only be limited to the case of homogeneous
boundary conditions, as our approach is based on an analysis on the space of
test functions, rather than trial functions, the extension to non-homogeneous
boundary conditions is straightforward as the test function space remains
unchanged.

###  2.1 Preliminaries

We start by defining relevant operators and function spaces. Herein, the
definitions are regarded as classical and can be found in [18]. First, we
denote $\mathcal{D}(\Omega)$ the space of smooth functions with compact
support in $\Omega$ and $\mathcal{D}^{\prime}(\Omega)$ is the space of
distributions. We let $L^{p}(\Omega)$ be the space of the $p$-integrable real-
valued functions equipped with the usual norm.

For smooth functions, we define the divergence, gradient and curl according to
their usual definitions (see [34]). In 2D, we define the curl and its adjoint
as

$\displaystyle\mathrm{curl}:$
$\displaystyle\quad[\mathcal{D}(\Omega)]^{2}\ni\bm{\phi}\quad\mapsto\quad\partial_{y}\bm{\phi}_{1}-\partial_{x}\bm{\phi}_{2}\in\mathcal{D}(\Omega),$
(2.3) $\displaystyle\mathrm{curl}^{*}:$
$\displaystyle\quad\mathcal{D}(\Omega)\ni\phi\quad\mapsto\quad\begin{pmatrix}-\partial_{y}\phi\\\
\partial_{x}\phi\end{pmatrix}\in[\mathcal{D}(\Omega)]^{2}.$

The curl operator in 2D and 3D can then be extended as mappings from
appropriate $L^{2}$ spaces to $\mathcal{D}^{\prime}(\Omega)$ by duality. We
then define the function space ${H}(\mathrm{curl},\Omega)$ consisting of
functions in $[L^{2}(\Omega)]^{n}$ whose curl, interpreted in the sense of
distributions, is in $[L^{2}(\Omega)]^{n^{\prime}}$, i.e.,

${H}(\mathrm{curl},\Omega):=\\{\bm{u}\in[L^{2}(\Omega)]^{n}:\mathrm{curl}(\bm{u})\in[L^{2}(\Omega)]^{n^{\prime}}\\},$

where $n^{\prime}=1$ if $n=2$ and $n^{\prime}=3$ if $n=3$. The space
${H}(\mathrm{curl},\Omega)$ is a Hilbert space with inner product given by

$(\bm{u},\bm{v})_{{H}(\mathrm{curl},\Omega)}=\int_{\Omega}\mathrm{curl}(\bm{u})\cdot\mathrm{curl}(\bm{v})+\bm{u}\cdot\bm{v}\,d{\bm{x}}\qquad\forall\bm{u},\bm{v}\in{H}(\mathrm{curl},\Omega).$
(2.4)

For any bounded Lipschitz domain $\Omega\subset{\mathbb{R}}^{n}$ with boundary
$\Gamma$ and outward normal $\bm{\mathrm{n}}$, the mapping
$\gamma_{t}:\mathcal{C}^{1}(\bar{\Omega})\to L^{2}(\Gamma)$ with
$\gamma_{t}(\bm{u})=\bm{u}|_{\Gamma}\times\bm{\mathrm{n}}$ can be uniquely
extended to the continuous tangential trace operator,
$\gamma_{t}:{H}(\mathrm{curl},\Omega)\to
H^{-\frac{1}{2}}(\Gamma,{\mathbb{R}}^{d^{\prime}})$ (see [5])111In two
dimensions, the cross product is interpreted as the scalar product
$\bm{u}\times\bm{v}=\bm{v}_{2}\bm{u}_{1}-\bm{u}_{2}\bm{v}_{1}$.. We then
define the space

$\begin{split}{H}_{0}(\mathrm{curl},\Omega):=&\\{\bm{u}\in{H}(\mathrm{curl},\Omega):\,\gamma_{t}(\bm{u})=0\\}.\end{split}$

###  2.2 Weak formulation

There are multiple weak formulations for the Maxwell system, all based on the
general idea of minimizing the functional that represents the energy of an
electromagnetic field. Notice that in (2.1) one can eliminate $\bm{H}$ or
$\bm{E}$ from each of the equations. Assuming $\mu$ and $\epsilon$ are real-
valued, bounded and non-zero functions, $\bm{J}\in[L^{2}(\Omega)]^{n}$, and
$\bm{E}^{I}=\mathbf{0}$, the weak formulation corresponding to the electric
field in the problem (2.1) is: Find $\bm{E}\in{H}_{0}(\mathrm{curl},\Omega)$
satisfying

$\int_{\Omega}\mu^{-1}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\bm{\phi})-\omega^{2}\epsilon\bm{E}\cdot\bm{\phi}\,d{\bm{x}}=\int_{\Omega}i\omega\bm{J}\cdot\bm{\phi}\,d{\bm{x}}\qquad\forall\bm{\phi}\in{H}_{0}(\mathrm{curl},\Omega).$
(2.5)

An analogous weak form exists for the magnetic field $\bm{H}$. Notice that,
Gauss’ Laws in (2.2) are satisfied weakly, by considering test functions
$\bm{\phi}=\nabla u$ for $u\in H^{1}_{0}(\Omega)$ in (2.5).

###  2.3 Residual minimization

The residual operator corresponding with the weak form (2.5) is
$\mathcal{R}:H\to H^{\prime}$ with $H={H}_{0}(\mathrm{curl},\Omega)$ and
$H^{\prime}$ being its dual. This weak residual operator may be expressed in
the general form

$\langle\mathcal{R}(\bm{E}),\bm{\phi}\rangle_{H^{\prime}\times
H}=b(\bm{E},\bm{\phi})-\ell(\bm{\phi}),$ (2.6)

where $\ell\in H^{\prime}$ is

$\ell(\bm{\phi})=\int_{\Omega}i\omega\bm{J}\cdot\bm{\phi}\,d{\bm{x}}\qquad\forall\bm{\phi}\in
H,$

and $b:H\times H\to{\mathbb{R}}$ is the following bilinear form

$b(\bm{E},\bm{\phi})=\int_{\Omega}\mu^{-1}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\bm{\phi})-\omega^{2}\epsilon\bm{E}\cdot\bm{\phi}\,d{\bm{x}}\qquad\forall\bm{\phi}\in
H.$ (2.7)

The existence and uniqueness of a solution of (2.5), for $\omega$ outside of a
countable set of resonant frequencies, is proved in, e.g. [29, Chapter 4] and
[26, Theorem 4.32]. Moreover, using the reasoning in [14, Section 25.3], we
know that the solution for the variational problem exists and is unique if and
only if the following bounds apply

$\displaystyle\|\mathcal{R}(\bm{E})\|_{H^{\prime}}^{2}$
$\displaystyle=\sup_{\bm{\phi}\in
H\backslash\\{0\\}}\frac{|\langle\mathcal{R}(\bm{E}),\bm{\phi}\rangle_{H^{\prime}\times
H}|}{\|\bm{\phi}\|_{H}}$ $\displaystyle\leq M\|\bm{E}-\bm{E}^{*}\|_{H},$ (2.8)
$\displaystyle\|\mathcal{R}(\bm{E})\|_{H^{\prime}}^{2}$
$\displaystyle=\sup_{\bm{\phi}\in
H\backslash\\{0\\}}\frac{|\langle\mathcal{R}(\bm{E}),\bm{\phi}\rangle_{H^{\prime}\times
H}|}{\|\bm{\phi}\|_{H}}$ $\displaystyle\geq\gamma\|\bm{E}-\bm{E}^{*}\|_{H},$

where $\gamma$ and $M$ are positive constants depending on $\mu$, $\omega$ and
$\epsilon$, and $\bm{E}^{*}$ denotes the exact solution of (2.5). We emphasize
that the coercive case of $\epsilon<0$, is mathematically interesting, as it
implies that the bilinear form (2.7) becomes equivalent to the inner product
on ${H}(\mathrm{curl},\Omega)$.

On the other hand, for any $\mathcal{F}\in H^{\prime}$, by the Riesz
representation theorem, there exists some $\bm{u}_{\mathcal{F}}\in H$ with
$\mathcal{F}(\bm{v})=(\bm{u}_{\mathcal{F}},\bm{v})_{H}$ for all $\bm{v}\in H$
and $\|{\mathcal{F}}\|_{H^{\prime}}=\|\bm{u}_{\mathcal{F}}\|_{H}$.
Furthermore, if $(\Phi_{k})_{k\in\mathcal{I}}$ is an orthonormal basis of $H$,
with $\mathcal{I}$ denoting a set of indices, then, by using the generalized
Parseval’s identity, we have that the dual norm of any ${\mathcal{F}}\in
H^{\prime}$ can be expressed as

$\|{\mathcal{F}}\|_{H^{\prime}}^{2}=\|\bm{u}_{\mathcal{F}}\|_{H}^{2}\stackrel{{\scriptstyle\text{Parseval}}}{{=}}\sum\limits_{k\in\mathcal{I}}(\bm{u}_{\mathcal{F}},\Phi_{k})_{H}^{2}=\sum\limits_{k\in\mathcal{I}}{\mathcal{F}}(\Phi_{k})^{2}.$
(2.9)

From (2.6) and (2.9) we obtain an expression for the dual norm of the residual
as

$\|\mathcal{R}(\bm{E})\|_{H^{\prime}}^{2}=\sum\limits_{k\in\mathcal{I}}\langle\mathcal{R}(\bm{E}),\Phi_{k}\rangle_{H^{\prime}\times
H}^{2},$ (2.10)

According to (2.10), determining a set of orthonormal basis functions
$\Phi_{k}$ is all that is required to calculate the dual norm of the residual.
Whilst this is generally a non-trivial task, in the following section, we will
construct such a set of basis functions in simplified geometries.

##  3 A set of basis functions for ${H}_{0}(\mathrm{curl},\Omega)$

In order to find basis functions for $H_{0}(\mathrm{curl},\Omega)$, we seek an
orthonormal eigenbasis for the differential operator corresponding to the
inner product (2.4). More details of these ideas can be found in [26, Chapter
4]. The inner product (2.4) is naturally associated with the differential
operator $(1+\mathrm{curl}\text{-}\mathrm{curl})$. So, we consider the weak
eigenpairs
$(\lambda_{k},\Phi_{k})\in{\mathbb{R}}\times{H}_{0}(\mathrm{curl},\Omega)$
solving

$\displaystyle\int_{\Omega}\mathrm{curl}(\Phi_{k})\cdot\mathrm{curl}(\bm{v})+\Phi_{k}\cdot\bm{v}\,d{\bm{x}}$
$\displaystyle=\lambda_{k}\int_{\Omega}\Phi_{k}\cdot\bm{v}\,d{\bm{x}}\qquad\forall\bm{v}\in{H}_{0}(\mathrm{curl},\Omega),$
(3.1) $\displaystyle\|\Phi_{k}\|_{{H}(\mathrm{curl},\Omega)}$
$\displaystyle=1.$

In strong form, we may write

$\displaystyle(1+\mathrm{curl}\text{-}\mathrm{curl})\Phi_{k}$
$\displaystyle=\lambda_{k}\Phi_{k}$ $\displaystyle\text{ in }\Omega,$ (3.2)
$\displaystyle\Phi_{k}\times\bm{\mathrm{n}}$ $\displaystyle=0$
$\displaystyle\text{ on }\Gamma.$

If the inverse of $(1+\mathrm{curl}\text{-}\mathrm{curl})$ were to be compact
and self-adjoint on ${H}_{0}(\mathrm{curl},\Omega)$, then, the application of
the Hilbert–Schmidt theorem to the inverse would provide the existence of an
eigenbasis for the operator itself and a complete eigenbasis of
${H}_{0}(\mathrm{curl},\Omega)$. Nevertheless, we notice that, when restricted
to functions of the form $\bm{v}=\nabla u$ for $u\in H_{0}^{1}(\Omega)$, the
inverse of the operator $(1+\mathrm{curl}\text{-}\mathrm{curl})$ is the
identity and thus not compact. However, this technical issue can be resolved
using the Helmholtz decomposition of the space
${H}_{0}(\mathrm{curl},\Omega)$, i.e., we decompose the problem (3.1) into two
sub-problems. We omit the specifics of this reasoning and refer to [3, Chapter
6] and [11] for more details.

From (3.1), one obtains that the eigenvalue $\lambda_{k}=1$ has an infinite-
dimensional eigenspace, that is, the null space of the $\mathrm{curl}$
operator. In simply connected domains, we have that $\mathrm{curl}(\bm{v})=0$
implies $\bm{v}=\nabla u$ for some $u\in H^{1}(\Omega)$. Moreover, if $u\in
H_{0}^{1}(\Omega)$, then $\nabla u$ is parallel to the unit normal vector
$\bm{\mathrm{n}}$ on $\partial\Omega$, which means that $\nabla
u\in{H}_{0}(\mathrm{curl},\Omega)$. Thus, we identify a large space of
eigenvectors with eigenvalue $\lambda_{k}=1$, $\nabla
H_{0}^{1}(\Omega)\subset{H}_{0}(\mathrm{curl},\Omega)$ defined as

$\nabla H^{1}_{0}(\Omega):=\\{\nabla u:u\in H^{1}_{0}(\Omega)\\}.$

We note that equipping $H^{1}_{0}(\Omega)$ with the inner product
$(u,v)_{H^{1}_{0}(\Omega)}=\int_{\Omega}\nabla u\cdot\nabla v\,d{\bm{x}}$, the
space $\nabla H_{0}^{1}(\Omega)$, as a subspace of
${H}_{0}(\mathrm{curl},\Omega)$, is isometric to $H_{0}^{1}(\Omega)$, i.e.,
$(\nabla u,\nabla v)_{{H}(\mathrm{curl},\Omega)}=(u,v)_{H_{0}^{1}(\Omega)}$.
Consequently, the space $\nabla H^{1}_{0}(\Omega)$ forms a closed subspace of
${H}_{0}(\mathrm{curl},\Omega)$ (see [26, Lemma 4.20]), and we can employ the
following orthogonal decomposition of ${H}_{0}(\mathrm{curl},\Omega)$

${H}_{0}(\mathrm{curl},\Omega)=X_{0}(\Omega)\oplus\nabla H^{1}_{0}(\Omega),$
(3.3)

where $X_{0}(\Omega):=(\nabla H^{1}_{0}(\Omega))^{\perp}$. Notice that for a
function $\bm{v}\in X_{0}(\Omega)$, one necessarily has that

$(\bm{v},\nabla u)_{{H}(\mathrm{curl},\Omega)}=0\qquad\forall u\in
H^{1}_{0}(\Omega),$

meaning $\mathrm{div}(\bm{v})=0$ weakly. That is to say that vector fields in
${H}_{0}(\mathrm{curl},\Omega)$ can be decomposed into two parts: a curl-free
component and a divergence-free component. Finding an orthonormal basis for
$H_{0}(\mathrm{curl},\Omega)$ reduces to finding an orthonormal basis for each
component. We consider each of these in turn.

###  3.1 A set of basis functions for $\nabla H^{1}_{0}(\Omega)$

To find a basis for the space $\nabla H^{1}_{0}(\Omega)$, we use the fact that
the differential operator $\nabla$ defines an isometry between
$H^{1}_{0}(\Omega)$ and $\nabla H^{1}_{0}(\Omega)$, viewed as a subset of
$H_{0}(\mathrm{curl},\Omega)$. The gradients of any orthonormal basis of
$H^{1}_{0}(\Omega)$ thus define an orthonormal basis of $\nabla
H^{1}_{0}(\Omega)\subset{H}_{0}(\mathrm{curl},\Omega)$.

From classical spectral theory, the following proposition defines an
orthonormal basis for $\nabla H^{1}_{0}(\Omega)$ in terms of the homogeneous-
Dirichlet eigenvectors of $-\Delta$ in $\Omega$.

######  Proposition 3.1.

Let $\Omega$ be a bounded, simply connected, and Lipschitz domain in
$\mathbb{R}^{n}$, where $n=2$ or $3$. There exists an orthonormal basis of
$H^{1}(\Omega)$, consisting of non-zero homogeneous-Dirichlet eigenvectors of
$-\Delta$ in $\Omega$, $(\phi_{k})_{k\in\mathcal{I}}$, for a countable index
set $\mathcal{I}$. Then, the sequence $(\bm{\phi}_{k})_{k\in\mathcal{I}}$,
defined as

$\bm{\phi}_{k}:=\frac{\nabla\phi_{k}}{\|\nabla\phi_{k}\|_{[L^{2}(\Omega)]^{n}}},$

forms an orthonormal basis for $\nabla
H^{1}_{0}(\Omega)\subset{H}_{0}(\mathrm{curl},\Omega)$.

Appendix A contains the proof of Proposition 3.1.

###  3.2 A set of basis functions for $X_{0}(\Omega)$

For constructing the eigenbasis for the space $X_{0}(\Omega)$, we resume the
prior discussion on the eigenvectors of the operator
$(1+\mathrm{curl}\text{-}\mathrm{curl})$. We recall that $X_{0}(\Omega)$ is
compactly embedded into $L^{2}(\Omega)$ (see [26, Theorem 4.23]), which
ensures that the $(1+\mathrm{curl}\text{-}\mathrm{curl})$ operator admits a
compact and self-adjoint inverse when restricted to $X_{0}(\Omega)$. Then, the
eigenbasis of the operator forms a complete and orthonormal basis of the space
$X_{0}(\Omega)$. Moreover, when restricted to divergence-free vector fields,
the $(1+\mathrm{curl}\text{-}\mathrm{curl})$ operator reduces to $(1-\Delta)$,
which is suggestive of the fact that one may construct eigenvectors of
$(1+\mathrm{curl}\text{-}\mathrm{curl})$ via eigenvectors of the negative
Laplacian, which we perform in the following propositions. Here, we divide the
construction of this eigenbasis into two scenarios based on the dimensionality
of $\Omega$.

######  Proposition 3.2.

Let $\Omega$ be a bounded, simply connected, and Lipschitz domain in
$\mathbb{R}^{2}$. There exists an orthonormal basis of $H^{1}(\Omega)$,
consisting of non-constant homogeneous-Neumann eigenvectors of $-\Delta$ in
$\Omega$, $(\bm{\psi}_{k})_{k\in\mathcal{I}}$, for a countable index set
$\mathcal{I}$. Then, the sequence $(\bm{\psi}_{k})_{k\in\mathcal{I}}$, defined
as

$\bm{\psi}_{k}:=\frac{\mathrm{curl}^{*}(\phi_{k})}{\|\mathrm{curl}^{*}(\phi_{k})\|_{{H}(\mathrm{curl},\Omega)}},$

where $\mathrm{curl}^{*}$ is as defined in (2.3), forms an orthonormal basis
for $X_{0}(\Omega)\subset{H}_{0}(\mathrm{curl},\Omega)$.

The proof of Proposition 3.2 is detailed in Appendix A.

Finding an orthonormal basis for the space $X_{0}(\Omega)$ in 3D is more
complex than in the 2D case. Here, we restrict the construction of basis
functions for the space $X_{0}(\Omega)$ to Cartesian product domains, i.e.,
domains defined as the Cartesian product of a simply connected domain in
${\mathbb{R}}^{2}$ and a closed interval in ${\mathbb{R}}$. In this particular
case, we can construct eigenvectors of the
$(1+\mathrm{curl}\text{-}\mathrm{curl})$ operator in a similar fashion,
according to differential operators acting upon scalar-valued eigenfunctions
of the Laplacian with appropriate boundary conditions.

In such geometries, the basis functions of $X_{0}(\Omega)$ in three dimensions
come in two distinct modes, called TM (Transverse Magnetic) and TE (Transverse
Electric) modes [26]. The construction of such an eigenbasis was developed in
[10]. There, the authors demonstrate that the TE and TM modes defined below
constitute a complete and orthonormal basis of $X_{0}(\Omega)$.

Without loss of generality, we state that a Cartesian product domain
$\Omega\in{\mathbb{R}}^{3}$ is of the form $\Omega^{*}\times I$, with
$\Omega^{*}\subset{\mathbb{R}}^{2}$ simply connected and
$I\subset{\mathbb{R}}$ an interval. We refer to the coordinate direction
corresponding to the interval $I$ as the distinguished direction. For this
specific definition of $\Omega$, the $z$-coordinate is the distinguished axis
of the Cartesian product domain, but it is worth noting that the following
derivation is direction-independent.

Considering the scope of this work, we refer to [10] for a more detailed
explanation and in Appendix A we give further details of the properties of the
basis functions for $X_{0}(\Omega)$ in 3D.

######  Proposition 3.3.

Let $\Omega\subset{\mathbb{R}}^{2}$ be a Cartesian product domain such that
$\Omega=\Omega^{*}\times I$, where $\Omega^{*}\subset\mathbb{R}^{2}$ is a
simply connected domain, and $I\subset{\mathbb{R}}$ is an interval. Given two
sets of indices $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$, consider the
following sets of functions:

  * •

The non-zero functions $(p_{k})_{k\in\mathcal{I}_{1}}$ forming a complete set
of eigenvectors of $-\Delta$ in $\Omega$ with homogeneous-Dirichlet boundary
conditions on $\overline{\Omega^{*}}\times\partial I$ and homogeneous-Neumann
boundary conditions on $\overline{I}\times\partial\Omega^{*}$.

  * •

The non-zero functions $(q_{h})_{h\in\mathcal{I}_{2}}$ forming a complete set
of eigenvectors of $-\Delta$ in $\Omega$ with homogeneous-Neumann boundary
conditions on $\overline{\Omega^{*}}\times\partial I$ and homogeneous-
Dirichlet boundary conditions on $\overline{I}\times\partial\Omega^{*}$.

Then, we define the sequences of vector fields
$(\bm{\psi}_{k}^{\mathrm{TM}})_{{k}\in\mathcal{I}_{1}}$ and
$(\bm{\psi}_{h}^{\mathrm{TE}})_{{h}\in\mathcal{I}_{2}}\in X_{0}(\Omega)$, via

$\bm{\psi}_{k}^{\mathrm{TM}}:=\frac{\mathrm{curl}(p_{k}\bm{e})}{\|\mathrm{curl}(p_{k}\bm{e})\|_{{H}(\mathrm{curl},\Omega)}},\text{
and
}\quad\bm{\psi}_{h}^{\mathrm{TE}}:=\frac{\mathrm{curl}(\mathrm{curl}(q_{h}\bm{e}))}{\|\mathrm{curl}(\mathrm{curl}(q_{h}\bm{e}))\|_{{H}(\mathrm{curl},\Omega)}},$

for each ${k}\in\mathcal{I}_{1}$ and ${h}\in\mathcal{I}_{2}$, and where
$\bm{e}$ is the unit vector in the distinguished direction. Then, the union of
the two sequences $(\bm{\psi}_{k}^{\mathrm{TM}})_{{k}\in\mathcal{I}_{1}}$ and
$(\bm{\psi}_{h}^{\mathrm{TE}})_{{h}\in\mathcal{I}_{2}}$ forms an orthonormal
basis for $X_{0}(\Omega)\subset{H}_{0}(\mathrm{curl},\Omega)$.

Note that the building of basis functions for the space $X_{0}(\Omega)$ is
valid, for instance, in three-dimensional rectangular, cubic, or cylindrical
domains, provided the eigenbasis of the Laplacian in $\Omega^{*}$. The
calculations are similar in all of these cases. In Tables 1 and 2 we show the
resulting basis functions for the spaces $\nabla H_{0}^{1}(\Omega)$ and
$X_{0}(\Omega)$ on $n$-dimensional cubes $\Omega=[0,\pi]^{n}$ for $n=2,3$. For
clarity, we detail the construction of the eigenbasis on $n$-dimensional cubes
in Appendix B. An extension to more general rectangular domains with distinct
side lengths is trivial, requiring only straightforward yet tedious
calculations and it is therefore omitted.

- | Basis functions | Index $k=(k_{1},k_{2})$
---|---|---
$\nabla H_{0}^{1}(\Omega)$ | $\bm{\phi}_{k}=\frac{2}{\pi|k|}\begin{pmatrix}k_{1}\cos(k_{1}x)\sin(k_{2}y)\\\ k_{2}\sin(k_{1}x)\cos(k_{2}y)\end{pmatrix}$ | $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$
$X_{0}(\Omega)$ | $\bm{\psi}_{k}=\frac{c_{k}}{\sqrt{|k|^{4}+|k|^{2}}}\begin{pmatrix}k_{2}\cos(k_{1}x)\sin(k_{2}y)\\\ -k_{1}\sin(k_{1}x)\cos(k_{2}y)\end{pmatrix}$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $k_{1}=0$ xor $k_{2}=0$
Table 1: The basis functions for the spaces $\nabla H_{0}^{1}(\Omega)$ and $X_{0}(\Omega)$ in 2D with $\Omega=[0,\pi]^{2}$. Here, $c_{k}=\frac{2}{\pi}$ if $k_{\mathrm{i}}>0$ and $c_{k}=\frac{\sqrt{2}}{\pi}$ if $k_{1}=0$ or $k_{2}=0$. - | Basis functions |  Index $k=(k_{1},k_{2})$
---|---|---
$\nabla H_{0}^{1}(\Omega)$ | $\bm{\phi}_{k}=\frac{2\sqrt{2}}{\pi^{3/2}|k|}\begin{pmatrix}k_{1}\cos(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)\\\ k_{2}\sin(k_{1}x)\cos(k_{2}y)\sin(k_{3}z)\\\ k_{3}\sin(k_{1}x)\sin(k_{2}y)\cos(k_{3}z)\end{pmatrix}$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$
$X_{0}(\Omega)$ | $\bm{\psi}_{k}^{\mathrm{TM}}=\frac{2\sqrt{2}}{\pi^{3/2}\sqrt{(1+|k|^{2})(k_{2}^{2}+k_{3}^{2})}}\begin{pmatrix}0\\\ -k_{3}\sin(k_{1}x)\cos(k_{2}y)\sin(k_{3}z)\\\ k_{2}\sin(k_{1}x)\sin(k_{2}y)\cos(k_{3}z)\end{pmatrix}$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $k_{1}>0$ and $k_{2}=0$ xor $k_{3}=0$
| $\bm{\psi}_{k}^{\mathrm{TE}}=\frac{1}{\sqrt{c^{\prime}_{k}(1+|k|^{2})|k|^{2}\left(|k|^{2}-k_{1}^{2}\right)}}\begin{pmatrix}(k_{2}^{2}+k_{3}^{2})\cos(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)\\\ -k_{1}k_{2}\sin(k_{1}x)\cos(k_{2}y)\sin(k_{3}z)\\\ -k_{1}k_{3}\sin(k_{1}x)\sin(k_{2}y)\cos(k_{3}z)\end{pmatrix}$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $k_{2}>0$ and $k_{3}>0$
Table 2: The basis functions for the spaces $\nabla H_{0}^{1}(\Omega)$ and
$X_{0}(\Omega)$ in 3D with $\Omega=[0,\pi]^{3}$. Here,
$c^{\prime}_{k}=\frac{\pi^{3}}{4}$ if $k_{1}=0$ and
$c^{\prime}_{k}=\frac{\pi^{3}}{8}$ if $k_{1}\neq 0$.

##  4 The DFR method

In this section we outline the structure of the Neural Networks that we will
employ and the fundamental principle underlying the DFR method, which is the
construction of a discretized and computable loss function.

###  4.1 Neural Networks

A neural network is a mathematical model comprising multiple compositions of
simple functions, called layers. Specifically, a NN is a non-linear function
$\mathcal{N}({\bm{x}};\bm{W},\bm{b})$ parametrized by a set of weights
$\bm{W}$ and biases $\bm{b}$, with ${\bm{x}}\in{\mathbb{R}}^{d}$ being the
input vector. We restrict our attention to fully-connected feed-forward NNs.
In this architecture, the weights $\bm{W}$ are represented as a collection of
dense matrices ${\bm{W}_{1},\bm{W}_{2},...,\bm{W}_{L}}$, where
$\bm{W}_{\mathrm{j}}\in{\mathbb{R}}^{d_{\mathrm{j}}\times d_{\mathrm{j}-1}}$
is the weight matrix for the layer $\mathrm{j}$ and $d_{\mathrm{j}}$ denotes
the number of nodes on each layer. And the set of biases $\bm{b}$ is a
collection of vectors ${\bm{b}_{1},\bm{b}_{2},...,\bm{b}_{L}}$, where
$\bm{b}_{\mathrm{j}}\in{\mathbb{R}}^{d_{\mathrm{j}}}$ is the bias vector for
the layer $\mathrm{j}$.

The output of each layer of the network is computed as

$\bm{A}_{\mathrm{j}}=\sigma_{\mathrm{j}}(\bm{W}_{\mathrm{j}}\bm{A}_{\mathrm{j}-1}+\bm{b}_{\mathrm{j}})$

where $\sigma_{\mathrm{j}}$ is a nonlinear activation function,
$\bm{A}_{0}={\bm{x}}$ and $\sigma_{L}$ is the identity function. The weights
and biases are obtained via a gradient-based optimization algorithm applied to
a loss function, whose gradients are efficiently calculated via
backpropagation. This process adjusts the parameters to minimize a specified
loss function, which in our case is a discretized and computable approximation
to the dual norm of the PDE residual, as defined in the following section.

We impose homogeneous Dirichlet-type boundary conditions on our candidate
solutions by using a cutoff function
$\xi:\overline{\Omega}\to{\mathbb{R}}^{d\times d}$. Specifically, we define
the approximation of the solution of (2.1) as $\bm{E}=\xi\tilde{\bm{E}}$,
where $\tilde{\bm{E}}$ is the output of the fully-connected feed-forward
neural network $\mathcal{N}({\bm{x}};\bm{W},\bm{b})$ and $\xi$ is smooth and
non-trainable. The function $\xi$ gives matrices that are positive definite in
$\Omega$, and enforces the tangential operator to be zero on the boundary,
i.e., the constraint $\gamma_{t}(\bm{E})=0$ on $\Gamma$.

###  4.2 The discretized loss

In this section, we construct a discretized and computable loss function
$\mathcal{L}$ as an approximation of the dual norm of the residual (2.10). For
the sake of simplicity and exposition, we restrict the details in this section
to squares, i.e., $\Omega=[0,\pi]^{2}$ using the basis functions in Table 1.
The upcoming calculations are similar for other Cartesian product domains in
2D and 3D.

Notice that calculating the integrals in
$\langle\mathcal{R}(\bm{E}),\Phi_{k}\rangle_{H^{\prime}\times H}$ may be
costly since it requires the evaluation of many functions and their
derivatives. However, we underline that the basis functions for the space
${H}_{0}(\mathrm{curl},\Omega)$ are specified in terms of Laplacian’s
eigenbasis and thus both the basis functions and their derivatives are
described in terms of sines and cosines. This feature motivates the use of
efficient methods for approximating such integrals, which are naturally based
on the Fast Fourier Transform (FFT). In this regard, applying Discrete
Sine/Cosine transformations (DST/DCT) is a well-known strategy that
significantly reduces computational complexity when applying a $N$ point
midpoint rule for $N$ basis functions. Specifically, using DST/DCT as a
quadrature rule reduces the amount of required calculations from $O(N^{2})$ to
$O(N\log(N))$. In our case, the Discrete Sine/Cosine transforms appear
naturally when one applies the mid-point integration rule to the integrals
appearing in the residual. Here, we use the type II Sine/Cosine transforms
defined in [4, Section 4.2]. Each transform is represented by an $N\times N$
matrix as

$\displaystyle(S^{II}_{N})_{\mathrm{i}\mathrm{j}}$
$\displaystyle:=\sqrt{\frac{2}{N}}\sigma_{\mathrm{i}}\sin\left(\frac{\pi}{N}\left(\mathrm{j}+\frac{1}{2}\right)(\mathrm{i}+1)\right)\text{
with }\sigma_{\mathrm{i}}=\left\\{\begin{array}[]{c
c}\frac{1}{\sqrt{2}}&\mathrm{i}=N-1,\\\ 1&\mathrm{i}\neq
N-1,\end{array}\right.$ $\displaystyle(C^{II}_{N})_{\mathrm{i}\mathrm{j}}$
$\displaystyle:=\sqrt{\frac{2}{N}}\sigma^{\prime}_{\mathrm{i}}\cos\left(\frac{\pi}{N}\left(\mathrm{j}+\frac{1}{2}\right)\mathrm{i}\right)\text{
with }\sigma^{\prime}_{\mathrm{i}}=\left\\{\begin{array}[]{c
c}\frac{1}{\sqrt{2}}&\mathrm{i}=0,\\\ 1&\mathrm{i}\neq 0,\end{array}\right.$

where $\mathrm{i},\mathrm{j}=0,\dots,N-1$.

Given cutoff frequencies (Fourier modes) $N>0$ and $M>0$, we aim to
approximate $\langle\mathcal{R}(\bm{E}),\Phi_{k}\rangle_{H^{\prime}\times H}$,
where $k=(k_{1},k_{2})$ is such that $0\leq k_{1}\leq N$ and $0\leq k_{2}\leq
M$. Moreover, the indices $k$ satisfy the conditions in Table 1, and we let
$\mathcal{I}$ denote the set of all appropriate indices. For each
$k\in\mathcal{I}$, the basis functions $\Phi_{k}$ are
$\Phi_{k}=\bm{\phi}_{k}+\bm{\psi}_{k}$, with $\bm{\phi}_{k}$ and
$\bm{\psi}_{k}$ as defined in Table 1. Next, we consider the mid-point
integration rule, that is,

$\int_{\Omega}f\,d{\bm{x}}\approx\frac{\pi^{2}}{NM}\sum_{\mathrm{i}=0}^{N-1}\sum_{\mathrm{j}=0}^{M-1}f(x_{\mathrm{i}},y_{\mathrm{j}}),$
(4.1)

where the integration points are $x_{\mathrm{i}}=\frac{2\mathrm{i}+1}{2N}\pi$
and $y_{\mathrm{j}}=\frac{2\mathrm{j}+1}{2M}\pi$. Applying (4.1) in (2.6), we
approximate the integrals appearing in the residual as

$\displaystyle\int_{\Omega}\mu^{-1}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\Phi_{k})\,d{\bm{x}}$
$\displaystyle\approx\sum_{\mathrm{i}=0}^{N-1}\sum_{\mathrm{j}=0}^{M-1}\frac{\pi^{2}}{\sqrt{NM}}\mu^{-1}\mathrm{curl}(\bm{E})\\!\\!\left(x_{\mathrm{i}},y_{\mathrm{j}}\right)\alpha_{k}\left(C^{II}_{N}\right)_{k_{1}\\!-\\!1\mathrm{i}}\left(C^{II}_{M}\right)_{k_{2}\\!-\\!1\mathrm{j}}$
(4.2) $\displaystyle=:\mathcal{R}_{1k}(\mu^{-1}\mathrm{curl}(\bm{E}))$
$\displaystyle\int_{\Omega}\omega^{2}\epsilon\bm{E}\cdot\Phi_{k}\,d{\bm{x}}$
$\displaystyle\approx\sum_{\mathrm{i}=0}^{N-1}\sum_{\mathrm{j}=0}^{M-1}\frac{\pi^{2}}{\sqrt{NM}}\omega^{2}\epsilon\,\bm{E}\\!\left(x_{\mathrm{i}},y_{\mathrm{j}}\right)\cdot{\bm{\alpha}}^{\prime}_{k}\mathbf{C}^{k}_{\mathrm{i}\mathrm{j}}\mathbf{S}^{k}_{\mathrm{i}\mathrm{j}}$
$\displaystyle=:\mathcal{R}_{2k}(\omega^{2}\epsilon\,\bm{E})$
$\displaystyle\int_{\Omega}i\omega\bm{J}\cdot\Phi_{k}\,d{\bm{x}}$
$\displaystyle\approx\sum_{\mathrm{i}=0}^{N-1}\sum_{\mathrm{j}=0}^{M-1}\frac{\pi^{2}}{\sqrt{NM}}i\omega\bm{J}\\!\left(x_{\mathrm{i}},y_{\mathrm{j}}\right)\cdot{\bm{\alpha}}^{\prime}_{k}\mathbf{C}^{k}_{\mathrm{i}\mathrm{j}}\mathbf{S}^{k}_{\mathrm{i}\mathrm{j}}$
$\displaystyle=\mathcal{R}_{2k}(i\omega\bm{J})$

where $\alpha_{k}=\frac{c_{k}(k_{1}k_{2}-k_{1}^{2})}{\sqrt{|k|^{4}+|k|^{2}}}$
and
${\bm{\alpha}}^{\prime}_{k}=\left(\frac{2k_{1}}{\pi|k|}+\frac{c_{k}k_{2}}{\sqrt{|k|^{4}+|k|^{2}}},\frac{2k_{2}}{\pi|k|}-\frac{c_{k}k_{1}}{\sqrt{|k|^{4}+|k|^{2}}}\right)^{t}\\!\\!\mathbb{I}$,
with $\mathbb{I}$ being the identity matrix. The matrices
$\mathbf{C}^{k}_{\mathrm{i}\mathrm{j}}$ and
$\mathbf{S}^{k}_{\mathrm{i}\mathrm{j}}$ contain the cosine and sine
transformations as follows

$\mathbf{C}^{k}_{\mathrm{i}\mathrm{j}}=\begin{pmatrix}\left(C^{II}_{N}\right)_{k_{1}\\!-\\!1\mathrm{i}}&0\\\
0&(C^{II}_{M})_{k_{2}\\!-\\!1\mathrm{j}}\end{pmatrix}\,\text{ and
}\,\mathbf{S}^{k}_{\mathrm{i}\mathrm{j}}=\begin{pmatrix}(S^{II}_{M})_{k_{2}\\!-\\!1\mathrm{j}}&0\\\
0&\left(S^{II}_{N}\right)_{k_{1}\\!-\\!1\mathrm{i}}\end{pmatrix}.$

Notice that, to evaluate the $\mathrm{curl}$ of the candidate solution
$\bm{E}$ one needs to use automatic differentiation as described in [1].
Finally, our discretized loss is

$\mathcal{L}(\bm{E}):=\sqrt{\sum_{k\in\mathcal{I}}|\mathcal{R}_{1k}(\mu^{-1}\mathrm{curl}(\bm{E}))-\mathcal{R}_{2k}(\omega^{2}\epsilon\,\bm{E})-\mathcal{R}_{2k}(i\omega\bm{J})|}.$
(4.3)

We point out that it is also possible to employ a different number of
integration points and basis functions (modes), but we omit the details here.
Notice that the use of less integration points that Fourier modes is expected
to reduce integration errors.

##  5 Numerical experiments

In this section, we present four numerical experiments that illustrate the
main capabilities and some of the limitations of the DFR method. We use
Tensorflow 2.8 and implement a feed-forward fully-connected NN that consists
of five hidden layers, each with $20$ neurons with a $\tanh$ activation
function. In our implementation, Adam serves as the optimizer. Through the use
of a Callback, we allow the optimiser to dynamically modify the learning rate
based on the decay of the loss and reject iteration steps that result in an
increase in loss (see [38]). We choose an starting learning rate of $10^{-4}$.
In addition, we use a validation set and compare the training and validation
losses every iteration. This method aids in detecting overfitting, as
mentioned in [31].

###  5.1 Case 1. Smooth solution in 2D

Let $\Omega=[0,\pi]^{2}$ and consider the variational form: find
$\bm{E}\in{H}_{0}(\mathrm{curl},\Omega)$ satisfying

$\int_{\Omega}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\bm{\phi})+\bm{E}\cdot\bm{\phi}\,d{\bm{x}}=\int_{\Omega}\tilde{\bm{J}}\cdot\bm{\phi}\,d{\bm{x}}\qquad\forall\bm{\phi}\in{H}_{0}(\mathrm{curl},\Omega).$
(5.1)

Here, $\tilde{\bm{J}}$ is chosen such that the exact solution is
$\bm{E}^{*}(x,y)=(xy(y-\pi),xy(x-\pi))^{t}$. Notice that, the bilinear form in
(5.1) is precisely the inner product on ${H}(\mathrm{curl},\Omega)$ and thus
the dual norm of the PDE residual and the $H$(curl)-norm of the error are
equal, in the sense that the constants $M$ and $\gamma$ in (2.8) are equal to
$1$. Whilst this physically does not correspond to the standard Maxwell’s
equations, we include this as a test of the mathematical accuracy of the
method.

We take $N=M=100$ integration points for training the NN and $117$ integration
points on each direction for validation. Moreover, we also select $100$ modes
in both training and validation.

In Figure 1(a) we show the evolution of the loss on the training and
validation data sets. After $10^{4}$ iterations, both losses stabilize and
converge to a limiting value. On the other hand, we show in Figure 1(b) the
contribution of the spaces $\nabla H_{0}^{1}(\Omega)$ and $X_{0}(\Omega)$ to
the squared loss on the training set. We let $\mathcal{L}_{\nabla H_{0}^{1}}$
and $\mathcal{L}_{X_{0}}$ denote the parts of $\mathcal{L}(\bm{E})^{2}$
computed by using only the basis functions of the corresponding subspaces.

(a) The evolution of the loss $\mathcal{L}(\bm{E})$ on both the training and
validation data sets.

(b) The contribution of the spaces $\nabla H_{0}^{1}(\Omega)$ and
$X_{0}(\Omega)$ to the loss $\mathcal{L}(\bm{E})^{2}$ on the training set.

Figure 1: The evolution of the loss in Case 1.

In Figure 2 we show the relationship between the losses and the relative error
of the solution during training and validation. We define the relative error
in terms of the $H$(curl)-norm of the error as:

$\mathcal{E}(\bm{E}):=\frac{\|\bm{E}-\bm{E}^{*}\|_{{H}(\mathrm{curl},\Omega)}}{\|\bm{E}^{*}\|_{{H}(\mathrm{curl},\Omega)}}$

and we always measure this error on the validation set. In Figure 2 we include
an straight line with slope one and highlight the linear relationship between
the loss and the relative error, as expected in view of (2.8). We conclude
that the proposed loss is an accurate approximation of the $H$(curl)-norm of
the error.

Figure 2: The correlation between the loss $\mathcal{L}(\bm{E})$ and the
relative error of the solution $\mathcal{E}(\bm{E})$ during training and
validation in Case 1.

Figure 3 shows the obtained solution, the curl of the solution, and the
corresponding errors calculated pointwise.

(a) Approximate solution.

(b) Curl of the approximate solution.

(c) Error in the solution.

(d) Error in the curl of the solution.

Figure 3: Solution and errors for the model Case 1.

###  5.2 Case 2. Discontinuous parameters in 2D

We take $\Omega=[0,\pi]^{2}$. We define
$\Omega_{0}=\left\\{(x,y):f_{+}(x,y)<1\right\\}$ with
$f_{+}(x,y)=\left(x-\frac{\pi}{2}\right)^{2}+\left(y-\frac{\pi}{2}\right)^{2}$,
and select $\mu$ and $\epsilon$ to be piecewise continuous, such that

$\begin{split}\mu(x,y)=&\mu_{1}{\bf 1}_{\Omega_{0}}+\mu_{2}(1-{\bf
1}_{\Omega_{0}}),\\\ \epsilon(x,y)=&\epsilon_{1}{\bf
1}_{\Omega_{0}}+\epsilon_{2}(1-{\bf 1}_{\Omega_{0}}),\end{split}$

with $\mu_{1}=3,\mu_{2}=1$, $\epsilon_{1}=1$ and $\epsilon_{2}=3$. With these
particular parameters, we will consider two PDEs, one corresponding the
coercive variational form, and a particular case of Maxwell’s equations. In
both cases, we seek for an exact solution
$\bm{E}^{*}=(E_{1}^{*},E_{2}^{*})^{t}$ defined as:

$\begin{split}E_{1}^{*}(x,y)=&\left\\{\begin{array}[]{c
c}-\mu_{1}k_{1}\left(1-f_{-}(x,y)\right)\left(y-\frac{\pi}{2}\right)&(x,y)\in\Omega_{0}\\\
-\mu_{2}k_{2}\left(1-f_{-}(x,y)\right)\left(r^{2}-f_{-}(x,y)\right)\left(y-\frac{\pi}{2}\right)&\text{else}\end{array}\right.\\\
E_{2}^{*}(x,y)=&\left\\{\begin{array}[]{c
c}-\mu_{1}k_{1}\left(1-f_{-}(x,y)\right)\left(x-\frac{\pi}{2}\right)&(x,y)\in\Omega_{0}\\\
-\mu_{2}k_{2}\left(1-f_{-}(x,y)\right)\left(r^{2}-f_{-}(x,y)\right)\left(x-\frac{\pi}{2}\right)&\text{else}\end{array}\right.\\\
\end{split}$ (5.2)

where
$f_{-}(x,y)=\left(x-\frac{\pi}{2}\right)^{2}-\left(y-\frac{\pi}{2}\right)^{2}$,
$k_{1}=1,k_{2}=35$ and $r=6$. In particular, we note that the solution admits
discontinuities in both the vector field itself and the curl across the
$\partial\Omega_{0}$, and thus the solution is in
${H}(\mathrm{curl},\Omega)\setminus H^{1}(\Omega)$.

###  5.3 Case 2.1. Coercive variational form

First, we consider the variational form: find
$\bm{E}\in{H}_{0}(\mathrm{curl},\Omega)$ satisfying

$\int_{\Omega}\mu^{-1}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\bm{\phi})+\epsilon\bm{E}\cdot\bm{\phi}\,d{\bm{x}}=\int_{\Omega}\tilde{\bm{J}}\cdot\bm{\phi}\,d{\bm{x}}\qquad\forall\bm{\phi}\in{H}_{0}(\mathrm{curl},\Omega).$
(5.3)

Here, $\tilde{\bm{J}}$ is chosen such that the exact solution is (5.2). Notice
that the exact solution does not satisfy the strong form of the equation due
to interface conditions at the discontinuities.

In this case, as the loss requires the integration of discontinuous functions,
we expect more significant integration errors during training. Thus, we employ
a larger number of integration points than modes to mitigate this problem.
Specifically, we take $N=M=200$ integration points for training the NN and
$234$ points on each direction for validation. In both cases, we use $150$
modes.

(a) The evolution of the loss $\mathcal{L}(\bm{E})$ on both the training and
validation data sets.

(b) The contribution of the spaces $\nabla H_{0}^{1}(\Omega)$ and
$X_{0}(\Omega)$ to the loss $\mathcal{L}(\bm{E})^{2}$ on the training set.

Figure 4: The evolution of the loss in Case 2.1.

The evolution of the loss on the training and validation data sets is shown in
Figure 4(a). Figure 4(b) shows the contribution of the spaces $\nabla
H_{0}^{1}(\Omega)$ and $X_{0}(\Omega)$ to the squared loss on the training
set, denoted $\mathcal{L}_{\nabla H_{0}^{1}}$ and $\mathcal{L}_{X_{0}}$,
respectively. The largest loss is attributed to the space ${X_{0}(\Omega)}$.
This means that any enhancements to the basis functions in this domain will
have a major impact on lowering the overall loss. Further research in this
area could include parameterizing the number of basis functions related to
each space separately.

Figure 5 shows the obtained solution, the curl of the solution and the
corresponding errors. As expected, the maximum errors are located close to the
discontinuities of the parameters, which is inevitable in our implementation
as we are using smooth NNs to approximate discontinuous functions.

(a) Approximate solution.

(b) Curl of the approximate solution.

(c) Error in the solution.

(d) Error in the curl of the solution.

Figure 5: Solution and errors for the model Case 2.1.

The relationship between losses and the relative error in the solution at each
iteration is shown in Figure 6. In the asymptotic regime, we obtain a linear
relationship between the loss and the error. Since (5.3) relates with the
inner product of ${H}(\mathrm{curl},\Omega)$, in this case we can use the
Riesz representation theorem to estimate the equivalence constants $M$ and
$\gamma$ in (2.8). A simple calculation give us $\frac{1}{M}=\frac{1}{3}$ and
$\frac{1}{\gamma}=3$. Figure 6 exhibits the expected oscillation between these
two parallel lines.

Figure 6: The correlation between the loss $\mathcal{L}(\bm{E})$ and the
relative error in the solution $\mathcal{E}(\bm{E})$ during training and
validation in Case 2.1.

Now, we remark the importance of accurately calculating the integrals in
(4.3). Figure 7 shows the loss evolution and the correlation between the
$H$(curl)-norm of the error and the loss when the number of integration points
is equal to the number of integration points in the Case 2.1. There we use
$N=M=100$ integration points for training the NN and $120$ points on each
direction for validation. In Figure 7(a) we observe a divergence in the
training and validation losses after roughly $2000$ iterations, owing to
integration errors, highlighting the need for accurate integration when
solutions are of lower regularity.

(a) The evolution of the loss $\mathcal{L}(\bm{E})$ on both the training and
validation data sets.

(b) The correlation between the loss and the relative error.

Figure 7: The evolution of the loss and the correlation between the loss
$\mathcal{L}(\bm{E})$ and the relative error in the solution
$\mathcal{E}(\bm{E})$ during training and validation in Case 2.1 when using an
inaccurate integration rule.

###  5.4 Case 2.2. The physical variational form

Now, we modify Case 2 and seek for $\bm{E}\in{H}_{0}(\mathrm{curl},\Omega)$
satisfying

$\int_{\Omega}\mu^{-1}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\bm{\phi})-\omega^{2}\epsilon\bm{E}\cdot\bm{\phi}\,d{\bm{x}}=\int_{\Omega}\tilde{\bm{J}}\cdot\bm{\phi}\,d{\bm{x}},$

with $\omega=1.25$. This modification corresponds to the weak formulation of
the time-harmonic Maxwell’s equations in (2.5) with positive discontinuous
parameters $\mu$ and $\epsilon$.

(a) The evolution of the loss $\mathcal{L}(\bm{E})$ on both the training and
validation data sets.

(b) The correlation between the loss and the relative error.

Figure 8: The evolution of the loss and the correlation between the loss
$\mathcal{L}(\bm{E})$ and the relative error in the solution
$\mathcal{E}(\bm{E})$ during training and validation in Case 2.2.

In contrast to the Case 2.1, in Figure 8 one notices that the evolution of the
loss decreases slower and the relation between the loss and the error is not
linear. This effect is a clear consequence of the inclusion of the frequency
$\omega$ that modifies the bounds of the error in (2.8).

###  5.5 Case 3. Smooth solution in 3D

Take $\Omega=[0,\pi]^{3}$. We consider the variational form: find
$\bm{E}\in{H}_{0}(\mathrm{curl},\Omega)$ satisfying

$\int_{\Omega}\mu^{-1}\mathrm{curl}(\bm{E})\cdot\mathrm{curl}(\bm{\phi})-\omega^{2}\epsilon\bm{E}\cdot\bm{\phi}\,d{\bm{x}}=\int_{\Omega}\tilde{\bm{J}}\cdot\bm{\phi}\,d{\bm{x}}\qquad\forall\bm{\phi}\in{H}_{0}(\mathrm{curl},\Omega).$
(5.4)

We choose $\mu$ and $\epsilon$ to be constant equal to $1$ and $\omega=1.5$.
Here, $\tilde{\bm{J}}$ is such that the exact solution is

$\bm{E}^{*}(x,y,z)=\begin{bmatrix}\sin(y)\sin(z)\sin(\omega x)\\\
\sin(x)\sin(z)\sin(\omega y)\\\ \sin(x)\sin(y)\sin(\omega z)\end{bmatrix}.$
(5.5)

Here, we impose homogeneous-Dirichlet boundary conditions on $\partial\Omega$.
We take a partition of $50$ integration points on each direction for training
the NN and $60$ points for validation. Moreover, we use $50$ modes in both
training and validation.

Similarly to the 2D cases, Figure 9(a) shows the evolution of the loss on the
training and validation data sets. After $10^{5}$ iterations, the losses
stabilize and reach values upto $10^{-3}$. Figure 9(b) shows the contribution
of the spaces $\nabla H_{0}^{1}(\Omega)$ and $X_{0}(\Omega)$ to the squared
loss on the training set.

(a) The evolution of the loss $\mathcal{L}(\bm{E})$ on both the training and
validation data sets.

(b) The contribution of the spaces $\nabla H_{0}^{1}(\Omega)$ and
$X_{0}(\Omega)$ to the loss $\mathcal{L}(\bm{E})^{2}$ on the training set.

Figure 9: The evolution of the loss in Case 3. Figure 10: The correlation
between the loss $\mathcal{L}(\bm{E})$ and the relative error in the solution
$\mathcal{E}(\bm{E})$ during training and validation in Case 3.

Finally, we illustrate the correlation between the loss and the relative error
in the solution during training and validation in Figure 10. In the asymptotic
regime, the relationship between the loss and the error is linear.

##  6 Conclusions

We extended the principles in [37] and implemented the DFR method for solving
Maxwell’s equations. In this case, we select the dual norm of the residual as
the loss function of a NN solving Maxwell’s problem. We rely on the weak
formulation of Maxwell’s electric field and the Helmholtz decomposition of the
space ${H}_{0}(\mathrm{curl},\Omega)$. We proposed orthonormal basis functions
for each sub-space and used them to construct a computable loss function. To
lower the computational cost of estimating integrals, we apply DST/DCT
transformations in our discretized loss function. In two and three spatial
dimensions, the numerical examples show a linear association between the loss
and the $H$(curl)-norm of the error. We included examples with discontinuous
parameters an noticed the importance of avoiding overfitting and integration
errors.

We note that the DFR discussed here suffers the curse of dimensionality and
that the implementation differs in 2D and 3D. Additionally, the DFR in 2D can
only be applied to general domains if the eigenbasis of the Laplacian on those
geometries is known. The 3D version is restricted to Cartesian product
domains, which are rectangular or cylindrical domains based on the Laplacian’s
eigenbasis in 2D, which again must be known in order to implement the method.
Future work will look at the method’s scalability and investigate the use of
subdomain-based local test functions. By doing so, we will consider more
general geometries.

##  7 Acknowledgements

Jamie M. Taylor is supported by the Basque Government through the BERC
2018-2021 program and by the Spanish State Agency of Research through “BCAM
Severo Ochoa” accreditation of excellence SEV-2017-0718 and through the
project (PID2020-114189RB-I00 / AEI / 10.13039 / 501100011033). David Pardo
has received funding from: the Spanish Ministry of Science and Innovation
projects with references TED2021-132783B-I00, PID2019-108111RB-I00 (FEDER/AEI)
and PDC2021-121093-I00 (MCIN / AEI / 10.13039/501100011033/Next Generation
EU), the “BCAM Severo Ochoa” accreditation of excellence CEX2021-001142-S /
MICIN / AEI / 10.13039/ 501100011033; the Spanish Ministry of Economic and
Digital Transformation with Misiones Project IA4TES (MIA.2021.M04.008 /
NextGenerationEU PRTR); and the Basque Government through the BERC 2022-2025
program, the Elkartek project SIGZE (KK-2021/00095), and the Consolidated
Research Group MATHMODE (IT1456-22) given by the Department of Education.
Ignacio Muga is supported by the Chilean National Agency for Research &
Development through the Fondecyt Project #1230091.

## References

  * [1] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of Marchine Learning Research, 18 (2018), pp. 1–43.
  * [2] I. Brevis, I. Muga, D. Pardo, O. Rodríguez, and K. G. van der Zee, Learning quantities of interest from parametric PDEs: An efficient neural-weighted Minimal Residual approach, arXiv preprint arXiv:2304.01722, (2023).
  * [3] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations, Springer New York, NY, 2010.
  * [4] V. Britanak, P. C. Yip, and K. R. Rao, Discrete cosine and sine transforms: general properties, fast algorithms and integer approximations, Elsevier, 2010.
  * [5] A. Buffa, M. Costabel, and D. Sheen, On traces for H(curl, Ω) in Lipschitz domains, Journal of Mathematical Analysis and Applications, 276 (2002), pp. 845–867.
  * [6] F. Calabrò, S. Cuomo, D. di Serafino, G. Izzo, and E. Messina, The effect of time discretization on the solution of parabolic PDEs with ANNs, arXiv preprint arXiv:2206.00452, (2022).
  * [7] F. Calabrò, G. Fabiani, and C. Siettos, Extreme learning machine collocation for the numerical solution of elliptic PDEs with sharp gradients, Computer Methods in Applied Mechanics and Engineering, 387 (2021), p. 114188.
  * [8] L. Camargo, B. López-Rodríguez, M. Osorio, and M. Solano, An HDG method for Maxwell’s equations in heterogeneous media, Computer Methods in Applied Mechanics and Engineering, 368 (2020), p. 113178.
  * [9] T. Chaumont-Frelet and P. Vega, Frequency-explicit a posteriori error estimates for discontinuous Galerkin discretizations of Maxwell’s equations, arXiv preprint arXiv:2208.01475, (2022).
  * [10] M. Costabel and M. Dauge, Maxwell eigenmodes in product domains, Maxwell’s equations: analysis and numerics, Berlin, Boston: De Gruyter, (2019), pp. 171–198.
  * [11] E. B. Davies, Spectral Theory and Differential Operators, Cambridge Studies in Advanced Mathematics, Cambridge University Press, 1995.
  * [12] L. Demkowicz and L. Vardapetyan, Modeling of electromagnetic absorption/scattering problems using hp-adaptive finite elements, Computer Methods in Applied Mechanics and Engineering, 152 (1998), pp. 103–124.
  * [13] K. Doległo, M. Paszyński, and L. Demkowicz, Deep Neural Networks and Smooth Approximation of PDEs, in Computational Science–ICCS 2022: 22nd International Conference, London, UK, June 21–23, 2022, Proceedings, Part II, Springer, 2022, pp. 328–332.
  * [14] A. Ern and J.-L. Guermond, Finite Elements II: Galerkin Approximation, Elliptic and Mixed PDEs, Springer Cham, 2021.
  * [15] L. Espath, P. Behnoudfar, and R. Tempone, Physics-informed Spectral Learning: the Discrete Helmholtz–Hodge Decomposition, arXiv preprint arXiv:2302.11061, (2023).
  * [16] L. Espath, D. Kabanov, J. Kiessling, and R. Tempone, Statistical learning for fluid flows: Sparse Fourier divergence-free approximations, Physics of Fluids, 33 (2021), p. 097108.
  * [17] L. C. Evans, Partial Differential Equations, vol. 19 of Graduate Studies in Mathematics, American Mathematical Society, 2nd ed., 1998.
  * [18] V. Girault and P.-A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, vol. 5, Springer Science & Business Media, 2012.
  * [19] B. Gralak, Exact Modal Methods, in Gratings: Theory and Numeric Applications, E. Popov, ed., AMU (PUP), Dec. 2012, pp. 10.1–10.21.
  * [20] J. Han, A. Jentzen, and W. E, Solving high-dimensional partial differential equations using deep learning, Proceedings of the National Academy of Sciences, 115 (2018), pp. 8505–8510.
  * [21] D. S. Jones, Methods in electromagnetic wave propagation, Oxford, (1979).
  * [22] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nature Reviews Physics, 3 (2021), pp. 422–440.
  * [23] E. Kharazmi, Z. Zhang, and G. Karniadakis, VPINNs: variational physics-informed neural networks for solving partial differential equations, arXiv preprint arXiv:1912.00873, (2019).
  * [24] E. Kharazmi, Z. Zhang, and G. E. Karniadakis, hp-VPINNs: Variational physics-informed neural networks with domain decomposition, Computer Methods in Applied Mechanics and Engineering, 374 (2021), p. 113547.
  * [25] Y. Khoo, J. Lu, and L. Ying, Solving parametric PDE problems with artificial neural networks, European Journal of Applied Mathematics, 32 (2021), pp. 421–435.
  * [26] A. Kirsch and F. Hettlich, The mathematical theory of time-harmonic Maxwell’s equations, volume 190 of Applied Mathematical Sciences, Springer, Cham, 2014.
  * [27] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature, 521 (2015), pp. 436–444.
  * [28] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, DeepXDE: A deep learning library for solving differential equations, SIAM review, 63 (2021), pp. 208–228.
  * [29] P. Monk and Y. Zhang, Finite element methods for Maxwell’s equations, Oxford University Press, 2003.
  * [30] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics, 378 (2019), pp. 686–707.
  * [31] J. A. Rivera, J. M. Taylor, Á. J. Omella, and D. Pardo, On quadrature rules for solving partial differential equations using neural networks, Computer Methods in Applied Mechanics and Engineering, 393 (2022), p. 114710.
  * [32] V. A. Rukavishnikov and A. O. Mosolapov, New numerical method for solving time-harmonic Maxwell equations with strong singularity, Journal of Computational Physics, 231 (2012), pp. 2438–2448.
  * [33] L. Ruthotto and E. Haber, Deep neural networks motivated by partial differential equations, Journal of Mathematical Imaging and Vision, 62 (2020), pp. 352–364.
  * [34] H. M. Schey and H. M. Schey, Div, grad, curl, and all that: an informal text on vector calculus, WW Norton & Company, New York, USA, 2004\.
  * [35] T. Służalec, R. Grzeszczuk, S. Rojas, W. Dzwinel, and M. Paszyński, Quasi-optimal hp-finite element refinements towards singularities via deep neural network prediction, Computers & Mathematics with Applications, 142 (2023), pp. 157–174.
  * [36] A. Taflove and S. C. Hagness, Computational Electromagnetics: The Finite-Difference Time-Domain Method, Artech House, 2005.
  * [37] J. M. Taylor, D. Pardo, and I. Muga, A Deep Fourier Residual method for solving PDEs using Neural Networks, Computer Methods in Applied Mechanics and Engineering, 405 (2023), p. 115850.
  * [38] C. Uriarte, D. Pardo, and Á. J. Omella, A finite element based deep learning solver for parametric PDEs, Computer Methods in Applied Mechanics and Engineering, 391 (2022), p. 114562.

##  Appendix A Basis for the space ${H}_{0}(\mathrm{curl},\Omega)$

This appendix contains classical results available, for example, in [11, 17].
First, in the following proposition, we construct a set of auxiliary
eigenvectors of the operator $-\Delta$ in $H^{1}(\Omega)$.

######  Proposition A.1.

Let $\Omega$ be a bounded, simply connected, and Lipschitz domain in
$\mathbb{R}^{n}$, with $n=2$ or $3$ and boundary $\partial\Omega$. Assume that
$\partial\Omega$ consists of two disjoint parts $\Gamma_{D}$ and $\Gamma_{N}$,
such that $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$. Then, there exists a
countable set of eigenfunctions $(\phi_{k})_{k\in\mathcal{I}}$, with
$\mathcal{I}$ being a set of indices, and corresponding eigenvalues
$\lambda_{k}>0$ such that the following results hold, for all
$k\in\mathcal{I}$

$\displaystyle-\Delta\phi_{k}$ $\displaystyle=\lambda_{k}\phi_{k}$
$\displaystyle\text{ in }\Omega,$ (A.1) $\displaystyle\phi_{k}$
$\displaystyle=0$ $\displaystyle\text{ on }\Gamma_{D},$
$\displaystyle\nabla\phi_{k}\cdot\bm{\mathrm{n}}$ $\displaystyle=0$
$\displaystyle\text{ on }\Gamma_{N}.$

Moreover, the sequence $(\phi_{k})_{k\in\mathcal{I}}$ forms an orthonormal
basis for $L^{2}(\Omega)$ and an orthogonal basis for $H^{1}(\Omega)$.

The proof of Proposition A.1 is a classical result available, for example, in
[17, Theorem D.6.7] and [11, Theorem 1.8]. Here, we only prove Propositions
3.1 and 3.2.

### Proof of Proposition 3.1

###### Proof.

First we note that, since the space $\nabla H_{0}^{1}(\Omega)$ as a subspace
of ${H}(\mathrm{curl},\Omega)$ is isometric to $H_{0}^{1}(\Omega)$, then the
gradients of any orthogonal basis of $H^{1}_{0}(\Omega)$ define an orthogonal
basis of $\nabla H^{1}_{0}(\Omega)$.

Take $\partial\Omega=\Gamma_{D}$ in Proposition A.1. Then, there exists a set
of eigenfunctions $(\phi_{k})_{k\in\mathcal{I}}$ satisfying (A.1), which are
called homogeneous-Dirichlet eigenvectors of $-\Delta$ in $\Omega$. For each
$k\in\mathcal{I}$, the function $\phi_{k}\in H^{1}_{0}(\Omega)$, so the
sequence $(\phi_{k})_{k\in\mathcal{I}}$ forms an orthonormal basis for
$L^{2}(\Omega)$ and an orthogonal basis for $H^{1}_{0}(\Omega)$. Therefore,
the sequence $(\nabla\phi_{k})_{k\in\mathcal{I}}$ forms an orthogonal basis of
$\nabla H^{1}_{0}(\Omega)$ due to the isometry between the spaces. The result
in Proposition 3.1 follows from the normalization of the functions
$\nabla\phi_{k}$. ∎

### Proof of Proposition 3.2

###### Proof.

Take $\partial\Omega=\Gamma_{N}$ in Proposition A.1. Then, there exists a set
of eigenfunctions $(\phi_{k})_{k\in\mathcal{I}}$ satisfying (A.1). These are
called homogeneous-Neumann eigenvectors of $-\Delta$ in $\Omega$. We disregard
the constant functions and define
$\hat{\bm{\psi}}_{k}:=\mathrm{curl}^{*}(\phi_{k})$ for all $k$ in the index
set $\mathcal{I}$. Each function $\phi_{k}$ is smooth, so one certainly has
that $\hat{\bm{\psi}}_{k}\in{H}(\mathrm{curl},\Omega)$ for all
$k\in\mathcal{I}$.

From (A.1), we obtain that
$\mathrm{curl}^{*}(\phi_{k})\times\bm{\mathrm{n}}=\nabla\phi_{k}\cdot\bm{\mathrm{n}}=0$
on $\Gamma_{N}$, so that
$\hat{\bm{\psi}}_{k}\in{H}_{0}(\mathrm{curl},\Omega)$. Moreover, for all $u\in
H_{0}^{1}(\Omega)$ we have

$(\hat{\bm{\psi}}_{k},\nabla
u)_{{H}(\mathrm{curl},\Omega)}=\int_{\Omega}\mathrm{curl}^{*}(\phi_{k})\cdot\nabla
u\,d{\bm{x}}=\int_{\Omega}\phi_{k}\,\mathrm{curl}(\nabla u)\,d{\bm{x}}=0,$

meaning $\mathrm{div}(\hat{\bm{\psi}}_{k})=0$ weakly, so
$\hat{\bm{\psi}}_{k}\in X_{0}(\Omega)$ for all $k\in\mathcal{I}$.

To show that the sequence $(\hat{\bm{\psi}}_{k})_{k\in\mathcal{I}}$ forms an
orthogonal basis of $X_{0}(\Omega)$, we turn first to the orthogonality of the
functions. Clearly,
$\mathrm{curl}(\hat{\bm{\psi}}_{k})=-\Delta\phi_{k}=\lambda_{k}\phi_{k}$.
Therefore, if $k\neq k^{\prime}$, one obtains

$(\hat{\bm{\psi}}_{k},\hat{\bm{\psi}}_{k^{\prime}})_{{H}(\mathrm{curl},\Omega)}=\lambda_{k^{\prime}}(\lambda_{k}+1)\int_{\Omega}\phi_{k}\phi_{k^{\prime}}\,d{\bm{x}}=0,$

as $\phi_{k}$ and $\phi_{k^{\prime}}$ are orthogonal in $L^{2}(\Omega)$.

Now, to show that they form a complete basis, we assume otherwise for the sake
of contradiction, that there exists some non-zero $\bm{v}\in X_{0}(\Omega)$
such that $(\hat{\bm{\psi}}_{k},\bm{v})_{{H}(\mathrm{curl},\Omega)}=0$ for all
$k\in\mathcal{I}$. In this case, we can write that

$\begin{split}0=&(\hat{\bm{\psi}}_{k},\bm{v})_{{H}(\mathrm{curl},\Omega)}=\int_{\Omega}(-\Delta\phi_{k}+\phi_{k})\mathrm{curl}(\bm{v})\,d{\bm{x}}=\left(\lambda_{k}+1\right)\int_{\Omega}\phi_{k}\mathrm{curl}(\bm{v})\,d{\bm{x}}.\end{split}$

As $\lambda_{k}>0$ and $(\phi_{k})_{k\in\mathcal{I}}$ is a complete orthogonal
basis for $L^{2}(\Omega)$, disregarding constant functions, this implies that
$\mathrm{curl}(\bm{v})$ is constant. We write
$\mathrm{curl}(\bm{v})=c\in{\mathbb{R}}$. Accordingly,
$\mathrm{curl}^{*}(\mathrm{curl}(\bm{v}))=0$, so simple calculations lead to

$c^{2}|\Omega|=\int_{\Omega}(\mathrm{curl}(\bm{v}))^{2}\,dx=\int_{\Omega}\bm{v}\cdot\mathrm{curl}^{*}(\mathrm{curl}(\bm{v}))\,dx=0.$

Thus, $\mathrm{curl}(\bm{v})=0$ and $\mathrm{div}(\bm{v})=0$ and this implies
that $\bm{v}=0$, contradicting our original assumption.

We conclude that $(\hat{\bm{\psi}}_{k})_{k\in\mathcal{I}}$ is a complete
orthogonal basis of $X_{0}(\Omega)$. The result in Proposition 3.2 follows
from the normalization of the functions $\hat{\bm{\psi}}_{k}$. ∎

As discussed in Section 3, to find an orthogonal basis for the space
$X_{0}(\Omega)$ in 3D is not as straightforward as in the 2D case.
Nonetheless, in [10] the authors demonstrate that the TE and TM modes
constitute a complete and orthogonal basis of $X_{0}(\Omega)$. Here, we only
show that the functions $\bm{\psi}_{k}^{\mathrm{TE}}\in X_{0}(\Omega)$ and
$\bm{\psi}_{k}^{\mathrm{TM}}\in X_{0}(\Omega)$ for all the appropriate indices
$k$.

#### TM modes

Take $\Gamma_{D}=\overline{\Omega^{*}}\times\partial I$ and
$\Gamma_{N}=\partial\Omega^{*}\times\bar{I}$ in Proposition A.1. Then, there
exists a set of eigenfunctions $(p_{k})_{k\in\mathcal{I}}$ satisfying (A.1).
Disregarding the zero function, for all $k$ in the index set $\mathcal{I}$, we
define $\hat{\bm{\psi}}_{k}^{\mathrm{TM}}:=\mathrm{curl}(p_{k}\bm{e})$, where
$\bm{e}$ is the unit vector in the $z$-direction. The functions
$\hat{\bm{\psi}}_{k}^{\mathrm{TM}}$ are the $\mathrm{curl}$ of vector fields
with a single non-zero component in the distinguished direction and they are
divergence free.

Using the vector identities of the curl, the operator
$\mathrm{curl}\text{-}\mathrm{curl}$ reduces to the negative Laplacian when
acting on the vector field $\hat{\bm{\psi}}_{k}^{\mathrm{TM}}$, i.e.,

$\mathrm{curl}(\mathrm{curl}(\hat{\bm{\psi}}_{k}^{\mathrm{TM}}))=\nabla(\mathrm{div}(\hat{\bm{\psi}}_{k}^{\mathrm{TM}}))-\Delta\hat{\bm{\psi}}_{k}^{\mathrm{TM}}=-\mathrm{curl}(\Delta
p_{k}\bm{e})=\lambda_{k}\bm{\psi}_{k}^{\mathrm{TM}},$ (A.2)

then $\hat{\bm{\psi}}_{k}^{\mathrm{TM}}$ are eigenvectors of the curl-curl
operator with corresponding eigenvalues $\lambda_{k}>0$. Direct calculations
verify that $\hat{\bm{\psi}}_{k}^{\mathrm{TM}}$ satisfies the boundary
condition $\hat{\bm{\psi}}_{k}^{\mathrm{TM}}\times\bm{\mathrm{n}}=0$, implying
that $\hat{\bm{\psi}}_{k}^{\mathrm{TM}}\in X_{0}(\Omega)$ for all
$k\in\mathcal{I}$.

#### TE Modes

Similarly to the TM modes, take
$\Gamma_{N}=\overline{\Omega^{*}}\times\partial I$ and
$\Gamma_{D}=\partial\Omega^{*}\times\bar{I}$ in Proposition A.1. Then, there
exists a set of eigenfunctions $(\phi_{k})_{k\in\mathcal{I}}$ satisfying
(A.1). Disregarding the zero function, for all $k$ in the index set
$\mathcal{I}$, we define
$\hat{\bm{\psi}}_{k}^{\mathrm{TE}}:=\mathrm{curl}(\mathrm{curl}(\phi_{k}\bm{e}))$,
where $\bm{e}$ is the unit vector in the $z$-direction. The functions
$\hat{\bm{\psi}}_{k}^{\mathrm{TE}}$ are curl of smooth functions and
divergence free. Moreover, $\hat{\bm{\psi}}_{k}^{\mathrm{TE}}$ are
eigenvectors of the curl-curl operator, which is inherited from the fact that
$\phi_{k}$ are eigenvectors of the Laplacian. A direct calculation verifies
that these satisfy the boundary condition
$\hat{\bm{\psi}}_{k}^{\mathrm{TE}}\times\bm{\mathrm{n}}=0$, implying that
$\hat{\bm{\psi}}_{k}^{\mathrm{TE}}\in X_{0}(\Omega)$ for all
$k\in\mathcal{I}$.

##  Appendix B Basis functions on $\Omega=[0,\pi]^{n}$

Here, we consider the domain $\Omega=[0,\pi]^{n}$ with $n=2$ or $3$. Table 3
shows the eigenvectors of the Laplacian in two and three dimensions for
different types of boundary conditions. In 3D, we take
$\Omega=[0,\pi]^{3}=[0,\pi]\times\Omega^{*}$, where $\Omega^{*}=[0,\pi]^{2}$.

We note that the eigenvectors for the Laplacian on rectangular domains with
distinct side lengths are readily obtained by rescaling those in Table 3. The
constructions of basis functions for ${H}_{0}(\mathrm{curl},\Omega)$ on
rectangular domains can be obtained via the same constructions outlined in
Sections B.I-B.IV. The appearance of numerous multiplicative constants makes
the arithmetic unwieldy, albeit simple, thus for simplicity we only include
the calculations for $n$-dimensional cubes.

In Table 3, the mixed boundary conditions of type 1 are homogeneous-Dirichlet
on $\\{0,\pi\\}\times\overline{\Omega^{*}}$ and homogeneous-Neumann on
$[0,\pi]\times\partial\Omega^{*}$, and, the type 2 are the reversed. We let
$\mathcal{I}$ be the set of indices $k$ such that $k=(k_{1},k_{2})$ in 2D,
$k=(k_{1},k_{2},k_{3})$ in 3D, and each $\phi_{k}$ satisfies the problem
(A.1). We disregard some indices $k$, as shown in Table 3. Notice that the
eigenvalues corresponding to all the eigenvectors $\phi_{k}$ in Table 3 are
$\lambda_{k}=|k|^{2}$. More details about the properties of these functions
can be found in [37].

|  Boundary conditions |  $\phi_{k}$ |  Index $k$
---|---|---|---
2D |  Dirichlet |  $\frac{2}{\pi}\sin(k_{1}x)\sin(k_{2}y)$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$
2D |  Neumann |  $\frac{2}{\pi}\cos(k_{1}x)\cos(k_{2}y)$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $k_{1}=0$ xor $k_{2}=0$
3D |  Dirichlet |  $\frac{2\sqrt{2}}{\pi^{3/2}}\sin(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$
3D |  Mixed type 1 |  $\sin(k_{1}x)\cos(k_{2}y)\cos(k_{3}z)$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $k_{1}>0$ and $k_{2}=0$ xor $k_{3}=0$
3D |  Mixed type 2 |  $\cos(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)$ |  $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $k_{2}>0$ and $k_{3}>0$
Table 3: The eigenvectors of $-\Delta$ on $n$-dimensional cubes (2D and 3D)
with different boundary conditions.

###  B.I Basis for $\nabla H^{1}_{0}(\Omega)$ in 2D

Let $\Omega=[0,\pi]^{2}$. From Table 3, we take
$\phi_{k}(x,y)=\frac{2}{\pi}\sin(k_{1}x)\sin(k_{2}y)$ satisfying (A.1) with
$\Gamma_{D}=\partial\Omega$ and $k=(k_{1},k_{2})$. Here, $k_{1}$, $k_{2}$ and
$k_{3}$ are integers greater or equal to zero, and we let $\mathcal{I}$ be the
set of all such indices $k$. Notice that
$\|\nabla\phi_{k}\|_{[L^{2}(\Omega)]^{2}}=|k|$ for all $k\in\mathcal{I}$.
Using Proposition 3.1, we conclude that the set of functions
$(\bm{\phi}_{k})_{k\in\mathcal{I}}$, where each function $\bm{\phi}_{k}$ is
defined as

$\bm{\phi}_{k}:=\frac{1}{|k|}\nabla\phi_{k}=\frac{2}{\pi|k|}\begin{pmatrix}k_{1}\cos(k_{1}x)\sin(k_{2}y)\\\
k_{2}\sin(k_{1}x)\cos(k_{2}y)\end{pmatrix},$

forms an orthonormal basis for $\nabla
H^{1}_{0}(\Omega)\subset{H}_{0}(\mathrm{curl},\Omega)$ in 2D.

###  B.II Basis for $X_{0}(\Omega)$ in 2D

Let $\Omega=[0,\pi]^{2}$. From Table 3, we take
$\phi_{k}(x,y)=\frac{2}{\pi}\cos(k_{1}x)\cos(k_{2}y)$ satisfying (A.1) with
$\Gamma_{N}=\partial\Omega$ and $k=(k_{1},k_{2},k_{3})$. Here,
$k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $\mathrm{i}=1$ or $2$, and we
let $\mathcal{I}$ be the set of all such indices $k$. We disregard the case
when $k_{1}=k_{2}=0$. Let $\hat{\bm{\psi}}_{k}$ be the functions obtained by
applying the operator $\mathrm{curl}^{*}$ to $\phi_{k}$, i.e.,
$\hat{\bm{\psi}}_{k}:=\mathrm{curl}^{*}(\phi_{k})$. A direct calculations lead
to

$\|\hat{\bm{\psi}}_{k}\|_{H(\mathrm{curl},\Omega)}^{2}=\int_{\Omega}|k|^{4}|\phi_{k}|^{2}+\left(\mathrm{curl}(\mathrm{curl}^{*}(\phi_{k}))\right)\phi_{k}\,d{\bm{x}}=|k|^{4}+|k|^{2},$

as $\phi_{k}$ have unit $L^{2}$-norm. For all ${k\in\mathcal{I}}$, we obtain

$\bm{\psi}_{k}:=\frac{\hat{\bm{\psi}}_{k}}{\sqrt{|k|^{4}+|k|^{2}}}=\frac{c_{k}}{\sqrt{|k|^{4}+|k|^{2}}}\begin{pmatrix}k_{2}\cos(k_{1}x)\sin(k_{2}y)\\\
-k_{1}\sin(k_{1}x)\cos(k_{2}y)\end{pmatrix},$

where

$c_{k}=\left\\{\begin{array}[]{c c}\frac{2}{\pi}&k_{1}>0\text{ and
}k_{2}>0,\\\ \frac{\sqrt{2}}{\pi}&k_{1}=0\text{ xor
}k_{2}=0.\end{array}\right.$

Using Proposition 3.2, we conclude that the set
$(\bm{\psi}_{k})_{k\in\mathcal{I}}$ forms an orthonormal basis of the space
$X_{0}(\Omega)$ in 2D.

###  B.III Basis for $\nabla H^{1}_{0}(\Omega)$ in 3D

Let $\Omega=[0,\pi]^{3}$. From Table 3, we take
$\phi_{k}(x,y,z)=\frac{2\sqrt{2}}{\pi^{3/2}}\sin(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)$
satisfying (A.1) with $\Gamma_{D}=\partial\Omega$ and $k=(k_{1},k_{2},k_{3})$.
Here, $k_{\mathrm{i}}$ are integers greater or equal to zero, and we let
$\mathcal{I}$ be the set of all such indices $k$. Notice that
$\|\nabla\phi_{k}\|_{[L^{2}(\Omega)]^{2}}=|k|$ for all $k\in\mathcal{I}$.
Using Proposition 3.1, we conclude that the set of functions
$(\bm{\phi}_{k})_{k\in\mathcal{I}}$, where each function $\bm{\phi}_{k}$ is
defined as

$\bm{\phi}_{k}:=\frac{1}{|k|}\nabla\phi_{k}=\frac{2\sqrt{2}}{\pi^{3/2}|k|}\begin{pmatrix}k_{1}\cos(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)\\\
k_{2}\sin(k_{1}x)\cos(k_{2}y)\sin(k_{3}z)\\\
k_{3}\sin(k_{1}x)\sin(k_{2}y)\cos(k_{3}z)\end{pmatrix},$

forms an orthonormal basis for $\nabla H^{1}_{0}(\Omega)$ in 3D.

###  B.IV Basis for $X_{0}(\Omega)$ in 3D

In cubic domains, i.e., $\Omega=[0,\pi]^{3}$, the construction of the TM and
TE modes is independent of the direction. Here, we take the $x$ direction
(with unit vector $\bm{e}_{1}$) to be the distinguished direction for the TM
and TE modes.

  * •

From Table 3, we take $\phi_{k}(x,y,z)=\sin(k_{1}x)\cos(k_{2}y)\cos(k_{3}z)$
satisfying (A.1) with $\Gamma_{D}=\\{0,\pi\\}\times\overline{\Omega^{*}}$,
$\Gamma_{N}=[0,\pi]\times\partial\Omega^{*}$ and $k=(k_{1},k_{2},k_{3})$.
Here, $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $\mathrm{i}=1,2$ or $3$,
and we let $\mathcal{I}$ be the set of all such indices $k$. We disregard the
cases where $k_{1}=0$ or $k_{2}=k_{3}=0$.

For all $k\in\mathcal{I}$, we define
$\hat{\bm{\psi}}_{k}^{\mathrm{TM}}:=\mathrm{curl}(\phi_{k}\bm{e}_{1})$. So, a
direct calculation shows that

$\|\hat{\bm{\psi}}_{k}^{\mathrm{TM}}\|_{{H}(\mathrm{curl},\Omega)}^{2}=c^{\prime}_{k}(1+|k|^{2})(k_{2}^{2}+k_{3}^{2}),$

where $c^{\prime}_{k}=\frac{\pi^{3}}{8}$ if $k_{2}>0$ and $k_{3}>0$, and
$c^{\prime}_{k}=\frac{\pi^{3}}{4}$ if $k_{2}=0$ or $k_{3}=0$. Thus, for all
$k\in\mathcal{I}$, the normalized TM modes are

$\bm{\psi}_{k}^{\mathrm{TM}}:=\frac{1}{\sqrt{c^{\prime}_{k}(1+|k|^{2})(k_{2}^{2}+k_{3}^{2})}}\begin{pmatrix}0\\\
-k_{3}\sin(k_{1}x)\cos(k_{2}y)\sin(k_{3}z)\\\
k_{2}\sin(k_{1}x)\sin(k_{2}y)\cos(k_{3}z)\end{pmatrix},\forall
k\in\mathcal{I}.$ (B.1)

  * •

From Table 3, we take $\phi_{k}(x,y,z)=\cos(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)$
satisfying (A.1) with $\Gamma_{N}=\\{0,\pi\\}\times\overline{\Omega^{*}}$,
$\Gamma_{D}=[0,\pi]\times\partial\Omega^{*}$ and $k=(k_{1},k_{2},k_{3})$.
Here, $k_{\mathrm{i}}\in{\mathbb{Z}}_{\geq 0}$ with $\mathrm{i}=1,2$ or $3$,
and we let $\mathcal{I}$ be the set of all such indices $k$. We disregard the
cases where $k_{2}=0$ or $k_{3}=0$.

For all $k\in\mathcal{I}$, we define
$\hat{\bm{\psi}}_{k}^{\mathrm{TE}}=\mathrm{curl}(\mathrm{curl}(\phi_{k}\bm{e}_{1}))$.
The norm of $\hat{\bm{\psi}}_{k}^{\mathrm{TE}}$ in $H(\mathrm{curl},\Omega)$
is

$\|\hat{\bm{\psi}}_{k}^{\mathrm{TE}}\|_{H(\mathrm{curl},\Omega)}^{2}=\frac{\pi^{3}}{8}|k|^{2}(1+|k|^{2})\left(|k|^{2}-k_{1}^{2}\right).$

This gives the TE modes as

$\bm{\psi}_{k}^{\mathrm{TE}}=\frac{2\sqrt{2}}{\pi^{3/2}\sqrt{(1+|k|^{2})|k|^{2}\left(|k|^{2}-k_{1}^{2}\right)}}\begin{pmatrix}(k_{2}^{2}+k_{3}^{2})\cos(k_{1}x)\sin(k_{2}y)\sin(k_{3}z)\\\
-k_{1}k_{2}\sin(k_{1}x)\cos(k_{2}y)\sin(k_{3}z)\\\
-k_{1}k_{3}\sin(k_{1}x)\sin(k_{2}y)\cos(k_{3}z)\end{pmatrix},\forall
k\in\mathcal{I}.$ (B.2)