remarkRemark hypothesisHypothesis claimClaim Hierarchies of semiclassical
Jacobi polynomialsI. P. A. Papadopoulos, T. S. Gutleb, R. M. Slevinsky, S.
Olver

# Building hierarchies of semiclassical Jacobi polynomials for spectral
methods in annuli ††thanks: Submitted DATE. This work was completed with the
support of the EPSRC grant EP/T022132/1 “Spectral element methods for
fractional differential equations, with applications in applied analysis and
medical imaging” and the Leverhulme Trust Research Project Grant RPG-2019-144
“Constructive approximation theory on and inside algebraic curves and
surfaces”.

Ioannis P. A. Papadopoulos Department of Mathematics, Imperial College London,
London, UK
<EMAIL_ADDRESS>Timon S. Gutleb Mathematical
Institute, University of Oxford, UK<EMAIL_ADDRESS>Richard
M. Slevinsky Department of Mathematics, University of Manitoba, Canada,
<EMAIL_ADDRESS>Sheehan Olver Department of Mathematics,
Imperial College London, UK<EMAIL_ADDRESS>
###### Abstract

We discuss computing with hierarchies of families of (potentially weighted)
semiclassical Jacobi polynomials which arise in the construction of
multivariate orthogonal polynomials. In particular, we outline how to build
connection and differentiation matrices with optimal complexity and compute
analysis and synthesis operations in quasi-optimal complexity. We investigate
a particular application of these results to constructing orthogonal
polynomials in annuli, called the generalised Zernike annular polynomials,
which lead to sparse discretisations of partial differential equations. We
compare against a scaled-and-shifted Chebyshev–Fourier series showing that in
general the annular polynomials converge faster when approximating smooth
functions and have better conditioning. We also construct a sparse spectral
element method by combining disk and annulus cells, which is highly effective
for solving PDEs with radially discontinuous variable coefficients and data.

###### keywords:

semiclassical orthogonal polynomials, multivariate orthogonal polynomials,
spectral methods, disk, annulus

33C45, 33C50, 65D05, 65N35

##  1 Introduction

Semiclassical Jacobi polynomials are univariate polynomials orthogonal with
respect to the weight $x^{a}(1-x)^{b}(t-x)^{c}$ on $[0,1]$ where $t>1$ and
$a,b>-1$. The semiclassical Jacobi polynomials are used to give explicit
expressions for the so-called generalised Zernike annular polynomials. These
are multivariate orthogonal polynomials ( in $x$ and $y$) on the annulus
$\Omega_{\rho}=\\{(x,y)\in\mathbb{R}^{2}:0<\rho\leq r\leq 1\\}$, where
$r^{2}=x^{2}+y^{2}$, orthogonal with respect to a Jacobi-like inner product:

(1)
$\displaystyle\left\langle{f,g}\right\rangle_{\rho,(a,b)}=\iint_{\Omega_{\rho}}f(x,y)g(x,y)(1-r^{2})^{a}(r^{2}-\rho^{2})^{b}\,{\rm
d}x{\rm d}y,$

where $a,b>-1$. We shall utilise the connection between semiclassical Jacobi
polynomials and Zernike annular polynomials to introduce quasi-optimal
complexity means for discretising and  solving partial differential equations
in annuli. This requires the construction of operators associated with a
hierarchy of semiclassical Jacobi polynomials where $c\in\\{0,1,\dots,C\\}$.

Tatian [35] and Mahajan [23] were the first to introduce the (non-generalised)
Zernike annular polynomials orthogonal with respect to the unweighted
$L^{2}$-inner product, $a=b=0$, on the annulus. Sometimes referred to as
Tatian–Zernike or fringe-Tatian polynomials, these have proven very popular in
the optics and other communities, cf. [11, 12, 22, 31]. Zernike annular
polynomials also underpin the recently introduced gyroscopic orthogonal
polynomials [14], which are used for solving PDEs in cylinders of varying
heights. In all the aforementioned works, the synthesis, analysis, connection
and differentiation operators are constructed via the Christoffel–Darboux
formula combined with quadrature, which does not achieve optimal complexity.
More specifically, using Christoffel–Darboux, which is equivalent to the
Cholesky factorisation technique described in Section 3, they calculate $N$
Jacobi matrices in $\mathcal{O}(N^{2})$ flops. For our applications this is
numerically unstable when $N\gg 0$ and alternative methods via QR
factorisations are preferable cf. Section 3. They also calculate the Laplacian
and identity operators in $\mathcal{O}(N^{3})$ flops as opposed to our methods
which require $\mathcal{O}(N^{2})$ flops.

By constructing multivariate orthogonal polynomials with respect to the non-
uniform weight in Eq. 1, one may build sparse spectral methods to solve
partial differential equations (PDEs) on such regimes. In Section 7, we focus
on using these multivariate orthogonal polynomials to solve PDEs on disks and
annuli via a _sparse_ spectral method and compare against a method based on
the Chebyshev–Fourier series, though the results extend to the setting of
gyroscopic orthogonal polynomials. We also construct a spectral element method
for problems with radial direction discontinuities in the variable
coefficients and right-hand side, as exemplified in Fig. 1. The boundary
conditions and continuity across cells are enforced via a tau-method [18, 10].

(a) Right-hand side $f_{1}$ Eq. 64.

(b) Right-hand side $f_{2}$ Eq. 65.

(c) Solution $u$ Eq. 61.

Figure 1: Plots of the right-hand sides with discontinuities in the radial
direction and the corresponding solution of Eq. 62–Eq. 63. The setup of the
problem is given in Section 7.3. These functions are resolved to machine
precision utilising the spectral element method described in Section 6 based
on Zernike and Zernike annular polynomials.

Spectral methods for the disk/ball [9, 37, 39, 8, 24, 20, 6, 15] and sphere
[38, 19, 36] have been thoroughly studied. On the annulus many methods utilise
variations of a Chebyshev–Fourier series, e.g. [25]. Barakat [7] constructs a
basis for the annulus but it is not polynomial in Cartesian coordinates. In
2011, Boyd and Yu [9] compared seven spectral methods for a disk. Out of the
seven, they found that Zernike polynomials and Chebyshev–Fourier series had
the best approximation properties where, for non-trivial examples, Zernike
polynomials required half the degrees of freedom of the next best method. This
motivates the construction of multivariate orthogonal polynomials for domains
of interest. Until recently, a common critique of Zernike polynomials was that
no fast transform for the radial direction was known [9, Sec. 6.1]. However,
this was resolved by Slevinsky [34, 32], see also [27]. This technique was
extended to generalised Zernike annular polynomials in [17, Sec. 4.4].

Note that as Zernike (annular) polynomials and the Chebyshev–Fourier series
feature a Fourier mode component, any operator that commutes with rotations
(such as the Laplacian or identity operator) decouples across the Fourier
modes. That is, we can decompose a two-dimensional PDE solve on the annulus
into $2M+1$ one-dimensional solves where $M$ is the highest Fourier mode
considered. For each one-dimensional solve, the scaled-and-shifted
Chebyshev–Fourier series results in an almost-banded matrix. Aside from the
two dense rows associated with the boundary conditions, one recovers matrices
with a bandwidth of five and nine for the Poisson and Helmholtz equations,
respectively. By contrast the Zernike annular polynomials result in a
tridiagonal and a pentadiagonal matrix per Fourier mode for the Poisson and
Helmholtz equations, respectively, as depicted in Fig. 2: a smaller bandwidth.
Furthermore, we note that the Laplacian matrix for the Zernike annular
polynomials is significantly better conditioned than the Laplacian matrix of
the Chebyshev–Fourier series, for increasing Fourier mode $m$, as plotted in
Fig. 3.

Due to the decoupling of the considered PDEs across Fourier modes, we require
$\mathcal{O}(N^{2})$ for the linear system solve. Thus the overall complexity,
from approximation of the right-hand side through to evaluating the
approximated solution on a grid is $\mathcal{O}(N^{2}\log(N))$.

(a) All modes. (b) $0^{\mathrm{th}}$-Fourier mode. (c) All modes. (d)
$0^{\mathrm{th}}$-Fourier mode.

Figure 2: Spy plots of the matrices after discretising the Laplacian,
$\Delta$, ((a) and (b)) and the Helmholtz operator $(\Delta+\mathcal{I})$ ((c)
and (d)) on the 2D annulus $\Omega_{1/2}$ truncated at polynomial degree 20
with generalised Zernike annular polynomials. (a) and (c) show all the Fourier
modes, whereas (b) and (d) show the matrix when one decouples the Fourier
modes and focuses on the $m=0$ mode in isolation.

(a) Condition number. (b) Matrix size.

Figure 3: The condition number and matrix size of the Laplacian matrix, for
increasing Fourier mode $m$, induced by the Chebyshev–Fourier series and the
Zernike annular polynomials on the domain $\Omega_{1/2}$ with truncation
degree $N=100$. In both cases the Laplacian matrices are preconditioned with
the inverse of the diagonal of the original matrix. We observe that the
Laplacian matrices for the Zernike annular polynomials are significantly
better conditioned and smaller in size.

We describe the scaled-and-shifted Chebyshev–Fourier series in Section 2. In
Section 3 we introduce the semiclassical Jacobi polynomials and discuss
factorisation techniques for deriving the relationships between different
family parameters. In Section 4, we introduce the generalised Zernike annular
polynomials and derive various relationships between families of parameters
that may be computed with optimal complexity. We outline the sparse spectral
methods for solving the Helmholtz equation in Section 5. In Section 6 we
construct a spectral element method where the cells are an inner disk and
concentric annuli. Finally, in Section 7 we consider various PDE examples and
compare the spectral methods.

######  Remark 1.1 (Weak formulation).

A major benefit of the generalised Zernike (annular) polynomials is that they
are easily amendable to discretising PDEs posed in weak form. This allows one
to construct sparse spectral element methods for variational problems. A
thorough construction of such a method is beyond the scope of this work,
however, we make some comments in Section 4.

######  Remark 1.2.

Zernike annular polynomials also allow one to construct random functions on
the annulus via the methods of Filip et al. [16].

######  Remark 1.3 (Generalisation to 3D).

Just as generalised Zernike polynomials on the disk can be extended to the
ball [13, Prop. 5.2.1], orthogonal polynomials on annuli naturally extend to
higher dimensional spherical shells with spherical harmonics in place of
Fourier modes. For simplicity we restrict our attention to 2D.

##  2 Scaled-and-shifted Chebyshev–Fourier series

The Chebyshev–Fourier series is constructed via a tensor product of Chebyshev
polynomials of the first kind (denoted $T_{n}(x)$, $n\geq 0$) [26, Sec. 18.3]
in the radial direction with the Fourier series in the angular direction.

Following [27], we use quasimatrix notation, i.e.

(2) $\displaystyle{\bf T}(x)\coloneqq{\begin{pmatrix}\,T_{0}(x)\ |\ T_{1}(x)\
|\ T_{2}(x)\ |\ \cdots\,\end{pmatrix}}.$

This is convenient for expressing recurrence relationships. For example, the
three term recurrence is expressed as

(3) $\displaystyle x{\bf T}(x)={\bf T}(x)\begin{pmatrix}0&1/2&&&\\\
1&0&1/2&&\\\ &1/2&0&1/2&\\\ &&\ddots&\ddots&\ddots\end{pmatrix}.$

Consider the annulus
$\Omega_{\rho}=\\{(x,y)\in\mathbb{R}^{2}:0<\rho\leq|(x,y)|\leq 1\\}$. Define

(4) $\displaystyle r_{\rho}=\frac{2}{1-\rho}\left(r-\frac{1+\rho}{2}\right).$

Let $F_{m,j}$ denote the Fourier series at mode $m$ where $j\in\\{0,1\\}$ if
$m\in\\{1,2,\dots\\}$, else $j=1$ if $m=0$. For clarity, $F_{0,1}(\theta)=1$
and for $m\geq 1$,

(5) $\displaystyle
F_{m,j}(\theta)\coloneqq\begin{cases}\sin(m\theta)&\text{if}\;j=0,\\\
\cos(m\theta)&\text{if}\;j=1.\end{cases}$

Then the scaled-and-shifted Chebyshev–Fourier series on the annulus is given
by the tensor product:

(6) $\displaystyle\begin{split}&{\bf T}(r_{\rho})\otimes{\bf F}(\theta)\\\
&\coloneqq\left(T_{0}(r_{\rho})F_{0,1}(\theta)\,|\,T_{1}(r_{\rho})F_{0,1}(\theta)\,|\,\cdots\,|\,T_{0}(r_{\rho})F_{1,0}(\theta)\,|\,\cdots\,|\,T_{2n}(r_{\rho})F_{n,0}(\theta)\,|\,\cdots\right).\end{split}$

An expansion in this basis takes the form:

$\displaystyle u(r,\theta)$
$\displaystyle=\sum_{n=0}^{\infty}u_{n,0,0}T_{n}(r_{\rho})+\sum_{n=0}^{\infty}\sum_{m=1}^{\infty}\left[u_{n,m,0}T_{n}(r_{\rho})\sin(m\theta)+u_{n,m,1}T_{n}(r_{\rho})\cos(m\theta)\right]$
$\displaystyle=\left[{\bf T}(r_{\rho})\otimes{\bf
F}(\theta)\right]\begin{pmatrix}u_{0,0,0}&u_{1,0,0}&\cdots&u_{2n,m,j}&\cdots\end{pmatrix}^{\top}=\left[{\bf
T}(r_{\rho})\otimes{\bf F}(\theta)\right]{\bf u}.$

If we rearrange the coefficients in ${\bf u}$ into a matrix such that

(7) $\displaystyle
U=\begin{pmatrix}u_{0,0,1}&u_{0,1,0}&u_{0,1,1}&u_{0,2,0}&\cdots\\\
u_{1,0,1}&u_{1,1,0}&u_{1,1,1}&u_{1,2,0}&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{pmatrix},$

then $\left[{\bf T}(r_{\rho})\otimes{\bf F}(\theta)\right]{\bf u}={\bf
T}(r_{\rho})U{\bf F}(\theta)^{\top}$. The choice of truncation degree $N$ and
Fourier mode $M$ are independent and may be custom chosen for each problem
individually. In this work we pick $M=N$ resulting in the truncated
coefficient matrix $U$ of size $(N+1)\times(2N+1)$.

##  3 Semiclassical Jacobi polynomials

We denote the orthonormalised Jacobi polynomials by $P^{(a,b)}_{n}(x)$, $n\geq
0$ which are orthonormal with respect to the weight $(1-x)^{a}(1+x)^{b}$ such
that $a,b>-1$ [26, Sec. 18.3]. The term _weighted Jacobi polynomials_ refers
to the Jacobi polynomial $P_{n}^{(a,b)}(x)$ multiplied by its weight, i.e.
$W^{(a,b)}_{n}(x)\coloneqq(1-x)^{a}(1+x)^{b}P_{n}^{(a,b)}(x)$.

The building blocks for the Zernike annular polynomials are the so-called
semiclassical Jacobi polynomials. Recall that these are univariate orthogonal
polynomials on the interval $(0,1)$ with respect to the inner product

(8) $\displaystyle\int_{0}^{1}f(x)g(x)x^{a}(1-x)^{b}(t-x)^{c}{\rm d}x,$

where $t>1$ and $a,b>-1$. First introduced by Magnus111The weight considered
by Magnus is the slightly different $(1-x)^{a}x^{b}(t-x)^{c}$, which has the
unfortunate property that the singularities are neither listed left-to-right
or right-to-left. We have changed the order to be left-to-right. in [21, Sec.
5], we denote the orthonormal semiclassical Jacobi polynomials as

$Q^{t,(a,b,c)}_{n}(x)=k_{n}x^{n}+O(x^{n-1}),$

where, for concreteness, $k_{n}>0$. Note that, when $c=0$, these become
scaled-and-shifted orthonormalised Jacobi polynomials and we drop the $t$
dependence. That is, we have for any $t$,

(9) $\displaystyle
Q_{n}^{t,(a,b,0)}(x)=(-1)^{n}2^{(a+b+1)/2}P_{n}^{(a,b)}(1-2x)=2^{(a+b+1)/2}P_{n}^{(b,a)}(2x-1).$

Many beautiful connections between semiclassical Jacobi polynomials and
Painlevé equations are known, and we refer the readers to [21].

######  Remark 3.1.

The semiclassical Jacobi polynomials are a special case of the generalised
Jacobi polynomials found in [14, Sec. 3].

We require recurrence relationships between families of the Zernike annular
polynomials in the construction of spectral methods. These recurrence
relationships are derived from recurrence relationships that relate
hierarchies of semiclassical Jacobi polynomial families. We denote the
connection matrix mapping the parameter family from $(\alpha,\beta,\gamma)$ to
$(a,b,c)$ by $R^{t,(\alpha,\beta,\gamma)}_{(a,b,c)}$. To indicate that maps
are from a (partially) weighted family to another (partially) weighted family,
we utilise the Roman script $\mathrm{a}$, $\mathrm{b}$, and $\mathrm{c}$ or
combinations thereof. For instance:

(10) $\displaystyle\begin{split}{\bf Q}^{t,(a,b,c)}(x)&={\bf
Q}^{t,(\alpha,\beta,\gamma)}(x)R^{t,(\alpha,\beta,\gamma)}_{(a,b,c)},\\\
x^{a}{\bf Q}^{t,(a,b,c)}&=x^{\alpha}{\bf
Q}^{t,(\alpha,\beta,\gamma)}R^{t,(\alpha,\beta,\gamma)}_{\mathrm{a},(a,b,c)},\\\
x^{a}(1-x)^{b}(t-x)^{c}{\bf
Q}^{t,(a,b,c)}(x)&=x^{\alpha}(1-x)^{\beta}(t-x)^{\gamma}{\bf
Q}^{t,(\alpha,\beta,\gamma)}(x)R^{t,(\alpha,\beta,\gamma)}_{\mathrm{abc},(a,b,c)}.\end{split}$

Similarly, we denote the differentiation connection matrix between families of
semiclassical Jacobi polynomials by $D^{t,(\alpha,\beta,\gamma)}_{(a,b,c)}$,
e.g.

(11) $\displaystyle\begin{split}\frac{\mathrm{d}}{{\rm d}x}{\bf
Q}^{t,(a,b,c)}(x)&={\bf
Q}^{t,(\alpha,\beta,\gamma)}(x)D^{t,(\alpha,\beta,\gamma)}_{(a,b,c)},\\\
\frac{\mathrm{d}}{{\rm d}x}x^{a}{\bf Q}^{t,(a,b,c)}&=x^{\alpha}{\bf
Q}^{t,(\alpha,\beta,\gamma)}D^{t,(\alpha,\beta,\gamma)}_{\mathrm{a},(a,b,c)},\\\
\frac{\mathrm{d}}{{\rm d}x}[x^{a}(1-x)^{b}(t-x)^{c}{\bf
Q}^{t,(a,b,c)}(x)]&=x^{\alpha}(1-x)^{\beta}(t-x)^{\gamma}{\bf
Q}^{t,(\alpha,\beta,\gamma)}(x)D^{t,(\alpha,\beta,\gamma)}_{\mathrm{abc},(a,b,c)}.\end{split}$

There exist choices of $(\alpha,\beta,\gamma)$ and $(a,b,c)$ such that the
matrices $R^{t,(\alpha,\beta,\gamma)}_{\boldsymbol{\cdot},(a,b,c)}$ and
$D^{t,(\alpha,\beta,\gamma)}_{\boldsymbol{\cdot},(a,b,c)}$ are sparse. In the
next three subsections, we introduce the factorisation techniques that allow
one to apply analysis and synthesis operators with quasi-optimal complexity as
well as compute Jacobi, connection and differentiation matrices.

###  3.1 Connection matrices

Given the quasimatrix of a semiclassical Jacobi family evaluated at a point,
${\bf Q}^{t,(a,b,c)}(x)$, the goal is to find the (truncated) infinite-
dimensional matrix of coefficients in Eq. 10 that relates the evaluation to a
different semiclassical Jacobi family evaluated at the same point. In recent
work by Gutleb et al. [17], the authors analyse six infinite-dimensional
matrix factorisations for this purpose. The main contribution of their work is
summarised in [17, Tab. 1]. We replicate the relevant results in Theorem 3.2.

######  Theorem 3.2 (Tab. 1 in [17]).

Let $w(x)$ denote a nonnegative and bounded weight on $\mathbb{R}$ and
consider the measure $\mathrm{d}\mu$ such that both $\mathrm{d}\mu$ and
$w(x)\mathrm{d}\mu$ are positive Borel measures on the real line whose support
contains an infinite number of points and has finite moments. Suppose that
${\bf P}(x)$ and ${\bf Q}(x)$ are the quasimatrices of two families of
orthonormal polynomials with respect to $\mathrm{d}\mu$ and
$w(x)\mathrm{d}\mu$, respectively. Let $X_{P}$ denote the Jacobi matrix of
${\bf P}(x)$, i.e. $x{\bf P}(x)={\bf P}(x)X_{P}$. Then

(12) $\displaystyle{\bf P}(x)={\bf Q}(x)R\;\;\text{and}\;\;w(x){\bf Q}(x)={\bf
P}(x)R^{\top}\;\;\text{where}\;\;w(X_{P})=R^{\top}R.$

Here $R$ is upper-triangular, i.e. $R^{\top}R$ is a Cholesky factorisation.
Moreover, if $w(x)$ is a polynomial, then $w(X_{P})$ and $R$ are banded.
Furthermore, if $\sqrt{w}(x)$ is also a polynomial, then

(13) $\displaystyle{\bf P}(x)={\bf Q}(x)R\;\;\text{and}\;\;\sqrt{w}(x){\bf
Q}(x)={\bf P}(x)Q\;\;\text{where}\;\;\sqrt{w}(X_{P})=QR.$

Here $Q$ is orthogonal and $R$ is upper-triangular, i.e. they are the factors
of a positive-phase QR-factorisation (and are unique).

Suppose that $\alpha-a\geq 0$, $\beta-b\geq 0$, $\gamma-c\geq 0$ are integers.
Let $X_{t,(a,b,c)}$ denote the Jacobi matrix of the orthormalised
semiclassical Jacobi quasimatrix ${\bf Q}^{t,(a,b,c)}(x)$, i.e.

$\displaystyle x{\bf Q}^{t,(a,b,c)}(x)={\bf Q}^{t,(a,b,c)}(x)X_{t,(a,b,c)}.$

Suppose one computes the following Cholesky factorisation

(14) $\displaystyle X_{t,(a,b,c)}^{\alpha-a}(I-X_{t,(a,b,c)})^{\beta-b}(tI-
X_{t,(a,b,c)})^{\gamma-c}=R^{\top}R,$

where $I$ is the identity matrix of conforming size. Then Eq. 12 in Theorem
3.2 reveals that $R^{t,(\alpha,\beta,\gamma)}_{(a,b,c)}=R$ and
$R^{t,(a,b,c)}_{\mathrm{abc},(\alpha,\beta,\gamma)}=R^{\top}$. Furthermore, we
note that $R$ is upper-triangular with bandwidth
$(\alpha-a)+(\beta-b)+(\gamma-c)+1$. For more general choices of
$\alpha,\beta,a,b>-1$, and $\gamma,c\in\mathbb{R}$, we refer the interested
reader to [17].

When $\alpha\gg a$, $\beta\gg b$, or $\gamma\gg c$, the matrix on the left-
hand side of Eq. 14 becomes increasingly less positive-definite and eventually
a Cholesky factorisation fails due to numerical roundoff. In such a case, we
recommend computing the connection matrix via an intermediate semiclassical
Jacobi family. Alternatively, one may consider Eq. 13 in Theorem 3.2:

(15) $\displaystyle{\bf Q}^{t,(a,b,c)}(x)={\bf
Q}^{t,(a,b,c+2)}(x)R,\;\;(t-x){\bf Q}^{t,(a,b,c+2)}(x)={\bf
Q}^{t,(a,b,c)}(x)Q.$

In summary ${\bf Q}^{t,(a,b,0)}(x)={\bf Q}^{t,(a,b,1)}(x)R_{1}$ where
$R_{1}^{\top}R_{1}=tI-X_{t,(a,b,0)}$ and ${\bf Q}^{t,(a,b,0)}(x)={\bf
Q}^{t,(a,b,2)}(x)R_{2}$ where $QR_{2}=tI-X_{t,(a,b,0)}$. Connection matrices
for higher values of $c$ may be computed by chaining together these
factorisations.

###  3.2 Jacobi matrices

The goal is to compute the entries of the Jacobi matrix $X_{t,(a,b,c)}$ with
$\mathcal{O}(cN)$ complexity, where $c\in\mathbb{N}$ and $N$ is the truncation
degree. We first note that the main difficulty arises in the parameter $c$. If
$c=0$ then explicit formulae for the Jacobi matrix exist.

######  Remark 3.3.

Given the base Jacobi matrix $X_{t,(a,b,c)}$ for any $c\in\mathbb{R}$, $c\geq
0$, then one may compute the Jacobi matrix $X_{t,(a,b,c+m)}$, for any
$m\in\mathbb{N}$, with the following techniques. Thus we are not restricted to
the case where $c\in\mathbb{N}$.

The identities in Eq. 14 and Eq. 15 provide three alternatives for computing
the Jacobi matrix when $c\neq 0$, $c\in\mathbb{N}$:

  1. 1.

via the Cholesky factorisation;

  2. 2.

via the $R$ in a QR factorisation;

  3. 3.

via the $Q$ in a QR factorisation.

Via Cholesky. One may find the connection matrix between the case $c=0$ and
any other choice of $c\in\mathbb{N}$ utilising Eq. 14. If $c\in\mathbb{N}$
then the connection matrix is banded. We note that in principle one may
compute the connection matrix to any $c\in\mathbb{N}$ in one step. However,
one must find the Cholesky factorisation of $(tI-X_{t,(a,b,0)})^{c}$ which is
numerically indefinite for even moderate values of $c$. Thus it is _highly_
recommended that one proceeds in steps of one for $c$. If $R^{\top}R=tI-
X_{t,(a,b,0)}$, then

(16) $\displaystyle\begin{split}x{\bf Q}^{t,(a,b,1)}(x)&=x{\bf
Q}^{t,(a,b,0)}(x)R^{-1}={\bf Q}^{t,(a,b,0)}(x)X_{t,(a,b,0)}R^{-1}\\\ &={\bf
Q}^{t,(a,b,1)}(x)RX_{t,(a,b,0)}R^{-1}.\end{split}$

Hence $X_{t,(a,b,1)}=RX_{t,(a,b,0)}R^{-1}$. When $c>1$, we increment the $c$
parameter in steps of one until we reach the desired Jacobi matrix.

Via QR. If one utilises a QR factorisation via Eq. 15 then it is possible to
increment the $c$ parameter in steps of two. If $QR=tI-X_{t,(a,b,0)}$, then

(17) $\displaystyle\begin{split}x(t-x){\bf Q}^{t,(a,b,2)}(x)&=x{\bf
Q}^{t,(a,b,0)}(x)Q={\bf Q}^{t,(a,b,0)}(x)X_{t,(a,b,0)}Q\\\ &=(t-x){\bf
Q}^{t,(a,b,2)}(x)Q^{\top}X_{t,(a,b,0)}Q.\end{split}$

Hence, $X_{t,(a,b,2)}=Q^{\top}X_{t,(a,b,0)}Q$. Alternatively, note that

(18) $\displaystyle\begin{split}{\bf Q}^{t,(a,b,0)}(x)X_{t,(a,b,0)}&=x{\bf
Q}^{t,(a,b,0)}(x)=x{\bf Q}^{t,(a,b,2)}(x)R\\\ &={\bf
Q}^{t,(a,b,2)}(x)X_{t,(a,b,2)}R={\bf
Q}^{t,(a,b,0)}(x)R^{-1}X_{t,(a,b,2)}R.\end{split}$

Thus $X_{t,(a,b,2)}=RX_{t,(a,b,0)}R^{-1}$.

For either of the three routes, one computes the triple-matrix-product in
$\mathcal{O}(N)$ complexity. See Section 3.5 for more details. Thus to compute
the Jacobi matrix for an arbitrary parameter $c\in\mathbb{N}$ requires
$\mathcal{O}(cN)$ complexity. As mentioned before, when building multivariate
orthogonal polynomials, one typically requires a hierarchy of Jacobi matrices,
one for each value of $c\in\\{c_{0},c_{0}+1,c_{0}+2,\dots\\}$, $c_{0}\geq 0$.
Hence, in practice, these intermediate Jacobi matrices are required.

Out of these three approaches, the most theoretically stable method is the
$Q$-factor variant of the QR factorisation as $\|Q\|_{2}=1$, whereas the upper
triangular factors (and their inverses) in Cholesky and QR do not have unit
$2$-norm.

###  3.3 Analysis & synthesis

The analysis and synthesis operators map a function to the coefficient vector
of an expansion and vice versa. More precisely, let $\mathrm{d}\mu$ denote a
positive Borel measure on the real line whose support contains an infinite
number of points and has finite moments. Let $L^{2}(\mathbb{R},\mathrm{d}\mu)$
denote the space of square integrable functions with respect to the measure
$\mathrm{d}\mu$ and let $\ell^{2}$ denote the space of square summable
sequences. Suppose that ${\bf P}(x)$ denotes a quasimatrix that is orthogonal
with respect to $\mathrm{d}\mu$. Then the analysis operator $\mathcal{A}$ and
synthesis operator $\mathcal{S}$ satisfy:

(19) $\displaystyle\mathcal{A}:L^{2}(\mathbb{R},\mathrm{d}\mu)\to\ell^{2},$
$\displaystyle f\mapsto\int_{\mathbb{R}}{\bf
P}(x)^{\top}f(x)\mathrm{d}\mu(x)={\bf f},$ (20)
$\displaystyle\mathcal{S}:\ell^{2}\to L^{2}(\mathbb{R},\mathrm{d}\mu),$
$\displaystyle{\bf f}\mapsto{\bf P}{\bf f}=f.$

For a given basis, it is favourable if there exists a fast transform between
function evaluations and coefficient vector expansions, i.e. an ability to
compute the analysis and synthesis operators in a fast manner.

Suppose that the goal is to apply the synthesis operator to $(t-x)^{c/2}{\bf
Q}^{t,(a,b,c)}(x){\bf f}$. We consider a “square-root” weighted expansion as
it directly connects to the synthesis operator for Zernike annular polynomials
in Section 4.2. Suppose that $c$ is even. Then, via Eq. 15,

(21) $\displaystyle\begin{split}(t-x)^{c/2}{\bf Q}^{t,(a,b,c)}(x){\bf
f}&=(t-x)^{c/2-1}{\bf Q}^{t,(a,b,c-2)}(x)Q_{c-2}{\bf f}\\\ &={\bf
Q}^{t,(a,b,0)}(x)Q_{0}Q_{2}\cdots Q_{c-2}{\bf f}\\\ &={\bf
T}(1-2x)R_{T}S^{-1}_{(a,b)}Q_{0}Q_{2}\cdots Q_{c-2}{\bf f}.\end{split}$

Here $R_{T}$ is the upper-triangular conversion matrix as computed via the
Chebyshev–Jacobi transform [33], $S_{(a,b)}$ is the diagonal scaling as
defined by Eq. 9, and $Q_{\gamma}R_{\gamma}=tI-X_{t,(a,b,\gamma)}$. When $c$
is odd, one converts down to $(t-x)^{1/2}{\bf Q}^{t,(a,b,1)}(x)$ with QR and
then a final Cholesky factorisation to recover $(t-x)^{1/2}{\bf
Q}^{t,(a,b,0)}(x)$. We then convert ${\bf Q}^{t,(a,b,0)}(x)$ to ${\bf
T}(1-2x)$ and there exists a fast transform for the synthesis operator of
$(t-x)^{1/2}{\bf T}(1-2x)$. Thus the complexity for computing the Chebyshev
coefficients is $\mathcal{O}(cN)$ where $N$ is the truncation degree. The
synthesis operator may then be applied in $\mathcal{O}(N\log N)$ complexity
via the DCT. We note that the analysis operator is the reverse process of the
synthesis operator.

For an unweighted expansion, we instead utilise the $R$ factors in the QR
factorisation. If $c$ is even, we find that

(22) $\displaystyle\begin{split}{\bf Q}^{t,(a,b,c)}(x){\bf f}={\bf
T}(1-2x)R_{T}S^{-1}_{(a,b)}R_{0}^{-1}R_{2}^{-1}\cdots R_{c-2}^{-1}{\bf
f}.\end{split}$

###  3.4 Differentiation

The goal is to compute the entries of the banded matrix
$D^{t,(a+1,b+1,c+1)}_{(a,b,c)}$ and weighted counterparts. A core concept for
computing the differentiation matrices with low complexity is the following
sparsity result.

######  Theorem 3.4 (Theorem 2.20 in [17]).

The differentiation matrix in the right-hand side of

$\displaystyle\frac{\mathrm{d}}{{\rm d}x}{\bf Q}^{t,(a,b,c)}(x)={\bf
Q}^{t,(a+1,b+1,c+1)}(x)D^{t,(a+1,b+1,c+1)}_{(a,b,c)}$

only has nonzero entries on the first two super-diagonals.

Similar results holds for the (partially) weighted counterparts.

$D^{t,(a-1,b+1,c+1)}_{\mathrm{a},(a,b,c)},D^{t,(a-1,b-1,c+1)}_{\mathrm{ab},(a,b,c)}\hbox{,
and }D^{t,(a-1,b-1,c-1)}_{\mathrm{abc},(a,b,c)}$

all have a bandwidth of two and similarly for the other combination of
weights, b, c, bc, and ac.

Given the differentiation and Jacobi matrix for ${\bf Q}^{t,(a,b,c)}(x)$ we
note that one obtains the differentiation matrix for ${\bf
Q}^{t,(a,b,c+1)}(x)$ in $\mathcal{O}(N)$ complexity as follows:

(23) $\displaystyle\begin{split}\frac{\mathrm{d}}{{\rm d}x}{\bf
Q}^{t,(a,b,c+1)}(x)&=\frac{\mathrm{d}}{{\rm d}x}{\bf
Q}^{t,(a,b,c)}(x)(R^{t,(a,b,c+1)}_{(a,b,c)})^{-1}\\\ &={\bf
Q}^{t,(a+1,b+1,c+1)}(x)D^{t,(a+1,b+1,c+1)}_{(a,b,c)}(R^{t,(a,b,c+1)}_{(a,b,c)})^{-1}\\\
&={\bf
Q}^{t,(a+1,b+1,c+2)}(x)R^{t,(a+1,b+1,c+2)}_{(a+1,b+1,c+1)}D^{t,(a+1,b+1,c+1)}_{(a,b,c)}(R^{t,(a,b,c+1)}_{(a,b,c)})^{-1}.\end{split}$

The connection coefficients are Cholesky factors [17]:

(24)
$\displaystyle(R^{t,(a,b,c+1)}_{(a,b,c)})^{\top}R^{t,(a,b,c+1)}_{(a,b,c)}$
$\displaystyle=tI-X_{t,(a,b,c)},$ (25)
$\displaystyle(R^{t,(a+1,b+1,c+2)}_{(a+1,b+1,c+1)})^{\top}R^{t,(a+1,b+1,c+2)}_{(a+1,b+1,c+1)}$
$\displaystyle=tI-X_{t,(a+1,b+1,c+1)}.$

The unknown Jacobi matrix $X_{t,(a+1,b+1,c+1)}$ is obtained by noting

(26)
$\displaystyle(R^{t,(a+1,b+1,c+1)}_{(a,b,c)})^{\top}R^{t,(a+1,b+1,c+1)}_{(a,b,c)}$
$\displaystyle=X_{t,(a,b,c)}(I-X_{t,(a,b,c)})(tI-X_{t,(a,b,c)})$

and

(27) $\displaystyle
X_{t,(a+1,b+1,c+1)}=R^{t,(a+1,b+1,c+1)}_{(a,b,c)}X_{t,(a,b,c)}(R^{t,(a+1,b+1,c+1)}_{(a,b,c)})^{-1}.$

Thanks to Theorem 3.4 the triple matrix product on the right-hand side of Eq.
23 has a bandwidth of two. Hence, the entries of the differentiation matrix
may be computed in linear complexity as one knows the upper bandwidth and thus
we are not required to compute every entry of each column. A similar
computation is conducted for (partially) weighted quasimatrices except the
parameter associated with the weight is incremented down instead of up.

If one does not have the differentiation and Jacobi matrix of ${\bf
Q}^{t,(a,b,c)}(x)$, then one recursively decrements the $c$ parameter until
they hit a case where the the matrices are known. In the worst case, they
reach the case where $c=0$ which is a scaling of the classical Jacobi
polynomial case. In particular

(28) $\displaystyle\frac{\mathrm{d}}{{\rm d}x}{\bf Q}^{t,(a,b,0)}(x)={\bf
Q}^{t,(a+1,b+1,0)}(x)D^{t,(a+1,b+1,0)}_{(a,b,0)}$

where $D^{t,(a+1,b+1,0)}_{(a,b,0)}$ has an explicit expression and has a
bandwidth of one [26, Eq. 18.9.15]. Thus the worst case complexity to compute
the differentiation matrix for an arbitrary ${\bf Q}^{t,(a,b,c)}(x)$
quasimatrix, $c\in\mathbb{N}$, with zero pre-computation is $\mathcal{O}(cN)$.

###  3.5 Computing the triple matrix products in $\mathcal{O}(N)$ complexity

In this section we briefly discuss how to compute the triple matrix products
in (16), (17) and (18) in $\mathcal{O}(N)$ complexity and provide a comparison
of the Cholesky method and the two methods based on QR decomposition. For the
$R$-based methods of Cholesky or QR, the object of concern is naturally the
right-application of the inverse $R^{-1}$ which, in general, will not be
banded. For the $Q$-based method via QR, the inverse computation is trivial
via $Q^{\top}$, but in exchange neither $Q$ nor $Q^{\top}$ are banded.

As usual in computational contexts, the main idea for optimal complexity for
the $R$ methods is to never explicitly construct the inverse matrix $R^{-1}$.
Instead one notes that the bandwidth of the matrix product $RX$ to which we
wish to right-apply $R^{-1}$ is known from the bandwidths of $R$ and $X$,
where $X$ is a Jacobi matrix. Since the resulting matrix $RXR^{-1}$ is a
Jacobi matrix and must thus be tridiagonal, one only needs to compute three
entries per column with $N$-independent complexity, i.e. in $\mathcal{O}(1)$.
Computing the entries required for an $N\times N$ principal subblock of the
new Jacobi matrix thus has complexity $\mathcal{O}(N)$.

Obtaining optimal complexity for the computation of $Q^{\top}XQ$ in the $Q$
method is slightly more work, in particular, since $Q^{\top}$ generally has
infinitely many subdiagonals. Applying $Q^{\top}$ naïvely to the tridiagonal
$X$ would thus require $O(N^{2})$ operations. To avoid this we require access
to the Householder matrices that comprise the orthogonal matrix $Q$.
Fortunately, competitive QR factorisations of symmetric tridiagonal (and
banded) matrices are already stored in such a way that $Q$ is not directly
constructed and instead the information to compute the individual Householder
matrices is retained. Using this information, the action of each Householder
triple matrix product $H_{j}^{\top}XH_{j}$ only affects a small block
(determined by the size of the _lower_ bandwidth of the $Q$ factor) on the
band and elements in the top left of $H_{j}^{\top}XH_{j}$ are progressively
finalised and no longer modified by applications of the remaining Householder
matrices. This means that the algorithm only has to update a small block of
$X$ along the band starting in the top left by multiplying with the active
subblock of the Householder matrix on both sides. The size of the blocks to be
updated does not change with $N$ if we further use the knowledge that the
resulting matrix, as a Jacobi matrix of univariate orthogonal polynomials,
must be tridiagonal and thus the modification to $X$ from each triple
Householder product $H_{j}^{\top}XH_{j}$ can be computed in $\mathcal{O}(1)$
complexity.

######  Remark 3.5.

The $Q$ and $R$ variants correspond to the two similarity transformations
implied by the QR algorithm with a constant shift for symmetric tridiagonal
matrices. Similarly, the $R$ method from Cholesky corresponds to the LR
algorithm with the same constant shift.

In Fig. 4a we show a plot demonstrating the $\mathcal{O}(N)$ computational
complexity of the above-described algorithms and compare the performance of
the three available approaches: Cholesky as in (16), $Q$-based QR as in (17)
and $R$-based QR as in (18). In Fig. 4b we show the relative 2-norm error
comparing single and double precision computations of computed Jacobi matrix
principal subblocks. The stability of all three methods is self-evident. In
our experiments, we have not observed more than a constant factor between all
three approaches; certain parameter values lead to one method outperforming
the others. By analogy to the QR and LR eigenvalue algorithms, our default
preference is the $Q$-method.

(a) Computational cost

(b) 2-norm relative error (single precision)

Figure 4: (a) shows the $\mathcal{O}(N)$ complexity of the Cholesky as well as
$Q$ and $R$ based QR methods for computing the $N\times N$ principal Jacobi
matrix subblock of the semiclassical Jacobi polynomials
$Q_{n}^{t,(1,1,22)}(x)$ from $Q_{n}^{t,(1,1,20)}(x)$. The Cholesky method
requires twice the amount of work in this example as it proceeds in two steps.
CPU timings were obtained on a Lenovo ThinkPad X1 Carbon laptop with a 12th
Gen Intel i5-1250P CPU. (b) shows the relative 2-norm error for $N\times N$
principal subblocks of the Jacobi matrices comparing single to double
precision computations.

##  4 Orthogonal polynomials on the annulus

We now have the ingredients to define the orthogonal polynomials on the
annulus as well as compute the connections between various families. Consider
an inner radius $0<\rho<1$ and without loss of generality take the outer
radius to be one. Throughout this section
$R^{t,(\alpha,\beta,\gamma)}_{\boldsymbol{\cdot},(a,b,c)}$ and
$D^{t,(\alpha,\beta,\gamma)}_{\boldsymbol{\cdot},(a,b,c)}$ denote the lowering
and differentiation matrix from the (potentially weighted) semiclassical
Jacobi quasimatrix with parameters $(\alpha,\beta,\gamma)$ to $(a,b,c)$ as
they were defined in Section 3.

On the two-dimensional disk ($\rho=0$) we have the (complex-valued)
generalised Zernike polynomials, written in terms of Jacobi polynomials as

(29) $\displaystyle Z^{(b)}_{n,m}(x,y)\coloneqq r^{m}{\rm e}^{{\rm
i}m\theta}P_{(n-m)/2}^{(b,m)}(2r^{2}-1),$

where $r^{2}=x^{2}+y^{2}$, $0\leq m\leq n$ and

(30) $\displaystyle\begin{cases}m=0&\text{if}\;n=0\\\ m=1,3,\dots,n&\text{if
$n$ is odd},\\\ m=0,2,\ldots,n&\text{if $n$ is even}.\end{cases}$

These are polynomials in Cartesian coordinates $x$ and $y$ since $r^{m}{\rm
e}^{\pm{\rm i}m\theta}=(x\pm{\rm i}y)^{m}$. Moreover, they are orthogonal on
the unit disk with respect to the weight $(1-r^{2})^{b}$.

Note the real analogues are deduced by replacing the complex exponentials with
sines and cosines. It is convenient for generalising to other dimensions to
write these in terms of harmonic polynomials

(31) $\displaystyle Y_{m,0}(x,y)$ $\displaystyle={\rm Im}\,(x+{\rm
i}y)^{m}=r^{m}\sin(m\theta),$ (32) $\displaystyle Y_{m,1}(x,y)$
$\displaystyle={\rm Re}\,(x+{\rm i}y)^{m}=r^{m}\cos(m\theta).$

Thus the real-valued Zernike polynomials are:

(33) $\displaystyle
Z^{(b)}_{n,m,j}(x,y)=Y_{m,j}(x,y)P_{(n-m)/2}^{(b,m)}(2r^{2}-1).$

For $n$ odd, $m=1,3,\ldots,n$, and for $n$ even, $m=0,2,\ldots,n$. If $m=0$
then $j=1$, otherwise $j\in\\{0,1\\}$.

######  Definition 4.1 (Generalised Zernike annular polynomials).

For $t=(1-\rho^{2})^{-1}$, $r^{2}=x^{2}+y^{2}$, we define the _generalised
Zernike annular polynomials_ , $Z^{\rho,(a,b)}_{n,m,j}$, as

(34) $\displaystyle Z^{\rho,(a,b)}_{n,m,j}(x,y)\coloneqq
Y_{m,j}(x,y)Q_{(n-m)/2}^{t,(a,b,m)}\\!\left({1-r^{2}\over 1-\rho^{2}}\right),$

where for $n$ odd, $m=1,3,\ldots,n$, and for $n$ even, $m=0,2,\ldots,n$. If
$m=0$ then $j=1$, otherwise $j\in\\{0,1\\}$.

First note the functions as defined in Eq. 34 are all polynomials.

######  Proposition 4.2 (Orthogonality of generalised Zernike annular
polynomials, cf. [14, Eq. 69]).

Let $r^{2}=x^{2}+y^{2}$ and $\Omega_{\rho}=\\{(x,y):0<\rho\leq r\leq 1\\}$.
Then $Z^{\rho,(a,b)}_{n,m,j}$ are orthogonal polynomials with respect to

(35)
$\displaystyle\left\langle{f,g}\right\rangle_{\rho,(a,b)}=\iint_{\Omega_{\rho}}f(x,y)g(x,y)(1-r^{2})^{a}(r^{2}-\rho^{2})^{b}{\rm
d}x{\rm d}y.$

###  4.1 Connection and differentiation matrices

In this subsection, our goal is to find the matrices $D$ and $L$ such that

(36) $\displaystyle\Delta{\bf W}^{\rho,(1,1)}(x,y)$ $\displaystyle={\bf
Z}^{\rho,(1,1)}(x,y)D,$ (37) $\displaystyle{\bf W}^{\rho,(1,1)}(x,y)$
$\displaystyle={\bf Z}^{\rho,(1,1)}(x,y)R,$

where ${\bf W}^{\rho,(a,b)}(x,y)$ denotes the weighted version of ${\bf
Z}^{\rho,(a,b)}(x,y)$. We discover that $D$ and $R$ are block tridiagonal and
pentadiagonal, respectively, where each block is diagonal.

Let ${\bf Z}^{\rho,(a,b)}_{m,j}(x,y)$ denote the quasimatrix containing only
the “$m$, $j$”-mode 2D annulus orthogonal polynomials with increasing $n$. For
instance:

(38) $\displaystyle{\bf Z}^{\rho,(a,b)}_{m,j}=\left(Z^{\rho,(a,b)}_{m,m,j}\ |\
Z^{\rho,(a,b)}_{m+2,m,j}\ |\ Z^{\rho,(a,b)}_{m+4,m,j}\ |\ \cdots\right).$

${\bf W}^{\rho,(1,1)}_{m,j}$ is defined analogously.

######  Proposition 4.3 (Lowering of weighted generalised Zernike annular
polynomials).

(39) $\displaystyle{\bf W}^{\rho,(1,1)}_{m,j}(x,y)=(1-\rho^{2})^{2}{\bf
Z}^{\rho,(1,1)}_{m,j}(x,y)R^{t,(1,1,m)}_{(0,0,m)}R^{t,(0,0,m)}_{\mathrm{ab},(1,1,m)}.$

$R^{t,(1,1,m)}_{(0,0,m)}R^{t,(0,0,m)}_{\mathrm{ab},(1,1,m)}$ is a
pentadiagonal matrix. Again the different $m$ and $j$ modes do not interact.
When interlacing the different Fourier modes, we result in a block
pentadiagonal matrix where each block is diagonal.

######  Proof 4.4.

Let $z=x+{\rm i}y$.

(40) $\displaystyle\begin{split}{\bf
W}^{\rho,(1,1)}_{m,1}(x,y)&=(1-r^{2})(r^{2}-\rho^{2}){\rm Re}\,[z^{m}]{\bf
Q}^{t,(1,1,m)}(\tau)\\\ &=t^{-2}{\rm Re}\,[z^{m}]\tau(1-\tau){\bf
Q}^{t,(1,1,m)}(\tau)\\\ &=t^{-2}{\rm Re}\,[z^{m}]{\bf
Q}^{t,(0,0,m)}(\tau)R^{t,(0,0,m)}_{\mathrm{ab},(1,1,m)}\\\ &=t^{-2}{\rm
Re}\,[z^{m}]{\bf
Q}^{t,(1,1,m)}(\tau)R^{t,(1,1,m)}_{(0,0,m)}R^{t,(0,0,m)}_{\mathrm{ab},(1,1,m)}\\\
&=(1-\rho^{2})^{2}{\bf
Z}^{\rho,(1,1)}_{m,1}(x,y)R^{t,(1,1,m)}_{(0,0,m)}R^{t,(0,0,m)}_{\mathrm{ab},(1,1,m)}.\end{split}$

The proof is similar for $j=0$.

In order to compute the Laplacian of the generalised Zernike annular
polynomials, we follow [37] and use the decomposition

$\Delta=4\partial\bar{\partial}=4\bar{\partial}\partial$

for the two linear operators222We have used slightly different operators from
[37], which used $D_{\pm}:={\partial\over\partial r}\pm{{\rm i}\over
r}{\partial\over\partial\theta}$. This is to ensure we mapped from polynomials
to polynomials.

$\displaystyle 2\partial=\partial_{x}-{\rm i}\partial_{y}={\rm e}^{-{\rm
i}\theta}\left({\partial_{r}-{{\rm i}\over
r}\partial_{\theta}}\right)\;\;\text{and}\;\;2\bar{\partial}=\partial_{x}+{\rm
i}\partial_{y}={\rm e}^{-{\rm i}\theta}\left({\partial_{r}+{{\rm i}\over
r}\partial_{\theta}}\right).$

######  Lemma 4.5.

For $\tau={1-r^{2}\over 1-\rho^{2}}$, $t=(1-\rho^{2})^{-1}$, $z=r{\rm e}^{{\rm
i}\theta}$, and a differentiable function $f$, we have

$\displaystyle\partial z^{m}f(\tau)\;$
$\displaystyle=z^{m-1}(mf(\tau)-tr^{2}f^{\prime}(\tau))=z^{m-1}(\tau-t)^{1-m}{\mathrm{d}\over\mathrm{d}\tau}[(\tau-t)^{m}f(\tau)],$
$\displaystyle\partial\bar{z}^{m}f(\tau)\;$
$\displaystyle=-t\bar{z}^{m+1}f^{\prime}(\tau),$
$\displaystyle\bar{\partial}z^{m}f(\tau)\;$
$\displaystyle=-tz^{m+1}f^{\prime}(\tau),$
$\displaystyle\bar{\partial}\bar{z}^{m}f(\tau)\;$
$\displaystyle=\bar{z}^{m-1}(mf(r^{2})-tr^{2}f^{\prime}(r^{2}))=\bar{z}^{m-1}(\tau-t)^{1-m}{\mathrm{d}\over\mathrm{d}\tau}[(\tau-t)^{m}f(\tau)].$

######  Proposition 4.6 (Laplacian of generalised Zernike annular
polynomials).

Let $t=(1-\rho^{2})^{-1}$, then

(41) $\displaystyle\Delta{\bf
Z}^{\rho,(a,b)}_{m,j}(x,y)=4(1-\rho^{2})^{-1}{\bf
Z}^{\rho,(a+2,b+2)}_{m,j}(x,y)D^{t,(a+2,b+2,m)}_{\mathrm{c},(a+1,b+1,m+1)}D^{t,(a+1,b+1,m+1)}_{(a,b,m)}.$

######  Proof 4.7.

Consider $j=1$ and $z=x+{\rm i}y$. Note that $\Delta=4\partial\bar{\partial}$
is a real operator and, therefore, commutes with the real operator ${\rm
Re}\,$. Hence,

(42) $\displaystyle\begin{split}&\Delta{\bf Z}^{\rho,(a,b)}_{m,1}(x,y)\\\
&\indent=4\partial\bar{\partial}{\rm Re}\,[z^{m}]{\bf Q}^{t,(a,b,m)}(\tau)\\\
&\indent=4{\rm Re}\,\partial\left(-tz^{m+1}{\bf
Q}^{t,(a+1,b+1,m+1)}(\tau)D^{t,(a+1,b+1,m+1)}_{(a,b,m)}\right)\\\
&\indent=-4t{\rm
Re}\,z^{m}(\tau-t)^{-m}{\mathrm{d}\over\mathrm{d}\tau}\left[(\tau-t)^{m+1}{\bf
Q}^{t,(a+1,b+1,m+1)}(\tau)\right]D^{t,(a+1,b+1,m+1)}_{(a,b,m)}\\\
&\indent=4t{\rm Re}\,[z^{m}]{\bf
Q}^{t,(a+2,b+2,m)}(\tau)D^{t,(a+2,b+2,m)}_{\mathrm{c},(a+1,b+1,m+1)}D^{t,(a+1,b+1,m+1)}_{(a,b,m)}\\\
&\indent=4t{\bf
Z}^{\rho,(a+2,b+2)}_{m,1}(x,y)D^{t,(a+2,b+2,m)}_{\mathrm{c},(a+1,b+1,m+1)}D^{t,(a+1,b+1,m+1)}_{(a,b,m)}.\end{split}$

The second and third lines in Eq. 42 hold thanks to Lemma 4.5. By replacing
${\rm Re}\,[z^{m}]$ with ${\rm Im}\,[z^{m}]$, the exact same steps hold for
$j=0$.

$D^{t,(a+2,b+2,m)}_{\mathrm{c},(a+1,b+1,m+1)}D^{t,(a+1,b+1,m+1)}_{(a,b,m)}$ is
an upper-triangular matrix with a bandwidth of three. We derive a similar
relationship for the weighted version $W^{\rho,(1,1)}_{n,m,j}(x,y)$.

######  Proposition 4.8 (Laplacian of weighted generalised Zernike annular
polynomials).

(43) $\displaystyle\Delta{\bf W}^{\rho,(1,1)}_{m,j}(x,y)=4(1-\rho^{2}){\bf
Z}^{\rho,(1,1)}_{m,j}(x,y)D^{t,(1,1,m)}_{\mathrm{c},(0,0,m+1)}D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}.$

######  Proof 4.9.

Consider $j=1$ and $z=x+{\rm i}y$. Then,

(44) $\displaystyle\begin{split}&\Delta{\bf W}^{\rho,(1,1)}_{m,1}(x,y)\\\
&\indent=4\partial\bar{\partial}\left[(1-r^{2})(r^{2}-\rho^{2}){\rm
Re}\,[z^{m}]{\bf Q}^{t,(1,1,m)}(\tau)\right]\\\ &\indent=4t^{-2}{\rm
Re}\,\partial\bar{\partial}\left[z^{m}\tau(1-\tau){\bf
Q}^{t,(1,1,m)}(\tau)\right]\\\ &\indent=-4t^{-1}{\rm
Re}\,\partial\left[z^{m+1}{\bf
Q}^{t,(0,0,m+1)}(\tau)D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}\right]\\\
&\indent=-4t^{-1}{\rm
Re}\,z^{m}(\tau-t)^{-m}{\mathrm{d}\over\mathrm{d}\tau}\left[(\tau-t)^{m+1}{\bf
Q}^{t,(0,0,m+1)}(\tau)\right]D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}\\\
&\indent=4t^{-1}{\rm Re}\,[z^{m}]{\bf
Q}^{t,(1,1,m)}(\tau)D^{t,(1,1,m)}_{\mathrm{c},(0,0,m+1)}D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}\\\
&\indent=4t^{-1}{\bf
Z}^{\rho,(1,1)}_{m,1}(x,y)D^{t,(1,1,m)}_{\mathrm{c},(0,0,m+1)}D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}.\end{split}$

The third and fourth lines in Eq. 44 hold thanks to Lemma 4.5. The case $j=0$
follows similarly.

The matrix
$D^{t,(1,1,m)}_{\mathrm{c},(0,0,m+1)}D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}$ is
tridiagonal. An important feature is the different $m$ and $j$ modes do not
interact. Thus we obtain tridiagonal Laplacians for each Fourier mode.

###  4.2 Analysis & synthesis

By computing the connection, Jacobi, and differentiation matrices for
hierarchies of semiclassical Jacobi families as described in Section 3, we may
compute the connection and differentiation matrices for the Zernike annular
polynomials in optimal complexity. In particular the connection and
differentiation matrices for all the Fourier modes up to the truncation mode
and degree $M=N$ in Propositions 4.3 and 4.6 are computed in
$\mathcal{O}(N^{2})$ complexity.

We now discuss the technique for $\mathcal{O}(N^{2}\log(N))$ complexity
analysis and synthesis operators. Suppose that the goal is to apply the
synthesis operator to the coefficient vector ${\bf f}$ to evaluate ${\bf
Z}^{\rho,(a,b)}(x,y){\bf f}$ on the grid
$(x,y)\in\\{(r_{k}\cos(l\theta),r_{k}\sin(l\theta))\\}$ where

(45) $\displaystyle r_{k}$
$\displaystyle=\sqrt{\cos^{2}\left(\frac{2k+1}{4K}\pi\right)+\rho^{2}\sin^{2}\left(\frac{2k+1}{4K}\pi\right)}\;\;\text{for}\;\;0\leq
k\leq K-1,$ (46) $\displaystyle\theta_{l}$
$\displaystyle=\frac{2l}{L}\pi\;\;\text{for}\;\;0\leq l\leq L-1.$

Here $K=\left\lfloor(N+1)/2\right\rfloor+1$, $L=4K-3$, and $N$ is the
truncation degree. Given the even truncation degree $N$, we first rearrange
the coefficients $f_{n,m,j}$ in ${\bf f}$ into the $N/2\times(N+1)$
rectangular matrix:

(47)
$\displaystyle\begin{pmatrix}f_{0,0,1}&f_{1,1,0}&f_{1,1,1}&f_{2,2,0}&f_{2,2,1}&\cdots&f_{N,N,0}&f_{N,N,1}\\\
f_{2,0,1}&f_{3,1,0}&f_{3,1,1}&f_{4,2,0}&f_{4,2,1}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
f_{N-2,0,1}&f_{N-1,1,0}&f_{N-1,1,1}&f_{N,2,0}&f_{N,2,1}&\cdots&0&0\\\
f_{N,0,1}&0&0&0&0&\cdots&0&0\end{pmatrix}.$

Each column corresponds to a separate Fourier mode $(m,j)$. The matrix of
coefficients Eq. 47 is converted to a matrix of Chebyshev–Fourier series
coefficients for which a fast transform exists via a DCT and FFT. The
conversion is best understood column-wise. Consider the last column and let
${\bf f}_{N,N,1}=(f_{N,N,1}\;0\;\cdots\;0)^{\top}\in\mathbb{R}^{N/2}$. Note
that

(48) $\displaystyle\begin{split}{\bf Z}^{\rho,(a,b)}_{N,1}(x,y){\bf
f}_{N,N,1}&=r^{N}\cos(N\theta){\bf Q}^{t,(a,b,N)}(\tau){\bf f}_{N,N,1}\\\
&=t^{-N/2}(t-\tau)^{N/2}\cos(N\theta){\bf Q}^{t,(a,b,N)}(\tau){\bf
f}_{N,N,1}.\end{split}$

Thus, as derived in Section 3.3, the final column in the matrix for the
equivalent Chebyshev–Fourier coefficients is given by
$t^{-N/2}R_{T}S^{-1}_{(a,b)}Q_{0}Q_{2}\cdots Q_{N-2}{\bf f}_{N,N,1}$.
Therefore, the final column is a product of $\mathcal{O}(N)$ $Q$-factors, each
of which costs $\mathcal{O}(N)$ to apply, culminating in $\mathcal{O}(N^{2})$
work. Naïvely, the complexity for computing the Chebyshev–Fourier coefficients
for all the columns is $\mathcal{O}(N^{3})$. However, this conversion has
considerable speedups if one utilises Givens rotations and a butterfly
application resulting in $\mathcal{O}(N^{2}\log N)$ complexity. We refer the
interested reader to [34] for a thorough description and numerical analysis of
these algorithms for spherical harmonic polynomials.

##  5 Sparse spectral method for PDEs

In this section we detail how one may use the Chebyshev–Fourier series and the
generalised Zernike annular polynomials to construct a _sparse_ spectral
method for the Helmholtz equation. Consider the Helmholtz equation, for some
$\lambda\in\mathbb{R}$ (potentially $\lambda\in L^{\infty}(\Omega_{\rho})$),

(49) $\displaystyle\begin{split}\Delta u(x,y)+\lambda
u(x,y)&=f(x,y)\;\;\text{in}\;\;\Omega_{\rho},\quad
u=0\;\;\text{on}\;\;\partial\Omega_{\rho}.\end{split}$

Boundary conditions. We consider two techniques to enforce boundary
conditions:

  1. 1.

Incorporate the boundary conditions into the basis. A homogeneous Dirichlet
boundary condition is automatically satisfied on the annulus if all the basis
polynomials contain the factor $(1-r^{2})(r^{2}-\rho^{2})$.

  2. 2.

Provided the equation decouples across Fourier modes, then one may add the
boundary conditions as two additional constraints. In general, these
additional constraints result in two dense rows.

An advantage of the second technique is that it is amendable to enforcing
inhomogeneous Dirichlet boundary conditions as well as Neumann and Robin
boundary conditions with only small modifications [28]. This results in an
almost banded matrix which can be solved optimally [28, Sec. 5]. However, it
may cause severe ill-conditioning if constructed naïvely. Tau methods are a
robust technique for enforcing boundary conditions in this way. Moreover if
they are combined with a preconditioner (via a Schur complement) then one may
recover the bandedness, cf. [18, 29] and [10, Sec. B]. In the first technique
the boundary conditions are automatically satisfied, however, incorporating
nontrivial boundary conditions can be difficult.

Chebyshev–Fourier series. Consider $x=r\cos(\theta)$ and $y=r\sin(\theta)$.
Then a change of variables for Eq. 49 gives:

(50) $\displaystyle\left(r^{2}\frac{\partial^{2}}{\partial
r^{2}}+r\frac{\partial}{\partial
r}+\frac{\partial^{2}}{\partial\theta^{2}}+\lambda
r^{2}\right)u(r,\theta)=r^{2}f(r,\theta).$

Let ${\bf C}^{(2)}$ denote the quasimatrix of the ultraspherical polynomials
[26, Sec. 18.3]. Recall the definition of $r_{\rho}$ from Eq. 4. Consider the
matrices $D$, $R$, and $X$ such that

$\displaystyle\left[r^{2}\frac{\partial^{2}}{\partial
r^{2}}+r\frac{\partial}{\partial r}\right]{\bf C}^{(2)}(r_{\rho})$
$\displaystyle={\bf T}(r_{\rho})D,\;\;{\bf C}^{(2)}(r_{\rho})R={\bf
T}(r_{\rho}),\;r{\bf C}^{(2)}(r_{\rho})={\bf C}^{(2)}(r_{\rho})X.$

It is known that $D$ and $R$ have a bandwidth of five [28] and $X$ is
tridiagonal. The algorithm for solving the Helmholtz equation on the annulus
via a scaled-and-shifted Chebyshev–Fourier series is provided in Algorithm 1.
Boundary conditions are enforced via the second technique.

Algorithm 1 Scaled-and-shifted Chebyshev–Fourier series

1:EXPAND. Expand the right-hand side $f(r,\theta)$ in the ${\bf
T}(r_{\rho})\otimes{\bf F}(\theta)$-basis. Truncate at polynomial degree $N$
and $2N-1$ Fourier modes for $(N+1)(2N+1)$ coefficients. The solve is then
decomposed into a parallelisable sequence of one-dimensional solves:

2:ASSEMBLE. (Lazily) assemble the matrices for $m\in\\{0,1,\dots,N\\}$:  (51)
$\displaystyle\begin{pmatrix}{\bf T}(-1)\\\ {\bf T}(1)\\\ D-m^{2}R+\lambda
X^{2}R\end{pmatrix}.$ and truncate to form the matrix
$A_{m}\in\mathbb{R}^{(N+2)\times(N+2)}$.

3:EXTRACT. Extract and truncate the coefficients from ${\bf f}$ associated
with the $(m,j)$-mode to form the coefficient vector ${\bf
f}_{m,j}\in\mathbb{R}^{N}$.

4:SOLVE. Solve the $2N+1$ linear systems of size $(N+2)\times(N+2)$,
$A_{m}{\bf u}_{m,j}=\left(0,\;0,\;(X^{2}R{\bf f}_{m,j})^{\top}\right)^{\top}$,
where $m\in\\{0,1,\dots,N\\}$, $j\in\\{0,1\\}$, $(m,j)\neq(0,0)$.

5:INTERLACE. Interlace ${\bf u}_{m,j}\in\mathbb{R}^{N+2}$ to form the
coefficient vector ${\bf u}$.

6:$u(r,\theta)\approx\left[{\bf T}(r_{\rho})\otimes{\bf F}(\theta)\right]{\bf
u}$.

The expansion of the right-hand side $f(r,\theta)$ in step 1 of Algorithm 1
may be achieved via a DCT and FFT for the Chebyshev and Fourier parts,
respectively, resulting in a complexity of $\mathcal{O}(N^{2}\log(N))$. In
step 4, the first two rows in the right-hand side vector enforce the
homogeneous Dirichlet boundary condition. The coefficient vector ${\bf
f}_{m,j}$ is an expansion of $f(r,\theta)$ in the ${\bf
T}(r_{\rho})\otimes{\bf F}(\theta)$ basis and must be mapped to an expansion
of $r^{2}f(r,\theta)$ in the ${\bf C}^{(2)}(r_{\rho})\otimes{\bf F}(\theta)$
basis. This is precisely $X^{2}R{\bf f}_{m,j}$.

Generalised Zernike annular polynomials. An algorithm for solving Eq. 49 with
a generalised Zernike annular polynomial discretisation is provided in
Algorithm 2. The Dirichlet boundary condition is directly incorporated into
the basis by expanding the solution in weighted generalised Zernike annular
polynomials.

Algorithm 2 Generalised Zernike annular polynomials

1:EXPAND. Expand $f$ in the ${\bf Z}^{\rho,(1,1)}$-basis and obtain the
coefficient vector $\mathbf{f}$. The solve is then decomposed into a
parallelisable sequence of one-dimensional solves:

2:ASSEMBLE. Per Propositions 4.3 and 4.8, (lazily) assemble the matrices
$4t^{-1}D^{t,(1,1,m)}_{\mathrm{c},(0,0,m+1)}D^{t,(0,0,m+1)}_{\mathrm{ab},(1,1,m)}+\lambda
t^{-2}R^{t,(1,1,m)}_{(0,0,m)}R^{t,(0,0,m)}_{\mathrm{ab},(1,1,m)}$ and truncate
at degree $N$ to form the matrix $A_{m}\in\mathbb{R}^{k\times k}$ where
$k=\left\lceil\frac{N+1-m}{2}\right\rceil$, $m\in\\{0,1,\dots,N\\}$.

3:EXTRACT. Extract the coefficients of ${\bf f}$ for each Fourier mode and
truncate to find the coefficient vectors ${\bf f}_{m,j}\in\mathbb{R}^{k}$,
$m\in\\{0,1,\dots,N\\}$, $j\in\\{0,1\\}$, $(m,j)\neq(0,0)$.

4:SOLVE. Solve (in parallel) the $2N+1$ sparse linear systems of size $k\times
k$, $A_{m}\mathbf{u}_{m,j}=\mathbf{f}_{m,j}$.

5:INTERLACE. Interlace the coefficient vectors $\mathbf{u}_{m,j}$ to construct
$\mathbf{u}$.

6:$u(x,y)\approx{\bf W}^{\rho,(1,1)}(x,y)\mathbf{u}$.

######  Remark 5.1.

We note that in Algorithms 1 and 2, one solves $2N+1$ linear systems where $N$
is the truncation degree. With Chebyshev–Fourier series, the size of all the
$2N+1$ linear systems is $(N+2)\times(N+2)$ with two dense rows and a leading
nonadiagonal band. Whereas the Zernike annular polynomial induced linear
systems are tridiagonal with decreasing size with increasing $m$, starting at
the size $\left\lceil(N+1)/2\right\rceil\times\left\lceil(N+1)/2\right\rceil$
and ending at the size $1\times 1$ as also remarked in Fig. 3.

##  6 Spectral element method

In this section we construct a spectral element method for the strong
formulation of Eq. 49. The disk mesh is constructed such that the first cell
is the innermost disk, and the subsequent cells are concentric annuli stacked
around the disk as visualised in Fig. 5. A mesh for the annulus omits the
innermost disk cell.

Figure 5: The first three cells of a mesh for the spectral element method on a
disk. The disk and annuli cells may vary in thickness.

For ease of notation, we consider a two-element method although the
methodology extends to an arbitrary number of elements. In the disk and
annulus cell, we utilise the Zernike and Zernike annular polynomials,
respectively. Consider a mesh with an inner disk cell of radius $\rho$ and the
surrounding annulus with inner radius $\rho$ and outer radius 1. Let
$\tilde{\bf Z}^{(0)}(r,\theta)\coloneqq{\bf Z}^{(0)}(r/\rho,\theta)$. The goal
is to approximate the solution $u(x,y)$ of the Helmholtz equation Eq. 49 with
the expansion

(52) $\displaystyle u(x,y)=\begin{pmatrix}\tilde{\bf Z}^{(0)}(x,y)&{\bf
Z}^{\rho,(0,0)}(x,y)\end{pmatrix}\begin{pmatrix}\mathbf{u}_{0}\\\
\mathbf{u}_{1}\end{pmatrix}.$

For any Fourier mode $(m,j)$ we have the lowering relationships

$\displaystyle\tilde{\bf Z}^{(0)}_{m,j}$ $\displaystyle=\tilde{\bf
Z}^{(2)}_{m,j}R^{(2)}_{m,(0)},\;\;\Delta\tilde{\bf Z}^{(0)}_{m,j}=\tilde{\bf
Z}^{(2)}_{m,j}\tilde{D}^{(2)}_{m,(0)}$ $\displaystyle{\bf
Z}^{\rho,(0,0)}_{m,j}$ $\displaystyle={\bf
Z}^{\rho,(2,2)}_{m,j}R^{\rho,(2,2)}_{m,(0,0)},\;\;\Delta{\bf
Z}^{\rho,(0,0)}_{m,j}={\bf Z}^{\rho,(2,2)}_{m,j}D^{\rho,(2,2)}_{m,(0,0)}.$

Note that $\tilde{D}^{(2)}_{m,(0)}$ has a bandwidth of one,
$R^{(2,2)}_{m,(0,0)}$ and $D^{\rho,(2,2)}_{m,(0,0)}$ have a bandwidth of
three, and $R^{\rho,(2,2)}_{m,(0,0)}$ has a bandwidth of five.

The Helmholtz operator with Dirichlet boundary conditions takes the following
form for each Fourier mode, $m\in\mathbb{N}_{0}$, $j\in\\{0,1\\}$,
$(m,j)\neq(0,0)$;

(53) $\displaystyle A_{m}=\begin{pmatrix}0&{\bf
Z}_{m,0}^{\rho,(0,0)}(1,0)&0\\\ \tilde{\bf Z}_{m,0}^{(0)}(\rho,0)&-{\bf
Z}_{m,0}^{\rho,(0,0)}(\rho,0)&1-\mathrm{e}^{-(m+1)}\\\
\frac{\mathrm{d}}{\mathrm{d}r}\tilde{\bf
Z}_{m,0}^{(0)}(\rho,0)&-\frac{\mathrm{d}}{\mathrm{d}r}{\bf
Z}_{m,0}^{\rho,(0,0)}(\rho,0)&0\\\ \tilde{D}^{(2)}_{m,(0)}+\lambda
R^{t,(2)}_{m,(0)}&0&E_{m}\\\ 0&D^{\rho,(2,2)}_{m,(0,0)}+\lambda
R^{\rho,(2,2)}_{m,(0,0)}&0\end{pmatrix}.$

The first row in $A_{m}$ enforces the boundary condition at $r=1$. As we are
discretising the equation in strong form, we are required to enforce
continuity of the expansion at $r=\rho$ between cells as well as the outward
normal of the derivative, $\partial_{n}=\frac{\mathrm{d}}{\mathrm{d}r}$, of
the expansion. This is achieved by the second and third rows, respectively.
The bottom left $2\times 2$ block matrix contains the recurrence relations
between the coefficients and has bandwidth three (one if $\lambda=0$) in the
top left block and bandwidth five (three if $\lambda=0$) in the bottom right
block. The final column is derived from a tau-method and ensures that, after a
square truncation, $A_{m}$ is well-conditioned.

Let $\tilde{A}_{m}$ denote the matrix $A_{m}$ without the last column. Given
an even truncation polynomial degree $N$, a discretisation truncates
$\tilde{A}_{m}$ to a rectangular matrix of size $(2k+3)\times(2k+2)$ where
$k=\left\lceil\frac{N+1-m}{2}\right\rceil$. The $(1,2)$, $(2,1)$, $(2,2)$,
$(3,1)$ and $(3,2)$ blocks are truncated at the size $1\times(k+1)$. The
$(4,1)$ and $(5,2)$ blocks are truncated at the size $k\times(k+1)$.
Unfortunately, this truncation leads to a rectangular overdetermined system.
This is a well-understood phenomenon when dealing with continuity and boundary
conditions in spectral discretisations of problems formulated in strong form.
Removing a row from either the $(4,1)$ or $(5,2)$ block leads to severe ill-
conditioning for increasing $N$ and $m$ and artificial numerical pollution in
the solutions. The final column in $A_{m}$ resolves this issue and may be
understood via the tau-method, c.f. [18, 29] and [10, Sec. B]. For each
Fourier mode, we augment the equation with the three tau-functions

(54) $\displaystyle\tau_{m,j,i}(x,y)$
$\displaystyle=c_{m,j,i}Z^{\rho,(0,0)}_{N_{m}+2i,m,j}(x,y),\;\;i\in\\{0,1,2\\}.$

where $N_{m}=N$ if $m$ is even, otherwise $N_{m}=N-1$. Thus for each Fourier
mode we wish to solve

(55) $\displaystyle\begin{split}&\Delta\tilde{\bf
Z}^{(0)}_{m,j}(x,y)\mathbf{u}_{0,m,j}+\Delta{\bf
Z}^{\rho,(0,0)}_{m,j}(x,y)\mathbf{u}_{1,m,j}+\sum_{i=0}^{2}c_{m,j,i}\tau_{m,j,i}(x,y)=f(x,y),\end{split}$

with the boundary condition, where $x^{2}+y^{2}=1$,

(56) $\displaystyle{\bf Z}^{\rho,(0,0)}_{m,j}(x,y)\mathbf{u}_{1,m,j}$
$\displaystyle=0,$

and the continuity conditions, where $x^{2}+y^{2}=\rho^{2}$,

(57) $\displaystyle\tilde{\bf Z}^{(0)}_{m,j}(x,y)\mathbf{u}_{0,m,j}-{\bf
Z}^{\rho,(0,0)}_{m,j}(x,y)\mathbf{u}_{1,m,j}+(1-\mathrm{e}^{-(m+1)})c_{m,j,0}\tau_{m,j,0}(x,y)$
$\displaystyle=0,$ (58) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}r}\tilde{\bf
Z}^{(0)}_{m,j}(x,y)\mathbf{u}_{0,m,j}-\frac{\mathrm{d}}{\mathrm{d}r}{\bf
Z}^{\rho,(0,0)}_{m,j}(x,y)\mathbf{u}_{1,m,j}$ $\displaystyle=0.$

The final two tau-functions, $i\in\\{1,2\\}$, are those chosen in a typical
ultraspherical method and result in two columns and associated rows with only
a one on the diagonal. Thus these may be immediately row eliminated. However,
$\tau_{m,j,0}(x,y)$ must be included as part of the linear system, resulting
in the (truncated) column vector
$E_{m}=(0,\cdots,\;\;0,\;\;1)^{\top}\in\mathbb{R}^{k}$. The inclusion of the
first tau-function in the continuity condition is atypical. Here it
significantly improves the conditioning of $A_{m}$ as $m\to\infty$. This is
because, at $x^{2}+y^{2}=\rho^{2}$,
$\tilde{Z}^{(0)}_{n,m,j}(x,y)\approx\mathcal{O}(1)$, whenever
$(n-m)\approx\mathcal{O}(1)$, whereas
${Z}^{\rho,(0,0)}_{n,m,j}(x,y)\approx\mathcal{O}(\rho^{m})$. Left unchecked
this introduces ill-conditioning for increasing $m$.

######  Remark 6.1.

A similar tau-method is used for the spectral element method where the basis
is the Chebyshev–Fourier series on the annuli cells. In this case the column
vector $E_{m}$ remains the same but the $1-\mathrm{e}^{-m}$ in the continuity
condition row is omitted.

Algorithm 3 Spectral element method

1:EXPAND. Restrict $f(x,y)$ to the disk and annulus cells and expand in the
bases $\tilde{\bf Z}^{(2)}(x,y)$ and ${\bf Z}^{\rho,(2,2)}(x,y)$ to find the
coefficient vectors $\mathbf{f}_{0}$ and $\mathbf{f}_{1}$, respectively. The
solve is then decomposed into a parallelisable sequence of one-dimensional
solves:

2:ASSEMBLE. Assemble and truncate $A_{m}$ as defined in Eq. 53 to form a
square matrix of size $(2k+3)\times(2k+3)$ where
$k=\left\lceil\frac{N+1-m}{2}\right\rceil$, $m\in\\{0,1,\dots,N\\}$, where $N$
is the truncation degree.

3:EXTRACT. Extract the coefficients from ${\bf f}_{0}$ and ${\bf f}_{1}$
associated with the $(m,j)$-mode to form the coefficient vector ${\bf
f}_{m,j}\in\mathbb{R}^{2k}$, $m\in\\{0,1,\dots,N\\}$, $j\in\\{0,1\\}$,
$(m,j)\neq(0,0)$.

4:SOLVE. Solve (in parallel) for $m\in\\{0,1,\dots,N\\}$, $j\in\\{0,1\\}$,
$(m,j)\neq(0,0)$,  (59) $\displaystyle A_{m}\begin{pmatrix}{\bf u}_{m,j}\\\
c_{m,j,0}\end{pmatrix}=\begin{pmatrix}0&0&0&{\bf
f}_{m,j}^{\top}\end{pmatrix}^{\top}.$

5:INTERLACE. Split and interlace ${\bf u}_{m,j}$ to form the coefficient
vectors ${\bf u}_{0}$ and ${\bf u}_{1}$.

##  7 Numerical examples

Code availability: For reproducibility, an implementation of the optimal
complexity algorithms for semiclassical Jacobi polynomials may be found in
ClassicalOrthogonalPolynomials.jl [3] and
SemiclassicalOrthogonalPolynomials.jl [5]. The Zernike annular polynomials are
implemented in the package AnnuliOrthogonalPolynomials.jl [1]. Their fast
analysis and synthesis operators are implemented in FastTransforms.jl [4].
Scripts to generate the plots and solutions of the numerical examples in this
section may be found in AnnuliPDEs.jl [2] which has been archived on Zenodo
[30].

In the examples we consider the problem of solving PDEs in the two-dimensional
annulus or disk, in particular the Poisson and Helmholtz equations with a
homogeneous Dirichlet boundary condition. We consider the sparse spectral
methods: Chebyshev–Fourier series (via Algorithm 1), Zernike annular
polynomials (via Algorithm 2) and spectral element methods (via Algorithm 3).
We emphasise that, for all the methods, the blocks decouple and one solves
over each Fourier mode separately, reducing the two-dimensional solve to
$2N+1$ one-dimensional solves, when truncating at the polynomial degree $N$.

Error measurement: In all examples we measure the $\ell^{\infty}$-norm error
of the approximation on each cell as evaluated on a heavily oversampled
generalised Zernike annular grid Eq. 46.

###  7.1 Forced Helmholtz equation

In this example we consider the forced Helmholtz equation Eq. 49 on the annuli
domains with inner radii $\rho=0.2$, $0.5$, and $0.8$, respectively, with
coefficient $\lambda=80^{2}$, and the data $f(x,y)=\sin(100x)$. Due to the
large positive choice of $\lambda$, the solution supports oscillations that
are, in general, difficult to resolve. Moreover, the right-hand side features
a high Fourier mode component. As the explicit solution is unavailable, we
measure the errors against over-resolved approximations with coefficients that
have long decayed to Float64 machine precision.

We provide plots of the solutions in Fig. 6 as well as the corresponding
convergence plots of the Chebyshev–Fourier series (via Algorithm 1) and
weighted Zernike annular discretisation (via Algorithm 2) in Fig. 7.

(a) $\rho=0.2$

(b) $\rho=0.5$

(c) $\rho=0.8$

Figure 6: Plots of the solutions of Eq. 49 with $\lambda=80^{2}$ and the
right-hand side $f(x,y)=\sin(100x)$ on the annuli domains with inradii
$\rho=0.2$, $0.5$, and $0.8$, respectively (Section 7.1).

(a) $\rho=0.2$ (b) $\rho=0.5$ (c) $\rho=0.8$

Figure 7: Convergence of the spectral methods for solving Eq. 49 with
$\lambda=80^{2}$ and the right-hand side $f(x,y)=\sin(100x)$ on the annuli
domains with inradii $\rho=0.2$, $0.5$, and $0.8$, respectively. The Zernike
annular polynomial approximation of the solution converges slower than the
Chebyshev–Fourier series analogue for $\rho=0.2$ but considerably faster when
$\rho=0.5$ and $0.8$.

The weighted Zernike annular discretisation converges much quicker than the
Chebyshev–Fourier series when $\rho=0.5$ and $0.8$ and considerably slower
when $\rho=0.2$. In fact, the convergence profile of the Chebyshev–Fourier
series is roughly independent of $\rho$, reaching Float64 machine precision
with 42,194 coefficients (truncation degree $N=140$) on all three domains. In
contrast, the convergence of the weighted Zernike annular method largely
depends on $\rho$. When the the inradius is large, the Zernike annular
discretisation requires fewer basis functions to resolve functions with large
Fourier mode components and fluctuating behaviour near the outer boundary. The
degradation in performance as $\rho\to 0$ is likely caused by the proximity of
the inner boundary to the logarithmic singularity of the solutions at the
origin (outside of the domain). The cause of this is due to the quadratic
variable transformation in the radial direction for the Zernike annular
discretisation, while the Chebyshev–Fourier series uses an affine variable
transformation. The former brings the logarithmic singularity at the origin
much closer in the Bernstein ellipse metric than the latter, cf. [9, Fig. 4].

We conclude by noting that, in practical examples, the Zernike annular
polynomials are likely to be used as part of a spectral element method (such
as in Sections 7.3 and 7.4). Thus typically $0\ll\rho<1$ as it is highly
unlikely one would utilise annuli cells with a very small inner radius
relative to the outer boundary.

###  7.2 Poisson equation with a modified Gaussian bump

Consider Eq. 49 with $\lambda=0$, $\rho=0.2$, and, for some coefficients
$a\in\mathbb{R}$, $b,c\in[\rho,1]$,

(60) $\displaystyle f(x,y)=-4a{\rm
e}^{-a((x-b)^{2}+(y-c)^{2})}(-a(b^{2}-2bx+c^{2}-2cy+x^{2}+y^{2})+1).$

The right-hand side $f$ is a smooth function and emulates a modified Gaussian
bump. For sufficiently positive $a\gg 0$, the solution is approximately
$u(x,y)\approx{\rm e}^{-a((x-b)^{2}+(y-c)^{2})}$. In Fig. 8 we plot the
solution as well as the convergence plot of the two spectral methods against
the number of basis functions utilised in the discretisation. We see that the
weighted Zernike annular polynomial discretisation converges the quickest.

Figure 8: Convergence of the spectral methods for solving Eq. 49 with
$\lambda=0$, $\rho=0.2$, and the right-hand side Eq. 60 with $a=250$, $b=0$,
and $c=0.6$ (left). Plot of the solution (right). The Zernike annular
polynomial approximation of the solution, with a smooth right-hand side,
converges faster than the Chebyshev–Fourier series analogue.

###  7.3 Discontinuous variable coefficients and data on a disk

Consider a variable-coefficient Helmholtz equation Eq. 49 such that the
coefficient $\lambda(r)$ and the right-hand side have discontinuities in the
radial direction. In this example we show that, although a traditional
spectral method struggles to resolve the jumps, the spectral element method of
Section 6 performs particularly well. Let the domain be the unit disk
$\Omega_{0}$. Define the jump function $\kappa(r)$,
$\kappa_{0},\kappa_{1}\in\mathbb{R}$,

$\kappa(r)\coloneqq\begin{cases}\kappa_{0}&\text{if}\;0\leq r<\rho,\\\
\kappa_{1}&\text{if}\;\rho\leq r\leq 1.\end{cases}$

Consider the continuous solution:

(61) $\displaystyle\begin{split}&u(x,y)=\left(\sum_{i=1}^{5}d_{i}{\rm
e}^{-a_{i}((x-b_{i})^{2}+(y-c_{i})^{2})}\right)\\\
&\indent\times\begin{cases}\kappa_{0}r^{2}/4+(\kappa_{1}-\kappa_{0})\rho^{2}/4-\kappa_{1}/4+(\kappa_{0}-\kappa_{1})\rho^{2}\log(\rho)/2&\text{if}\;0\leq
r\leq\rho,\\\
\kappa_{1}r^{2}/4-\kappa_{1}/4+(\kappa_{0}-\kappa_{1})\rho^{2}\log(r)/2&\text{if}\;\rho<r\leq
1.\end{cases}\end{split}$

Throughout this example we fix $d_{i}=1$ and $a_{i}=-10i$ for all $i$.
Moreover, for $i\in\\{1,2,3,4\\}$,
$(b_{i},c_{i})=(\rho\cos(\theta_{i}),\rho\sin(\theta_{i}))$ where
$\theta_{i}=0,\pi/2,\pi/3,5\pi/4$, respectively, and
$(b_{5},c_{5})=(0.9\cos(3\pi/4),0.9\sin(3\pi/4))$. Furthermore,
$\kappa(r)=\kappa_{0}=100$ if $0\leq r\leq\rho$ and $\kappa(r)=\kappa_{1}=1$
if $\rho<r\leq 1$ with $\rho=1/2$. We consider two equations, with the
homogeneous Dirichlet boundary condition $u|_{\partial\Omega_{0}}=0$,

(62) $\displaystyle\Delta u(x,y)$ $\displaystyle=f_{1}(x,y),$ (63)
$\displaystyle(\Delta+\kappa(r))u(x,y)$ $\displaystyle=f_{2}(x,y).$

The first equation Eq. 62 is the Poisson equation posed on the unit disk with
the right-hand side

(64) $\displaystyle\begin{split}f_{1}(x,y)&=\sum_{i=1}^{5}d_{i}{\rm
e}^{-a_{i}((x-b_{i})^{2}+(y-c_{i})^{2})}\\\
&\indent\times\left(-4a_{i}(-a_{i}(b_{i}^{2}-2b_{i}x+c_{i}^{2}-2c_{i}y+x^{2}+y^{2})+1)+\kappa(r)\right).\end{split}$

Note that $f_{1}(x,y)$ has discontinuities in the radial direction at $r=1/2$.
The second equation Eq. 63 is a variable-coefficient Helmholtz equation where
the Helmholtz coefficient is $\lambda(r)=\kappa(r)$. Hence, the coefficients
of the equation also contain jumps in the radial direction. The corresponding
right-hand side is

(65) $\displaystyle f_{2}(x,y)=f_{1}(x,y)+\kappa(r)u(x,y).$

In Fig. 9 we plot the conditioning and size of the Laplacian matrices for
increasing Fourier mode $m$ with a truncation degree $N=100$. We observe that
the two-element method of Zernike and Zernike annular polynomials has, on
average, better conditioning and smaller matrix size. The discontinuous right-
hand sides $f_{1}$ and $f_{2}$, together with the continuous solution Eq. 61,
are those provided in Section 1 in Fig. 1. In Fig. 10 we plot the convergence
of various spectral methods. For the Poisson equation Eq. 62 we compare the
spectral methods: (1) a one-element weighted Zernike spectral method
(Algorithm 2), (2) a two-element method Zernike (annular) method (Algorithm 3)
and (3) a two-element method via a Zernike discretisation on the disk cell and
a Chebyshev–Fourier series on the outer annulus cell. In the two-element
methods, the inner cell is the disk $\\{0\leq r\leq 1/2\\}$ and the outer cell
is the annulus $\\{1/2\leq r\leq 1\\}$. For the variable-coefficient Helmholtz
equation we only plot the convergence of (2) and (3). A one-element method
would require resolving the function $\kappa(r)$ over the whole disk,
requiring high-degree polynomials, and would result in systems with large
bandwidths leading to impractical methods. Whereas in the spectral element
method, $\kappa(r)$ is a constant in each cell and the problem effectively
reduces to a non-variable Helmholtz problem on each cell. We avoid any
strategy that utilises a Chebyshev–Fourier series discretisation in a disk
cell due to well-studied issues at the origin [9, 37, 39].

The two-element Zernike (annular) method of Algorithm 3 performs the best.
This is followed by the two-element Zernike–Chebyshev–Fourier series method.
The one-element method struggles to reduce the accuracy below
$\mathcal{O}(10^{-4})$.

(a) Condition number. (b) Matrix size.

Figure 9: The condition number and matrix size of the Laplacian matrix, for
increasing Fourier mode $m$, of the two-element methods with a Zernike
discretisation on the disk cell and either a Chebyshev–Fourier series or a
Zernike annular discretisation on the annulus cell (inner radius $\rho=1/2$)
on the domain $\Omega_{0}$ with truncation degree $N=100$. We observe that the
Laplacian matrices for the two-element Zernike (annular) polynomials is better
conditioned, on average, and smaller in size.

(a) Convergence of Eq. 62. (b) Convergence of Eq. 63.

Figure 10: Convergence of the spectral methods for solving Eq. 62 and Eq. 63.
For both right-hand sides, the two-element methods with the boundary of the
cells at the radial discontinuity performs significantly better than the one-
element analogue. Moreover, the Zernike (annular) two-element methods performs
the best overall.

###  7.4 Discontinuous variable coefficients and data on an annulus

We study a similar problem to the one described in Section 7.3 but on the
annulus domain $\Omega_{1/10}$. We consider the Helmholtz equation Eq. 49 with
the right-hand side Eq. 65. The parameters are $\rho=1/10$, $\alpha=1/2$, and
for $i\in\\{1,2,3,4\\}$, $d_{i}=1$, $a_{i}=-10i$,
$(b_{i},c_{i})=(\alpha\cos(\theta_{i}),\alpha\sin(\theta_{i}))$ where
$\theta_{i}=0,\pi/2,\pi/3,5\pi/4$, respectively. For $i=5$ we pick
$(d_{5},a_{5},b_{5},c_{5})=(10,80,-0.95,0)$. Finally $\lambda(r)=\kappa(r)$
where $\kappa(r)=\kappa_{0}=100$ if $\rho\leq r\leq\alpha$ and
$\kappa(r)=\kappa_{1}=1$ if $\alpha<r\leq 1$. The right-hand side has a radial
discontinuity at $r=1/2$.

We compare three spectral methods: (1) a two-element method Zernike annular
method, (2) a two-element method via a Chebyshev–Fourier series on the inner
annulus cell and a Zernike annular discretisation on the outer annulus cell
and (3) a two-element Chebyshev–Fourier series. The inner cell is the annulus
$\\{1/10\leq r\leq 1/2\\}$ and the outer cell is the annulus $\\{1/2\leq r\leq
1\\}$. The discontinuous right-hand, the continuous solution and the
convergence plots are depicted in Fig. 11. The quickest convergence is
achieved by (1) followed by (2) and then (3).

(a) Right-hand side.

(b) Solution.

(c) Convergence.

Figure 11: Convergence of the spectral methods for solving the example in
Section 7.4. The domain is the annulus with inner radius $\rho=1/10$. The
right-hand side has a radial discontinuity at $r=1/2$. The mesh consists of
two annuli cells which meet at $r=1/2$. The Zernike annular two-element method
performs the best overall, followed by the Chebyshev–Fourier/Zernike annular
method and lastly the Chebyshev–Fourier discretisation on both cells.

##  8 Conclusions

In this work we detailed an optimal complexity algorithm for computing
connection and differentiation matrices as well a quasi-optimal complexity
algorithm for the analysis and synthesis operators of hierarchies of
semiclassical Jacobi polynomials. These allowed us to develop similar optimal
complexity computations for the generalised Zernike annular polynomials. With
the generalised polynomials, we construct several sparse spectral methods for
solving PDEs posed on the annulus and disk. In particular we focused on the
scaled-and-shifted Chebyshev–Fourier series and the Zernike annular
polynomials. Akin to similar observations by Boyd and Yu [9] we observed that,
for problems that feature a high Fourier mode, the Zernike annular polynomials
often converge faster. A key note is that the Helmholtz operator discretised
with generalised Zernike annular polynomial results in matrices of bandwidth
five as opposed to the bandwidth of nine given by the scaled-and-shifted
Chebyshev–Fourier series. We also constructed a spectral element method for
the disk and annulus where the cells are an inner disk (omitted if the domain
is an annulus) and subsequent concentric annuli of varying thickness. We used
this spectral element method to solve the Helmholtz equation with a right-hand
side and a variable coefficient with discontinuities in the radial direction.
The spectral element method converged quickly to machine precision whereas one
cell counterparts did not reduce the error below $\mathcal{O}(10^{-4})$ for
the truncation degrees considered. Moreover, the Laplacian has better
conditioning when a Zernike annular basis is used.

Unlike the Chebyshev–Fourier series, generalised Zernike annular polynomials
may be used to discretise the Helmholtz equation in weak form resulting in a
banded and sparse discretisation that preserves symmetry. Constructing a
sparse spectral element method that utilises this approach will the subject of
future work. We conclude by noting that the definition of the Zernike annular
polynomials as given in this work naturally extends to three dimensions. This
is achieved by simply considering the three-dimensional spherical harmonics in
Definition 4.1. Hence, it is possible to construct sparse spectral methods for
spherical shells in three dimensions.

## Acknowledgments

We are grateful to Tom H. Koornwinder for providing us with previous
literature on Zernike annular polynomials.

## References

  * [1] AnnuliOrthogonalPolynomials.jl, 2023, https://github.com/JuliaApproximation/AnnuliOrthogonalPolynomials.jl.
  * [2] AnnuliPDEs.jl, 2023, https://github.com/ioannisPApapadopoulos/AnnuliPDEs.jl.
  * [3] ClassicalOrthogonalPolynomials.jl, 2023, https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl.
  * [4] FastTransforms.jl, 2023, https://github.com/JuliaApproximation/FastTransforms.jl.
  * [5] SemiclassicalOrthogonalPolynomials.jl, 2023, https://github.com/JuliaApproximation/SemiclassicalOrthogonalPolynomials.jl.
  * [6] K. Atkinson, D. Chien, and O. Hansen, Spectral methods using multivariate polynomials on the unit ball, CRC Press, 2019, https://doi.org/10.1201/9780429344374.
  * [7] R. Barakat, Optimum balanced wave-front aberrations for radially symmetric amplitude distributions: generalizations of Zernike polynomials, JOSA, 70 (1980), pp. 739–742, https://doi.org/10.1364/JOSA.70.000739.
  * [8] N. Boullé and A. Townsend, Computing with functions in the ball, SIAM Journal on Scientific Computing, 42 (2020), pp. C169–C191, https://doi.org/10.1137/19M1297063.
  * [9] J. P. Boyd and F. Yu, Comparing seven spectral methods for interpolation and for solving the Poisson equation in a disk: Zernike polynomials, Logan–Shepp ridge polynomials, Chebyshev–Fourier series, cylindrical Robert functions, Bessel–Fourier expansions, square-to-disk conformal mapping and radial basis functions, Journal of Computational Physics, 230 (2011), pp. 1408–1438, https://doi.org/10.1016/j.jcp.2010.11.011.
  * [10] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Physical Review Research, 2 (2020), p. 023068, https://doi.org/10.1103/PhysRevResearch.2.023068.
  * [11] L. de Winter, T. Tudorovskiy, J. van Schoot, K. Troost, E. Stinstra, S. Hsu, T. Gruner, J. Mueller, R. Mack, B. Bilski, et al., High NA EUV scanner: obscuration and wavefront description, in Extreme Ultraviolet Lithography 2020, vol. 11517, SPIE, 2020, pp. 91–114, https://doi.org/10.1117/12.2572878.
  * [12] L. de Winter, T. Tudorovskiy, J. van Schoot, K. Troost, E. Stinstra, S. Hsu, T. Gruner, J. Mueller, R. Mack, B. Bilski, et al., Extreme ultraviolet scanner with high numerical aperture: obscuration and wavefront description, Journal of Micro/Nanopatterning, Materials, and Metrology, 21 (2022), pp. 023801–023801, https://doi.org/10.1117/1.JMM.21.2.023801.
  * [13] C. F. Dunkl and Y. Xu, Orthogonal polynomials of several variables, no. 155, Cambridge University Press, 2014.
  * [14] A. C. Ellison and K. Julien, Gyroscopic polynomials, Journal of Computational Physics, (2023), p. 112268, https://doi.org/10.1016/j.jcp.2023.112268.
  * [15] A. C. Ellison, K. Julien, and G. M. Vasil, A gyroscopic polynomial basis in the sphere, Journal of Computational Physics, 460 (2022), p. 111170, https://doi.org/10.1016/j.jcp.2022.111170.
  * [16] S. Filip, A. Javeed, and L. N. Trefethen, Smooth random functions, random ODEs, and Gaussian processes, SIAM Review, 61 (2019), pp. 185–205, https://doi.org/10.1137/17M1161853.
  * [17] T. S. Gutleb, S. Olver, and R. M. Slevinsky, Polynomial and rational measure modifications of orthogonal polynomials via infinite-dimensional banded matrix factorizations, arXiv preprint arXiv:2302.08448, (2023).
  * [18] C. Lanczos, Trigonometric interpolation of empirical and analytical functions, Journal of Mathematics and Physics, 17 (1938), pp. 123–199, https://doi.org/10.1002/sapm1938171123.
  * [19] D. Lecoanet, G. M. Vasil, K. J. Burns, B. P. Brown, and J. S. Oishi, Tensor calculus in spherical coordinates using Jacobi polynomials. Part-II: implementation and examples, Journal of Computational Physics: X, 3 (2019), p. 100012, https://doi.org/10.1016/j.jcpx.2019.100012.
  * [20] H. Li and Y. Xu, Spectral approximation on the unit ball, SIAM Journal on Numerical Analysis, 52 (2014), pp. 2647–2675, https://doi.org/10.1137/130940591.
  * [21] A. P. Magnus, Painlevé-type differential equations for the recurrence coefficients of semi-classical orthogonal polynomials, Journal of Computational and Applied Mathematics, 57 (1995), pp. 215–237, https://doi.org/10.1016/0377-0427(93)E0247-J.
  * [22] W. E. Maguire, R. Sejpal, and B. W. Smith, Defining Tatian–Zernike polynomials for use in a lithography simulator, in Optical and EUV Nanolithography XXXVI, vol. 12494, SPIE, 2023, pp. 222–232, https://doi.org/10.1117/12.2660861.
  * [23] V. N. Mahajan, Zernike annular polynomials for imaging systems with annular pupils, JOSA, 71 (1981), pp. 75–85, https://doi.org/10.1364/JOSA.71.000075.
  * [24] M. M. Meyer and F. R. P. Medina, Polar differentiation matrices for the Laplace equation in the disk under nonhomogeneous Dirichlet, Neumann and Robin boundary conditions and the biharmonic equation under nonhomogeneous Dirichlet conditions, Computers & Mathematics with Applications, 89 (2021), pp. 1–19, https://doi.org/10.1016/j.camwa.2021.02.005.
  * [25] M. Molina-Meyer and F. R. P. Medina, A collocation-spectral method to solve the bi-dimensional degenerate diffusive logistic equation with spatial heterogeneities in circular domains, Rendiconti dell’Istituto di matematica dell’Università di Trieste, 52 (2020), pp. 311–344, https://doi.org/10.13137/2464-8728/30917.
  * [26] F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, NIST Digital Library of Mathematical Functions.  http://dlmf.nist.gov/, Release 1.1.4 of 2022-01-15, 2022, http://dlmf.nist.gov/.
  * [27] S. Olver, R. M. Slevinsky, and A. Townsend, Fast algorithms using orthogonal polynomials, Acta Numerica, 29 (2020), pp. 573–699, https://doi.org/10.1017/S0962492920000045.
  * [28] S. Olver and A. Townsend, A fast and well-conditioned spectral method, SIAM Review, 55 (2013), pp. 462–489, https://doi.org/10.1137/120865458.
  * [29] E. L. Ortiz, The tau method, SIAM Journal on Numerical Analysis, 6 (1969), pp. 480–492.
  * [30] I. P. A. Papadopoulos, ioannispapapadopoulos/annulipdes.jl: v0.0.1, Oct. 2023, https://doi.org/10.5281/zenodo.8430270.
  * [31] J. P. Rolland, M. A. Davies, T. J. Suleski, C. Evans, A. Bauer, J. C. Lambropoulos, and K. Falaggis, Freeform optics for imaging, Optica, 8 (2021), pp. 161–176, https://doi.org/10.1364/OPTICA.413762.
  * [32] R. M. Slevinsky, FastTransforms, 2018, https://github.com/MikaelSlevinsky/FastTransforms.
  * [33] R. M. Slevinsky, On the use of Hahn’s asymptotic formula and stabilized recurrence for a fast, simple and stable Chebyshev–Jacobi transform, IMA Journal of Numerical Analysis, 38 (2018), pp. 102–124, https://doi.org/10.1093/imanum/drw070.
  * [34] R. M. Slevinsky, Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series, Applied and Computational Harmonic Analysis, 47 (2019), pp. 585–606, https://doi.org/10.1016/j.acha.2017.11.001.
  * [35] B. Tatian, Aberration balancing in rotationally symmetric lenses, JOSA, 64 (1974), pp. 1083–1091, https://doi.org/10.1364/JOSA.64.001083.
  * [36] A. Townsend, H. Wilber, and G. B. Wright, Computing with functions in spherical and polar geometries I. The sphere, SIAM Journal on Scientific Computing, 38 (2016), pp. C403–C425, https://doi.org/10.1137/15M1045855.
  * [37] G. M. Vasil, K. J. Burns, D. Lecoanet, S. Olver, B. P. Brown, and J. S. Oishi, Tensor calculus in polar coordinates using Jacobi polynomials, Journal of Computational Physics, 325 (2016), pp. 53–73, https://doi.org/10.1016/j.jcp.2016.08.013.
  * [38] G. M. Vasil, D. Lecoanet, K. J. Burns, J. S. Oishi, and B. P. Brown, Tensor calculus in spherical coordinates using Jacobi polynomials. Part-I: Mathematical analysis and derivations, Journal of Computational Physics: X, 3 (2019), p. 100013, https://doi.org/10.1016/j.jcpx.2019.100013.
  * [39] H. Wilber, A. Townsend, and G. B. Wright, Computing with functions in spherical and polar geometries II. The disk, SIAM Journal on Scientific Computing, 39 (2017), pp. C238–C262, https://doi.org/10.1137/16M1070207.