remRemark propProposition lmLemma thmTheorem remarkRemark hypothesisHypothesis
claimClaim Divergence-free spectral algorithm for curl-curl equationL. Qin, C.
Sheng, and Z. Yang [][nocite]ex_supplement

# A highly efficient and accurate divergence-free spectral method for curl-
curl equation in two and three dimensions††thanks: Submitted to the editors
DATE.

Lechang Qin School of Mathematical Sciences, Shanghai Jiao Tong University,
Shanghai, 200240, China. ( (L. Qin<EMAIL_ADDRESS>Changtao
Sheng School of Mathematics, Shanghai University of Finance and Economics,
Shanghai 200433, China. The research of this author is partially supported by
by the National Natural Science Foundation of China (Nos. 12201385 and
12271365), Shanghai Pujiang Program 21PJ1403500, and the Fundamental Research
Funds for the Central Universities 2021110474. ( (C. Sheng)).
<EMAIL_ADDRESS>Zhiguo Yang Corresponding author. School of Mathematical
Sciences, MOE-LSC and CMA-Shanghai, Shanghai Jiao Tong University, Shanghai
200240, China. The research of this author is partially supported by the
National Natural Science Foundation of China (No. 12101399), the Shanghai
Sailing Program (No. 21YF1421000), the Strategic Priority Research Program of
Chinese Academy of Sciences (No. XDA25010402) and the Key Laboratory of
Scientific and Engineering Computing (Ministry of Education). ( (Z. Yang)).
<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a fast divergence-free spectral algorithm (FDSA) for
the curl-curl problem. Divergence-free bases in two and three dimensions are
constructed by using the generalized Jacobi polynomials. An accurate spectral
method with exact preservation of the divergence-free constraint point-wisely
is then proposed, and its corresponding error estimate is established. We then
present a highly efficient solution algorithm based on a combination of
matrix-free preconditioned Krylov subspace iterative method and a fully
diagonalizable auxiliary problem, which is derived from the spectral
discretisations of generalized eigenvalue problems of Laplace and biharmonic
operators. We rigorously prove that the dimensions of the invariant subspace
of the preconditioned linear system resulting from the divergence-free
spectral method with respect to the dominate eigenvalue $1$, are $(N-3)^{2}$
and $2(N-3)^{3}$ for two- and three-dimensional problems with $(N-1)^{2}$ and
$2(N-1)^{3}$ unknowns, respectively. Thus, the proposed method usually takes
only several iterations to converge, and astonishingly, as the problem size
(polynomial order) increases, the number of iterations will decrease, even for
highly indefinite system and oscillatory solutions. As a result, the
computational cost of the solution algorithm is only a small multiple of
$N^{3}$ and $N^{4}$ floating number operations for 2D and 3D problems,
respectively. Plenty of numerical examples for solving the curl-curl problem
with both constant and variable coefficients in two and three dimensions are
presented to demonstrate the accuracy and efficiency of the proposed method.

###### keywords:

Spectral method, curl-curl problem, divergence-free condition, property-
preserving discretization, matrix diagonalisation method

65N35, 65N22, 65F05, 35J05

##  1 Introduction

The curl-curl problem with a divergence-free constraint arises naturally in
many real-world applications related to electromagnetics, such as
semiconductor manufacturing, astrophysics, inertial confinement fusion,
photonic crystals, and more [7, 11, 22, 27, 42]. When dealing with problems of
this type, which involve variable coefficients, inhomogeneities, or nonlinear
couplings, obtaining analytical solutions can often be challenging. In such
scenarios, efficient and accurate simulations of these problems are of great
importance in gaining a deep understanding of the underlying physical
mechanisms and in addressing practical engineering problems.

Nevertheless, the aforementioned problem is notoriously difficult to solve
numerically. It poses challenges and difficulties that can be attributed to at
least three different aspects. Firstly, even a small violation of the
divergence-free constraint can lead to a significant deviation of the
numerical solution from the actual solution [3, 4]. Secondly, the indefinite
nature of the system can cause the solution algorithm to converge slowly [2,
28]. Finally, an excessively large number of degrees of freedom may be
required to accurately approximate highly oscillatory solutions [1, 14]. On
the one hand, it is essential to construct approximation spaces that ensure
well-posedness of the discrete problems and adherence to the divergence-free
constraint (either strongly or weakly). On the other hand, it is equally
important to devise efficient algorithms for solving the resulting indefinite
linear system, particularly for solutions with high oscillation.

There has been a long standing interest in developing numerical algorithms for
the curl-curl problem that preserve the divergence-free constraint at the
discrete level. The common methods of choice include, but not limited to,
staggered-grid finite difference Yee algorithm [37, 40], the vector potential
method [10, 21, 31], the divergence-regularzation method [8, 12, 18], and the
Lagrange multiplier method [28]. These methods require either the solution of
an extra equation to project the numerical solution to certain divergence-free
space, or augmenting/reformulating the original system with new variables,
leading to equations with higher spatial order or large coupled system to be
solved. In contrast to the aforementioned methods, the method of constructing
divergence-free basis can not only preserve the divergence-free constraint in
a straightforward way, but also significantly reduces the degrees of freedom
to be solved. In [39], a weakly divergence-free finite element basis was
constructed for the Stokes equation. In [15], Guo and Jiao proposed a
divergence-free basis functions for $n$-dimensional Navier-Stokes equation.
Cai et. al in [6] constructed a divergence-free $\boldsymbol{H}({\rm div})$
hierarchical basis for the MHD system. Local divergence-free methods based on
discontinuous-Garlerkin (DG) or weak-Galerkin (WG) frameworks were also
proposed in [19, 25, 29]. However, to the best of the authors’ knowledge,
there are no fast divergence-free spectral algorithm available for the curl-
curl problem.

The main contributions of this paper are twofold. Firstly, we derive a novel
divergence-free spectral approximation bases using generalized Jacobi
polynomials for the curl-curl problem in both two and three dimensions. We
conduct a rigorous error estimate to justify the spectral accuracy of this
method. Secondly, we propose a highly efficient solution algorithm, which
draws inspirations from spectral discretizations of the eigenvalue problems of
Laplace and biharmonic operators and the fast diagonalisation method for
Poisson-type equations. Lemma 2.3 demonstrates that the mass and stiffness
matrices resulted from the proposed method are closely related, inspiring the
construction of a fast diagonalizable auxiliary problem. A highly efficient
solution algorithm is then proposed by a proper combination of the auxiliary
problem and matrix-free preconditioned Krylov subspace iterative method. We
prove rigorously that the dimensions of the invariant subspace of the
preconditioned system with respect to eigenvalue 1 are precisely $(N-3)^{2}$
and $2(N-3)^{3}$ for two- and three-dimensional problems with $(N-1)^{2}$ and
$2(N-1)^{3}$ unknowns, respectively. As a result, the proposed fast
divergence-free spectral algorithm usually takes only several iterations to
converge, and as the problem size increases, the number of iterations will
decrease, even for highly indefinite system and oscillatory solutions. Thanks
to these properties, the presented divergence-free spectral method and
corresponding fast solution algorithms are very competitive and
computationally attractive.

The remainder of the paper is structured as follows: In Section 2, we
introduce the divergence-free spectral method and conduct the corresponding
error estimate for the curl-curl problem in two dimensions. Then, we present a
highly efficient solution algorithm that combines a fast diagonalizable
auxiliary problem and matrix-free preconditioned Krylov subspace iterative
method. Theoretical results on the dimensions of the invariant subspace of the
preconditioned system with respect to the dominate eigenvalue 1 are also
presented. In Section 3, we expand upon the proposed divergence-free spectral
method and fast solution algorithm to the three-dimensional case. In Section
4, we present ample numerical results to demonstrate the efficiency and
accuracy of the proposed method, particularly for indefinite systems and
highly oscillatory solutions. Finally, in Section 5, we conclude the
discussion with closing remarks.

##  2 A fast divergence-free spectral algorithm for 2D curl-curl problem

We first introduce some notation and review some relevant properties of the
generalised Jacobi polynomials (cf. [16, 33]), which will be an indispensable
building block for the construction of divergence-free basis.

###  2.1 Generalised Jacobi polynomials

For $\alpha,\beta>-1$, the classical Jacobi polynomials
$P_{n}^{(\alpha,\beta)}(\xi)$ are mutually orthogonal with respect to the
weight function $\omega^{\alpha,\beta}(\xi)=(1-\xi)^{\alpha}(1+\xi)^{\beta}$
on $\Lambda=(-1,1)$, that is

(1)
$\int_{-1}^{1}P_{n}^{(\alpha,\beta)}P_{m}^{(\alpha,\beta)}\omega^{\alpha,\beta}\,\mathrm{d}\xi=\frac{2^{\alpha+\beta+1}\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{(2n+\alpha+\beta+1)\Gamma(n+1)\Gamma(n+\alpha+\beta+1)}\delta_{nm},$

where $\delta_{nm}$ is the Kronecker symbol. The classical Jacobi polynomials
can be generalized to cases with general $\alpha,\beta\in\mathbb{R}$, that is,
the generalised Jacobi functions (cf. [16, 33]). The flexibility afforded by
the parameters $\alpha,\beta$ allows us to design a suitable basis for exact
preservation of the divergence-free constraint. Hereafter, we adopt the
following generalized Jacobi polynomials:

(2) $\displaystyle P_{n}^{(0,0)}(\xi)=L_{n}(\xi),\;\;\xi\in\Lambda,\;n\geq 0,$
(3) $\displaystyle
P_{n}^{(-1,-1)}(\xi)=\frac{\xi^{2}-1}{4}P_{n-2}^{(1,1)}(\xi)=\frac{n-1}{2(2n-1)}\big{(}L_{n}(\xi)-L_{n-2}(\xi)\big{)},\;\;n\geq
2,$

where $L_{n}(\xi)$ also known as the Legendre polynomial of order $n$. In view
of (1), these polynomials satisfy the following orthogonality properties:

(4) $\displaystyle\int_{-1}^{1}L_{m}(\xi)L_{n}(\xi)\,{\rm
d}\xi=\frac{2}{2n+1}\delta_{mn},$ (5)
$\displaystyle\int_{-1}^{1}P_{m}^{(-1,-1)}(\xi)P_{n}^{(-1,-1)}(\xi)(1-\xi^{2})^{-1}\,{\rm
d}\xi=\frac{n-1}{2n(2n-1)}\delta_{mn}.$

One prominent feature of the generalized Jacobi polynomials is the derivative
property

(6)
$\big{(}P_{n}^{(-1,-1)}\big{)}^{\prime}(\xi)=\frac{n-1}{2}L_{n-1}(\xi),\;\;n\geq
2,$

which plays a significant role in constructing the high-order divergence-free
spectral bases.

###  2.2 The divergence-free approximation scheme and its error estimate

We are concerned with the following curl-curl problem in two and three
dimensions

(7a)
$\displaystyle\nabla\times\nabla\times{\boldsymbol{u}}+\kappa\,{\boldsymbol{u}}={\boldsymbol{f}}\;\;{\rm
in}\;\;\Omega\subset\mathbb{R}^{d},$ (7b)
$\displaystyle{\boldsymbol{n}}\cdot\boldsymbol{u}=0,\quad\boldsymbol{n}\times(\nabla\times{\boldsymbol{u}})=0\;\;{\rm
on}\;\;\partial\Omega,$

where $d=2,3$, $\boldsymbol{u}(\boldsymbol{x})$ usually represents the
magnetic field, $\kappa$ is a constant, $\boldsymbol{f}(\boldsymbol{x})$ is
the source term satisfying $\nabla\cdot\boldsymbol{f}=0$, and $\boldsymbol{n}$
is the unit outward vector normal to the boundary $\partial\Omega$.
$\nabla\times$ and $\nabla\cdot$ are the usual curl and divergence operators,
respectively. Define $\mathbb{X}^{{\rm
div}0}(\Omega):=\big{\\{}\boldsymbol{u}\in H({\rm curl};\Omega)\cap H_{0}({\rm
div},\Omega),\;\nabla\cdot\boldsymbol{u}=0\big{\\}}$, where the Sobolev space
$H({\rm curl};\Omega)$ and $H_{0}({\rm div};\Omega)$ are

$\displaystyle H({\rm
curl};\Omega)=\big{\\{}\boldsymbol{u}\in(L^{2}(\Omega))^{d},\nabla\times\boldsymbol{u}\in(L^{2}(\Omega))^{d}\big{\\}},\;\;$
$\displaystyle H_{0}({\rm
div};\Omega)=\big{\\{}\boldsymbol{u}\in(L^{2}(\Omega))^{d},\nabla\cdot\boldsymbol{u}\in
L^{2}(\Omega),\;\;\boldsymbol{n}\cdot\boldsymbol{u}=0\;{\rm
at}\;\partial\Omega\big{\\}}.$

The weak formulation of equation (7) reads: find
$\boldsymbol{u}\in\mathbb{X}^{{\rm div}0}(\Omega)$ such that

(8)
$(\nabla\times\boldsymbol{u},\nabla\times\boldsymbol{v})+\kappa(\boldsymbol{u},\boldsymbol{v})=(\boldsymbol{f},\boldsymbol{v}),\quad\forall\boldsymbol{v}\in\mathbb{X}^{{\rm
div}0}(\Omega).$

In what follows, we let $\Omega:=\Lambda^{2}=(-1,1)^{2}$. We propose the
divergence-free approximation space in two dimensions as follows.

######  Proposition 1.

$\mathbb{X}_{N,2d}^{{\rm div}0}(\Lambda^{2})$ is a divergence-free
approximation space for $\mathbb{X}^{{\rm div}0}(\Lambda^{2})$ taking the form

(9) $\mathbb{X}_{N,2d}^{{\rm div}0}(\Lambda^{2})={\rm
span}\big{\\{}\boldsymbol{\Phi}_{m,n},\;1\leq m,n\leq N-1\big{\\}},$

where
${\boldsymbol{\Phi}}_{m,n}(\boldsymbol{x})=\big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2}),-\phi_{m}(x_{1})\psi_{n+1}(x_{2})\big{)}^{\intercal}$
and

(10)
$\displaystyle\psi_{m+1}(\xi)=\frac{\sqrt{2(2m+1)}}{m}P_{m+1}^{(-1,-1)}(\xi),\quad\phi_{m}(\xi)=\sqrt{\frac{2m+1}{2}}L_{m}(\xi).$

######  Proof 2.1.

With the help of the derivative property (6), one readily verifies that

(11) $\psi_{m+1}^{\prime}(\xi)=\phi_{m}(\xi),$

from where we obtain $\nabla\cdot\boldsymbol{\Phi}_{m,n}(\boldsymbol{x})=0.$

Correspondingly, the approximation scheme for the weak formulation (8) in two
dimensions takes the form: find $\boldsymbol{u}_{N}\in\mathbb{X}_{N,2d}^{{\rm
div}0}(\Lambda^{2})$ such that

(12)
$\mathcal{A}(\boldsymbol{u}_{N},\boldsymbol{v}):=(\nabla\times{\boldsymbol{u}}_{N},\nabla\times{\boldsymbol{v}})+\kappa({\boldsymbol{u}}_{N},\boldsymbol{v})=({\boldsymbol{f}},\boldsymbol{v}),\quad\forall\boldsymbol{v}\in\mathbb{X}_{N,2d}^{{\rm
div}0}(\Lambda^{2}).$

In order to conduct the error estimate for the proposed scheme, let us
introduce some notation. Let $\mathbb{N}$ be the set of all positive integers,
and let $\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$. We denote the $l^{1}$-norm of
$\boldsymbol{\xi}$ in $\mathbb{R}^{d}$ by $|\boldsymbol{\xi}|_{1}$, and the
set of all algebraic polynomials of degree at most $N$ by $\mathcal{P}_{N}$.
For the sake of simplicity, we also assume $\kappa>0$. Wavenumber-explicit
priori estimate and error analysis for curl-curl problem with large negative
$\kappa$ can be conducted by Rellich identities, but much more involved. We
refer the readers to relevant works for Helmholtz equation ([9, 34, 35]) and
time-harmonic Maxwell’s equations ([20, 26]).

The bilinear form (8) satisfies the following continuity and coercivity:

(13) $\begin{split}&\mathcal{A}(\boldsymbol{u},\boldsymbol{v})\leq
C\|\boldsymbol{u}\|_{H({\rm curl};\Omega)}\|\boldsymbol{v}\|_{H({\rm
curl};\Omega)},\quad{\rm
for}\quad\boldsymbol{u},\boldsymbol{v}\in\mathbb{X}^{{\rm div}0}(\Omega),\\\
&\mathcal{A}(\boldsymbol{u},\boldsymbol{u})\geq\alpha\|\boldsymbol{u}\|^{2}_{H({\rm
curl};\Omega)},\qquad\qquad{\rm for}\quad\boldsymbol{u}\in\mathbb{X}^{{\rm
div}0}(\Omega),\end{split}$

where the positive constants $C$ and $\alpha$ are independent of
$\boldsymbol{u}.$

By the Helmholtz-Hodge decomposition theory, $\boldsymbol{u}$ can be
decomposed into two components as

(14) $\boldsymbol{u}=\nabla\times\phi+\nabla\psi,$

which leads to

(15)
$\nabla\cdot\boldsymbol{u}=\Delta\psi,\;\;\;\nabla\times\boldsymbol{u}=\nabla\times\nabla\times\phi.$

By the divergence-free condition, we obtain that
$\boldsymbol{u}=\nabla\times\phi.$ We find from (11) that

(16)
$\begin{split}&\nabla\times\big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2}),-\phi_{m}(x_{1})\psi_{n+1}(x_{2})\big{)}^{\intercal}\\\
&=-\phi_{m}^{\prime}(x_{1})\psi_{n+1}(x_{2})-\psi_{m+1}(x_{1})\phi_{n}^{\prime}(x_{2})\\\
&=-\psi_{m+1}^{\prime\prime}(x_{1})\psi_{n+1}(x_{2})-\psi_{m+1}(x_{1})\psi_{n+1}^{\prime\prime}(x_{2})\\\
&=-\Delta(\psi_{m+1}(x_{1})\psi_{n+1}(x_{2})),\end{split}$

which implies $\phi|_{\partial\Omega}=0.$

To measure the error, we further define the $d$-dimensional Jacobi-weighted
Sobolev space below

$B^{m}(\Omega):=\left\\{\phi:\partial_{\boldsymbol{x}}^{\boldsymbol{k}}\phi\in
L^{2}_{\boldsymbol{\omega}^{\boldsymbol{k},\boldsymbol{k}}}(\Omega),0\leq|\boldsymbol{k}|_{1}\leq
m\right\\},\quad\forall m\in\mathbb{N}_{0},$

equipped with the norm and semi-norm

$\displaystyle\|\phi\|_{B^{m}(\Omega)}=\bigg{(}\sum_{0\leq|\boldsymbol{k}|_{1}\leq
m}\left\|\partial_{\boldsymbol{x}}^{k}\phi\right\|_{\boldsymbol{\omega}^{k,k}}^{2}\bigg{)}^{1/2},\;\;\;|\phi|_{B^{m}(\Omega)}=\bigg{(}\sum_{j=1}^{d}\left\|\partial_{x_{j}}^{m}\phi\right\|_{\boldsymbol{\omega}^{me_{j},me_{j}}}\bigg{)}^{1/2}.$

We consider the $H^{2}$-orthogonal projection
$\pi^{2}_{N}:H^{1}_{0}(\Lambda^{2})\cap
H^{2}(\Lambda^{2})\rightarrow(\mathcal{P}_{N})^{2}\cap
H^{1}_{0}(\Lambda^{2})\cap H^{2}(\Lambda^{2})$, defined by

(17)
$(\Delta(\phi-\pi^{2}_{N}\phi),\Delta\varphi)=0,\;\;\;\varphi\in(\mathcal{P}_{N})^{2}\cap
H^{1}_{0}(\Lambda^{2})\cap H^{2}(\Lambda^{2}).$

Then, we can derive the following approximation result of $\pi^{2}_{N}$.

######  Lemma 2.2.

If $\phi\in H^{1}_{0}(\Lambda^{2})\cap H^{2}(\Lambda^{2})$ and $\Delta\phi\in
B^{r}(\Lambda^{2})$, there holds

(18) $\|\Delta(\pi_{N}^{2}\phi-\phi)\|_{\Lambda^{2}}\leq
cN^{-r}|\Delta\phi|_{B^{r}(\Lambda^{2})}.$

######  Proof 2.3.

In view of (17), we have

(19) $\pi_{N}^{2}=\pi_{N}^{(1)}\circ\pi_{N}^{(2)},$

where the one dimensional orthogonal projection
$\pi^{(j)}_{N}:H^{2}(\Lambda)\rightarrow\mathcal{P}_{N}$, $j=1,2$ such that
$\pi^{(j)}_{N}\phi(\pm 1)=\phi(\pm 1)=0$ and $(\pi^{(j)}_{N}\phi)^{\prime}(\pm
1)=\phi^{\prime}(\pm 1)$. According to [33, (4.82)], for any $\phi\in
H^{2}(\Lambda)$ and $\partial^{2}_{x}\phi\in
B^{r-2}_{\ast}(\Lambda):=\\{\phi:\partial^{k}_{x}\phi\in
L^{2}_{\omega^{k,k}}(\Lambda),0\leq k\leq r-2\\}$, we have

(20) $\|\pi_{N}^{(j)}\phi-\phi\|_{H^{\mu}(\Lambda)}\leq
cN^{\mu-r}(\|\phi\|_{H^{2}}+\|\partial_{x}^{r}\phi\|_{\omega^{r-2,r-2}}),\;\;\;\mu\in[0,2].$

Using integration by parts and the Cauchy-Schwarz inequality, we obtain that

(21)
$\begin{split}\|\Delta(\pi_{N}^{2}\phi-\phi)\|^{2}_{\Lambda^{2}}&=\|\partial^{2}_{x_{1}}(\pi_{N}^{2}\phi-\phi)\|^{2}_{\Lambda^{2}}+\|\partial^{2}_{x_{2}}(\pi_{N}^{2}u-u)\|^{2}_{\Lambda^{2}}\\\
&\quad+2\int_{\Lambda^{2}}\partial^{2}_{x_{1}}(\pi_{N}^{2}\phi-\phi)\partial^{2}_{x_{2}}(\pi_{N}^{2}\phi-\phi){\rm
d}\boldsymbol{x}\\\ &\leq
c\|\partial^{2}_{x_{1}}(\pi_{N}^{2}\phi-\phi)\|^{2}_{\Lambda^{2}}+c\|\partial^{2}_{x_{2}}(\pi_{N}^{2}u-u)\|^{2}_{\Lambda^{2}}:=I_{1}+I_{2}.\end{split}$

We only consider $I_{1}$, since the other terms in $I_{2}$ can be derived in a
similar fashion. By virtue of (20) with $\mu=2$ and (19), we derive that

$\begin{split}I_{1}&\leq
2\Big{(}\|\partial^{2}_{x_{1}}(\pi_{N}^{(1)}\phi-\phi)\|^{2}_{\Lambda^{2}}+\|\partial_{x_{1}}^{2}\pi_{N}^{(1)}\circ(\pi_{N}^{(2)}\phi-\phi)\|^{2}_{\Lambda^{2}}\\\
&\quad\quad+\|\partial^{2}_{x_{2}}(\pi_{N}^{(2)}\phi-\phi)\|^{2}_{\Lambda^{2}}+\|\partial_{x_{2}}^{2}\pi_{N}^{(2)}\circ(\pi_{N}^{(1)}\phi-\phi)\|^{2}_{\Lambda^{2}}\Big{)}\\\
&\leq
cN^{-2r}|\Delta\phi|^{2}_{B^{r}(\Lambda^{2})}+2\|\pi_{N}^{(2)}(\partial_{x_{1}}^{2}\phi)-(\partial_{x_{1}}^{2}\phi)\|^{2}_{\Lambda^{2}}+\|\pi_{N}^{(1)}(\partial_{x_{2}}^{2}\phi)-(\partial_{x_{2}}^{2}\phi)\|^{2}_{\Lambda^{2}}\\\
&\leq cN^{-2r}|\Delta\phi|^{2}_{B^{r}(\Lambda^{2})}.\end{split}$

A combination of the above facts leads to the desired result.

n order to characterize the error bound of the numerical solution, we need to
introduce two one-dimensional orthogonal projections. More specifically, we
define the $L^{2}(\Lambda)$-orthogonal projection
$\pi_{N}:L^{2}(\Lambda)\rightarrow\mathcal{P}_{N}$ such that

(22) $\left(\pi_{N}u-u,v\right)_{L^{2}(\Lambda)}=0,\quad\forall
v\in\mathcal{P}_{N},$

or equivalently,

(23) $\left(\pi_{N}u\right)(x)=\sum_{n=0}^{N}\hat{u}_{n}\phi_{n}(x).$

Consider the orthogonal projection
${\pi}_{N}^{1,0}:H^{1}_{0}(\Lambda)\rightarrow\mathcal{P}_{N}^{0}:=\\{u\in\mathcal{P}_{N}:u(\pm
1)=0\\}$, defined by

(24)
$\big{(}\partial_{x}(\pi^{1,0}_{N}u-u),\partial_{x}v\big{)}_{L^{2}(\Lambda)}=0,\quad\forall
v\in\mathcal{P}_{N},$

or equivalently,

(25)
$\big{(}\pi^{1,0}_{N}u\big{)}(x)=\sum_{n=1}^{N-1}\tilde{u}_{n}\psi_{n+1}(x).$

Then the basic approximation results for (31) and (33) are stated as follows
(cf. [33]). If $u\in H^{1}_{0}(\Lambda)$ and $\partial_{x}^{r}u\in
L^{2}_{\omega^{r-1}}(\Lambda)$, then for $1\leq r\leq N+1,$

(26) $\begin{split}&\big{\|}u-\pi_{N}^{1,0}u\big{\|}_{\mu}\leq
cN^{\mu-r}\|\partial_{x}^{r}u\|_{\omega^{r-1,r-1}},\;\;\;{\rm
with}\;\;\mu=0,1.\end{split}$

Moreover, if $\partial_{x}^{r}u\in L^{2}_{\omega^{r}}(\Lambda)$, then for
$1\leq r\leq N+1,$

(27) $\begin{split}&\big{\|}u-\pi_{N}u\big{\|}\leq
cN^{-r}\|\partial_{x}^{r}u\|_{\omega^{r,r}},\end{split}$

and

(28) $\begin{split}&\big{\|}\partial_{x}(u-\pi_{N}u)\big{\|}\leq
cN^{\frac{3}{2}-r}\|\partial_{x}^{r}u\|_{\omega^{r,r}}\,.\end{split}$

In view of (LABEL:eq:_unt), (32), and (34), for
$\boldsymbol{u}\in\mathbb{X}^{{\rm div}0}(\Lambda^{2})$, we define

(29)
$\begin{split}&\Pi_{N}\boldsymbol{u}(\boldsymbol{x}):=(\pi^{1,0}_{N_{1}}\pi_{N_{2}-1}u_{1},\pi_{N_{1}-1}\pi^{1,0}_{N_{2}}u_{2})^{\intercal}\\\
&=\sum_{m,n=1}^{N-1}\tilde{u}_{mn}\big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2}),-\phi_{m}(x_{1})\psi_{n+1}(x_{2})\big{)}^{\intercal}=\sum_{m,n=1}^{N-1}\tilde{u}_{mn}\boldsymbol{\Phi}_{mn}(\boldsymbol{x}),\end{split}$

where $\pi_{N_{1}}$ and $\pi_{N_{2}}$ denote the one dimensional orthogonal
projections in $x_{1}$\- and $x_{2}$-direction, respectively. One verifies
readily that $\nabla\cdot\Pi_{N}\boldsymbol{u}(\boldsymbol{x})=0.$ Note that
the subscript of $N_{1}$ and $N_{2}$ are only used to distinguish the degree
in $x_{1}$\- and $x_{2}$-direction. In what follows, we take $N=N_{1}=N_{2}$
for simplification.

######  Theorem 2.

Let $\boldsymbol{u},\boldsymbol{u}_{N}$ be respectively the solutions of (8)
and (12). If $\boldsymbol{u}\in\mathbb{X}^{{\rm div}0}(\Lambda^{2})$ and
$\nabla\times\boldsymbol{u}\in B^{r}(\Lambda^{2})$ with $r\geq 0$, then we
have

(30)
$\|\nabla\times(\boldsymbol{u}-\boldsymbol{u}_{N})\|_{\Lambda^{2}}+\kappa\|\boldsymbol{u}-\boldsymbol{u}_{N}\|_{\Lambda^{2}}\leq
cN^{-r}|\nabla\times\boldsymbol{u}|_{B^{r}(\Lambda^{2})}.$

where $c$ is a positive constant independent of $r$, $N$, and
$\boldsymbol{u}$.

######  Proof 2.4.

Applying the first Strang Lemma (see, e.g., [36]) to (12), we can obtain

$\begin{split}\|\nabla\times(\boldsymbol{u}-\boldsymbol{u}_{N})\|_{\Lambda^{2}}+\kappa\|\boldsymbol{u}-\boldsymbol{u}_{N}\|_{\Lambda^{2}}&\lesssim\inf_{\boldsymbol{w}_{N}\in\mathbb{X}^{{\rm
div}0}_{N,2d}}\|\nabla\times(\boldsymbol{u}-\boldsymbol{w}_{N})\|_{\Lambda^{2}}+\kappa\|\boldsymbol{u}-\boldsymbol{w}_{N}\|_{\Lambda^{2}}.\end{split}$

Taking
$\boldsymbol{w}_{N}=\Pi_{N}\boldsymbol{u}(\boldsymbol{x})=\nabla\times\pi^{2}_{N}\phi$
in above, we obtain from (15), (18), and Poincaré inequality that

$\begin{split}&\|\nabla\times(\boldsymbol{u}-\Pi_{N}\boldsymbol{u})\|_{\Lambda^{2}}+\kappa\|\boldsymbol{u}-\Pi_{N}\boldsymbol{u}\|_{\Lambda^{2}}\\\
&=\|\nabla\times\nabla\times(\phi-\pi^{2}_{N}\phi)\|_{\Lambda^{2}}+\kappa\|\nabla\times(\phi-\pi^{2}_{N}\phi)\|_{\Lambda^{2}}\\\
&\leq\|\nabla\times\nabla\times(\phi-\pi^{2}_{N}\phi)\|_{\Lambda^{2}}=\|\Delta(\phi-\pi^{2}_{N}\phi)\|_{\Lambda^{2}}\\\
&\leq
cN^{-r}|\Delta\phi|_{B^{r}(\Lambda^{2})}=cN^{-r}|\nabla\times\boldsymbol{u}|_{B^{r}(\Lambda^{2})}.\end{split}$

This ends the proof.

n order to characterize the error bound of the numerical solution, we need to
introduce two one-dimensional orthogonal projections. More specifically, we
define the $L^{2}(\Lambda)$-orthogonal projection
$\pi_{N}:L^{2}(\Lambda)\rightarrow\mathcal{P}_{N}$ such that

(31) $\left(\pi_{N}u-u,v\right)_{L^{2}(\Lambda)}=0,\quad\forall
v\in\mathcal{P}_{N},$

or equivalently,

(32) $\left(\pi_{N}u\right)(x)=\sum_{n=0}^{N}\hat{u}_{n}\phi_{n}(x).$

Consider the orthogonal projection
${\pi}_{N}^{1,0}:H^{1}_{0}(\Lambda)\rightarrow\mathcal{P}_{N}^{0}:=\\{u\in\mathcal{P}_{N}:u(\pm
1)=0\\}$, defined by

(33)
$\big{(}\partial_{x}(\pi^{1,0}_{N}u-u),\partial_{x}v\big{)}_{L^{2}(\Lambda)}=0,\quad\forall
v\in\mathcal{P}_{N},$

or equivalently,

(34)
$\big{(}\pi^{1,0}_{N}u\big{)}(x)=\sum_{n=1}^{N-1}\tilde{u}_{n}\psi_{n+1}(x).$

Then the basic approximation results for (31) and (33) are stated as follows
(cf. [33]). If $u\in H^{1}_{0}(\Lambda)$ and $\partial_{x}^{r}u\in
L^{2}_{\omega^{r-1}}(\Lambda)$, then for $1\leq r\leq N+1,$

(35) $\begin{split}&\big{\|}u-\pi_{N}^{1,0}u\big{\|}_{\mu}\leq
cN^{\mu-r}\|\partial_{x}^{r}u\|_{\omega^{r-1,r-1}},\;\;\;{\rm
with}\;\;\mu=0,1.\end{split}$

Moreover, if $\partial_{x}^{r}u\in L^{2}_{\omega^{r}}(\Lambda)$, then for
$1\leq r\leq N+1,$

(36) $\begin{split}&\big{\|}u-\pi_{N}u\big{\|}\leq
cN^{-r}\|\partial_{x}^{r}u\|_{\omega^{r,r}},\end{split}$

and

(37) $\begin{split}&\big{\|}\partial_{x}(u-\pi_{N}u)\big{\|}\leq
cN^{\frac{3}{2}-r}\|\partial_{x}^{r}u\|_{\omega^{r,r}}\,.\end{split}$

To measure the error, we further define the $d$-dimensional Jacobi-weighted
Sobolev space below

$B^{m}(\Omega):=\left\\{u:\partial_{\boldsymbol{x}}^{\boldsymbol{k}}u\in
L^{2}_{\boldsymbol{\omega}^{\boldsymbol{k},\boldsymbol{k}}}(\Omega),0\leq|\boldsymbol{k}|_{1}\leq
m\right\\},\quad\forall m\in\mathbb{N}_{0},$

equipped with the norm and semi-norm

$\displaystyle\|u\|_{B^{m}(\Omega)}=\bigg{(}\sum_{0\leq|\boldsymbol{k}|_{1}\leq
m}\left\|\partial_{\boldsymbol{x}}^{k}u\right\|_{\boldsymbol{\omega}^{k,k}}^{2}\bigg{)}^{1/2},\;\;\;|u|_{B^{m}(\Omega)}=\bigg{(}\sum_{j=1}^{d}\left\|\partial_{x_{j}}^{m}u\right\|_{\boldsymbol{\omega}^{me_{j},me_{j}}}\bigg{)}^{1/2}.$

In view of (LABEL:eq:_unt), (32), and (34), for
$\boldsymbol{u}\in\mathbb{X}^{{\rm div}0}(\Lambda^{2})$, we define

(38)
$\begin{split}&\Pi_{N}\boldsymbol{u}(\boldsymbol{x}):=(\pi^{1,0}_{N_{1}}\pi_{N_{2}-1}u_{1},\pi_{N_{1}-1}\pi^{1,0}_{N_{2}}u_{2})^{\intercal}\\\
&=\sum_{m,n=1}^{N-1}\tilde{u}_{mn}\big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2}),-\phi_{m}(x_{1})\psi_{n+1}(x_{2})\big{)}^{\intercal}=\sum_{m,n=1}^{N-1}\tilde{u}_{mn}\boldsymbol{\Phi}_{mn}(\boldsymbol{x}),\end{split}$

where $\pi_{N_{1}}$ and $\pi_{N_{2}}$ denote the one dimensional orthogonal
projections in $x_{1}$\- and $x_{2}$-direction, respectively. One verifies
readily that $\nabla\cdot\Pi_{N}\boldsymbol{u}(\boldsymbol{x})=0.$ Note that
the subscript of $N_{1}$ and $N_{2}$ are only used to distinguish the degree
in $x_{1}$\- and $x_{2}$-direction. In what follows, we take $N=N_{1}=N_{2}$
for simplification.

######  Theorem 3.

Let $\boldsymbol{u},\boldsymbol{u}_{N}$ be respectively the solutions of (8)
and (12). If $\boldsymbol{u}\in\mathbb{X}^{{\rm
div}0}(\Lambda^{2})\cup(B^{r}(\Lambda^{2}))^{2}$ with $r\geq 0$, then we have

(39) $\|\boldsymbol{u}-\boldsymbol{u}_{N}\|_{H({\rm curl};\Lambda^{2})}\leq
cN^{\frac{3}{2}-r}\|\boldsymbol{u}\|_{(B^{r}(\Lambda^{2}))^{2}}.$

where $c$ is a positive constant independent of $r$, $N$, and
$\boldsymbol{u}$.

######  Proof 2.5.

Applying the first Strang Lemma (see, e.g., [36]) to (12), we can obtain

$\begin{split}\|\boldsymbol{u}-\boldsymbol{u}_{N}\|_{H({\rm
curl};\Lambda^{2})}&\lesssim\inf_{\boldsymbol{w}_{N}\in\mathbb{X}^{{\rm
div}0}_{N,2d}}\|\boldsymbol{u}-\boldsymbol{w}_{N}\|_{H({\rm
curl};\Lambda^{2})}.\end{split}$

Taking $\boldsymbol{w}_{N}=\Pi_{N}\boldsymbol{u}(\boldsymbol{x})$ in above, we
obtain from (LABEL:eq:_curl), (38), and triangle inequality that

$\begin{split}&\|\boldsymbol{u}-\Pi_{N}\boldsymbol{u}\|_{H({\rm
curl};\Lambda^{2})}=\|\nabla\times(\boldsymbol{u}-\Pi_{N}\boldsymbol{u})\|_{L^{2}(\Lambda^{2})}+\|\boldsymbol{u}-\Pi_{N}\boldsymbol{u}\|_{L^{2}(\Lambda^{2})}\\\
&\leq\big{\|}\partial_{x_{1}}(u_{2}-\pi_{N_{1}-1}\pi^{1,0}_{N_{2}}u_{2})\big{\|}_{L^{2}(\Lambda^{2})}+\big{\|}\partial_{x_{2}}(\pi^{1,0}_{N_{1}}\pi_{N_{2}-1}u_{1}-u_{1})\big{\|}_{L^{2}(\Lambda^{2})}\\\
&\quad+\|u_{1}-\pi^{1,0}_{N_{1}}\pi_{N_{2}-1}u_{1}\|_{L^{2}(\Lambda^{2})}+\|u_{2}-\pi_{N_{1}-1}\pi^{1,0}_{N_{2}}u_{2}\|_{L^{2}(\Lambda^{2})}\\\
&:=I_{1}+I_{2}+I_{3}+I_{4}.\end{split}$

We find from the above and (LABEL:pnerr1)-(LABEL:pnerr3) that

$\begin{split}I_{1}&\leq\|\partial_{x_{1}}(u_{2}-\pi_{N_{1}-1}u_{2})\|_{L^{2}(\Lambda^{2})}+\|\partial_{x_{1}}(\pi_{N_{1}-1}u_{2}-\pi_{N_{1}-1}\pi^{1,0}_{N_{2}}u_{2})\|_{L^{2}(\Lambda^{2})}\\\
&\leq
cN^{3/2-r}\|u_{2}\|_{B^{r}(\Lambda^{2})}+\|\partial_{x_{1}}(\pi_{N_{1}-1}-\mathcal{I})(u_{2}-\pi^{1,0}_{N_{2}}u_{2})\|_{L^{2}(\Lambda^{2})}\\\
&\quad+\|\partial_{x_{1}}(u_{2}-\pi^{1,0}_{N_{2}}u_{2})\|_{L^{2}(\Lambda^{2})}\\\
&\leq cN^{3/2-r}\|u_{2}\|_{B^{r}(\Lambda^{2})},\end{split}$

and

$\begin{split}I_{2}&\leq\|\partial_{x_{2}}(u_{1}-\pi_{N_{2}-1}u_{1})\|_{L^{2}(\Lambda^{2})}+\|\partial_{x_{2}}(\pi_{N_{2}-1}u_{1}-\pi^{1,0}_{N_{1}}\pi_{N_{2}-1}u_{1})\|_{L^{2}(\Lambda^{2})}\\\
&\leq
cN^{3/2-r}\|u_{1}\|_{B^{r}(\Lambda^{2})}+\|\partial_{x_{2}}(\pi_{N_{2}-1}-\mathcal{I})(u_{1}-\pi^{1,0}_{N_{1}}u_{1})\|_{L^{2}(\Lambda^{2})}\\\
&\quad+\|\partial_{x_{2}}(u_{1}-\pi^{1,0}_{N_{1}}u_{1})\|_{L^{2}(\Lambda^{2})}\\\
&\leq cN^{3/2-r}\|u_{1}\|_{B^{r}(\Lambda^{2})}.\end{split}$

Similarly, we have

$\begin{split}I_{3}&\leq\|u_{1}-\pi_{N_{2}-1}u_{1}\|_{L^{2}(\Lambda^{2})}+\|\pi_{N_{2}-1}u_{1}-\pi^{1,0}_{N_{1}}\pi_{N_{2}-1}u_{1}\|_{L^{2}(\Lambda^{2})}\\\
&\leq
cN^{-r}\|u_{1}\|_{B^{r}(\Lambda^{2})}+\|u_{1}-\pi^{1,0}_{N_{1}}u_{1}\|_{L^{2}(\Lambda^{2})}\leq
cN^{-r}\|u_{1}\|_{B^{r}(\Lambda^{2})},\end{split}$

and

$\begin{split}I_{4}&\leq\|u_{2}-\pi_{N_{1}-1}u_{2}\|_{L^{2}(\Lambda^{2})}+\|\pi_{N_{1}-1}u_{2}-\pi_{N_{1}-1}\pi^{1,0}_{N_{2}}u_{2}\|_{L^{2}(\Lambda^{2})}\\\
&\leq
cN^{-r}\|u_{2}\|_{B^{r}(\Lambda^{2})}+\|u_{2}-\pi^{1,0}_{N_{2}}u_{2}\|_{L^{2}(\Lambda^{2})}\leq
cN^{-r}\|u_{2}\|_{B^{r}(\Lambda^{2})}.\end{split}$

A combination of the above facts leads to the desired result.

###  2.3 Solution algorithm based on fast diagonalizable auxiliary problem

In this subsection, we propose a highly efficient solution algorithm for the
divergence-free spectral approximation scheme (12), which draws inspirations
from the fast diagonalisation method (cf. [5, 17, 32]) for Poisson-type
equations and the matrix-free preconditioned Krylov subspace iterative method.

######  Proposition 4.

Denote symmetric matrices $M=(M_{mn})$, $I=(I_{mn})$ and $S=(S_{mn})$, where
${1\leq m,n\leq N-1}$ and the entries are given by

(40) $\displaystyle
M_{mn}=\big{(}\psi_{n+1},\psi_{m+1}\big{)}=\begin{cases}\dfrac{1}{2n+1}\big{(}\dfrac{1}{2n-1}+\dfrac{1}{2n+3}\big{)},\quad&m=n,\\\\[5.0pt]
-\dfrac{1}{\sqrt{2n+1}}\dfrac{1}{\sqrt{2n+5}}\dfrac{1}{2n+3},\quad&m=n+2,\\\\[5.0pt]
0,\quad&{\rm otherwise};\end{cases}$ (41) $\displaystyle
I_{mn}=\big{(}\phi_{n},\phi_{m}\big{)}=\delta_{mn};$ (42) $\displaystyle
S_{mn}=\big{(}\phi_{n}^{\prime},\phi_{m}^{\prime}\big{)}=\frac{1}{4}{\sqrt{(2m+1)(2n+1)}}\,n(n+1)\big{(}1+(-1)^{m+n}\big{)},\;\;m\geq
n.$

######  Proof 2.6.

The analytic expressions of $M_{mn}$ and $I_{mn}$ in equations (40)-(41) have
been given in [32] and can be easily derived with the help of the
orthogonality of Legendre polynomials in equation (4) and the recurrence
relation

(43) $(2n+1)L_{n}(\xi)=L_{n+1}^{\prime}(\xi)-L_{n-1}^{\prime}(\xi).$

Next, we look into the derivation of $S$. Since $S$ is symmetric, we consider
only $m\geq n$. By integration by parts, one obtains

(44)
$S_{mn}=\frac{\sqrt{(2m+1)(2n+1)}}{2}\Big{\\{}[L_{n}^{\prime}(\xi)L_{m}(\xi)]\big{|}_{-1}^{1}-\int_{-1}^{1}L_{n}^{\prime\prime}(\xi)L_{m}(\xi)\,{\rm
d}\xi\Big{\\}}.$

By the fact that when $m\geq n$,
$L_{n}^{\prime\prime}(\xi)\in\mathcal{P}_{m-2}$, thus
$L_{n}^{\prime\prime}(\xi)$ can be expanded by a linear combination of
$\\{L_{k}(\xi)\\}_{k=0}^{m-2}$. Thus, the last term in the above equation
vanishes due to orthogonality (4). Note that the boundary values of Legendre
polynomials satisfy

(45) $L_{n}(\pm 1)=(\pm 1)^{n},\quad L_{n}^{\prime}(\pm 1)=\frac{1}{2}(\pm
1)^{n-1}n(n+1).$

Inserting equation (45) into equation (44) directly leads to the analytic
expression (42).

One can expand the numerical solution $\boldsymbol{u}_{N}(\boldsymbol{x})$ of
the approximation scheme (12) by the divergence-free basis functions
$\boldsymbol{u}_{N}(\boldsymbol{x})=\sum_{m,n=1}^{N-1}u_{mn}\boldsymbol{\Phi}_{mn}(\boldsymbol{x}),$
and denote

$f_{mn}=\big{(}\boldsymbol{f},\boldsymbol{\Phi}_{mn}\big{)},\quad
F=(f_{mn})_{m,n=1,\cdots,N-1},\quad U=(u_{mn})_{m,n=1,\cdots,N-1}.$

Taking
$\boldsymbol{v}(\boldsymbol{x})=\boldsymbol{\Phi}_{\hat{m}\hat{n}}(\boldsymbol{x})$
in equation (8) for $\hat{m},\hat{n}=1,\cdots,N-1$, and noticing the
definitions of matrices $M$, $S$ and $I$ in equations (40)-(42), one arrives
at the following matrix equation

(46) $SUM+2IUI+MUS+\kappa\big{(}MUI+IUM\big{)}=F.$

One can remap the unknown matrix $U$ and source term matrix $F$ into vector
forms ${\rm vec}(U)$ and ${\rm vec}(F)$, following a natural index ordering as
follows:

$\\{({\rm vec}(U))_{p},({\rm
vec}(F))_{p}\\}_{p=1,\cdots,(N-1)^{2}}=\\{u_{mn},f_{mn}\\}_{m,n=1}^{N-1},\;\;p=m+(N-1)(n-1).$

The matrix equation (46) can be written equivalently into a linear system
$A\,{\rm vec}(U)$ $={\rm vec}(F)$, with the global matrix $A$ be assembled by

(47) $A:=M\otimes S+S\otimes M+2I\otimes I+\kappa(I\otimes M+M\otimes I).$

In order to design fast solution algorithm, it is necessary to look into the
sparsity of the global matrix $A$. Though $M$ is penta-diagonal, due to the
dense nature of the stiffness matrix $S$, the resultant global matrix $A$ is a
dense matrix. Thus, direct solvers based on fast factorization [13, 38] that
perform well for sparse matrix will no longer work. Generally speaking,
matrix-free Krylov subspace iteration methods [23, 30] would be the method of
choice for solving such large dense linear systems formed by summation of
Kronecker products of sub-matrices, due to significant reduction on the number
of floating point operations needed for operator evaluation and memory
comsumption, compared with solving the fully assembled system. Nevertheless,
these iterative methods would encounter severe convergence issue for ill-
conditioned system, arising from high-order spectral discretisation.

On the other hand, it is known that matrix diagonalisation method [5, 17, 32]
is an efficient solution algorithm for solving Poisson-type equations. The
essential idea of this technique resides in that a symmetric matrix and an
identity matrix can be simultaneously diagonalised by the same orthonormal
matrix. However, it is highly non-trivial and seems impoissble to extend this
idea to the current problem as it involves the simultaneous diagonalisation of
multiple matrices $M$, $S$ and $I$. In what follows, we show how to link these
matrices through two auxiliary eigenvalue problems and reveal the fact that
these three matrices are almost simultaneously diagonalizable.

To facilitate the development of the efficient algorithm for solving matrix
equation (46), one needs to resort to the following two auxiliary eigenvalue
problems:

(48) $y^{(4)}(x)=\lambda y(x),\;\;x\in\Lambda,\quad y(\pm
1)=0,\;\;y^{\prime\prime}(\pm 1)=0,$

and

(49) $-w^{(2)}(x)=\mu w(x),\;\;x\in\Lambda,\quad w(\pm 1)=0.$

Their weak formulations are as follows:

  * •

the biharmonic eigenvalue problem: find
$y(x)\in\mathring{H}^{2}(\Lambda):=\big{\\{}u\in H^{2}(\Lambda),$ $u(\pm
1)=0\big{\\}}$ and $\lambda\in\mathbb{R}$, such that

(50)
$(y^{\prime\prime},v^{\prime\prime})_{L^{2}(\Lambda)}=\lambda(y,v)_{L^{2}(\Lambda)},\;\;\forall\,v(x)\in\mathring{H}^{2}(\Lambda),$

  * •

the Laplace eigenvalue problem: find $w(x)\in H^{1}_{0}(\Lambda)$ and
$\mu\in\mathbb{R}$, such that

(51)
$(w^{\prime},v^{\prime})_{L^{2}(\Lambda)}=\mu(w,v)_{L^{2}(\Lambda)},\;\;\forall\,v(x)\in
H^{1}_{0}(\Lambda).$

Let us define the approximation space $\mathbb{Y}_{N}:={\rm
span}\\{\psi_{m+1}(x)\\}_{m=1}^{N-1}$ for equations (50) and (51), the related
approximation schemes for these two problems are to

  * •

find $y_{N}\in\mathbb{Y}_{N}$ and $\lambda_{h}\in\mathbb{R}$, such that

(52)
$(y_{N}^{\prime\prime},v^{\prime\prime})_{L^{2}(\Lambda)}=\lambda_{h}(y_{N},v)_{L^{2}(\Lambda)},\;\;\forall
v(x)\in\mathbb{Y}_{N},$

  * •

and find $w_{N}\in\mathbb{Y}_{N}$ and $\mu_{h}\in\mathbb{R}$, such that

(53)
$(w_{N}^{\prime},v^{\prime})_{L^{2}(\Lambda)}=\mu_{h}(w_{N},v)_{L^{2}(\Lambda)},\;\;\forall
v(x)\in\mathbb{Y}_{N}.$

Expanding $y_{N}(x)$ and $w_{N}(x)$ by

(54) $y_{N}(x)=\sum_{n=1}^{N-1}y_{n}\psi_{n+1}(x),\qquad
w_{N}(x)=\sum_{n=1}^{N-1}w_{n}\psi_{n+1}(x),$

and taking $v(x)=\psi_{m+1}(x)$ for $m=1,\cdots,N-1$ in (52) and (53),
respectively, one obtains

(55)
$\displaystyle\sum_{n=1}^{N-1}\big{(}\psi_{n+1}^{\prime\prime},\psi_{m+1}^{\prime\prime}\big{)}_{L^{2}(\Lambda)}y_{n,i}=\lambda_{h,i}\sum_{n=1}^{N-1}\big{(}\psi_{n+1},\psi_{m+1}\big{)}_{L^{2}(\Lambda)}y_{n,i},$
(56)
$\displaystyle\sum_{n=1}^{N-1}\big{(}\psi_{n+1}^{\prime},\psi_{m+1}^{\prime}\big{)}_{L^{2}(\Lambda)}w_{n,i}=\mu_{h,i}\sum_{n=1}^{N-1}\big{(}\psi_{n+1},\psi_{m+1}\big{)}_{L^{2}(\Lambda)}w_{n,i},$

where $\lambda_{h,i}$ and $\mu_{h,i}$ are the $i$-th eigenvalues of the
generalized matrix eigenvalue problems (55) and (56), respectively, and
$y_{n,i}$ and $w_{n,i}$ are respectively the $n$-th expansion coefficients of
the $i$-th eigenvector under the basis $\\{\psi_{n+1}(x)\\}_{n=1}^{N-1}$ of
equations (55) and (56). Aware of the derivative property (11) and the
definition of matrix $M$,$S$, and $I$ in equations (40)-(42), the above
equations (55)-(56) lead to the following generalized matrix eigenvalue
problems as follows:

(57) $IE=ME\Theta,\quad SE_{1}=ME_{1}\Theta_{1}.$

Here, $E,$ $E_{1}$ are orthonormal matrices and $\Theta$ and $\Theta_{1}$ are
diagonal matrices of size $(N-1)$-by-$(N-1)$, respectively, defined by
$(E)_{ni}=w_{n,i},$ $(E_{1})_{ni}=y_{n,i}$ and

(58) $\Theta={\rm diag}(\mu_{h,1},\cdots,\mu_{h,N-1}),\quad\Theta_{1}={\rm
diag}(\lambda_{h,1},\cdots,\lambda_{h,N-1}).$

It is crucial to observe that the eigenvalues and eigenvectors of equations
(48) and (49) have analytic expressions and are closely related to each other
as given below: for $i=1,\cdots,N-1$,

(59)
$y_{i}(x)=w_{i}(x)=\sin\big{(}{i}{\pi}(x+1)/{2}\big{)},\;\;\mu_{i}=\big{(}{i}{\pi}/{2}\big{)}^{2},\;\;\lambda_{i}=\mu_{i}^{2}.$

Thus, one can infer readily that $E_{1}\approx E$ and
$\Theta_{1}\approx\Theta^{2}$, and as a result, $MS\approx SM\approx I$. Note
that if $M$ and $S$ commute, i.e. $MS=SM$, there exists a set of common
orthonormal eigenvectors such that $M$ and $S$ are simultaneously
diagonalizable (see e.g. [24] p.74, Theorem 14). Next, we show how close $MS$
and $SM$ to the identity matrix in the lemma below. $MS$ and $SM$ are mutually
transposed matrices of size $(N-1)$-by-$(N-1)$, with the following analytic
expression:

(60) $\begin{split}&(SM)_{nm}=(MS)_{mn}\\\
&=\begin{cases}\delta_{mn},\quad&1\leq m\leq N-3,\\\\[5.0pt]
\delta_{(N-2)n}+\dfrac{1}{4}\sqrt{\dfrac{2n+1}{2N-3}}\dfrac{n(n+1)}{2N-1}\Big{(}1+(-1)^{n+N}\Big{)},\;\;&m=N-2,\\\\[5.0pt]
\delta_{(N-1)n}+\dfrac{1}{4}\sqrt{\dfrac{2n+1}{2N-1}}\dfrac{n(n+1)}{2N+1}\Big{(}1+(-1)^{n+N+1}\Big{)},\;\;&m=N-1.\end{cases}\end{split}$

######  Proof 2.7.

We postpone the detailed proof in Appendix A.

It can be observed from Lemma 2.3 that $MS$ is a lower triangular matrix, with
its first $(N-3)$-by-$(N-3)$ sub-matrix been identity. This inspires us to
approximate $S$ by $M^{-1}$. Note that from equation (57), one has

(61) $M=E\Theta^{-1}E^{\intercal}=EDE^{\intercal},$

where $D={\rm diag}(d_{1},\cdots,d_{N-1}):=\Theta^{-1}$ is a diagonal matrix
with its main diagonal entries the eigenvalues of $M$ and $E$ is the
corresponding orthonormal eigenvector matrix. Correspondingly, $S$ is
approximated by

(62) $S\approx M^{-1}=ED^{-1}E^{\intercal}.$

######  Remark 5.

The numerical eigenvalues/eigenvectors of $M$ and $S$ in equations (57), i.e.
$E,E_{1}$ and $\Theta,\Theta_{1}$, respectively, could also be approximated by
their analytic counterparts. To be more specific, $y_{n,i}$, $w_{n,i}$ and
$\lambda_{h,i}$ and $\mu_{h,i}$ can be approximated by
$\lambda_{h,i}\approx\lambda_{i}=\big{(}\frac{\pi}{2}i\big{)}^{4},$
$\mu_{h,i}\approx\mu_{i}=\big{(}\frac{\pi}{2}i\big{)}^{2},$ $y_{n,i}\approx
w_{n,i}\approx\hat{y}_{n,i}=\hat{w}_{n,i},$ and

(63)
$\sum_{n=1}^{N-1}\hat{y}_{n,i}\psi_{n+1}(x)=\sum_{n=1}^{N-1}\hat{w}_{n,i}\psi_{n+1}(x)\approx\sin\Big{(}\frac{\pi}{2}i(x+1)\Big{)},\;\;i=1,\cdots
N-1.$

By taking derivative on both sides of the above equation, using equation (11)
and the orthogonality of Legendre polynomials under the $L^{2}$ inner product,
one readily verifies

(64)
$\hat{y}_{n,j}=\hat{w}_{n,j}=\pi\sqrt{\frac{(2n+1)j}{2}}\cos\Big{(}\frac{\pi}{2}(j+n)\Big{)}J_{n+\frac{1}{2}}\Big{(}\frac{\pi}{2}j\Big{)}.$

Here, $J_{n}(r)$ is the Bessel function of the first kind and the following
identity has been used

(65) $\int_{-1}^{1}e^{{\rm i}rx}L_{n}(x)dx={\rm
i}^{n}\sqrt{\frac{2\pi}{r}}J_{n+1/2}(r),$

where ${\rm i}$ is the complex unit. Thus, define $\hat{D}={\rm
diag}(\hat{d}_{1},\cdots,\hat{d}_{N-1}),$
$\hat{d}_{j}=\big{(}\frac{\pi}{2}i\big{)}^{-2}$, and

(66)
$(\hat{E})_{ij}=\pi\sqrt{\frac{(2i+1)j}{2}}\cos\Big{(}\frac{\pi}{2}(j+i)\Big{)}J_{i+\frac{1}{2}}\Big{(}\frac{\pi}{2}j\Big{)}.$

$M$ and $S$ could also be approximated by

(67) $M\approx\hat{E}\hat{D}\hat{E}^{\intercal},\quad
S\approx\hat{E}\hat{D}^{-1}\hat{E}^{\intercal}.$

However, further numerical evidence illustrates the deficiency of using the
above approximation, compared with using the approximation in equation (62).

(a) Eigenvalues of $S$: exact vs approximation  (b) Eigenvectors of $S$: exact
vs approximation

Figure 1: Comparisons of the exact eigenvalues and eigenvectors of $S$ with
their approximations. (a) Eigenvalues of $S$ (in ascending order) as a
function of the eigenvalue index $i$ and compare them with the eigenvalues of
$M^{-1}$ and the analytic expression ${(\pi i/2)^{2}}$; (b) the differences in
$l^{2}$-norm between the associated eigenvectors of $S$ and $E$ (the
eigenvector matrix of $M^{-1}$) and the analytic expression $\hat{E}$ in
equation (66).

In order to quantify the accuracy of the approximations (62) and (67), let us
denote the eigenvalues of S by $\\{s_{i}\\}$ (in ascending order) and the
corresponding orthonormal eigenvector matrix by $E^{S}$. In Fig. 1 (a), we fix
the polynomial order $N=50$ and depict the spectrum of $S$ and $M^{-1}$, as a
function of the eigenvalue index $i$ in ascending order. The analytic formula
$(\pi i/2)^{2}$ obtained from equation (59) is also included for comparison.
One can observe that the differences of the top $70\%$ of these eigenvalue
pairs are negligible. While for the remaining $30\%$ of the eigenvalues, the
result obtained from $M^{-1}$ is obviously better than that from the analytic
formula $(\pi i/2)^{2}$. In Fig. 1 (b), the difference in $l^{2}$-norm of the
eigenvectors between $S$ and $M^{-1}$ and the analytic formula in (67), i.e.,
$\|E^{S}_{:,i}-E_{:,i}\|_{2}$ and $\|E^{S}_{:,i}-\hat{E}_{:,i}\|_{2}$ are
depicted as a function of the eigenvector index $i$. It can be seen that the
accuracy of both approximations to the top $30\%$ of the eigenvectors is up to
the machine roundoff error. While for the rest of the $70\%$ of the
eigenvectors, the approximation errors gradually increase from $10^{-14}$ to
$10^{-1}$ for both approximations. Moreover, it can be observed clearly that
in general, $\|E^{S}_{:,i}-E_{:,i}\|_{2}$ is smaller than
$\|E^{S}_{:,i}-\hat{E}_{:,i}\|_{2}$. These numerical evidences demonstrate
that $M^{-1}$ serves as a viable approximation of $S$, especially for the low
and medium frequency modes. This observation plays a significant role in the
development of the proposed method.

In what follows, we construct an auxiliary matrix equation by replacing $S$ by
$M^{-1}$ in equation (46) and show that this matrix equation is fully
diagonalizable and easy to compute.

######  Proposition 6.

Define an auxiliary matrix equation for the original system (46) as

(68) $M^{-1}UM+2IUI+MUM^{-1}+\kappa\big{(}MUI+IUM\big{)}=F,$

which can be written equivalently as the following linear system

(69) $\tilde{A}\,{\rm vec}(U)={\rm vec}(F),$

with $\tilde{A}:=M\otimes M^{-1}+M^{-1}\otimes M+2I\otimes I+\kappa(I\otimes
M+M\otimes I).$ Here, $M$ is the penta-diagonal matrix defined in equation
(40) and can be diagonalized by (61). The solution $U$ takes the form

(70)
$U=EVE^{\intercal},\quad\text{with}\,\,\,\,V_{ij}=\frac{(E^{\intercal}FE)_{ij}}{(d_{i}/d_{j}+d_{j}/d_{i}+2)+\kappa(d_{i}+d_{j})}.$

######  Proof 2.8.

Inserting equation (61) into equation (68) and introducing the auxiliary
matrix $V=E^{\intercal}UE$, one readily obtains

$D^{-1}VD+DVD^{-1}+2V+\kappa(DV+VD)=E^{\intercal}FE,$

which directly leads to the desired formula (70).

A comprehensive numerical investigation of the spectrum of ${\tilde{A}}^{-1}A$
leads to the following observations.

  1. (i)

The eigenvalues of $A$ spread widely in between $-10^{4}$ to $10^{7}$, while
that of ${\tilde{A}}^{-1}A$ is confined in between $1$ to $40$ and more than
$96\%$ of the eigenvalues of $\tilde{A}^{-1}A$ are clustered at $1$, computed
with a fixed $N=120$ and various $\kappa=-2000,-5000,-10000$, as can be seen
in Fig. 2.

  2. (ii)

The differences in $l^{\infty}$-norm between the first $13689$ eigenvalues
(i.e. $(120-3)^{2}$) and 1 is negligible, and the differences jump abruptly to
$O(1)$ for the rest of the eigenvalues, indicating that there is an extremely
large invariant subspace associated with the dominate eigenvalue $1$, as shown
in Fig. 2 (b).

(a) Spectrum of $A$ (b) Spectrum of $\tilde{A}^{-1}A$

Figure 2: Eigenvalue distributions of A and $\tilde{A}^{-1}A$ in 2D. (a)
Spectrum of $A$; (b) spectrum of $\tilde{A}^{-1}A$; (c) the differences in
$l^{\infty}$-norm between the eigenvalues of $\tilde{A}^{-1}A$ and $1$. These
results are computed for various $\kappa=-2000,-5000,-10000$ and a fixed
polynomial order $N=120$.

######  Theorem 7.

The dimension of the invariant subspace of the $(N-1)^{2}$-by-$(N-1)^{2}$
matrix $\tilde{A}^{-1}A$ with respect to the eigenvalue $1$ is $(N-3)^{2}$.

######  Proof 2.9.

It suffices to find $(N-3)^{2}$ linearly independent solutions ${\rm vec}(U)$
that satisfy

$A\,{\rm vec}(U)=\tilde{A}\,{\rm vec}(U),$

or more specificically, $(N-3)^{2}$ linearly independent $U$ that satisfy

$SUM+MUS=M^{-1}UM+MUM^{-1}.$

Let us introduce an auxiliary matrix $\bar{V}$ such that $U=M\bar{V}M$, the
above equation can be written equivalently as

(71) $(SM-I)\bar{V}M^{2}+M^{2}\bar{V}(MS-I)=0.$

With the help of the analytic expressions of $MS$ and $SM$ in Lemma 2.3, one
can obtain that all the elements in the first $(N-3)$ columns of $SM-I$ are
zeros. Thus, by choosing $\bar{V}$ be a $(N-1)$-by-$(N-1)$ matrix defined by

$\bar{V}=\begin{cases}\bar{v}_{ij},&i,j\leq N-3,\\\ 0,&{\rm
otherwise},\end{cases}$

where $\\{\bar{v}_{ij}\\}_{i,j=1}^{N-3}$ are $(N-3)^{2}$ degree of freedom,
one can readily verify that $\bar{V}$ satisfies equation (71). This ends the
proof.

These observations and Theorem 7 inspire us to combine the auxiliary problem
in Proposition 6 with a preconditioned Krylov subspace iterative method. In
what follows, we present a highly efficient solution algorithm, based on a
proper combination of the fast diagonalisation solver for the auxiliary
problem, the matrix-free operator evaluation approach without the assembly of
global matrix, and a preconditioned GMRES method. The algorithm is summarised
in Algorithm 1.

Algorithm 1 Fast divergence-free spectral algorithm for equation (7)

0: The polynomial order $N$, the matrices $M$ and $S$ in equation (40)-(42),
the matrices $E$ and $D$ in (61), the basis
$\big{\\{}\boldsymbol{\Phi}_{m,n}\big{\\}}_{m,n=1}^{N-1}$ in equation (9), the
positive integer $J_{\rm max}$ and the threshold $\varepsilon$

0: The divergence-free expansion coefficients $U$ of the numerical solution
$\boldsymbol{u}_{N}$

1: Compute $F=(f_{mn})_{m,n=1,\cdots,N-1},$ with
$f_{mn}=\big{(}\boldsymbol{f},\boldsymbol{\Phi}_{mn}\big{)}$

2: Evaluate $U_{0}=(u^{0}_{pq})_{p,q=1,\cdots,N-1}$ with
$U_{0}=E^{\intercal}FE$

3: Assign $u^{0}_{pq}\leftarrow
u^{0}_{pq}/\big{\\{}(d_{p}/d_{q}+d_{p}/d_{q}+2)+\kappa(d_{p}+d_{q})\big{\\}}$

4: Evaluate
$R=(r_{mn})_{m,n=1,\cdots,N-1}:=F-\big{(}SU_{0}M+2U_{0}+MU_{0}S+\kappa(MU_{0}+U_{0}M)\big{)}$

5: Assign $R\leftarrow E\tilde{R}E^{\intercal}$, $\boldsymbol{r}\leftarrow{\rm
vec}(R)$ and $j\leftarrow 0$, compute $\delta=\|\boldsymbol{r}\|_{2}$,
$\rho=\delta$, $\gamma=\|{\rm vec}(F)\|_{2}$ and
$\boldsymbol{v}^{(1)}=\boldsymbol{r}/\delta$

6: while $\rho/\gamma>\varepsilon$, and $j\leq J_{\rm max}$ do

7: Assign $j\leftarrow j+1$

8: Reshape $\boldsymbol{v}^{(j)}$ to a matrix $V$ using natural column
indexing and compute  $W=(w_{mn})_{m,n=1,\cdots,N-1}=SVM+2V+MVS+\kappa(MV+VM)$

9: Evaluate $\tilde{W}=(\tilde{w}_{pq})_{p,q=1,\cdots,N-1}$ with
$\tilde{W}=E^{\intercal}WE$

10: Assign
$\tilde{w}_{pq}\leftarrow\tilde{w}_{pq}/\big{\\{}(d_{p}/d_{q}+d_{q}/d_{p}+2)+\kappa(d_{p}+d_{q})\big{\\}}$

11: Assign $W=E\tilde{W}E^{\intercal}$ and reshape the matrix $W$ to a vector
$\boldsymbol{w}$ using natural column indexing

12: for $i=1,\cdots,j$ do

13: Evaluate $h_{i,j}=(\boldsymbol{w},\boldsymbol{v}^{(i)})$, assign
$\boldsymbol{w}\leftarrow\boldsymbol{w}-h_{i,j}\boldsymbol{v}^{(i)}$

14: end for

15: Evaluate $h_{j+1,j}=\|\boldsymbol{w}\|_{2}$ and
$\boldsymbol{v}^{(j+1)}=\boldsymbol{w}/h_{j+1,j}$

16: Define $H_{j}=\\{h_{p,q}\\}_{1\leq p\leq j+1;1\leq q\leq j}$ and
$\boldsymbol{e}_{1}=(1,0,\cdots,0)^{\intercal}\in\mathbb{R}^{j+1}$

17: Evaluate ${\boldsymbol{y}^{(j)}={\rm
argmin}_{\boldsymbol{y}}\|\delta\boldsymbol{e}_{1}-H_{j}\boldsymbol{y}\|_{2}}$,
assign $\rho=\|\delta\boldsymbol{e}_{1}-H_{j}\boldsymbol{y}^{(j)}\|_{2}$

18: end while

19: Define $V_{j}:=[\boldsymbol{v}^{(1)},\cdots,\boldsymbol{v}^{(j)}]$,
evaluate $\boldsymbol{x}=V_{j}\boldsymbol{y}^{(j)}$

20: Reshape $\boldsymbol{x}$ to matrix $X$ and compute $U=U_{0}+X$

We expect that $\tilde{A}^{-1}$ be not only an effective preconditioner for
$A$, but also the number of iterations will decrease, as the polynomial order
$N$ increases, since the portion of the invariant subspace of
$\tilde{A}^{-1}A$ w.r.t. the eigenvalue 1 increases as
$({N-3})^{2}/({N-1})^{2}$.

Moreover, the proposed method can be readily extended to solve the curl-curl
problem with variable coefficients

(72)
$\nabla\times\big{(}\alpha(\boldsymbol{x})\nabla\times(\beta(\boldsymbol{x}){\boldsymbol{u}}(\boldsymbol{x}))\big{)}+\kappa\,{\boldsymbol{u}}(\boldsymbol{x})={\boldsymbol{f}}(\boldsymbol{x}),\;\;\;\;\boldsymbol{x}\in\Lambda^{2},$

where $\alpha,\beta$ are functions depending on $\boldsymbol{x}$ and
$\nabla\cdot\boldsymbol{f}=0$. This can be accomplished with minor
modifications in Algorithm 1 as follows:

  1. 1.

Replace step 4 in Algorithm 1 by evaluating

(73) $\begin{split}R&=(r_{mn})_{m,n=1,\cdots,N-1}\\\
&:=\big{(}\boldsymbol{f},\boldsymbol{\Phi}_{mn}\big{)}-\big{(}\alpha(\boldsymbol{x})\nabla\times(\beta(\boldsymbol{x})\boldsymbol{u}_{0}),\nabla\boldsymbol{\Phi}_{mn}\big{)}+\kappa(\boldsymbol{u}_{0},\boldsymbol{\Phi}_{mn}),\end{split}$

and using Gaussian quadrature, where
${\boldsymbol{u}}_{0}(\boldsymbol{x})=\sum_{m,n=1}^{N-1}u^{0}_{mn}\boldsymbol{\Phi}_{mn}(\boldsymbol{x}).$

  2. 2.

Replace step 8 in Algorithm 1 by reshaping $\boldsymbol{v}^{(j)}$ to a matrix
$V$ using natural column indexing and computing

(74)
$W=(w_{mn})_{m,n=1,\cdots,N-1}=\big{(}\alpha(\boldsymbol{x})\nabla\times(\beta(\boldsymbol{x})\boldsymbol{v}),\nabla\boldsymbol{\Phi}_{mn}\big{)}+\kappa(\boldsymbol{v},\boldsymbol{\Phi}_{mn}),$

where
$\boldsymbol{v}(\boldsymbol{x})=\sum_{m,n=1}^{N-1}(V)_{mn}\boldsymbol{\Phi}_{mn}(\boldsymbol{x})$.

######  Remark 8.

Note that one can split the penta-diagonal matrix $M$ into two symmetric
tridiagonal sub-matrices and compute its associated eigenvalues and
eigenvectors in $O(N^{2})$ operations (see, e.g. a detailed description of the
solution procedure in [41]). The overall computational cost is dominated by
the matrix-matrix multiplication in step 2 and 9 in Algorithm 1, which is only
a small multiple of $N^{3}$ operations, as opposed to $O(N^{4})$ operations
for assembled global system.

##  3 A fast algorithm for the three-dimensional curl-curl equation

The proposed divergence-free spectral method and its fast solution algorithm
in two dimensions could be straightforwardly extended to the solution of
three-dimensional curl-curl problem with both constant and variable
coefficients. Let us propose the divergence-free approximation space first.

######  Proposition 9.

$\mathbb{X}_{N,3d}^{{\rm div}0}(\Lambda^{3})$ is a divergence-free
approximation space for $\mathbb{X}^{{\rm div}0}(\Lambda^{3})$ taking the form

$\mathbb{X}_{N,3d}^{{\rm div}0}={\rm
span}\Big{(}\Big{\\{}\boldsymbol{\Phi}_{m,n,l}^{1},\boldsymbol{\Phi}_{m,n,l}^{2}\Big{\\}}_{m,n,l=1}^{N-1},\,\Big{\\{}\boldsymbol{\Phi}^{x}_{n,l}\Big{\\}}_{n,l=1}^{N-1},\,\Big{\\{}\boldsymbol{\Phi}^{y}_{m,l}\Big{\\}}_{m,l=1}^{N-1},\,\Big{\\{}\boldsymbol{\Phi}^{z}_{m,n}\Big{\\}}_{m,n=1}^{N-1}\Big{)}.$

Here,
$\\{\boldsymbol{\Phi}_{m,n,l}^{1}(\boldsymbol{x}),\,\boldsymbol{\Phi}_{m,n,l}^{2}(\boldsymbol{x})\\}$
stands for the interior modal basis functions defined by

(75)
$\displaystyle\boldsymbol{\Phi}_{m,n,l}^{1}(\boldsymbol{x})=\Big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2})\phi_{l}(x_{3}),\;\;-\phi_{m}(x_{1})\psi_{n+1}(x_{2})\phi_{l}(x_{3}),\;\;0\Big{)}^{\intercal},$
$\displaystyle\boldsymbol{\Phi}_{m,n,l}^{2}(\boldsymbol{x})=\Big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2})\phi_{l}(x_{3}),\;\;0,\;\;-\phi_{m}(x_{1})\phi_{n}(x_{2})\psi_{l+1}(x_{3})\Big{)}^{\intercal},$

and
$\big{\\{}\boldsymbol{\Phi}^{x}_{n,l},\,\boldsymbol{\Phi}^{y}_{m,l},\,\boldsymbol{\Phi}^{z}_{m,n}\big{\\}}$
stands for the face modal basis functions defined by

(76)
$\displaystyle\boldsymbol{\Phi}^{x}_{n,l}(\boldsymbol{x})=\big{(}0,\;\;\psi_{n+1}(x_{2})\phi_{l}(x_{3}),\;\;-\phi_{n}(x_{2})\psi_{l+1}(x_{3})\big{)}^{\intercal},$
(77)
$\displaystyle\boldsymbol{\Phi}^{y}_{m,l}(\boldsymbol{x})=\big{(}\psi_{m+1}(x_{1})\phi_{l}(x_{3}),\;\;0,\;\;-\phi_{m}(x_{1})\psi_{l+1}(x_{3})\big{)}^{\intercal},$
(78)
$\displaystyle\boldsymbol{\Phi}^{z}_{m,n}(\boldsymbol{x})=\big{(}\psi_{m+1}(x_{1})\phi_{n}(x_{2}),\;\;-\phi_{m}(x_{1})\psi_{n+1}(x_{2}),\;\;0\big{)}^{\intercal},$

where $\\{\psi_{m+1}\\}_{m=1}^{N-1}$ and $\\{\phi_{m}\\}_{m=1}^{N-1}$ are
defined in equation (10).

######  Proof 3.1.

By the identity (11), one confirms that

(79)
$\nabla\cdot\boldsymbol{\chi}(\boldsymbol{x})=0,\quad\boldsymbol{\chi}=\boldsymbol{\Phi}_{m,n,l}^{1},\;\boldsymbol{\Phi}_{m,n,l}^{2},\;\boldsymbol{\Phi}^{x}_{n,l},\;\boldsymbol{\Phi}^{y}_{m,l},\;\boldsymbol{\Phi}^{z}_{m,n}.$

This ends the proof.

Correspondingly, the approximation scheme of the three-dimensional curl-curl
problem with constant coefficients in equation (8) takes the form: find
$\boldsymbol{u}_{N}\in\mathbb{X}_{N,3d}^{{\rm div}0}$ such that

(80)
$\mathcal{A}(\boldsymbol{u}_{N},\boldsymbol{v}):=(\nabla\times{\boldsymbol{u}}_{N},\nabla\times{\boldsymbol{v}})+\kappa({\boldsymbol{u}}_{N},\boldsymbol{v})=({\boldsymbol{f}},\boldsymbol{v}),\quad\forall\boldsymbol{v}\in\mathbb{X}_{N,3d}^{{\rm
div}0}.$

Let us expand the numerical solution $\boldsymbol{u}_{N}$ by

$\begin{split}\boldsymbol{u}_{N}(\boldsymbol{x})=&\sum_{m,n,l=1}^{N-1}\big{\\{}u_{mnl}^{1}\boldsymbol{\Phi}_{m,n,l}^{1}+u_{mnl}^{2}\boldsymbol{\Phi}_{m,n,l}^{2}\big{\\}}\\\
&+\sum_{n,l=1}^{N-1}u^{x}_{nl}\boldsymbol{\Phi}^{x}_{n,l}+\sum_{m,l=1}^{N-1}u^{y}_{ml}\boldsymbol{\Phi}^{y}_{m,l}+\sum_{m,n=1}^{N-1}u^{z}_{mn}\boldsymbol{\Phi}^{z}_{m,n},\end{split}$

where
$f^{i}_{\hat{m}\hat{n}\hat{l}}=\big{(}\boldsymbol{f},\boldsymbol{\Phi}^{i}_{\hat{m},\hat{n},\hat{l}}\big{)},\;\;i=1,2$,
and

$\displaystyle
f^{x}_{\hat{n}\hat{l}}=\big{(}\boldsymbol{f},\boldsymbol{\Phi}^{x}_{\hat{n},\hat{l}}\big{)},\;\;f^{y}_{\hat{m}\hat{l}}=\big{(}\boldsymbol{f},\boldsymbol{\Phi}^{y}_{\hat{m},\hat{l}}\big{)},\;\;f^{z}_{\hat{m}\hat{n}}=\big{(}\boldsymbol{f},\boldsymbol{\Phi}^{z}_{\hat{m},\hat{n}}\big{)}.$

Remapping the multidimensional variables
$\\{u^{1}_{mnl},u^{2}_{mnl},u^{x}_{nl},u^{y}_{ml},u^{z}_{mn}\\}$ and the
source terms $\\{f^{1}_{mnl},$ $f^{2}_{mnl},$ $f^{x}_{nl},$
$f^{y}_{ml},f^{z}_{mn}\\}$ into a vector form following a natural index
ordering as follows:

$\displaystyle\big{\\{}\big{(}{\rm vec}(u^{i})\big{)}_{p_{1}},({\rm
vec}(f^{i}))_{p_{1}}\big{\\}}_{p_{1}=1,\cdots,(N-1)^{3}}=\\{u^{i}_{mnl},f^{i}_{mnl}\\}_{m,n,l=1}^{N-1},\;\;i=1,2,$
$\displaystyle\big{\\{}\big{(}{\rm vec}(u^{x})\big{)}_{p},\big{(}{\rm
vec}(f^{x})\big{)}_{p}\big{\\}}_{p=1,\cdots,(N-1)^{2}}=\\{u^{x}_{nl},f^{x}_{nl}\\}_{n,l=1}^{N-1},\;\;p=n+(N-1)(l-1),$
$\displaystyle\big{\\{}\big{(}{\rm vec}(u^{y})\big{)}_{p},\big{(}{\rm
vec}(f^{y})\big{)}_{p}\big{\\}}_{p=1,\cdots,(N-1)^{2}}=\\{u^{y}_{ml},f^{y}_{ml}\\}_{m,l=1}^{N-1},\;\;p=m+(N-1)(l-1),$
$\displaystyle\big{\\{}\big{(}{\rm vec}(u^{z})\big{)}_{p},\big{(}{\rm
vec}(f^{z})\big{)}_{p}\big{\\}}_{p=1,\cdots,(N-1)^{2}}=\\{u^{z}_{mn},f^{z}_{mn}\\}_{m,n=1}^{N-1},\;\;p=m+(N-1)(n-1),$

where $p_{1}=m+(N-1)(n-1)+(N-1)^{2}(l-1)$ for the first identity in the above.
Let us take
$\boldsymbol{v}(\boldsymbol{x})=\Big{\\{}\boldsymbol{\Phi}_{\hat{m},\hat{n},\hat{l}}^{1},\;\boldsymbol{\Phi}_{\hat{m},\hat{n},\hat{l}}^{2},$
$\boldsymbol{\Phi}^{x}_{\hat{n},\hat{l}},\;\boldsymbol{\Phi}^{y}_{\hat{m},\hat{l}},\;\boldsymbol{\Phi}^{z}_{\hat{m},\hat{n}}\Big{\\}}$
in equation (80) for $\hat{m},\hat{n},\hat{l}=1,\cdots,N-1$, one obtains the
following system of linear equations:

(81) $\displaystyle A^{11}\,{\rm vec}(u^{1})+A^{12}\,{\rm vec}(u^{2})={\rm
vec}(f^{1}),\;\;\;\;A^{21}\,{\rm vec}(u^{1})+A^{22}\,{\rm vec}(u^{2})={\rm
vec}(f^{2}),$ $\displaystyle A^{xx}\,{\rm vec}(u^{x})={\rm
vec}(f^{x}),\;\;\;\;A^{yy}\,{\rm vec}(u^{y})={\rm
vec}(f^{y}),\;\;\;\;A^{zz}\,{\rm vec}(u^{z})={\rm vec}(f^{z}),$

where

(82) $\begin{split}A^{11}:=&S\otimes M\otimes I+S\otimes I\otimes M+I\otimes
M\otimes S+I\otimes S\otimes M\\\ &+2I\otimes I\otimes I+\kappa\big{(}I\otimes
I\otimes M+I\otimes M\otimes I\big{)},\\\ A^{22}:=&M\otimes S\otimes
I+S\otimes I\otimes M+M\otimes I\otimes S+I\otimes S\otimes M\\\ &+2I\otimes
I\otimes I+\kappa\big{(}I\otimes I\otimes M+M\otimes I\otimes I\big{)},\\\
A^{12}=&A^{21}:=S\otimes I\otimes M+I\otimes S\otimes M+I\otimes I\otimes
I+\kappa I\otimes I\otimes M,\\\ A^{xx}=&A^{yy}=A^{zz}:=2\Big{(}M\otimes
S+S\otimes M+2I\otimes I+\kappa\big{(}I\otimes M+M\otimes
I\big{)}\Big{)}.\end{split}$

Similar with the two-dimensional case, we construct an auxiliary equation by
replacing $S$ by $M^{-1}$ in equation (82) and show that the system of matrix
equation is also fully diagonalizable and the solution can be obtained via
analytic expressions.

######  Proposition 10.

Define an auxiliary system of matrix equations of the original system (81)

(83) $\displaystyle\tilde{A}^{11}\,{\rm vec}(u^{1})+\tilde{A}^{12}\,{\rm
vec}(u^{2})={\rm vec}(f^{1}),\;\;\;\tilde{A}^{21}\,{\rm
vec}(u^{1})+\tilde{A}^{22}\,{\rm vec}(u^{2})={\rm vec}(f^{2}),$
$\displaystyle\tilde{A}^{xx}\,{\rm vec}(u^{x})={\rm
vec}(f^{x}),\;\;\;\tilde{A}^{yy}\,{\rm vec}(u^{y})={\rm
vec}(f^{y}),\;\;\;\tilde{A}^{zz}\,{\rm vec}(u^{z})={\rm vec}(f^{z}),$

where $\tilde{A}^{xx}=\tilde{A}^{yy}=\tilde{A}^{zz}$ and

(84) $\begin{split}\tilde{A}^{xx}=&:=2\Big{(}M\otimes M^{-1}+M^{-1}\otimes
M+2I\otimes I+\kappa\big{(}I\otimes M+M\otimes I\big{)}\Big{)},\\\
\tilde{A}^{11}:=&M^{-1}\otimes M\otimes I+M^{-1}\otimes I\otimes M+I\otimes
M\otimes M^{-1}+I\otimes M^{-1}\otimes M\\\ &+2I\otimes I\otimes
I+\kappa\big{(}I\otimes I\otimes M+I\otimes M\otimes I\big{)},\\\
\tilde{A}^{22}:=&M\otimes M^{-1}\otimes I+M^{-1}\otimes I\otimes M+M\otimes
I\otimes M^{-1}+I\otimes M^{-1}\otimes M\\\ &+2I\otimes I\otimes
I+\kappa\big{(}I\otimes I\otimes M+M\otimes I\otimes I\big{)},\\\
\tilde{A}^{12}=&\tilde{A}^{21}:=M^{-1}\otimes I\otimes M+I\otimes
M^{-1}\otimes M+I\otimes I\otimes I+\kappa I\otimes I\otimes M.\end{split}$

Here, $M$ is the penta-diagonal matrix defined in equation (40) and can be
diagonalized by (61). Let us denote the auxiliary variables
$\\{v^{1}_{ijk},v^{2}_{ijk},v_{jk}^{x},v_{ik}^{y},v_{ij}^{z}\\}_{i,j,k=1}^{N-1}$
and define for $m,n,l=1,\cdots,N-1$,

(85)
$\displaystyle\tilde{f}^{1}_{mnl}=\sum_{i,j,k=1}^{N-1}E^{\intercal}_{mi}E^{\intercal}_{nj}E^{\intercal}_{lk}f^{1}_{ijk},\quad\tilde{f}^{2}_{mnl}=\sum_{i,j,k=1}^{N-1}E^{\intercal}_{mi}E^{\intercal}_{nj}E^{\intercal}_{lk}f^{2}_{ijk},$
$\displaystyle\tilde{f}^{x}_{nl}=\sum_{i,j=1}^{N-1}E_{ni}^{\intercal}E_{lj}^{\intercal}f^{x}_{ij},\quad\tilde{f}^{y}_{ml}=\sum_{i,j=1}^{N-1}E_{mi}^{\intercal}E_{lj}^{\intercal}f_{ij}^{y},\quad\tilde{f}^{z}_{mn}=\sum_{i,j=1}^{N-1}E_{mi}^{\intercal}E_{nj}^{\intercal}f^{z}_{ij}.$

The solution $\\{u_{mnl}^{1},u_{mnl}^{2},u_{nl}^{x},u_{ml}^{y},u_{mn}^{z}\\}$
for $m,n,l=1,\cdots,N-1$ takes the form

(86) $\displaystyle
u^{1}_{mnl}=\sum_{i,j,k=1}^{N-1}E_{mi}E_{nj}E_{lk}v^{1}_{ijk},\quad
u^{2}_{mnl}=\sum_{i,j,k=1}^{N-1}E_{mi}E_{nj}E_{lk}v^{2}_{ijk},$ $\displaystyle
u_{nl}^{x}=\sum_{j,k=1}^{N-1}E_{nj}E_{lk}v_{jk}^{x},\quad
u_{ml}^{y}=\sum_{i,k=1}^{N-1}E_{mi}E_{lk}v_{ik}^{y},\quad
u_{mn}^{z}=\sum_{i,j=1}^{N-1}E_{mi}E_{nj}v_{ij}^{z},$

where

(87) $\displaystyle
v_{mnl}^{1}=\dfrac{\gamma_{mnl}^{22}\tilde{f}_{mnl}^{1}-\gamma_{mnl}^{12}\tilde{f}_{mnl}^{2}}{\gamma_{mnl}^{22}\gamma_{mnl}^{11}-\gamma_{mnl}^{12}\gamma_{mnl}^{21}},\qquad
v_{mnl}^{2}=\dfrac{-\gamma_{mnl}^{21}\tilde{f}_{mnl}^{1}+\gamma_{mnl}^{11}\tilde{f}_{mnl}^{2}}{\gamma_{mnl}^{22}\gamma_{mnl}^{11}-\gamma_{mnl}^{12}\gamma_{mnl}^{21}},$
$\displaystyle
v_{nl}^{x}=\dfrac{\tilde{f}^{x}_{nl}}{d_{n}/d_{l}+d_{l}/d_{n}+2+\kappa(d_{n}+d_{l})},$
$\displaystyle
v_{ml}^{y}=\dfrac{\tilde{f}^{y}_{ml}}{d_{m}/d_{l}+d_{l}/d_{m}+2+\kappa(d_{m}+d_{l})},$
$\displaystyle
v_{mn}^{z}=\dfrac{\tilde{f}^{z}_{mn}}{d_{m}/d_{n}+d_{n}/d_{m}+2+\kappa(d_{m}+d_{n})},$

with

(88)
$\displaystyle\gamma_{mnl}^{11}=d_{n}/d_{l}+d_{m}/d_{l}+d_{n}/d_{m}+d_{m}/d_{n}+2+\kappa(d_{m}+d_{n}),$
$\displaystyle\gamma_{mnl}^{12}=d_{m}/d_{l}+d_{m}/d_{n}+1+\kappa
d_{m},\;\;\;\gamma_{mnl}^{21}=d_{m}/d_{l}+d_{m}/d_{n}+1+\kappa d_{m},$
$\displaystyle\gamma_{mnl}^{22}=d_{l}/d_{n}+d_{m}/d_{l}+d_{l}/d_{m}+d_{m}/d_{n}+2+\kappa(d_{m}+d_{l}).$

######  Proof 3.2.

The proof is similar with that of Proposition 6. For the sake of conciseness,
we skip the detail.

It can be observed that the face modes $\big{\\{}{\rm vec}(u^{x}),{\rm
vec}(u^{y}),{\rm vec}(u^{z})\big{\\}}$ are fully decoupled with the interior
modes $\big{\\{}{\rm vec}(u^{1}),{\rm vec}(u^{2})\big{\\}}$ and the governing
equations for the face modes are exactly the same with the two-dimensional
case in equation (46). Thus, we could focus on the solution of the interior
modes. With a slight abuse of notation, we denote $A$ and $\tilde{A}$,
respectively, by

(89) $A=\begin{bmatrix}A^{11}&A^{12}\\\
A^{21}&A^{22}\end{bmatrix},\quad\tilde{A}=\begin{bmatrix}\tilde{A}^{11}&\tilde{A}^{12}\\\
\tilde{A}^{21}&\tilde{A}^{22}\end{bmatrix}.$

Once again, we investigate theoretically the spectrum of $\tilde{A}^{-1}A$ and
the following result holds.

######  Theorem 11.

The dimension of the invariant subspace of the $2(N-1)^{3}$-by-$2(N-1)^{3}$
matrix $\tilde{A}^{-1}A$ with respect to the eigenvalue $1$ is $2(N-3)^{3}$.

######  Proof 3.3.

The proof is similar with Theorem 7. It suffices to find $2(N-3)^{2}$ linearly
independent solutions $\\{{\rm vec}(u_{1}),{\rm vec}(u_{2})\\}$ such that

(90) $\displaystyle A^{11}\,{\rm vec}(u^{1})+A^{12}\,{\rm
vec}(u^{2})=\tilde{A}^{11}\,{\rm vec}(u^{1})+\tilde{A}^{12}\,{\rm
vec}(u^{2}),$ $\displaystyle A^{21}\,{\rm vec}(u^{1})+A^{22}\,{\rm
vec}(u^{2})=\tilde{A}^{21}\,{\rm vec}(u^{1})+\tilde{A}^{22}\,{\rm
vec}(u^{2}).$

Let us introduce the auxiliary variables
$\\{\bar{v}^{1}_{mnl},\bar{v}^{2}_{mnl}\\}$ such that

$u^{1}_{mnl}=\sum_{i,j,k=1}^{N-1}M_{mi}M_{nj}M_{lk}\bar{v}^{1}_{ijk},\quad
u^{2}_{mnl}=\sum_{i,j,k=1}^{N-1}M_{mi}M_{nj}M_{lk}\bar{v}^{2}_{ijk}.$

Inserting the above relation into equation (90), using the fact that the
elements in the first $(N-3)$ columns of $SM-I$ and the first $(N-3)$ rows of
$MS-I$ are zeros and choosing $\\{\bar{v}^{1}_{mnl},\bar{v}^{2}_{mnl}\\}$ such
that

$\bar{v}^{i}_{mnl}=\begin{cases}x^{i}_{mnl},\;\;&1\leq m,n,l\leq N-3,\\\
0,&{\rm otherwise},\end{cases},\quad i=1,2.$

where $\\{x^{i}_{mnl}\\}^{i=1,2}_{1\leq m,n,l\leq N-3}$ are $2(N-3)^{3}$ free
variables. One can readily show that $\\{u^{1}_{mnl},u^{2}_{mnl}\\}$ satisfies
the linear system (90), which ends the proof.

Thus, one can expect that the auxiliary problem in Proposition 10, when
combined with a preconditioned Krylov subspace iterative method, gives rise to
an efficient solution algorithm for the three dimensional curl-curl problem
with constant or variable coefficients. The solution procedure resembles
Algorithm 1. For the sake of conciseness, we skip the details of this
algorithm.

######  Remark 12.

The overall computational cost for the three-dimensional problem is dominated
by the tensor-matrix multiplication, which is only a small multiple of $N^{4}$
operations, as opposed to $O(N^{6})$ operations for assembled global system.

##  4 Representative numerical examples

In this section, we provide ample numerical examples to verify the accuracy
and efficiency of our proposed fast divergence-free spectral algorithm. We
conduct convergence tests employing both two- and three-dimensional contrived
solutions, with both constant and variable coefficients. We present the
convergence rates in both $L^{2}$\- and $H({\rm curl})$-norms. Additionally,
we also subject the method to challenging numerical tests on indefinite
systems with highly oscillatory solutions induced by inhomogeneous point
sources. In what follows, we adopt a stopping threshold $\varepsilon=10^{-12}$
for the proposed iterative method, without stated otherwise.

######  Example 4.1.

(Convergence test for curl-curl problem with constant coefficients in 2D). We
first consider (7) with the following manufactured solution

(91)
$\boldsymbol{u}(x,y)=\nabla\times\Big{(}\big{(}1-x^{2}\big{)}^{3}\big{(}1-y^{2}\big{)}^{3}\big{(}1+\sin(\pi
x)\cos(\pi y)\big{)}\exp\big{(}\sin(\pi x)\cos(\pi y)\big{)}\Big{)}.$

Accordingly, the source term $\boldsymbol{f}(x,y)$ in equation (7) is chosen
such that the analytic expression (91) satisfies equation (7).

(a) Errors vs $N$ (b) Iteration vs no iteration

Figure 3: Example 4.1: convergence tests of the curl-curl problem with
constant coefficients in 2D. (a) $L^{2}$\- and $H({\rm curl})$-errors of
$\boldsymbol{u}$ versus polynomial order $N$ for $\kappa=\pm 100$; (b) compare
the errors between the solution obtained from the proposed method and those
obtained directly from the auxiliary problem in Proposition 6.

In Fig. 3 (a), we plot the discrete $L^{2}$\- and $H({\rm curl})$-errors on
the semi-log scale against various $N$ using the proposed fast divergence-free
spectral algorithm. We observe the expected exponential convergence of the
proposed method for both positive and negative parameters $\kappa$ as $N$
increases. Fig. 3 (b) compares the numerical errors of the current method with
those obtained from directly solving the auxiliary problem stated in
Proposition 6, This comparison indicates that the auxiliary problem provides
numerical solutions with almost the same accuracy as the iterative algorithm
for rather smooth solutions. Thus, the auxiliary problem can serve as a good
approximation for the original problem, which is in agreement with previous
discussions in Section 2.3.

Table 1: Example 4.1: iteration numbers and numerical errors with various $\kappa$. | $\kappa=1$ |  | $\kappa=100$ |  | $\kappa=-1$ |  | $\kappa=-100$
---|---|---|---|---|---|---|---
$N$ | iter | Error |  | iter | Error |  | iter | Error |  | iter | Error
12 | 13 | 5.27e-2 |  | 12 | 5.31e-2 |  | 13 | 5.27e-2 |  | 16 | 5.91e-2
20 | 9 | 5.61e-5 |  | 9 | 5.62e-5 |  | 9 | 5.61e-5 |  | 10 | 5.63e-5
28 | 4 | 9.10e-9 |  | 4 | 9.10e-9 |  | 4 | 9.10e-9 |  | 4 | 9.11e-9
36 | 1 | 9.21e-13 |  | 1 | 9.21e-13 |  | 1 | 9.20e-13 |  | 0 | 9.26e-13
40 | 0 | 6.72e-14 |  | 0 | 6.62e-14 |  | 0 | 6.57e-14 |  | 0 | 1.24e-13

Next, we provide numerical evidence for one prominent feature of the proposed
method, namely, the iteration count will decrease as the polynomial order
increases. We tabulate in Table 1 the $H({\rm curl};\Omega)$-errors and the
corresponding iteration numbers obtained by Algorithm 1 with various $\kappa$
and $N$. Notably, as $N$ increases, the iteration count decreases, and
remarkably, it drops to $0$ when $N=40.$ We then set $\kappa=1$ and further
increase the polynomial order, comparing the iteration numbers and CPU time of
the proposed algorithm with those obtained without any preconditioner or with
Jacobi preconditioner, see Table 2. It can be observed that the proposed
method significantly outperforms the other two methods in terms of both
iteration numbers and CPU time. Unlike the algorithms without a preconditioner
or with a Jacobi preconditioner, whose iteration numbers increase with the
polynomial order, the iteration numbers of the proposed method remain at 0.
This occurs because the proportion of linearly independent eigenvectors of
$\tilde{A}^{-1}A$ with respect to 0 is $(N-3)^{2}/(N-1)^{2}$, which approaches
$100\%$ as $N$ gets larger.

Table 2: Example 4.1: iteration numbers and CPU time for large $N$. | No preconditioner |  | Jacobi preconditioner |  | Current method
---|---|---|---|---|---
$N$ | iter | CPU time (s) |  | iter | CPU time (s) |  | iter | CPU time (s)
50 | 237 | 1.34 |  | 158 | 0.94 |  | 0 | 1.82e-3
100 | 523 | 30.62 |  | 248 | 7.42 |  | 0 | 6.76e-3
500 | Not conv. | N.A. |  | 622 | 1040.21 |  | 0 | 0.29
1000 | Not conv. | N.A. |  | Not conv. | N.A. |  | 0 | 1.58

######  Example 4.2.

(Curl-curl problem with point source in 2D). We consider equation (7) with
Gaussian point source

$\boldsymbol{f}=\nabla\times\Big{(}{\rm
exp}\Big{(}{-\frac{(x+0.5)^{2}+y^{2}}{\sigma^{2}}}\Big{)}+{\rm
exp}\Big{(}{-\frac{(x-0.5)^{2}+y^{2}}{\sigma^{2}}}\Big{)}\Big{)},$

where $\sigma=0.01.$

(a) $\kappa=-100$ (b) $\kappa=-400$ (c) $\kappa=-2500$ (d) $\kappa=-10000$

Figure 4: Example 4.2: The profiles of the $x$-component of the numerical
solution versus different $\kappa$. These results are obtained with a fixed
polynomial order $N=128.$

In this case, the analytic solution is unavailable and numerical solution is
obtained via the proposed fast divergence-free spectral algorithm with
polynomial order $N=128$. We challenge our algorithm via highly indefinite
system with highly oscillatory solutions by choosing the parameter $\kappa$ as
large negative values $-100,-400,-2500,$ $-10000$, see Figure 4. It can be
seen that the proposed method provides good approximations, even for solutions
with highly oscillation. We also tabulate in Table 3 the number of iterations
of the proposed method as a function of the polynomial order, where we take
$\kappa=-10000.$ Astonishingly, It can be observed that the number of
iterations remains small and decreases from $35$ to $15$ as the polynomial
order increases from $120$ to $520$, demonstrating the superiority of the
proposed method for curl-curl problem with highly oscillatory solutions.

Table 3: Example 4.2: iteration numbers for various $N$ with $\kappa=-10000$. $N$ | 120 | 160 | 200 | 240 | 280 | 320 | 360 | 400 | 440 | 480 | 520
---|---|---|---|---|---|---|---|---|---|---|---
iter | 35 | 32 | 35 | 34 | 37 | 34 | 32 | 26 | 22 | 18 | 15

######  Example 4.3.

(Time-dependent Maxwell’s equations with variable coefficients in 2D).
Consider the following time-dependent problem with variable coefficients in
$\Lambda^{2}$ for $t\in(0,T]:$

(92) $\begin{cases}\dfrac{\partial\boldsymbol{B}}{\partial
t}+\nabla\times\boldsymbol{E}=0,\\\\[5.0pt] \dfrac{\partial}{\partial
t}(\epsilon_{0}\epsilon_{r}\boldsymbol{E})-\frac{1}{\mu_{0}}\nabla\times\boldsymbol{B}=-\boldsymbol{J},\end{cases}$

where $\boldsymbol{E}$ is the electric field and $\boldsymbol{B}$ is the
magnetic field taking the form

$\boldsymbol{E}=(0,0,E_{3}(x,y,t))^{\intercal},\quad\boldsymbol{B}=(B_{1}(x,y,t),B_{2}(x,y,t),0)^{\intercal}.$

Here, $\epsilon_{0}$ and $\mu_{0}$ are the permittivity and permeability in
vaccum, respectively, and $\epsilon_{r}$ is the relative permittivity. The
electric current density $\boldsymbol{J}$ is chosen as a point source given by

$\boldsymbol{J}=(0,0,J_{3}(x,y))^{\intercal},\quad
J_{3}(x,y)=\exp\Big{(}-\frac{(x+0.8)^{2}+y^{2}}{\sigma^{2}}\Big{)}.$

Here, we choose $\epsilon_{0}=1,$ $\mu_{0}=1$, $\sigma=0.04,$ and

$\epsilon_{r}(x,y)=1.5+0.5\tanh\Big{(}\frac{0.16+0.07\sin\big{(}6(\theta+\pi/4)\big{)}-x^{2}-y^{2}}{\lambda}\Big{)},$

where $\theta(x,y)$ is the angle of point $\boldsymbol{x}$ and $\lambda$ is
set to be $0.05$. The boundary condition is prescribed as

$\boldsymbol{n}\cdot\boldsymbol{B}=0,\quad\boldsymbol{n}\times(\nabla\times\boldsymbol{B})=0\;\;{\rm
on}\;\partial\Lambda^{2}.$

Let us partition the time interval $(0,T]$ using uniform time step size $\tau$
by

$t_{n}=n\tau,\;\;n=0,1,\cdots,N_{t},\quad\tau=T/N_{t}.\ $

and denote $\chi^{n}$ the numerical approximation of an arbitrary varaibe
$\chi$ at time step $n$, corresponding to the time $t_{n}.$ In the temporal
direction, we adopt the Crank-Nicolson (CN) temporal discretization

(93)
$\begin{cases}\dfrac{\boldsymbol{B}^{n+1}-\boldsymbol{B}^{n}}{\tau}+\nabla\times(\dfrac{\boldsymbol{E}^{n+1}+\boldsymbol{E}^{n}}{2})=0,\\\\[5.0pt]
\epsilon_{0}\epsilon_{r}\dfrac{\boldsymbol{E}^{n+1}-\boldsymbol{E}^{n}}{\tau}-\dfrac{1}{\mu_{0}}\nabla\times(\dfrac{\boldsymbol{B}^{n+1}+\boldsymbol{B}^{n}}{2})=-\boldsymbol{J},\end{cases}$

which can be reformulated into a curl-curl problem of $\boldsymbol{B}^{n+1}$.
In the spatial direction, we use the proposed fast divergence-free spectral
algorithm to discretize $\boldsymbol{B}^{n+1}$.

We depict the temporal sequence of snapshots (see Fig. 5 (b)-(d)) of $B_{2}$
at time $t=0.31,0.71,0.81$. It can be seen that the proposed method can
accurately simulate the the interaction of the electromagnetic waves with the
heterogeneous medium (see Fig. 5 (a)) and there is no numerical oscillations
occur for long time simulations. Moreover, the average number of iterations is
only around 36, indicating that the auxiliary problem serves as an effective
preconditioner for Maxwell’s equations with variable coefficients.

(a) Field phase (b) t=0.31 (c) t=0.71 (d) t=0.81

Figure 5: Example 4.3: electromagnetic wave propagation in inhomogeneous
medium. (a) Profile of the relative permittivity $\epsilon_{r}$; (b)-(d)
snapshots of $B_{2}$ at $t=0.31,0.71,0.81$. The simulations are obtained with
$N=250$ and $\tau=0.01.$

######  Example 4.4.

(Convergence test for curl-curl problem with constant coefficients in 3D).
Next, we consider the curl-curl problem (7) in three dimensions with the
manufactured solution

$\boldsymbol{u}(x,y,z)=\nabla\times\Big{(}\big{(}1-x^{2}\big{)}^{3}\big{(}1-y^{2}\big{)}^{3}\big{(}1-z^{2}\big{)}^{3}\big{(}1+w(x,y,z)\big{)}\exp(w(x,y,z))\Big{)}(1,1,1)^{\intercal},$

where $w(x,y,z)=\sin(\pi x)\sin(\pi y)\sin(\pi z)$, and the source term
$\boldsymbol{f}(x,y,z)$ in equation (7) is chosen such that the analytic
expression (4.4) satisfies equation (7).

(a) Errors vs $N$ (b) Iteration numbers vs $N$

Figure 6: Example 4.4: convergence tests of the curl-curl problem with
constant coefficients in 3D. (a) $L^{2}$\- and $H({\rm curl})$-errors of
$\boldsymbol{u}$; (b) the number of iterations of the proposed method. The
simulations are obtained with $\kappa=\pm 100$.

In Fig. 6 (a), we plot the discrete $L^{2}$\- and $H({\rm curl})$-errors on
the semi-log scale of the proposed fast divergence-free spectral algorithm for
(7) with (4.4) versus $N$. Observe that the errors decay exponentially in
terms of $N$ increases for both positive and negative parameters $\kappa$. We
plot the iteration numbers of the proposed method versus $N$, see Fig. 6 (b).
Again, for a given smooth function, we can see that the iteration number
reduces to $0$ as $N$ increases, which is in agreement with the analysis in
Theorem 11.

######  Example 4.5.

(Time-dependent Maxwell’s equations with point source in 3D). Consider the
time-dependent Maxwell’s equations as in Example 4.3, where

$\boldsymbol{E}=(E_{1}(\boldsymbol{x},t),E_{2}(\boldsymbol{x},t),E_{3}(\boldsymbol{x},t))^{\intercal},\quad\boldsymbol{B}=(B_{1}(\boldsymbol{x},t),B_{2}(\boldsymbol{x},t),B_{3}(\boldsymbol{x},t))^{\intercal},$

and the electric current density $\boldsymbol{J}=(J_{1},0,0)^{\intercal}$ is
given by

(94)
$J_{1}(x,y,z)=\exp\Big{(}-\frac{x^{2}+y^{2}+(z+0.5)^{2}}{\sigma^{2}}\Big{)}+\exp\Big{(}-\frac{x^{2}+y^{2}+(z-0.5)^{2}}{\sigma^{2}}\Big{)}.$

Here, we choose $\sigma=0.05,$ $\epsilon_{0}=\mu_{0}=1,$ and $\epsilon_{r}=1.$

Once again, we resort to the Crank-Nicolson scheme for temporal discretization
and the proposed fast divergence-free spectral algorithm for spatial
discretization of the resultant curl-curl problem of $\boldsymbol{B}$ in three
dimensions. In these simulations, we prescribe the polynomial order $N=100$
and the time step size $\tau=0.02$. It can be observed that the proposed
method can accurately simulate the the interaction of the electromagnetic
waves induced by inhomogeneous point sources, see Fig. 7. Moreover, the
average number of iterations is merely $8.21$, demonstrating that the proposed
method is highly efficient for long time 3D simulations.

(a) $t=0.1$ (b) $t=0.3$ (c) $t=0.46$ (d) $t=0.8$

Figure 7: Example 4.5: electromagnetic wave propagation induced by point
source in 3D. (a)-(d) snapshots of $B_{2}$ at $t=0.1,0.3,0.46,0.8$. The
simulations are obtained with $N=100$ and $\tau=0.02.$

##  5 Concluding Remarks

In this paper, we have developed fast and accurate divergence-free spectral
methods using generalized Jacobi polynomials for curl-curl problems in two and
three dimensions. The fast solver is integrated with two critical components:
(i). the matrix-free preconditioned Krylov subspace iterative method; and (ii)
the efficient preconditioner driven by a fully diagonalizable auxiliary
problem. With these, the new iteration method finds remarkably efficient in
solving the curl-curl problem, where the number of iterations can even be
reduced to $0$ for (7) with the smooth solution. More importantly, we prove
rigorously that the dimensions of the invariant subspace of the resultant
preconditioned linear system of the divergence-free spectral method for the
dominate eigenvalue $1$ is $(N-3)^{2}$ and $2(N-3)^{3}$, for two- and three-
dimensional problems with $(N-1)^{2}$ and $2(N-1)^{3}$ unknowns, respectively.
Finally, numerical experiments illustrate the efficiency of the proposed
scheme. Future studies along this line will include the extension of the
current method to unstructured meshes using spectral-element method for two-
and three-dimensional problems.

##  Appendix A Proof of Lemma 2.3

The derivation can be divided into three cases with respect to the range of
$m$.

$\underline{\text{\bf Case\,1}:1\leq m\leq 2,\,\,1\leq n\leq N-1:}$ Let us fix
$m=1$ first. By the analytic expression of $M$ in equation (40) and the
definition of $S$ in equation (42), one obtains

$\begin{split}&M_{11}=\frac{2}{5},\;\;\;M_{13}=-\frac{1}{5\sqrt{21}},\\\
&S_{1n}=\sqrt{\frac{3}{2}}\sqrt{\frac{2n+1}{2}}\int_{-1}^{1}L_{1}^{\prime}L_{n}^{\prime}\,{\rm
d}\xi,\;\;\;S_{1n}=\sqrt{\frac{7}{2}}\sqrt{\frac{2n+1}{2}}\int_{-1}^{1}L_{3}^{\prime}L_{n}^{\prime}\,{\rm
d}\xi.\end{split}$

Due to the fact that $M$ is penta-diagonal, one has

(A.95)
$(MS)_{1n}=M_{11}S_{1n}+M_{13}S_{3n}=\frac{\sqrt{3(2n+1)}}{5}\int_{-1}^{1}\big{(}L_{1}^{\prime}(\xi)-L_{3}^{\prime}(\xi)/6\big{)}L_{n}^{\prime}(\xi)\,{\rm
d}\xi.$

By using integration by parts and orthogonality of Legendre polynomials,
equation (A.95) leads to

(A.96)
$\begin{split}(MS)_{1n}&=\frac{\sqrt{3(2n+1)}}{4}\Big{\\{}[(1-\xi^{2})L_{n}(\xi)]\big{|}_{-1}^{1}+2\int_{-1}^{1}L_{1}(\xi)L_{n}(\xi)\,{\rm
d}\xi\Big{\\}}\\\
&=\frac{\sqrt{3(2n+1)}}{3}\delta_{1n}=\delta_{1n}.\end{split}$

Following the same procedure, one can obtain that

(A.97) $\displaystyle(MS)_{2n}=$ $\displaystyle
M_{22}S_{2n}+M_{24}S_{4n}=\frac{1}{14}\sqrt{\frac{2n+1}{5}}\int_{-1}^{1}\Big{(}\frac{10}{3}L_{2}^{\prime}(\xi)-L_{4}^{\prime}(\xi)\Big{)}L_{n}^{\prime}(\xi)\,{\rm
d}\xi$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2n+1}{5}}\frac{5}{2}\int_{-1}^{1}L_{2}(\xi)L_{n}(\xi)\,{\rm
d}\xi=\sqrt{\frac{2n+1}{5}}\delta_{2n}=\delta_{2n}.$

$\underline{\text{\bf Case}\,2:3\leq m\leq N-3,\,\,1\leq n\leq N-1:}$ Since
$M$ is penta-diagonal, $(MS)_{mn}$ reduces to

(A.98) $(MS)_{mn}=M_{m(m-2)}S_{(m-2)n}+M_{mm}S_{mn}+M_{m(m+2)}S_{(m+2)n},$

among which

$\displaystyle
M_{m(m-2)}S_{(m-2)n}=\frac{1}{2}\sqrt{\frac{2n+1}{2m+1}}\int_{-1}^{1}\Big{(}\frac{-1}{2m-1}L_{m-2}^{\prime}(\xi)\Big{)}\,L_{n}^{\prime}(\xi)\,{\rm
d}\xi,$ $\displaystyle
M_{mm}S_{mn}=\frac{1}{2}\sqrt{\frac{2n+1}{2m+1}}\int_{-1}^{1}\Big{(}\frac{1}{2m-1}+\frac{1}{2m+3}\Big{)}L_{m}^{\prime}(\xi)\,L_{n}^{\prime}(\xi)\,{\rm
d}\xi,$ $\displaystyle
M_{m(m+2)}S_{(m+2)n}=\frac{1}{2}\sqrt{\frac{2n+1}{2m+1}}\int_{-1}^{1}\Big{(}\frac{-1}{2m+3}L_{m+2}^{\prime}(\xi)\Big{)}L_{n}^{\prime}(\xi)\,{\rm
d}\xi.$

By using the recurrence relation (43), we have that

$\frac{-1}{2m-1}L_{m-2}^{\prime}(\xi)+\Big{(}\frac{1}{2m-1}+\frac{1}{2m+3}\Big{)}L_{m}^{\prime}(\xi)+\frac{-1}{2m+3}L_{m+2}^{\prime}(\xi)=L_{m-1}(\xi)-L_{m+1}(\xi).$

Thus, it is direct to obtain

$\displaystyle(MS)_{mn}=\frac{1}{2}\sqrt{\frac{2n+1}{2m+1}}\int_{-1}^{1}\big{(}L_{m-1}(\xi)-L_{m+1}(\xi)\big{)}L_{n}^{\prime}(\xi)\,{\rm
d}\xi$
$\displaystyle=\frac{1}{2}\sqrt{\frac{2n+1}{2m+1}}\Big{\\{}\big{[}(L_{m-1}(\xi)-L_{m+1}(\xi))L_{n}(\xi)\big{]}_{-1}^{1}+\int_{-1}^{1}(L_{m+1}(\xi)-L_{m-1}(\xi))^{\prime}L_{n}(\xi)\,{\rm
d}\xi\Big{\\}}$
$\displaystyle=\frac{1}{2}\sqrt{\frac{2n+1}{2m+1}}\int_{-1}^{1}(2m+1)L_{m}(\xi)L_{n}(\xi)\,{\rm
d}\xi=\sqrt{\frac{2m+1}{2n+1}}\delta_{mn}=\delta_{mn},$

where integration by parts, the recurrence relation (43), the boundary values
of Legendre polynomials in (45), and the orthogonality relation (4) has been
applied.

$\underline{\text{\bf Case}\,3:N-2\leq m\leq N-1,\,\,1\leq n\leq N-1:}$ Fix
$m=N-2$ and recall that $M$ is penta-diagonal, one has

$\displaystyle(MS)_{(N-2)n}=M_{(N-2)(N-4)}S_{(N-4)n}+M_{(N-2)(N-2)}S_{(N-2)n}$
$\displaystyle=$
$\displaystyle\Big{(}M_{(N-2)(N-4)}S_{(N-4)n}+M_{(N-2)(N-2)}S_{(N-2)n}+M_{(N-2)N}S_{Nn}\Big{)}-M_{(N-2)N}S_{Nn}$
$\displaystyle=$
$\displaystyle\delta_{(N-2)n}+\frac{1}{4}\sqrt{\frac{2n+1}{2N-3}}\frac{n(n+1)}{2N-1}\Big{(}1+(-1)^{n+N}\Big{)},$

where equations (40)-(42) and the fact that the three-term formula in equation
(A.98) equals to the Dirac-delta function have been used. Again, by invoking
(40)-(42), one obtains that for $m=N-1$,

$\displaystyle(MS)_{(N-1)n}=M_{(N-1)(N-3)}S_{(N-3)n}+M_{(N-1)(N-1)}S_{(N-1)n}$
$\displaystyle=$
$\displaystyle\Big{(}M_{(N-1)(N-3)}S_{(N-3)n}+M_{(N-1)(N-1)}S_{(N-1)n}+M_{(N-1)(N+1)}S_{(N+1)n}\Big{)}$
$\displaystyle-M_{(N-1)(N+1)}S_{(N+1)n}$ $\displaystyle=$
$\displaystyle\delta_{(N-1)n}+\frac{1}{4}\sqrt{\frac{2n+1}{2N-1}}\frac{n(n+1)}{2N+1}\Big{(}1+(-1)^{n+N+1}\Big{)}.$

This ends the proof.

## Acknowledgment

Helpful discussions with Professor Huiyuan Li (Institute of Software, Chinese
Academy of Sciences) are gratefully acknowledged.

## References

  * [1] I.M. Babuska and S.A. Sauter.  Is the pollution effect of the fem avoidable for the helmholtz equation considering high wave numbers?  SIAM J. Numer. Anal., 34(6):2392–2423, 1997.
  * [2] D. Boffi, F. Brezzi, and M. Fortin.  Mixed Finite Element Methods and Applications, volume 44.  Springer, 2013.
  * [3] J. Brackbill.  Fluid modeling of magnetized plasmas.  Space Sci. Rev., 42(1-2):153–167, 1985.
  * [4] J. Brackbill and D. Barnes.  The effect of nonzero $\nabla\cdot b$ on the numerical solution of the magnetohydrodynamic equations.  J. Comp. Phys., 35(3):426–430, 1980.
  * [5] B. Buzbee, G. Golub, and C. Nielson.  On direct methods for solving Poisson’s equations.  SIAM J. Numer. Anal., 7(4):627–656, 1970.
  * [6] W. Cai, J. Wu, and J. Xin.  Divergence-free $\boldsymbol{H}({\rm div})$-conforming hierarchical bases for magnetohydrodynamics (mhd).  Commun. Math. Stat., 1(1):19–35, 2013.
  * [7] F. Chen.  Introduction to Plasma Physics and Controlled Fusion, volume 1.  Springer, 1984.
  * [8] M. Costabel and M. Dauge.  Weighted regularization of Maxwell equations in polyhedral domains: A rehabilitation of nodal finite elements.  Numer. Math., 93:239–277, 2002.
  * [9] P. Cummings and X. Feng.  Sharp regularity coefficient estimates for complex-valued acoustic and elastic Helmholtz equations.  Math. Mod. Methods Appl. Sci., 16(01):139–160, 2006.
  * [10] E. Cyr, J. Shadid, R. Tuminaro, R. Pawlowski, and L. Chacón.  A new approximate block factorization preconditioner for two-dimensional incompressible (reduced) resistive mhd.  SIAM J. Sci. Comput., 35(3):B701–B730, 2013.
  * [11] P.A. Davidson.  An introduction to magnetohydrodynamics, 2002.
  * [12] H. Duan, S. Li, R. Tan, and W. Zheng.  A delta-regularization finite element method for a double curl problem with divergence-free constraint.  SIAM J. Numer. Anal., 50(6):3208–3230, 2012.
  * [13] G. Golub and C. Van Loan.  Matrix Computations.  JHU press, 2013.
  * [14] D. Gottlieb and S.A. Orszag.  Numerical Analysis of Spectral Methods: Theory and Applications.  SIAM, 1977.
  * [15] B. Guo and Y. Jiao.  Spectral method for Navier–Stokes equations with non-slip boundary conditions by using divergence-free base functions.  J. Sci. Comput., 66(3):1077–1101, 2016.
  * [16] B. Guo, J. Shen, and L.-L. Wang.  Generalized Jacobi polynomials/functions and their applications.  Appl. Numer. Math., 59(5):1011–1028, 2009.
  * [17] D.B. Haidvogel and T. Zang.  The accurate solution of Poisson’s equation by expansion in chebyshev polynomials.  J. Comput. Phys., 30(2):167–180, 1979.
  * [18] C. Hazard and M. Lenoir.  On the solution of time-harmonic scattering problems for Maxwell’s equations.  SIAM J. Math. Anal., 27(6):1597–1630, 1996.
  * [19] R. Hiptmair, L. Li, S. Mao, and W. Zheng.  A fully divergence-free finite element method for magnetohydrodynamic equations.  Math. Mod. Methods Appl. Sci., 28(04):659–695, 2018.
  * [20] R. Hiptmair, A. Moiola, and I. Perugia.  Stability results for the time-harmonic Maxwell equations with impedance boundary conditions.  Math. Mod. Methods Appl. Sci., 21(11):2263–2287, 2011.
  * [21] S. Jardin.  Computational Methods in Plasma Physics.  CRC press, 2010.
  * [22] J.D. Joannopoulos, S.G. Johnson, J.N. Winn, and R.D. Meade.  Photonic Crystals: Molding the Flow of Light - Second Edition.  Princeton University Press, 2011.
  * [23] C.T. Kelley.  Iterative Methods for Linear and Nonlinear Equations.  SIAM, 1995.
  * [24] P. Lax.  Linear Algebra and Its Applications, volume 78.  John Wiley & Sons, 2007.
  * [25] F. Li and C.-W. Shu.  Locally divergence-free discontinuous Galerkin methods for MHD equations.  J. Sci. Comput., 22(1):413–442, 2005.
  * [26] L. Ma, J. Shen, L.-L. Wang, and Z. Yang.  Wavenumber explicit analysis for time-harmonic Maxwell equations in a spherical shell and spectral approximations.  IMA J. Numer. Anal., 38(2):810–851, 05 2017.
  * [27] P.A. Markowich, C.A. Ringhofer, and C. Schmeiser.  Semiconductor Equations.  Springer Science & Business Media, 2012.
  * [28] P. Monk.  Finite Element Methods for Maxwell’s Equations.  Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2003.
  * [29] L. Mu, J. Wang, X. Ye, and S. Zhang.  A discrete divergence free weak Galerkin finite element method for the stokes equations.  Appl. Numer. Math., 125:172–182, 2018.
  * [30] Y. Saad.  Iterative Methods for Sparse Linear Systems.  SIAM, 2003.
  * [31] J. Shadid, R. Pawlowski, J. Banks, L. Chacón, P. Lin, and R. Tuminaro.  Towards a scalable fully-implicit fully-coupled resistive MHD formulation with stabilized fe methods.  J. Comput. Phys., 229(20):7649–7671, 2010.
  * [32] J. Shen.  Efficient spectral-Galerkin method I. direct solvers for second- and fourth-order equations by using Legendre polynomials.  SIAM J. Sci. Comput., 15:1489–1505, 1994.
  * [33] J. Shen, T. Tang, and L.-L. Wang.  Spectral Methods: Algorithms, Analysis and Applications, volume 41 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, Heidelberg, 2011.
  * [34] J. Shen and L.-L. Wang.  Spectral approximation of the Helmholtz equation with high wave numbers.  SIAM J. Numer. Anal., 43(2):623–644, 2005.
  * [35] J. Shen and L.L. Wang.  Analysis of a spectral-Galerkin approximation to the Helmholtz equation in exterior domains.  SIAM J. Numer. Anal., 45(5):1954–1978, 2007.
  * [36] G. Strang and G.J. Fix.  An Analysis of The Finite Element Method.  Prentice-Hall Series in Automatic Computation. Prentice-Hall, Inc., Englewood Cliffs, N.J., 1973.
  * [37] A. Taflove.  Application of the finite-difference time-domain method to sinusoidal steady-state electromagnetic-penetration problems.  IEEE Trans. Elect. Compat., (3):191–202, 1980.
  * [38] L.N. Trefethen and D. Bau.  Numerical Linear Algebra, volume 181.  SIAM, 2022.
  * [39] X. Ye and C.A. Hall.  A discrete divergence-free basis for finite element methods.  Numer. Algor., 16(3):365–380, 1997.
  * [40] K. Yee.  Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media.  IEEE Trans. Antennas and Prop., 14(3):302–307, 1966.
  * [41] J. Zhang, L.-L. Wang, H. Li, and Z. Zhang.  Optimal spectral schemes based on generalized prolate spheroidal wave functions of order-1-1.  J. Sci. Comput., 70:451–477, 2017.
  * [42] J. Zhang, W. Wang, X. Yang, et al.  Double-cone ignition scheme for inertial confinement fusion.  Phil. Trans. Royal Soc. A, 378(2184):20200015, 2020.