# Spectral solver for Cauchy problems in polar coordinates using discrete
Hankel transforms

Rundong Zhou111Now at Department of Physics, Chalmers University of
Technology, SE-412 96 Göteborg, Sweden<EMAIL_ADDRESS>Nicolas
Grisouard<EMAIL_ADDRESS>Division of Engineering Science,
University of Toronto, 40 St. George Street, Toronto, Canada, M5S 2E4
Department of Physics, University of Toronto, 60 St. George Street, Toronto,
Canada, M5S 1A7

###### Abstract

We introduce a Fourier-Bessel-based spectral solver for Cauchy problems
featuring Laplacians in polar coordinates under homogeneous Dirichlet boundary
conditions. We use FFTs in the azimuthal direction to isolate angular modes,
then perform discrete Hankel transform (DHT) on each mode along the radial
direction to obtain spectral coefficients. The two transforms are connected
via numerical and cardinal interpolations. We analyze the boundary-dependent
error bound of DHT; the worst case is $\sim N^{-3/2}$, which governs the
method, and the best $\sim e^{-N}$, which then the numerical interpolation
governs. The complexity is $O[N^{3}]$. Taking advantage of Bessel functions
being the eigenfunctions of the Laplacian operator, we solve linear equations
for all times. For non-linear equations, we use a time-splitting method to
integrate the solutions. We show examples and validate the method on the two-
dimensional wave equation, which is linear, and on two non-linear problems: a
time-dependent Poiseuille flow and the flow of a Bose-Einstein condensate on a
disk.

###### keywords:

Spectral methods, Discrete Hankel transforms , Initial value problems ,
Nonlinear partial differential equations , Bose-Einstein condensates

††journal: Journal of Computational Physics

A simple spectral solving algorithm for initial value problems featuring
Laplacian in polar coordinates.

Reformulating discrete Hankel transforms as pseudospectral collocation methods
and analyzing their error bounds.

Numerical results agree with analytical solutions with high accuracy, and
error bounds agree with the theory.

Solving the Gross-Pitaevskii equation for quantum vortices in Bose-Einstein
condensates.

##  1 Introduction

The Laplacian operator is associated with many important physical problems,
including diffusion, the Schrödinger equation or the wave equation. It is more
complicated under polar and spherical coordinates than in Cartesian
coordinates, and spectral approaches are often used to avoid coordinate
singularities [1, 2]. In this paper, we introduce a novel spectral solver
suited for time-dependent problems featuring Laplacian in polar coordinates.

We consider a function $\psi$ on a unit disk $r\in[0,1]$ that satisfies the
homogeneous Dirichlet boundary condition $\psi(1,\theta)=0$. The periodicity
in the azimuthal direction allows for a decomposition of $\psi$ as a Fourier
series, namely,

$\psi(r,\theta)=\sum_{q=-\infty}^{\infty}f_{q}(r)\text{e}^{\text{i}q\theta},\quad\text{with}\quad
f_{q}(r)=\frac{1}{2\pi}\int_{0}^{2\pi}\psi(r,\theta)\text{e}^{-\text{i}q\theta}\,\textrm{d}\theta$
(1)

and where $q$ denotes the angular mode number. There are various choices of
basis functions to decompose the radial function $f_{q}(r)$. Boyd and Yu [3]
list five spectral methods that use decompositions of the radial function.
Among the listed basis functions, four are families of orthogonal polynomials,
such as Zernike [4], Logan-Shepp [5], modified Chebyshev [6, 7, 1] or modified
Jacobi polynomials [8]. Most of them focus on solving boundary value problems
and time-independent partial differential equations (PDEs) in polar
coordinates, while Cauchy problems involving time evolution receive
comparatively little attention.

To build a Laplacian solver suitable for time-dependent problems, instead of
choosing orthogonal polynomial bases, we decided to use Bessel functions of
the first kind as our basis to decompose the radial function $f_{q}(r)$, given
by the Fourier-Bessel series [9],

$f_{q}(r)=\sum_{j=1}^{\infty}a_{q,j}J_{q}\left(k_{q,j}r\right),\quad\text{with}\quad
a_{q,j}=\frac{2}{J_{q+1}^{2}(k_{q,j})}\int_{0}^{1}rf_{q}(r)J_{q}(k_{q,j}r)\,\textrm{d}r.$
(2)

$J_{q}$ is the $q^{\text{th}}$-order Bessel function of the first kind and
$k_{q,j}$ denotes its $j^{\text{th}}$ non-negative zero. Temme [10] provides a
fast and accurate algorithm to compute these roots.

The Fourier-Bessel basis (i.e., the combination of the azimuthal Fourier basis
of Eqn. 1 and radial Bessel basis of Eqn. 2) has an algebraic rate of
convergence $\lesssim 1/N^{5/2}$ for functions satisfying homogeneous
Dirichlet boundary conditions [3], where $N$ is the number of basis functions
used for approximation. This compares unfavourably with many orthogonal
polynomial bases which have an exponential rate. However, Fourier-Bessel modes
are eigenfunctions of the Laplacian operator, and the boundary conditions are
enforced by the basis functions themselves. Such virtues ensure that a pure
spectral scheme can be applied at each iteration. It makes the Fourier-Bessel
basis a competitive choice for solving time-dependent initial value problems
associated with Laplacians under homogeneous Dirichlet conditions in polar
coordinates.

To decompose the radial function $f_{q}(r)$ into Fourier-Bessel series, we use
the discrete Hankel transform (DHT). It was first introduced mathematically by
Johnson [11], then independently re-invented by Yu et al. [12] for the zeroth
order mode ($q=0$) and Guizar-Sicairos et al. [13] for all integer orders.
Finally, It was categorized by Baddour [14] as a discrete variation of general
Fourier transform.

This paper is organized as follows: we first reformulate the DHTs as
pseudospectral collocation method by introducing their discrete inner
products, quadrature weights, pseudospectral grid points, and cardinal
interpolations in Section 2.1 and 2.2. We analyze error and complexity and
discuss the boundary dependency of the error bound in Section 2.3. Then we
introduce a systematic way to apply the method to compute Laplacian and solve
nonlinear time-dependent equations in Section 2.4. In Section 3 we show three
examples: the linear 2-D wave equation where we test the convergence of the
method numerically (§ 3.1), a time-dependent Poiseuille flow equation (§ 3.2),
and the Gross-Pitaevskii equation (§ 3.3). In Section 4 we conclude the paper
and discuss potential improvements.

##  2 Method

###  2.1 Discrete Hankel transform and pseudospectral grid points

Our first step is to reformulate the discrete Hankel transforms as
pseudospectral collocation methods. The optimal pseudospectral grid points are
the zeros of the $(N+1)^{\text{th}}$ basis function $\phi_{N+1}(r)$ [15, §
4.3], where $N$ is the total number of basis functions, or radial modes, used
to approximate the radial function $f_{q}(r)$. The basis functions of a
$q^{\text{th}}$-order Hankel transform are then

$\phi_{q,i}(r)=J_{q}\left(k_{q,i}r\right).$ (3)

Thus, the $(N+1)^{\text{th}}$ basis function is $J_{q}(k_{q,N+1}r)$, and we
choose its pseudospectral grid points to be

$\\{r_{q,i}\\}_{i=1,2,...,N}=\frac{k_{q,i}}{k_{q,N+1}}.$ (4)

We can now define the discrete inner product as [16, § 3.1.4]

$(f,\>g)_{q}\equiv\sum_{i=1}^{N}w_{q,i}f(r_{q,i})g(r_{q,i}),\quad\text{where}\quad
w_{q,i}=\frac{2}{k_{q,N+1}^{2}J_{q+1}^{2}(k_{q,i})}$ (5)

are the quadrature weights.

Bessel functions are orthogonal under our discrete inner product, namely [11]

$\frac{2}{J_{q+1}^{2}(k_{q,m})}\Big{(}J_{q}(k_{q,m}r),\>J_{q}(k_{q,n}r)\Big{)}_{q}=\delta_{mn}.$
(6)

The discrete Hankel transform of the radial function $f_{q}(r)$ is now defined
through the discrete inner product as

$F_{q,j}=\Big{(}f_{q}(r),J_{q}(k_{q,j}r)\Big{)}_{q}.$ (7)

By substituting the definition of discrete inner product and pseudospectral
grid points $\\{r_{q,i}\\}$, the forward discrete Hankel transform formula is
then given by

$F_{q,j}=\frac{2}{k_{q,N+1}^{2}}\sum_{i=1}^{N}\frac{J_{q}\left(k_{q,i}k_{q,j}/k_{q,N+1}\right)}{J_{q+1}^{2}\left(k_{q,i}\right)}f_{q}\left(r_{q,i}\right).$
(8)

The value of the radial function $f_{q}(r)$ at each pseudospectral grid point
$r_{q,i}$ can be recovered by the Fourier-Bessel series Eqn. (2), namely,

$f_{q}\left(r_{q,i}\right)=2\sum_{j=1}^{N}\frac{J_{q}\left(k_{q,i}k_{q,j}/k_{q,N+1}\right)}{J_{q+1}^{2}\left(k_{q,j}\right)}F_{q,j}.$
(9)

The transform formulae Eqns. (8) and (9) demonstrate symmetry and can be
implemented as matrix multiplications. Introducing

$f_{q,i}\equiv f_{q}\left(r_{q,i}\right)\quad\text{and}\quad
M_{q,ij}\equiv\phi_{q,i}(r_{q,j})w_{q,i}=\frac{2J_{q}\left(k_{q,i}k_{q,j}/k_{q,N+1}\right)}{k_{q,N+1}^{2}J_{q+1}^{2}\left(k_{q,i}\right)},$
(10)

we obtain

$\vec{F}_{q}=\textbf{M}_{q}\cdot\vec{f}_{q}\>;\quad\vec{f}_{q}=k_{N+1}^{2}\textbf{M}_{q}\cdot\vec{F}_{q},$
(11)

where $\textbf{M}_{q}$ is unitary, i.e.,
$\textbf{M}_{q}\textbf{M}_{q}^{-1}=\textbf{I}$. The discrete Hankel transform
can be categorized as a Matrix Multiplication Transformation [15, § 10.4], and
it demonstrates the spectral methods’ virtue of one-to-one correspondence
between spectral coefficients $F_{q,j}$ and pseudospectral grid points
$\\{f_{q}(r_{q,i})\\}$.

###  2.2 Interpolation and cardinal functions

To obtain the spectral coefficients $a_{q,j}$ of Eqn. (2), we perform a Fast
Fourier Transform (FFT) along the azimuthal direction followed in the radial
direction by one DHT for each angular mode $q$. The pseudospectral grid points
$\\{r_{q,i}\\}$ in Eqn. (4) defined for the DHTs are not evenly spaced and
vary with the angular mode $q$. Instead of using Baddour and Yao’s unevenly
sampled grid [17, 18], which causes unwanted artifacts in the center of the
domain, we choose an evenly sampled grid that is finer in the radial direction
before performing the azimuthal FFT (hereafter referred to as the FFT grid),
as visualized in Figure 1(a). We then numerically interpolate each radial
function $f_{q}$ on the corresponding DHT grid $\\{r_{q,i}\\}_{i=1,\dots,N}$
to obtain $\\{f_{q}(r_{q,i})\\}_{i=1,\dots,N}$, as shown in Figure 1(b).

(a) FFT grid

(b) DHT grid

Figure 1: Difference between the evenly sampled FFT grid on $\psi(r,\theta)$
and the pseudospectral $\\{r_{q,i}\\}$ grid for the DHTs in the case of $N=8$.
In panel (b), the sampling points for $q=0$ (marked as purple squares) do not
contain the origin since $r_{0,1}>0$.

To achieve optimal numerical accuracy for a given $N$, the number of radial
sampling points of the finer FFT grid should be at least $2N$ (see § 2.3.1)
The number of azimuthal sampling points is twice that of the desired angular
modes $q$. Conversely, when we transform from the spectral domain $a_{q,j}$
back to the physical domain $\psi(r,\theta)$, we perform inverse DHTs followed
by an inverse FFT. We use the cardinal functions of Fourier-Bessel series to
interpolate from the pseudospectral grid $\\{r_{q,i}\\}$ to the finer FFT grid
[14], namely,

$f_{q}(r)=\sum_{i=1}^{\infty}f_{q}(r_{q,i})C_{q,i}(r),\quad\textrm{where}\quad
C_{q,i}(r)=\frac{2k_{q,i}}{k_{q,i}^{2}-r^{2}/k_{q,N+1}^{2}}\frac{J_{q}(r/k_{q,N+1})}{J_{q+1}(k_{q,i})}.$
(12)

We discuss the error introduced by these two interpolations in the following
section.

###  2.3 Error analysis and complexity

####  2.3.1 Numerical interpolation

The numerical interpolation connecting the finer FFT grid to the
pseudospectral grid $\\{r_{q,i}\\}$ introduces an error that largely depends
on the resolution of the finer FFT grid and the choice of numerical method. We
choose the cubic spline interpolation, whose global error $E_{I}$ to
approximate the radial function $f_{q}(r)$ is bounded by [19]

$E_{I}\leq\frac{5}{384}\left|f_{q}^{(4)}\right|h^{4},$

where $h$ is the size of the radial interval of the FFT grid. According to our
sampling strategy explained in § 2.2, there are $2N$ radial FFT sampling
points and hence $h=(2N)^{-1}$. By assuming that the radial function is well-
behaved, the numerical interpolation error goes as

$E_{I}(N)\sim N^{-4}.$ (13)

####  2.3.2 Cardinal interpolation, general error of DHT

This section discusses the error introduced by the cardinal interpolation Eqn.
(12), which is also the general error associated with using discrete Hankel
transforms. According to the convergence theorem of Fourier-Bessel series, an
arbitrary radial function $f_{q}(r)$ can be fully recovered by Eq. (2). We
denote its truncation up to the $N^{\text{th}}$ term by

$f_{q}^{N}(r)\equiv\sum_{j=1}^{N}a_{q,j}J_{q}\left(k_{q,j}r\right).$ (14)

The truncation error $E_{T}$ is then

$E_{T}(q,N)\equiv\left|f_{q}^{N}(r)-f_{q}(r)\right|\leq\sum_{j=N+1}^{\infty}|a_{q,j}|\sup_{x>0}\left\\{|J_{q}\left(x\right)|\right\\}.$
(15)

The $q^{\textrm{th}}$-order Bessel function of the first kind is characterized
by [20]

$|J_{q}(x)|<b_{q}q^{-1/3},$

where the $b_{q}$’s are such that $b_{1}\simeq 0.58$, increasing monotonically
to around 0.67 as $q$ tends to infinity. We can conclude that, for orders
$q\geq 1$, the truncation error is

$E_{T}(q,N)<b_{q}q^{-1/3}\sum_{j=N+1}^{\infty}|a_{q,j}|.$ (16)

For $q=0$, $\sup_{x}\left\\{|J_{0}\left(x\right)|\right\\}=1$ and therefore,

$E_{T}(0,N)\leq\sum_{j=N+1}^{\infty}|a_{0,j}|.$ (17)

For the algebraically-converging Fourier-Bessel series, the magnitude of the
truncation error is given by [15, § 2.13]

$E_{T}(q=0,N)\sim N|a_{0,N}|,\quad\textrm{and}\quad E_{T}(q\geq 1,N)\sim
b_{q}q^{-1/3}N|a_{q,N}|.$ (18)

Now we denote the interpolated function approximated by $N$ radial modes via
DHT’s cardinal functions, i.e., Eqn. (12), as $S_{q}^{N}(r)$. It agrees with
the original radial function $f_{q}(r)$ on the pseudospectral grid points
$\\{r_{q,i}\\}$ defined in Eqn. (4), namely,

$S_{q}^{N}(r)=\sum_{i=1}^{N}f_{q}(r_{q,i})C_{q,i}(r)=2\sum_{j=1}^{N}\frac{F_{q,j}}{J_{q+1}^{2}\left(k_{q,j}\right)}J_{q}\left(k_{q,j}r\right),\quad\\\
\textrm{and}\quad\left\\{S_{q}^{N}(r_{q,i})\right\\}_{i=1,\dots,N}=\left\\{f_{q}(r_{q,i})\right\\}_{i=1,\dots,N},$
(19)

where the $F_{q,j}$’s are the result of the DHT, i.e., Eqn. (8). Note that
because $J_{q}(k_{q,j})=0$ by definition of the $k_{q,j}$’s, the Dirichlet
boundary condition is automatically enforced on the radial function, i.e.,
$f_{q}(1)=0$. This solves the problem that the pseudospectral grid does not
sample the boundary $r=1\notin\\{r_{q,i}\\}$. Using the discrete orthogonality
relation Eqn. (6), we show in A that the discrete Hankel transform on $N$
pseudospectral grid points produces the exact first $N$ Fourier-Bessel
coefficients $a_{q,j}$. In other words,

$F_{q,j}=\frac{J_{q+1}^{2}(k_{q,j})}{2}a_{q,j}=\int_{0}^{1}rf_{q}(r)J_{q}(k_{q,j}r)\,\textrm{d}r,\quad\textrm{and
thus}\quad S_{q}^{N}(r)=f_{q}^{N}(r).$ (20)

Therefore, the discrete Hankel transform introduces the same amount of error
as truncating an infinite Fourier-Bessel series would (i.e., $E_{T}$). From
Eqn. (18) and the convergence rate of Fourier-Bessel series for functions
satisfying the homogeneous Dirichlet condition, $a_{q,N}\lesssim N^{-5/2}$, we
obtain the general error of DHTs,

$E_{D}(q=0,N)\lesssim N^{-3/2},\quad\textrm{and}\quad E_{D}(q\geq 1,N)\lesssim
b_{q}q^{-1/3}N^{-3/2}.$ (21)

Note that ‘$\lesssim$’ above denotes the worst case, whereas the actual error
may be smaller depending on the radial function $f_{q}(r)$. Indeed, the
convergence rate of Fourier-Bessel series depends mainly on the boundary
behavior of the function [3]. The series has a faster convergence rate and
hence a smaller error if the function satisfies the following condition at the
boundary up to its $p^{\textrm{th}}$ derivative:

$f_{q}(1)=f_{q}^{\prime}(1)=f_{q}^{\prime\prime}(1)=\dots=f_{q}^{(p)}(1)=0.$
(22)

The higher $p$, the faster the convergence rate. Following Boyd’s result, such
an error bound is given by $E_{D}(p,N)\sim N^{-2p-3/2}$. That is, for the
worst-case scenario in which $p=0$, i.e., in which the homogeneous Dirichlet
condition is satisfied for the function and not its first derivative, Eqn.
(21) gives the error. For the best case, $p=\infty$, we have an exponential
rate of convergence. The DHT error is then governed by
$E_{D}(N)\sim\text{e}^{-N}$, although this occurs very rarely.

####  2.3.3 Interpolation vs. DHT errors

We saw that there were two competing sources of error in our method, the
interpolation error $E_{I}$ (Eqn. 13) and the DHT error $E_{D}$ due to the
truncation of the infinite Fourier-Bessel series (Eqn. 21). Because $E_{I}\sim
N^{-4}$ vs. $E_{D}\sim N^{-3/2}$ in the worst-case scenario described above
(case $p=0$ in Eqn. 22), the truncation error dominates. However, if the
function is such that $p$ is large (cf. Eq. 22), numerical interpolation can
become the main source of error.

All the error analyzes conducted so far are in the radial direction,
associated with decomposing the radial function $f_{q}(r)$. Indeed, in the
azimuthal direction, FFTs have an exponential rate of convergence; thus, their
error is negligible compared to the largest of the interpolation or DHT
errors.

####  2.3.4 Runtime analysis

Since the main purpose of the method is to solve time-dependent initial value
problems, multiple transforms are executed for a large number of time
iterations. Therefore, the complexity of the operations is of great concern.
We now analyze how the number of operations scales with the spatial
resolution. Consider the function $\psi(r,\theta)$, sampled with $2N$
azimuthal and radial points, where $N$ is the number of radial modes. FFTs on
such grids have a complexity of $4N^{2}\log(2N)$.

After the Fourier transform, we perform a cubic spline interpolation
(complexity $N$ [21]) within each angular mode $q$, where the $q$’s are
integers ranging from $-N$ to $N$. The complexity of the numerical
interpolation with non-periodic boundary conditions is therefore $2N^{2}$.

The discrete Hankel transform is a multiplication of a $N$-vector with a
$N\times N$ dense matrix. The total complexity is then $2N^{3}$. Thus, the DHT
operation governs the forward transform from $\psi(r,\theta)$ to spectral
coefficients $\\{a_{q,j}\\}$, the complexity is given by

$T_{\textrm{forward}}(N)\sim 2N^{3}.$

Conversely, for the backward transform, the FFT and DHT terms remain the same.
We use the cardinal interpolation introduced in § 2.2, instead of cubic
splines. It is equivalent to a $2N\times N$ matrix with a $N$ vector
multiplication. Thus, its complexity is $2N^{3}$. Hence, both cardinal
interpolation and DHT terms govern the backward transform, the complexity is
then

$T_{\textrm{backward}}(N)\sim 4N^{3}.$

Thus, the total complexity associated with the method is $O[N^{3}]$. The
matrices for DHT and cardinal interpolation Eqn. (8), (9) and (12), can be
pre-computed outside temporal loops and do not affect the complexity when
solving initial value problems.

###  2.4 Applying the Laplacian operator and spectral time scheme

Bessel functions are the eigenfunctions of the Laplacian operator in polar
coordinates. That is, for the function $\psi(r,\theta)$ decomposed under our
method with $N$ azimuthal and radial modes, namely,

$\psi(r,\theta)=\sum_{q=-N}^{N}\sum_{j=1}^{N}a_{q,j}\text{e}^{\text{i}q\theta}J_{q}(k_{q,j}r),$
(23)

the Laplacian of this function becomes

$\nabla^{2}\psi(r,\theta)=\sum_{q=-N}^{N}\sum_{j=1}^{N}-k_{q,j}^{2}a_{q,j}\text{e}^{\text{i}q\theta}J_{q}(k_{q,j}r).$
(24)

The set of basis functions remains unchanged. A pure spectral time scheme can
be applied here and the propagators are the Bessel roots $k_{q,j}$. The
implementation is straightforward: we use the transform method introduced in
previous sections to obtain the spectral coefficients $a_{q,j}$, then apply a
time scheme (the operator splitting in our case), and finally transform back
to the physical domain. The homogeneous Dirichlet boundary condition is
enforced by the Bessel functions without extra operations.

To solve PDEs associated with operators other than the Laplacian, we use the
operator splitting method [22]. It can be applied to differential equations of
the form

$\frac{\partial\psi}{\partial t}=c\nabla^{2}\psi+\mathfrak{D}[\psi],$ (25)

where $c$ is a constant, by creating an intermediate state
$\psi_{\textrm{int}}$ that would solve
$\partial_{t}\psi_{\textrm{int}}=c\nabla^{2}\psi$. In our case,

$\psi_{\textrm{int}}=\sum_{q=-N}^{N}\sum_{j=1}^{N}\left\\{a_{q,j}\text{e}^{-ck_{q,j}^{2}\Delta
t}\right\\}\text{e}^{\text{i}q\theta}J_{q}(k_{q,j}r)$ (26)

where $\Delta t$ is the time step. The other operator $\mathfrak{D}$ now acts
on $\psi_{\textrm{int}}$, and can be treated with other numerical methods. Due
to the nature of this spectral time scheme, calculating $\psi_{\textrm{int}}$
is mathematically exact, the sole source of error is from approximating the
function via DHTs. For PDEs associated with second-order time derivatives,
readers can refer to the first example illustrated in § 3.1. The method is
fast, accurate, easy to program, and it has many applications in polar
coordinates associated with Laplacian operators, as we demonstrate next.

##  3 Examples

###  3.1 Drum beats

Vibrations of an elastic membrane are governed by the 2-D wave equation, that
is,

$\frac{\partial^{2}u}{\partial t^{2}}=c^{2}\nabla^{2}u,$ (27)

where $u(r,\theta,t)$ denotes the transverse displacement of the membrane and
$c$ is the wave speed. It has the following general solution in 2-D polar
coordinates under homogeneous Dirichlet conditions,

$u(r,\theta,t)=\sum_{q=-\infty}^{\infty}\sum_{j=1}^{\infty}\left[a_{q,j}\cos(ck_{q,j}t)+\frac{b_{q,j}}{ck_{q,j}}\sin(ck_{q,j}t)\right]\text{e}^{\text{i}q\theta}J_{q}(k_{q,j}r),$
(28)

where $a_{q,j}$ and $b_{q,j}$ are the Fourier-Bessel coefficients of the
initial displacement and velocity fields, namely,

$\begin{cases}\sum_{q}\sum_{j}a_{q,j}\text{e}^{\text{i}q\theta}J_{q}(k_{q,j}r)=u|_{t=0}\\\
\sum_{q}\sum_{j}b_{q,j}\text{e}^{\text{i}q\theta}J_{q}(k_{q,j}r)=\partial_{t}u|_{t=0},\end{cases}$
(29)

can be directly obtained by the method introduced in this paper. This example
demonstrates the great advantage of spectral methods with eigenfunction bases
to solve linear equations, because the state at any arbitrary future time $t$
can be determined directly from the initial conditions. The time scheme is
fully spectral and no time evolution is required.

To validate the method, we construct a circular membrane with initial
displacement consisting of two off-centered Gaussians with a prefactor term
$\alpha(r)$ such that we can control its boundary behavior, and no initial
velocity (in which case $b_{q,j}=0$ $\forall(q,j)$). In other words,

$u|_{t=0}=\alpha(r)\left(\text{e}^{-15[r^{2}+r_{0}^{2}-2rr_{0}\cos(\theta-\theta_{0})]}-\text{e}^{-15[r^{2}+r_{1}^{2}-2rr_{1}\cos(\theta-\theta_{1})]}\right)\quad\text{and}\quad\partial_{t}u|_{t=0}=0,$
(30)

where $r_{0},\theta_{0}$, $r_{1}$, and $\theta_{1}$ locate the centers of the
two Gaussian peaks, respectively,

$\begin{cases}r_{0}=r_{1}=0.3\\\ \theta_{0}=\theta_{1}-\pi=1.\end{cases}$

We perform pairs of forward and backward transforms without time propagation
with different $N$’s (number of azimuthal and radial modes) on Eqn. (30) to
check the convergence of the method. We define the total error as

$E_{\textrm{tot}}=\int\frac{|u_{\textrm{transformed}}-u|}{\pi\max(u)}\,\textrm{d}\mathbf{a},$
(31)

where $\,\textrm{d}\mathbf{a}=r\,\textrm{d}r\,\textrm{d}\theta$ is the
infinitesimal patch of area. Figure 2 shows $E_{\textrm{tot}}$ vs. $N$ (number
of radial and azimuthal modes) with three different prefactors $\alpha(r)$.

Figure 2: Comparison of the method’s convergence rates on Eqn. (30) with
different prefactors.

For $\alpha(r)=1-r$ (blue dots), the convergence rate starts at $\sim
N^{-3.57}$, then decreases to $\sim N^{-2.21}$ when $N$ reaches $10^{2}$. The
displacement field satisfies the Dirichlet conditions, but its radial
derivative is non-zero at the boundary, and thus $p=0$, as discussed in
§2.3.2. However, the derivative is small: $\partial_{r}u|_{r=1}\simeq
6.4\times 10^{-4}$ according to our calculation. This explains why the
convergence rate here is faster than the theoretical bound; that is, the
method is less sensitive to the function’s higher-order derivatives’ boundary
conditions when the resolution is low ($N$ is small). It is not practical to
check the convergence at very large $N$, but we can predict that it will
eventually reach the worst-case error bound $E_{D}\sim N^{-3/2}$ as
$N\rightarrow\infty$. This is a positive feature, since in practical usages we
generally choose $N\sim 10^{2}$, and for such $N$’s, the method converges
faster and thus has less error than the theory.

For $\alpha(r)=\exp\left[-(1-x^{2})^{-1}\right]$, namely, the bump function
(black triangles), the convergence rate is given by $\sim N^{-4.12}$, and
shows no sign of decreasing. The bump function has the property that all its
derivatives cancel at $r=1$, corresponding to $p=\infty$ in § 2.3.2 and the
Fourier-Bessel series of the resulting initial condition has an exponential
convergence rate. Thus, the error bound of the method is no longer governed by
the DHT, but by the cubic spline interpolation $E_{I}\sim N^{-4}$, § 2.3.1,
which agrees with our numerical result.

Finally, we let $\alpha(r)=1$, thus the displacement field no longer satisfies
the homogeneous Dirichlet conditions. The convergence rate (purple squares)
drops below the theoretical worst-case error bound to around $\sim N^{-1}$, as
the theory predicts.

We have verified our statements on the error bound and demonstrated the
boundary-dependent nature of the method. When applying the method to a
function, we should always consider the boundary behavior to achieve the
desired accuracy.

Figure 3 shows an example using the method to solve the evolution of (30) with
$\alpha(r)=1-r$ at six different times. For simplicity, we set the sound speed
to unity ($c=1$). A remarkable feature of this example is that the system has
an infinite period if we have more than one mode, since the Bessel roots
$k_{q,j}$’s are transcendental and have no multiple of one another [23].

Figure 3: Solving the wave equation (27) in polar coordinates. Top-left panel:
Initial condition of the membrane displacement $u(r,\theta,t=0)$, consisting
of two off-centred Gaussians with the prefactor $\alpha=1-r$ (cf. Eq. 30).
Other panels: Membrane displacements at $t=\\{0.5,1.0,1.5,2.0,2.5\\}$.

###  3.2 Impulse start of a Poiseuille pipe flow

Our first non-linear example is the time-dependent Poiseuille pipe flow
described by

$\partial_{t}u=G(r,\theta,t)+Re^{-1}\nabla^{2}u\quad\text{and}\quad
u|_{r=1}=0,$ (32)

which solves for the velocity field $u$ of an incompressible, homogeneous
laminar flow across a unit circle, where $G(r,\theta,t)$ is the normalized
pressure gradient and $Re$ is a constant called the Reynolds number of the
flow. The initial condition is $u|_{t=0}=0$ everywhere on the circle.

First, we consider the simplest case where the pressure gradient $G$ is
constant throughout the domain at positive times. The problem now has radial
symmetry and only the $q=0$ angular mode is involved, reducing to a one-
dimensional (1D) problem, whose analytical solution is [24, § 4.3]

$u(r,t)=\left[\frac{1-r^{2}}{4}-2\sum_{i=1}^{\infty}\frac{J_{0}(k_{0,i}r)}{k_{0,i}^{3}J_{1}(k_{0,i})}\exp\left(-\frac{k_{0,i}^{2}t}{Re}\right)\right]GRe.$
(33)

Note that as $t\rightarrow\infty$, we recover the classical stationary
Poiseuille flow $u_{\infty}=(1-r^{2})GRe/4$.

Our simulation time-steps the prognostic equation in (32) with the operator-
splitting method described in § 2.4. The constant pressure gradient is dealt
with the forward Euler method, i.e., $\Delta u=G\Delta t$. Figures 4(a) and
5(a) show the time evolution of a slice of the velocity field with $GRe=1$ and
$\Delta t=10^{-3}$. We use $N=16$ and $32$ radial modes, respectively. The
blue triangles are the simulated values at pseudospectral grid points
$\\{r_{0,i}\\}$. Figures 4(b) and 5(b) show the relative percentage errors
with respect to the analytical solution,
$|u_{\textrm{simulated}}-u_{\textrm{analytic}}|/u_{\textrm{analytic}}$. The
error is below 4% except in the boundary layer, where the velocity field in
the denominator is close to zero. Increasing the number of radial modes $N$
yields a smoother error plot and a diminished the error boundary layer.
However, the overall percentage is not significantly reduced. This is due to
the error introduced by the Euler method, which we verified by reducing the
time step (see Figures 4(c) and 5(c)).

(a) Time evolution.

(b) Relative percentage error, $\Delta t=10^{-3}$.

(c) Relative percentage error, $\Delta t=10^{-4}$.

Figure 4: Simulation of the impulse start of a radially-symmetric Poiseuille
pipe flow with $N=16$ radial modes.

(a) Time evolution.

(b) Relative percentage error, $\Delta t=10^{-3}$.

(c) Relative percentage error, $\Delta t=10^{-4}$.

Figure 5: Same as for Figure 4 for $N=32$ radial modes.

Figure 6 illustrates an example subjected to a nonconstant pressure gradient
$G=y=r\sin{\theta}$ (varies along the vertical axis) and $Re=1$. The problem
is now two-dimensional, and multiple angular modes are involved. Although no
analytical solution is available for such a pressure gradient, the resulting
steady velocity field satisfies the Dirichlet condition, and the solution
respects the symmetry of the applied pressure gradient. With $N=15$ azimuthal
and radial modes, and $\Delta t=10^{-3}$, the simulation took a desktop Intel
Core-i7 8700K processor 7.33 seconds to reach $t=2$, or 2000 iterations.

Figure 6: Steady solution to Eqn. (32) at $t=2$ for $GRe=y=r\sin{\theta}$.

###  3.3 Quantum vortices in 2-D Bose-Einstein condensate

This example was studied by Grisouard [25], who laid out the basic features of
our method. We consider a Bose-Einstein condensate (BEC) inside a circular
infinite potential well. The dynamics of BECs are governed by the time-
dependent Gross-Pitaevskii equation (GPE) [26, 27], whose nondimensional form
is given by

$\frac{\partial\psi}{\partial
t}=\frac{\text{i}}{2}\nabla^{2}\psi+\frac{\text{i}}{2\xi^{2}}(1-|\psi|^{2})\psi,\quad\text{where}\quad\xi=\frac{\hbar}{R\sqrt{2mU_{0}n_{0}}}$
(34)

is a dimensionless so-called “healing length” built from the system’s
parameters, with $\hbar$ the reduced Planck constant, $R$ the radius of the
potential well that serves as our unit length, $m$ the particle mass, $U_{0}$
a constant potential representing repulsion between particles, and $n_{0}$ the
normalization constant of the probability density function $|\psi|^{2}$. The
details of nondimensionalization are further explained by Guo & Bühler [28].
The external potential is set to be infinite (nullifying $\psi$) outside the
unit radius ($R$ in dimensional units), such that the wave function $\psi$
satisfies the homogeneous Dirichlet condition.

Vortex circulation in a BEC is quantized and thus vortex strength cannot be
continuously varied. We initialize a point vortex of circulation $\Gamma=2\pi
p$, $p\in\mathbb{Z}$, centered at $(r_{0},\theta_{0})$ by adding the wave
function

$\psi_{\textrm{v}}^{\\{r_{0};\theta_{0};p\\}}(r,\theta)=\frac{r\text{e}^{\text{i}p\theta}-r_{0}\text{e}^{\text{i}p\theta_{0}}}{\sqrt{\xi^{2}+\left|r\text{e}^{\text{i}p\theta}-r_{0}\text{e}^{\text{i}p\theta_{0}}\right|^{2}}}$
(35)

to an initial background wave function $\psi_{bg}$ that is meant to
approximate the steady solution to the GPE, such that

$\psi(t=0)=\psi_{\textrm{bg}}\psi_{\textrm{v}}^{\\{r_{0};\theta_{0};p\\}},\quad\textrm{where}\quad\psi_{\textrm{bg}}(r,\theta)=\tanh\left(\frac{1-r}{\sqrt{2}\xi}\right).$
(36)

The form, $\psi_{bg}$ takes, derives from the steady solution to the GPE in
the case of an infinite, straight wall located at $x=0$, namely,
$\psi_{\infty}=\tanh(x/2\xi)$ [29]. That is, a smooth transition between
$\psi_{\infty}=0$ at the wall and $\psi_{\infty}\approx 1$ a few units of
$\xi$ away from it. In our circular domain, this solution is only approximate:
having $\psi(t=0)=\psi_{bg}$ as our sole initial condition results in the
creation of waves, oscillating about the exact and steady solution. Here, we
use $\xi=0.1$ and therefore, $\psi_{bg}$ is a good approximation of the steady
solution to the GPE, i.e., radiation of transient waves is weak.

BEC dynamics share many similarities with classical shallow water systems, as
revealed by the Madelung transformation [30]. Thus, vortex motion can be
analyzed classically. For example, the time for a point vortex to loop around
a unit disk can be calculated by using the method of images. As explained by
Mason et al. [29], due to the boundary layer effect caused by the infinite
potential well, the actual free-slip boundary (virtual wall) for the method of
images to be applied is a distance $\sqrt{2}\xi$ away from the potential wall.
The loop time is then given by [31]

$\tau_{\textrm{loop}}=\frac{4\pi^{2}[(1-\sqrt{2}\xi)^{2}-r_{0}^{2}]}{\Gamma}.$
(37)

We initialize a simulation with an initial wave function described by Eqn.
(36). We set $\xi=0.1$ and the size of the time step is $\Delta t=5\times
10^{-5}$. We sample the wave function with $256$ angular and radial points, so
there are $N=128$ radial modes involved in the simulation, as discussed in §
2.2 and 2.3.4. The nondimensional equation (34) is treated with the operator
splitting method mentioned in § 2.4. For the nonlinear part, we have

$\frac{\partial\psi}{\partial
t}=\frac{\text{i}}{2\xi^{2}}\left(1-|\psi|^{2}\right)\psi,$ (38)

which can be solved analytically because this equation leaves $|\psi|$
unchanged. Indeed, the variations of $\psi$ go as $\text{i}K\psi$, where $K$
is real, meaning that only the phase of $\psi$ changes. Thus, taking
$\psi_{\textrm{int}}$ as the initial condition, we can integrate over one time
step as

$\psi(t+\Delta
t)=\psi_{\textrm{int}}\exp\left[\frac{\text{i}\left(1-|\psi|^{2}\right)}{2\xi^{2}}\Delta
t\right],$ (39)

where we recall that $\psi_{\textrm{int}}$ is computed by Eqn. (26). Figure 7
shows $|\psi(t=1)|^{2}$ in the presence of a point vortex with $r_{0}=0.6$.

Figure 7: Probability density $|\psi|^{2}$ for a vortex with $r_{0}=0.6$ at
$t=1$. The vortex is at the centre of the patch of lower values of
$|\psi|^{2}$ in the top left quadrant.

Figure 8 shows the invariants of the vortex with $r_{0}=0.6$ over the time
period $\Delta t=1.5$. The total mass $n_{\textrm{tot}}$, angular momentum
$L$, and energy $E$ of the BEC are, respectively, [25, 28]

$n_{\textrm{tot}}=\int|\psi|^{2}\,\textrm{d}\textbf{a};$ (40a)
$L=-\text{i}\int\psi^{*}\frac{\partial\psi}{\partial\theta}\,\textrm{d}\textbf{a},\quad\text{and}\quad$
(40b)
$E=\frac{1}{2}\int\left[|\nabla\psi|^{2}+\frac{1}{2\xi^{2}}\left(|\psi|^{2}+1\right)^{2}\right]\,\textrm{d}\textbf{a}.$
(40c)

The mass is constant within 0.4‱ (Figure 8(a)), the angular momentum
fluctuates within 0.5‱ (Figure 8(b)), while the energy fluctuates by about
2.7% (Figure 8(c)). The energy expression (Eqn. 40c) includes the most
derivatives, both azimuthal and radial. The former derivative is calculated
spectrally using FFTs, and the latter via a second-order finite-difference
scheme. The latter suggests that the fluctuations might be due to the lack of
precision of our postprocessing methods, rather than to that of the spectral
method itself.

(a) Total mass

(b) Angular momentum

(c) Energy

Figure 8: Evolution of the invariants for a BEC vortex initially centered at
$r_{0}=0.6$.

We then run a suite of single-vortex simulations, each with a different
initial radius for the vortex center: $r_{0}=0.3$, $0.4$, $0.5$, $0.6$, $0.7$,
$0.75$, $0.8$ and $0.9$. All vortices start from $\theta_{0}=0$. Figure 9
shows the trajectories of the vortex centers over a period of time of $\Delta
t=1.5$. The marker pairs show the start and end positions of the vortex
centers. As mentioned earlier, the vortex centre moves parallel to the
circular wall, describing circles around the domain center. However, those
circles are not perfect due to interactions with BEC waves.

Figure 9: Traces of vortices’ centers over a period of 1.5 unit time.

Figure 10 shows measurements of the time $\tau_{0}$ it takes for a vortex
center to describe a full circle around the center of the domain, versus the
theoretical values predicted by the image method (Eqn. 37). As expected, the
experimental results deviate from the theoretical values when the vortex
center is close to or inside the virtual wall at $r=1-\xi\sqrt{2}$.

Figure 10: Theoretical loop time $\tau_{0}$ v.s. experimental values.

By comparing the looping times and computing the invariants, we validated our
method on a nonlinear, pseudo-spectral simulation. This example also
demonstrates the method at its full potential. Each simulation requires 30,000
iterations, which takes around 40 minutes to run on the same system as in §
3.2. The method solves a PDE involving Laplacian on a complex field with
$128\times 128$ angular and radial modes below 0.1 s.

##  4 Conclusions and perspectives

This paper introduces a novel spectral method suited for Cauchy problems
featuring Laplacians in polar coordinates. We reformulated the discrete Hankel
transforms as pseudospectral methods, and we introduced DHTs’ quadrature
weights, cardinal functions, and pseudospectral grids. We discussed the choice
of sampling grids for both FFTs and DHTs, and analyzed the error bounds of the
numerical interpolation and DHTs. We demonstrated a systematic way to use the
method to compute Laplacians and apply spectral time schemes to time-dependent
partial differential equations. We then showed a range of time-dependent
examples including linear and nonlinear PDEs, fields with complex numbers, and
the 2-D wave equation involving the second order time derivative. We verified
the theoretical error bound on functions with different boundary conditions.
The method was validated by comparing the simulation results with analytical
solutions and by checking energy, mass, and angular momentum conservations.
Although all our equations and examples are on unit disks throughout this
paper for concision, the method can be rescaled to arbitrary radii [17, 13].

A way to greatly speed up the method is to implement the program for parallel
computing. After the FFT step, for which parallel libraries are widely
available, each angular mode is independent of the others, and the following
interpolations and DHTs can be distributed to different processes at the same
time. Due to the nature of Bessel functions, the method converges poorly if
the homogeneous Dirichlet condition is not satisfied. However, there are many
known ways to homogenize boundary conditions [32, 33, 34], which may be
combined with our method for more general problems. We are currently
considering a new transform based on Dini series, a variation of the Fourier-
Bessel series which works on Neumann and Robin boundary conditions [9, §
18.11]. A combination of these transforms has the potential to compute
Laplacians under arbitrary boundary conditions. Holman and Kunyansky suggest a
reduced polar sampling grid [35], which becomes more sparse towards the domain
center. Our current sampling grid has redundancy in the inner part at high
angular modes, and adopting their strategy may improve efficiency.

## References

  * [1] K. Mohseni, T. Colonius, Numerical Treatment of Polar Coordinate Singularities, Journal of Computational Physics 157 (2) (2000) 787–795.  doi:10.1006/jcph.1999.6382.
  * [2] B. Prochnow, O. O’Reilly, E. M. Dunham, N. A. Petersson, Treatment of the polar coordinate singularity in axisymmetric wave propagation using high-order summation-by-parts operators on a staggered grid, Computers and Fluids 149 (2017) 138–149.  doi:10.1016/j.compfluid.2017.03.015.
  * [3] J. P. Boyd, 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 confor, Journal of Computational Physics 230 (4) (2011) 1408–1438.  doi:10.1016/j.jcp.2010.11.011.
  * [4] A. B. Bhatia, E. Wolf, On the circle polynomials of Zernike and related orthogonal sets, Mathematical Proceedings of the Cambridge Philosophical Society 50 (1) (1954) 40–48.  doi:10.1017/S0305004100029066.
  * [5] K. Atkinson, D. Chien, O. Hansen, A spectral method for elliptic equations: The Dirichlet problem, Advances in Computational Mathematics 33 (2) (2010) 169–189.  arXiv:arXiv:0809.0021v1, doi:10.1007/s10444-009-9125-8.
  * [6] J. Shen, New fast Chebyshev-Fourier algorithm for Poisson-type equations in polar geometries, Applied Numerical Mathematics 33 (1) (2000) 183–190.  doi:10.1016/S0168-9274(99)00082-3.
  * [7] T. Matsushima, P. S. Marcus, A spectral method for polar coordinates, Journal of Computational Physics 120 (2) (1995) 365–374.  doi:10.1006/jcph.1995.1171.
  * [8] G. M. Vasil, K. J. Burns, D. Lecoanet, S. Olver, B. P. Brown, J. S. Oishi, Tensor calculus in polar coordinates using Jacobi polynomials, Journal of Computational Physics 325 (2016) 53–73.  arXiv:1509.07624, doi:10.1016/j.jcp.2016.08.013.
  * [9] G. N. Watson, A Treatise on the Theory of Bessel Functions, Cambridge Mathematical Library, Cambridge University Press, 1995.
  * [10] N. Temme, An algorithm with ALGOL 60 program for the computation of the zeros of ordinary bessel functions and those of their derivatives, Journal of Computational Physics 32 (2) (1979) 270–279.  doi:10.1016/0021-9991(79)90134-7.
  * [11] H. F. Johnson, An improved method for computing a discrete Hankel transform, Computer Physics Communications 43 (2) (1987) 181–202.  doi:10.1016/0010-4655(87)90204-9.
  * [12] L. Yu, M. Huang, M. Chen, W. Chen, W. Huang, Z. Zhu, Quasi-discrete Hankel transform, Opt. Lett. 23 (6) (1998) 409–411.  doi:10.1364/OL.23.000409.
  * [13] M. Guizar-Sicairos, J. C. Gutierrez-Vega, Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields, Journal of the Optical Society of America A 21 (1) (2004) 53.  doi:10.1364/JOSAA.21.000053.
  * [14] N. Baddour, The Discrete Hankel Transform, in: G. S. Nikolic, D. Z. Markovic-Nikolic (Eds.), Fourier Transforms - Century of Digitalization and Increasing Expectations, Vol. 32, IntechOpen, London, UK, 2019, pp. 137–144.  doi:10.5772/intechopen.84399.
  * [15] J. P. Boyd, Chebyshev and Fourier Spectral Methods: Second Revised Edition, Dover Books on Mathematics, Dover Publications, 2013.
  * [16] J. Shen, T. Tang, L.-L. Wang, Spectral Methods, Vol. 41 of Springer Series in Computational Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011\.  doi:10.1007/978-3-540-71041-7.
  * [17] N. Baddour, Discrete two-dimensional fourier transform in polar coordinates part I: Theory and operational rules, Mathematics 7 (8) (2019) 1–32.  doi:10.3390/math7080698.
  * [18] X. Yao, N. Baddour, Discrete two dimensional Fourier Transform in polar coordinates Part II: Numerical computation and approximation of the continuous transform, PeerJ Computer Science 2020 (3) (2020) 1–38.  doi:10.7717/peerj-cs.257.
  * [19] C. A. Hall, W. Meyer, Optimal error bounds for cubic spline interpolation, Journal of Approximation Theory 16 (2) (1976) 105–122.  doi:https://doi.org/10.1016/0021-9045(76)90040-X.
  * [20] L. J. Landau, Bessel functions: Monotonicity and bounds, Journal of the London Mathematical Society 61 (1) (2000) 197–215.  doi:10.1112/S0024610799008352.
  * [21] K. Toraichi, K. Katagishi, I. Sekita, R. Mori, Computational complexity of spline interpolation, International Journal of Systems Science 18 (5) (1987) 945–954.  doi:10.1080/00207728708964021.
  * [22] W. H., S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing, 3rd Edition, Cambridge University Press, USA, 2007.
  * [23] C. L. Siegel, Über einige anwendungen diophantischer approximationen, 1929, in: U. Zannier (Ed.), On Some Applications of Diophantine Approximations: a translation of Carl Ludwig Siegel’s Über einige Anwendungen diophantischer Approximationen by Clemens Fuchs, with a commentary and the article Integral points on curves: Siegel’s theorem after Siegel’s proof by Clemens Fuchs and Umberto Zannier, Scuola Normale Superiore, Pisa, 2014, pp. 81–138.  doi:10.1007/978-88-7642-520-2_2.
  * [24] G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge Mathematical Library, Cambridge University Press, 2000.  doi:10.1017/CBO9780511800955.
  * [25] N. Grisouard, Quantum vortices in a glass of bose-einstein condensate, in: O. Bühler, K. Helfrich (Eds.), 2009 program of studies: nonlinear waves, Woods Hole Oceanographic Institution, 2010, pp. 269–289.  doi:10.1575/1912/3213.
  * [26] E. P. Gross, Structure of a quantized vortex in boson systems, Nuovo Cimento 20 (1960) 454–477.  doi:10.1007/BF02731494.
  * [27] L. P. Pitaevskii, Vortex lines in an imperfect Bose gas, Sov. Phys. JETP 13 (2) (1961) 451–454.
  * [28] Y. Guo, O. Bühler, Wave-vortex interactions in the nonlinear schrödinger equation, Physics of Fluids 26 (2) (2014) 1–22.  doi:10.1063/1.4865837.
  * [29] P. Mason, N. G. Berloff, A. L. Fetter, Motion of a vortex line near the boundary of a semi-infinite uniform condensate, Physical Review A - Atomic, Molecular, and Optical Physics 74 (2006).  doi:10.1103/PhysRevA.74.043611.
  * [30] E. Madelung, Quantentheorie in hydrodynamischer Form, Zeitschrift für Physik 40 (3) (1927) 322–326.  doi:10.1007/BF01400372.
  * [31] P. K. Kundu, I. M. Cohen, D. R. Dowling, Chapter 5 - Vorticity Dynamics, in: P. K. Kundu, I. M. Cohen, D. R. Dowling (Eds.), Fluid Mechanics (Sixth Edition), sixth edition Edition, Academic Press, Boston, 2016, pp. 195–226.  doi:https://doi.org/10.1016/B978-0-12-405935-1.00005-8.
  * [32] P. Biswas, S. Singh, H. Bindra, Homogenization of time dependent boundary conditions for multi-layer heat conduction problem in cylindrical polar coordinates, International Journal of Heat and Mass Transfer 129 (2019) 721–734.  doi:https://doi.org/10.1016/j.ijheatmasstransfer.2018.10.036.
  * [33] A. Friedman, B. Hu, A non-stationary multi-scale oscillating free boundary for the laplace and heat equations, Journal of Differential Equations 137 (1) (1997) 119–165.  doi:10.1016/S0022-0396(06)80006-9.
  * [34] Y. Achdou, P. Le Tallec, F. Valentin, O. Pironneau, Constructing wall laws with Domain Decomposition or asymptotic expansion techniques, Computer Methods in Applied Mechanics and Engineering 151 (1-2) (1998) 215–232.  doi:10.1016/S0045-7825(97)00118-7.
  * [35] B. Holman, L. Kunyansky, DRAFT A Second-order Finite Difference Scheme For The Wave Equation on a Reduced Polar Grid, [Unpublished], Department of Mathematics, The University of Arizona (2015).
URL https://www.math.arizona.edu/~leonk/papers/polarFD7.pdf

##  Appendix A Proof of $S_{q}^{N}(r)=f_{q}^{N}(r)$

We prove the equality in Eqn. (20). By substituting the Fourier-Bessel
expansion Eqn. (2) into the forward transform formula Eqn. (8), we have

$F_{q,j}=\frac{2}{k_{q,N+1}^{2}}\sum_{i=1}^{N}\frac{J_{q}\left(k_{q,j}r_{q,i}\right)}{J_{q+1}^{2}\left(k_{q,i}\right)}\left[\sum_{m=1}^{\infty}a_{q,m}J_{q}\left(k_{q,m}r_{q,i}\right)\right].$
(41)

We rearrange this expression into

$F_{q,j}=\frac{2}{k_{q,N+1}^{2}}\sum_{i=1}^{N}\sum_{m=1}^{\infty}\frac{a_{q,m}}{J_{q+1}^{2}\left(k_{q,i}\right)}J_{q}\left(k_{q,j}r_{q,i}\right)J_{q}\left(k_{q,m}r_{q,i}\right).$
(42)

Summing the above equation over the pseudospectral grid points
$\\{r_{q,i}\\}$, from the orthogonality relation (6), we find the expression
we were looking for, namely,

$F_{q,j}=\frac{1}{2}\sum_{m=1}^{\infty}a_{q,m}J_{q+1}^{2}\left(k_{q,m}\right)\delta_{jm}=\frac{J_{q+1}^{2}(k_{q,j})}{2}a_{q,j}.$
(43)

Therefore, applying a DHT with $N$ radial modes is equivalent to truncating a
Fourier-Bessel series at the $N^{\textrm{th}}$ term.


# LPV sequential loop closing for high-precision motion systems*

Yorick Broens, Hans Butler and Roland Tóth *This work has received funding
from the ECSEL Joint Undertaking (JU) under grant agreement No 875999 and from
the Ministry of Innovation and Technology NRDI Office within the framework of
the Autonomous Systems National Laboratory Program.Y.Broens, H.Butler and
R.Tóth are with the Department of Electrical Engineering, Eindhoven University
of Technology, Eindhoven, The Netherlands. H.Butler is also affiliated with
ASML, Veldhoven, The Netherlands. R.Tóth is also affiliated with the Systems
and Control Laboratory, Institute for Computer Science and Control, Hungary,
(email: Y.L.C.Broens@tue.nl).

###### Abstract

Increasingly stringent throughput requirements in the industry necessitate the
need for lightweight design of high-precision motion systems to allow for high
accelerations, while still achieving accurate positioning of the moving-body.
The presence of position dependent dynamics in such motion systems severely
limits achievable position tracking performance using conventional sequential
loop closing (SLC) control design strategies. This paper presents a novel
extension of the conventional SLC design framework towards linear-parameter-
varying systems, which allows to circumvent limitations that are introduced by
position dependent effects in high-precision motion systems. Advantages of the
proposed control design approach are demonstrated in simulation using a high-
fidelity model of a moving-magnet planar actuator system, which exhibits
position dependency in both actuation and sensing.

##  I Introduction

Growing demands in the industry result in increasingly stringent requirements
on throughput and positioning accuracy of motion systems, such as
wafersteppers, scanners, pick-and-place machines and wire bonders, see [1, 2,
3, 4]. Traditionally, control design for multiple-input multiple-output (MIMO)
systems is simplified by excellent mechanical design which ensures high
stiffness and reproducibility of the design, see [5]. Moreover, highly stiff
mechanical designs result in motion behavior dominated by rigid body dynamics
which simplify the control design procedure and allow to use methods like
rigid body decoupling to efficiently handle MIMO systems, see [6]. In many
industrial applications, _sequential loop-closing_ (SLC) strategies are
applied for rigid body feedback control design due to several (practical)
advantages. First, these methods employ well-understood control design
strategies, such as loop-shaping techniques (see [7], [8]). Additionally, the
SLC framework allows for feedback control design using a _non-parameteric
model_ of the motion system, e.g. frequency response function (FRF)
measurements, thus circumventing the necessity of identifying an accurate
dynamic model of the system that can capture complicated high-frequent
resonance dynamics.

However, the ever increasing throughput demands in the industry necessitate
the design of lightweight mechanical structures to allow for ultra-high
accelerations, while keeping the power demands relatively low. Therefore,
increased throughput is obtained at the cost of introducing low-frequent
flexible dynamics due to the limited stiffness of the mechanical structure,
see [9].

Additionally, control design of high-precision motion systems is further
complicated by the presence of position dependent effects. For many high-
precision motion systems, position dependent resonance dynamics are introduced
by relative displacements of the moving-body with respect to the sensor frame
and/or actuation frame, resulting in the need for coordinate frame
transformations to relate the input forces and/or measurement signals to the
center of mass of the moving-body.

In general, the presence of flexible dynamics severely limits the achievable
feedback control bandwidth (see [10]), which is even more critical in case
they are position dependent, since position dependent effects manifest as
dynamic uncertainties. While previously SLC control design strategies handled
position dependency in terms of robustness of the resulting LTI controller,
the performance price of robustness in case of the novel lightweight designs
would be intolerable in practise.

In order to circumvent the aforementioned performance limitations, a novel
control design approach is presented, which extends the conventional _linear-
time-invariant_ SLC method (e.g. see [11]) to _linear-parameter varying_ (LPV)
systems. The proposed LPV SLC control design strategy preserves the advantages
of the traditional SLC approach, while allowing for increased closed-loop
performance for motion systems affected by position dependency. The presented
control design approach is based on local controller design with loop-shaping
strategies, thus allowing for an interpolation based parameterization of the
controller coefficients on the scheduling vector respectively. Additionally,
an efficient implementation strategy is presented, which allows for real-time
implementation of the designed position dependent feedback controllers.

The main contributions of this paper are:

  * (C1)

The development of an extension of the SLC control design approach, such that
position dependency is explicitly taken into consideration, thereby allowing
for additional degrees of freedom to the control design, such that increased
closed-loop performance is achieved.

  * (C2)

The development of an efficient implementation strategy for the position
dependent controllers, which allows for real-time implementation of the
proposed control concept.

This manuscript is organized as follows. First, the problem formulation is
presented in Section II. Next, Section III presents the conventional
sequential loop-closing framework for rigid-body feedback control design for
high-precision motion systems. In Section IV, the proposed LPV SLC framework
approach is introduced together with an efficient implementation strategy for
real-time application of the resulting position dependent controllers. Section
V provides a simulation study of the proposed LPV SLC control approach on a
high-fidelity model of a state-of-the-art moving-magnet planar actuator.
Finally, in Section VI, conclusions on the proposed methodology are drawn.

##  II Problem formulation

###  II-A Background

Many high-precision motion systems exhibit position dependent effects, see
Figure 1, which are introduced by relative sensing and actuation of the moving
body. Therefore, such systems are often represented in an LPV form, see [12],
where position dependency is expressed by a so-called scheduling variable.
Consider the equations of motion of a mechatronic system that is subject to
relative actuation and relative sensing of the moving-body:

$M\ddot{q}(t)+D\dot{q}(t)+Kq(t)=\Phi_{a}(p(t))u(t),$ (1)

where $M$, $D$ and $K$ are the real symmetric mass, damping and stiffness
matrices of dimension $n_{q}\times n_{q}$ and
$\Phi_{a}(p(t))\in\mathbb{R}^{n_{q}\times n_{u}}$ maps the control forces
$u(t)$ to the appropriate masses based on the scheduling vector
$p:\mathbb{R}\rightarrow\mathbb{P}\subseteq\mathbb{R}^{n_{p}}$. Alternatively,
(1) is represented in LPV state-space form as:

$P=\left[\begin{array}[]{cc|c}0&I&0\\\
-M^{-1}K&-M^{-1}D&M^{-1}\Phi_{a}(p(t))\\\
\hline\cr\Phi_{s}(p(t))&0&0\end{array}\right],$ (2)

where $\Phi_{s}(p(t))\in\mathbb{R}^{n_{y}\times n_{q}}$ corresponds to a
mapping of $q(t)$ to the output $y(t)$ based on the scheduling vector $p(t)$.
Note that the $B(p(t))$ and $C(p(t))$ matrices are position dependent due to
relative actuation and sensing of the moving-body. If the scheduling vector
$p(t)$ is constant, implying $p(t)={\tt p}\in\mathbb{P}$ for all
$t\in\mathbb{R}$, (2) becomes an LTI system, which is often referred to as
_frozen dynamics_ for a particular fixed position of the motion system and is
denoted by $P_{\tt p}$. Typically, such a family of _frozen dynamics_ is used
in _sequential loop closing_ (SLC) control design strategies by viewing the
position dependent effects as a dynamic uncertainty. Consequently, during
conventional SLC control design, closed-loop performance is traded off against
robustness in order to ensure local stability of the closed-loop system.
Moreover, conventional SLC control design strategies are no longer sufficient
to meet with the increasing performance requirements in the industry.
Therefore, there is a necessity for a (practical) control design approach that
is able to push performance beyond LTI control design solutions, such that
satisfactory closed-loop performance is achieved.

###  II-B Problem statement

The problem that is addressed in this paper is to construct a position
dependent control law that considers all local aspects of the system by
following a loop-shaping based approach, capable to handle the MIMO nature of
the plant. The objective is to design a controller $K$, such that the
following requirements are satisfied:

  * (R1)

The closed-loop system is locally stabilized by $K$ for all $\tt
p\in\mathbb{P}$.

  * (R2)

The control design approach utilizes non-parametric models, e.g. FRF
measurements, thus circumventing the necessity of identifying an accurate
dynamic model of the system that can capture complicated high-frequent
flexible dynamics.

  * (R3)

The control design approach allows for modular design of the resulting
feedback controller $K$ such that a desirable structure can be imposed on the
controller.

The goal of this paper is the development of an extension of the conventional
LTI SLC control design framework towards a position dependent approach, such
that the position dependent nature of the system is exploited during feedback
control design in order to maximize closed-loop performance.

Figure 1: Low-frequent position dependent resonance dynamics of a high-
precision moving-magnet planar actuator system, showing the rigid body
decoupled transfer in $z$-direction for several _frozen positions_ of the
mover.

##  III Conventional SLC control design framework

This section gives a brief overview of the conventional SLC control design
framework (e.g. see [11]), which is widely employed in industrial
applications. The key concept of the SLC framework is that the multi-variable
control design is decomposed into a number of equivalent SISO control designs,
for which well-known control design tools are available, e.g. loop-shaping
techniques. The equivalent plant for the $i$-th design step of the _frozen_
high-precision motion system dynamics, denoted by $g_{\tt p}^{i}$, is
described by:

$g_{\tt p}^{i}=\mathcal{F}(P_{\tt p}^{i},-K^{i}),$ (3)

where $K^{i}$ = diag($k_{j}$) with $j=\\{1,\ldots,n_{u}\\}$ and $j\neq i$, see
[13]. The diagonal feedback controller corresponds to
$K=\text{diag}(k_{1},\ldots,k_{n_{u}})$. Furthermore, $P_{\tt
p}^{i}=W_{i}P_{\tt p}W_{i}$, where $W_{i}$ corresponds to an identity matrix
which has its $1$-st and $i$-th row interchanged. Closed loop stability of the
multi-variable rigid body decoupled system is assessed by considering the
characteristic equation $\mathrm{det}(I+P_{\tt p}K)\ \forall\tt
p\in\mathbb{P}$, which is related to the equivalent plant as (see [14]):

$\mathrm{det}(I+P_{\tt p}K)=\prod_{i=1}^{n_{u}}(1+g_{\tt p}^{i}k_{i})\
\forall\tt p\in\mathbb{P}$ (4)

In terms of (4), stability of the multi-variable system is reduced to the
assessment of $n_{u}$ SISO Nyquist criteria.

Generally, the controllers $k_{i}$ are designed in an iterative manner using
the rigid body decoupled dynamics of the corresponding high-precision motion
system, such that (local) stability is achieved for all _frozen positions_
$\tt p\in\mathbb{P}$, thus sacrificing performance for robustness during
control design due to the presence of position dependent flexible dynamics.

Figure 2: Typical structure of a SLC based rigid body feedback controller,
where the controller is composed from a cascade interconnection of sub-
components, such as a proportional gain, an integrator, a lead filter and
notch filter.

Typically, a SLC based rigid body feedback controller is structured as a
cascade interconnection of various LTI filters as illustrated by Figure 2,
where $C_{\mathrm{K}}^{i}$ corresponds to a proportional gain, which ensures
that the open loop dynamics of the equivalent plant $g_{\tt p}^{i}k_{i}$ cross
the 0 dB line at the desired rigid body target bandwidth. To suppress low-
frequent disturbances, integral action is imposed on $k_{i}$ by means of the
filter $C_{\mathrm{I}}^{i}$, which corresponds to:

$\begin{split}\dot{x}_{\mathrm{I}}(t)&=u_{\mathrm{I}}(t)\\\
y_{\mathrm{I}}(t)&=x_{\mathrm{I}}(t)\end{split}$ (5)

In order to stabilize the closed-loop system, lead filters are added to the
SLC based controller. The time domain representation of a lead filter
$C_{\mathrm{D}}^{i}$ is given by:

$\begin{split}\dot{x}_{\mathrm{D}}(t)&=-2\pi\alpha
f_{\mathrm{bw}}x_{\mathrm{D}}(t)+2\pi\alpha
f_{\mathrm{bw}}u_{\mathrm{D}}(t)\\\
y_{D}(t)&=1-\alpha^{2}x_{\mathrm{D}}(t)+\alpha^{2}u_{\mathrm{D}}(t)\end{split},$
(6)

such that approximately 45 degrees phase margin is present at the rigid body
target bandwidth $f_{\mathrm{bw}}>0$, while being subject to integral action
$C_{\mathrm{I}}^{i}$. The parameter $\alpha>0$ is used to shift the cut-off
frequencies of the differential action and the integral action of the lead
filter, such that sufficient phase lead is present at the desired target
bandwidth. Typically, $\alpha$ is chosen to be 3. At last, a notch filter
$C_{\mathrm{N}}^{i}$ is imposed on $k_{i}$ to suppress the effects of position
dependent flexible dynamics. The time domain representation of a notch filter
$C_{\mathrm{N}}^{i}$ is denoted by:

$\begin{split}\dot{x}_{\mathrm{N}}(t)&=\begin{bmatrix}-4\pi\beta_{2}f_{2}&-4f_{2}^{2}\pi^{2}\\\
1&0\end{bmatrix}x_{\mathrm{N}}(t)+\begin{bmatrix}4f_{2}^{2}\pi^{2}\\\
0\end{bmatrix}u_{\mathrm{N}}(t)\\\
y_{\mathrm{N}}(t)&=\begin{bmatrix}\frac{\beta_{1}f_{1}-\beta_{2}f_{2}}{f_{1}^{2}\pi}&1-\frac{f_{2}^{2}}{f_{1}^{2}}\end{bmatrix}x_{\mathrm{N}}(t)+\frac{f_{2}^{2}}{f_{1}^{2}}u_{\mathrm{N}}(t)\end{split},$
(7)

where $\beta_{1}$, $\beta_{2}$, $f_{1}$ and $f_{2}$ correspond to the notch
filter coefficients respectively.

##  IV LPV SLC control design approach

This section presents a novel extension of the conventional SLC control design
framework towards a position dependent approach, which allows to push closed-
loop performance of the position dependent system beyond LTI control
strategies.

Similar to the SLC control design framework, the (local) equivalent plant for
a particular frozen position of the motion system, denoted by $g_{\tt p}^{i}$,
corresponds to a lower fractional transformation (LFT):

$g_{\tt p}^{i}=\mathcal{F}(P_{\tt p}^{i},-K_{\tt p}^{i}),$ (8)

where $K_{\tt p}^{i}$ = diag(${k_{{\tt p}_{j}}}$), $j=\\{1,\ldots,n\\}$,
$j\neq i$. Moreover, the local diagonal feedback controller, which corresponds
to the local plant dynamics $P_{\tt p}$, is defined as $K_{\tt
p}=\text{diag}(k_{{\tt p}_{1}},\ldots,k_{{\tt p}_{n_{u}}})$, where local
stability is assessed using the characteristic equation det($I+P_{\tt p}K_{\tt
p}$) $\forall\tt p\in\mathbb{P}$, which, in case of a position dependent
controller, is related to the equivalent plant as:

$\mathrm{det}(I+P_{\tt p}K_{\tt p})=\prod_{i=1}^{n_{\mathrm{u}}}(1+g_{\tt
p}^{i}{k_{\tt p}}_{i})\ \forall\tt p\in\mathbb{P}$ (9)

For the design of the position dependent controllers $k_{{\tt p}_{i}}$, the
structured controller is partitioned into a low-frequent LTI controller,
denoted by $\Gamma_{i}$, and an interconnection of LPV filters, denoted by
$\Psi_{i}$, as illustrated by Figure 3. Typically, $\Gamma_{i}$ consists of a
proportional gain $C_{\mathrm{k}}^{i}$, an integrator $C_{\mathrm{I}}^{i}$,
see (5), and lead filters $C_{\mathrm{D}}^{i}$, see (6). The position
dependent cascade filter interconnection $\Psi_{i}$ consists of LPV notch
filters, which actively combat the effects of position dependent flexible
dynamics in a position dependent manner, therefore allowing for increased
rigid body feedback control bandwidth. The time-domain representation of a LPV
notch filter corresponds to:

Figure 3: Schematic representation of the proposed LPV SLC control design
concept, where $\Gamma_{i}$ denotes a cascade interconnection of LTI filters
and $\Psi_{i}$ corresponds to a cascade interconnection of LPV filters.

$\left[\begin{array}[]{c}\dot{x}_{\mathrm{N}}(t)\\\
y_{\mathrm{N}}(t)\end{array}\right]=\left[\begin{array}[]{c|c}\mathcal{A}(p(t))&\mathcal{B}(p(t))\\\
\hline\cr\mathcal{C}(p(t))&\mathcal{D}(p(t))\end{array}\right]\left[\begin{array}[]{c}x_{\mathrm{N}}(t)\\\
u_{\mathrm{N}}(t)\end{array}\right],$ (10)

with:

$\begin{split}\mathcal{A}(p(t))&=\begin{bmatrix}-4\pi\beta_{2}(p(t))f_{2}(p(t))&-4f_{2}(p(t))^{2}\pi^{2}\\\
1&0\end{bmatrix}\\\
\mathcal{B}(p(t))&=\begin{bmatrix}4f_{2}(p(t))^{2}\pi^{2}\\\ 0\end{bmatrix}\\\
\mathcal{C}(p(t))&=\begin{bmatrix}\frac{\beta_{1}(p(t))f_{1}(p(t))-\beta_{2}(p(t))f_{2}(p(t))}{f_{1}(p(t))^{2}\pi}&1-\frac{f_{2}(p(t))^{2}}{f_{1}(p(t))^{2}}\end{bmatrix}\\\
\mathcal{D}(p(t))&=\frac{f_{2}(p(t))^{2}}{f_{1}(p(t))^{2}}\end{split}$ (11)

Similar to the SLC control design approach, a set of _frozen dynamics_ ,
denoted by $\\{P_{\tt p}^{l}\\}_{l=1}^{n_{g}}$, is used to construct a set of
_frozen controllers_ $\\{K_{\tt p}^{l}\\}_{l=1}^{n_{g}}$ by utilizing loop
shaping techniques, thus allowing for an interpolation based parameterization
of the position dependent filter coefficients using the corresponding _frozen
scheduling vectors_ $\tt p\in\mathbb{P}$ respectively, where $n_{g}$ denotes
the number of _frozen positions_ that are considered for controller design.
The parameterization of the position dependent filters can be accomplished
using two approaches. First, a global parameterization of the filter
coefficients can be considered with polynomial dependence on the parameters.
Therefore, the polynomial dependence can be tuned using SLC control design
approaches , thus taking all the local dynamics into account during the tuning
of the position dependent filters. As an alternative, local parameterization
strategies can be applied, for which each individual filter coefficient is
interpolated based on a polynomial relationship which is dependent on the
scheduling vector. In this paper, the latter strategy is investigated.

Consider the vector $\Phi(p(t))$, which contains the position dependent notch
filter coefficients of (10):

$\small\Phi(p(t))=\begin{bmatrix}\beta_{1}(p(t))&\beta_{2}(p(t))&f_{1}(p(t))&f_{2}(p(t))\end{bmatrix}^{\top}$
(12)

Using the set of local controller designs $\\{K_{\tt p}^{l}\\}_{l=1}^{n_{g}}$,
the position dependent filter coefficients can be individually parameterized
using a polynomial dependence on the scheduling vector as:

$\Phi_{m}(p(t))=\sum_{v=0}^{(i-1)}\sum_{w=0}^{(j-1)}\theta_{vw}^{\Phi_{m}}q_{x}(t)^{v}q_{y}(t)^{w},$
(13)

where $\Phi_{m}(p(t))$ denotes the $m^{\mathrm{th}}$ element of the vector
$\Phi(p(t))$ and $q_{x},q_{y}\subseteq\tt p$. The parameters $i$ and $j$
denote the order of the assumed dependency, governing the smoothness of the
allowed variation. Additionally, from (13), it is observed that each parameter
$\Phi_{m}(p(t))$ is expressed by spatial coordinates (e.g. dependency on the
$q_{x}(t),q_{y}(t)$ position) and weighting coefficients
$\theta_{vw}^{\Phi_{m}}$, thus allowing for reformulation of the
parameterization as:

$\Phi_{m}(p(t))=\chi(p(t))\begin{bmatrix}\theta_{11}^{\Phi_{m}}&\ldots&\theta_{ij}^{\Phi_{m}}\end{bmatrix}^{\top},$
(14)

where the vector $\chi(p(t))\in\mathbb{R}^{i\cdot j}$ corresponds to:

$\chi(p(t))=\begin{bmatrix}1&\ldots&q_{x}(t)^{i-1}\end{bmatrix}\otimes\begin{bmatrix}1&\ldots&q_{y}(t)^{j-1}\end{bmatrix}$
(15)

In order to obtain the weighting coefficient vectors, the polynomial
parameterization, expressed by (14), is reformulated into a least-squares
regression problem using the designed local filter coefficients of $\\{K_{\tt
p}^{l}\\}_{l=1}^{n_{g}}$ and their corresponding _frozen positions_
respectively. Moreover, the least-squares formulation for the approximation of
$\Phi_{m}(p(t))$ is given by:

$\underbrace{\begin{bmatrix}{\Phi_{m}}^{1}\\\ \vdots\\\
{\Phi_{m}}^{n_{g}}\end{bmatrix}}_{Y_{\Phi_{m}}}=\underbrace{\begin{bmatrix}\chi({\tt
p}_{1})\\\ \vdots\\\ \chi({\tt
p}_{n_{g}})\end{bmatrix}}_{A}\underbrace{\begin{bmatrix}\theta_{11}^{\Phi_{m}}\\\
\vdots\\\ \theta_{ij}^{\Phi_{m}}\end{bmatrix}}_{\Theta_{\Phi_{m}}},$ (16)

where the coefficient vectors $\Theta_{\Phi_{m}}$ are obtained by minimizing
$||A\Theta_{\Phi_{m}}-Y_{\Phi_{m}}||^{2}$. Therefore, the parameterization of
the position dependent filter coefficients is based on minimizing the $l_{2}$
loss between the local coefficient tuning of $\\{K_{\tt
p}^{l}\\}_{l=1}^{n_{g}}$ and the fit of the polynomial parameterization.

Real-time implementation of the position dependent controller blocks is
achieved by using the controller structure given by (10), which allows for
implementation of the position dependent notch filter using an integrator
chain with position dependent coefficients as illustrated in Figure 4.

Figure 4: Implementation concept of a position dependent notch filter using an
integrator chain.

In total, a novel control design approach has been proposed, which allows for
position dependent sequential loop closing based control design. In this way,
position dependent phenomena in high-precision motion systems can be
explicitly handled, allowing for increased closed-loop performance.
Additionally, the structure of the position dependent filters allows for
efficient implementation of the filters using an integrator chain, where the
total feedback controller $k_{i}$ is obtained by interconnecting the LTI
filters $\Gamma_{i}$ with the LPV filters $\Psi_{i}$.

##  V Simulation study

In this Section, a simulation study of the LPV SLC control design approach is
presented. First Subsection V-A presents a brief description of a moving-
magnet planar actuator system. Next, simulation results are presented in
Subsection V-B, where the proposed LPV SLC control design approach is compared
to conventional SLC control design strategies on a state-of-the-art moving-
magnet planar actuator model, which exhibits position dependency in both
actuation and sensing.

###  V-A System description

Magnetically levitated planar actuators are high-precision motion systems
which are used in lithographic machines for nanometer accurate positioning of
silicon wafers under projection optics, see [1]. A moving-magnet planar
actuator system, which is illustrated by Figure 5, is comprised of three main
components: the stator base, the translator and the metrology frame. The
stator base is a double layer coil array, consisting of 160 coils of which 40
coils are simultaneously activated at every time instant using 40 power
converters, depending on the relative position of the translator (see [15]).
Proper actuation of the coils offers the means of both levitation and
propulsion of the magnet plate in 6 Degrees of Freedom (DoF). The translator,
comprised of 281 permanent magnets structured in a Hallbach array, is
constructed to be a lightweight magnetic plate, thus enabling high
accelerations. However, due to its low mass, low-frequent flexible dynamics
are introduced (first flexible mode at $\approx$ 226.5 Hz), which severely
limit position tracking performance capabilities. The metrology frame, which
rests on air mounts to suppress the effects of floor disturbances, is used as
a global reference frame in order to assess positioning accuracy of the mover.
On the metrology frame, 9 laser interferometers (LIFMs) are mounted to measure
the relative displacement of the translator with respect to the metrology
frame. The relative actuation and sensing of the moving-body introduces the
presence of position dependent flexible dynamics, which severely limit the
achievable rigid body feedback control bandwidth, see Figure 1. A detailed
description of the moving-magnet planar actuator prototype is given in [16].

Figure 5: Photograph of a moving-magnet planar actuator system prototype.

###  V-B Simulation results

For validation of the proposed control design approach, both a conventional
SLC feedback controller and a LPV SLC controller are designed for the $z$-axis
of the moving-magnet planar actuator model, which is severely limited in
achieving sufficient rigid body feedback control bandwidth due to the presence
of position dependent flexible dynamics as illustrated in Figure 1. In order
to compare the performance of both controllers with respect to each other, a
performance criteria is introduced by means of a 6 dB upper-bound on the
sensitivity function of the closed-loop. Additionally, both controllers are
constructed from a proportional gain, an integrator, lead filters and
(position dependent) notch filters. Using the aforementioned filters, the
conventional SLC feedback controller achieves a rigid body feedback control
bandwidth of approximately 95 Hz. The LPV SLC feedback controller achieves a
feedback control bandwidth of approximately 175 Hz due to the additional
degrees of freedom during controller design. The introduction of position
dependent notch filters to the controller design has several advantageous
properties. First, the effects of position dependent flexible dynamics can be
explicitly handled in a position dependent manner, thus allowing to push
performance beyond robust LTI SLC control design. Secondly, the position
dependent notch filters allow for relaxation of the feedback control design
constraints, since phase lead can be introduced in a position dependent manner
by means of skewing the notch filter (shifting the notch filter frequency
$f_{2}$, see (10), to a higher frequency compared to $f_{1}$), thus allowing
for increased feedback control bandwidth, while still satisfying the 6 dB
upper-bound of the sensitivity function respectively. The constructed position
dependent notch filter is illustrated in Figure 6 for the positions $q_{x}$ =
0.1 $m$, $q_{y}\in[0m\quad 0.2m]$.

Figure 6: Visualization of a designed position dependent notch filter for the
positions $q_{x}$ = 0.1$m$ and $q_{y}\in[0m\quad 0.2m]$.

From Figure 6 it is observed that the constructed position dependent notch
filter allows for introduction of phase lead for critical positions of the
system due to the additional design freedom of the position dependent filters,
thus allowing for increased closed-loop performance compared to conventional
SLC control design approaches. Additionally, the designed notch filter is
suppressing the effects of the resonance dynamics. In order to evaluate the
performance improvement of the proposed LPV SLC control design approach, both
the conventional SLC controller and the LPV SLC controller are implemented in
simulation together with a mass feedforward. For simulation purposes, a 4th
order reference trajectory is generated based on the methodology proposed by
[17]. Figure 7 illustrates the motion profiles that are used for simulation.

Figure 7: Motion profiles that are considered for simulation of the
conventional SLC controller and the LPV SLC controller.

Figure 8: Position tracking error in $z$-direction, where the green area
denotes the (de-)acceleration interval, the yellow area denotes the settling
time interval and the red area denotes the constant velocity interval.

The resulting tracking error is illustrated in Figure 8, where the position
tracking error using the conventional SLC controller is depicted in blue and
the position tracking error using the LPV SLC controller is depicted in red.
Furthermore, the green area denotes the acceleration and de-acceleration
interval, the yellow area denotes the settling time interval and the red area
denotes the constant velocity interval.

In industry, accuracy of the lithographic exposure process, which takes place
during the constant velocity interval, is assessed by two values: the Moving
Average ($\mathrm{MA}(t)$) and the Moving Standard Deviation
($\mathrm{MSD}(t)$), see [1]. These values are calculated a posteriori from
the measured position tracking error using:

$\begin{split}\mathrm{MA}(t)&=\frac{1}{T}\int_{t-\frac{T}{2}}^{t+\frac{T}{2}}e(\tau)d\tau\\\
\mathrm{MSD}(t)&=\sqrt{\frac{1}{T}\int_{t-\frac{T}{2}}^{t+\frac{T}{2}}(e(\tau)-\mathrm{MA}(t))^{2}d\tau}\end{split},$
(17)

where $T$ corresponds to the exposure window. The mean values of the
$\mathrm{MA(t)}$ and $\mathrm{MSD(t)}$ filtered error signals during the
constant velocity interval are listed in Table I together with the relative
reduction of the performance values.

TABLE I: Mean values of $\mathrm{MA(t)}$ and $\mathrm{MSD(t)}$ for both controllers during the constant velocity interval | Conventional SLC | LPV SLC | Relative reduction
---|---|---|---
MA [m] | 5.3139e-10 | 4.4360e-11 | 91.65 %
MSD [m] | 4.8913e-08 | 4.0785e-08 | 16.62 %

Based on Table I, it is clear that the LPV SLC design based controller has a
relative reduction of the mean $\mathrm{MA(t)}$ error of 91.65 % and a
relative reduction of the mean $\mathrm{MSD}(t)$ error of 16.62 %, thus
highlighting the effectiveness of the proposed control design approach. From
Figure 8 it is difficult to observe the improvement in $\mathrm{MA}(t)$ error
as the low-frequency behavior is obscured by the resonance. Moreover, an
important observation is that the LPV SLC controller is better equipped to
handle the position dependent resonance, which leads to increased feedback
control bandwidth. Nonetheless, the resonant behavior ($\mathrm{MSD}$ error)
is hardly reduced as it is excited by the mass feedforward.

##  VI Conclusions

This paper presents a novel extension of the conventional SLC control design
framework towards a linear-parameter-varying approach, which allows to
circumvent limitations that are introduced by position dependent effects of
high-precision motion systems. Moreover, the proposed framework introduces
additional degrees of freedom to the controller design by means of exploiting
position dependent filter design, which allows for increased closed-loop
performance of high-precision motion systems. Simulation results based on a
state-of-the-art moving-magnet planar actuator model highlight the
effectiveness of the proposed control design approach.

## References

  * [1] H. Butler, “Position control in lithographic equipment [applications of control],” IEEE Control Systems Magazine, vol. 31, no. 5, pp. 28–47, 2011\.
  * [2] N. Tamer and M. Dahleh, “Feedback control of piezoelectric tube scanners,” in Proceedings of 1994 33rd IEEE Conference on Decision and Control, vol. 2, pp. 1826–1831 vol.2, 1994.
  * [3] M. Heertjes, “Data-based motion control of wafer scanners,” IFAC-PapersOnLine, vol. 49, no. 13, pp. 1–12, 2016.  12th IFAC Workshop on Adaptation and Learning in Control and Signal Processing ALCOSP 2016.
  * [4] X. Ye, Y. Zhang, and Y. Sun, “Robotic pick-place of nanowires for electromechanical characterization,” in 2012 IEEE International Conference on Robotics and Automation, pp. 2755–2760, 2012.
  * [5] T. Oomen, “Advanced motion control for precision mechatronics: control, identification, and learning of complex systems,” IEEJ Journal of Industry Applications, vol. 7, pp. 127–140, Jan. 2018.
  * [6] M. Steinbuch, “Design and control of high tech systems,” in 2013 IEEE International Conference on Mechatronics (ICM), pp. 13–17, 2013.
  * [7] M. Steinbuch, R. Merry, M. Boerlage, M. Ronde, and M. Molengraft, van de, Advanced Motion Control Design, pp. 27–1/25.  CRC Press, 2010.
  * [8] M. Hovd and S. Skogestad, “Sequential design of decentralized controllers,” Automatica, vol. 30, no. 10, pp. 1601–1607, 1994.
  * [9] P. Hughes, “Space structure vibration modes: How many exist? which ones are important?,” IEEE Control Systems Magazine, vol. 7, no. 1, pp. 22–28, 1987\.
  * [10] M. Steinbuch and M. Norg, “Advanced motion control: An industrial perspective,” European Journal of Control, vol. 4, no. 4, pp. 278–293, 1998.
  * [11] S. Skogestad and I. Postlethwaite, Multivariable feedback control: analysis and design, vol. 2.  Citeseer, 2007.
  * [12] R. Toth, H. S. Abbas, and H. Werner, “On the state-space realization of lpv input-output models: Practical approaches,” IEEE Transactions on Control Systems Technology, vol. 20, no. 1, pp. 139–153, 2012.
  * [13] T. Oomen and M. Steinbuch, “Model-based control for high-tech mechatronic systems,” in Mechatronics and Robotics, pp. 51–80, CRC Press, 2020.
  * [14] D. Q. Mayne, “Sequential design of linear multivariable systems,” in Proceedings of the Institution of Electrical Engineers, vol. 126, pp. 568–572, IET, 1979.
  * [15] C. M. M. van Lierop, Magnetically levitated planar actuator with moving magnets: Electromechanical analysis and design II.  PhD thesis, Eindhoven University of Technology, 2007.
  * [16] C. Custers, I. Proimadis, J. W. Jansen, H. Butler, R. Tóth, E. Lomonova, and P. van den Hof, “Active compensation of the deformation of a magnetically levitated mover of a planar motor,” in 2019 IEEE International Electric Machines Drives Conference (IEMDC), pp. 854–861, 2019.
  * [17] P. Lambrechts, M. Boerlage, and M. Steinbuch, “Trajectory planning and feedforward design for electromechanical motion systems,” Control Engineering Practice, vol. 13, no. 2, pp. 145–157, 2005.