# The Adaptive Spectral Koopman Method for Dynamical Systems

Bian Li Department of Industrial and Systems Engineering, Lehigh University,
Bethlehem, PA<EMAIL_ADDRESS><EMAIL_ADDRESS>Yi-An Ma Halicioğlu Data
Science Institute & Department of Computer Science and Engineering, University
of California San Diego, San Diego, CA<EMAIL_ADDRESS>J. Nathan Kutz
Department of Applied Mathematics, University of Washington, Seattle, WA
<EMAIL_ADDRESS>Xiu Yang11footnotemark: 1

###### Abstract

Dynamical systems have a wide range of applications in mechanics, electrical
engineering, chemistry, and so on. In this work, we propose the adaptive
spectral Koopman (ASK) method to solve nonlinear autonomous dynamical systems.
This novel numerical method leverages the spectral-collocation (i.e., pseudo-
spectral) method and properties of the Koopman operator to obtain the solution
of a dynamical system. Specifically, this solution is represented as a linear
combination of the multiplication of Koopman operator’s eigenfunctions and
eigenvalues, and these eigenpairs are approximated using the spectral method.
Unlike conventional time evolution algorithms such as Euler’s scheme and the
Runge-Kutta scheme, ASK is mesh-free, and hence is more flexible when
evaluating the solution. Numerical experiments demonstrate high accuracy of
ASK for solving one-, two- and three-dimensional dynamical systems.

Keywords:  dynamical systems, Koopman operator, spectral-collocation method

##  1 Introduction

The Koopman operator, introduced in 1931 by B. O. Koopman [11], is an
infinite-dimensional _linear_ operator that describes the evolution of a set
of observables rather than the system state itself. The Koopman operator
approach to _nonlinear_ dynamical systems has attracted considerable attention
recently, as it provides a rigorous method for globally linearizing the system
dynamics. Specifically, because it is a linear operator, one can define its
eigenvalues, eigenfunctions, and modes, and use them to represent dynamically
interpretable low-dimensional embeddings of high-dimensional state spaces,
which helps to understand the behavior of the underlying system and construct
solutions through linear superposition [4]. In this procedure, the system
dynamics is typically decomposed into linearly independent Koopman modes even
if the system is nonlinear. In particular, as pointed out in [18, 12, 19], if
the dynamics is ergodic but non-chaotic, the spectrum of the Koopman operator
in properly defined spaces does not contain continuous spectra, and the
observable of the system can be represented as a linear combination of
eigenfunctions associated with discrete eigenvalues of the Koopman operator.

The Koopman operator provides powerful analytic tools to understand behaviors
of dynamical systems. For example, dynamical evolution of a finite-dimensional
system described by ordinary differential equations (ODEs) can be studied by
conducting Koopman mode analysis. Such analysis starts with a choice of a set
of linearly independent observables, and the Koopman operator is then analyzed
through its action on the subspace spanned by the chosen observables [17].
Moreover, it is also shown that the Koopman operator approach can be formally
generalized to infinite-dimensional dynamical systems described by partial
differential equations (PDEs), providing new perspectives on the analysis and
control of these nonlinear spatiotemporal dynamics [35, 20, 22, 19]. In
addition, ergodic quotients and eigenquotients allow the Koopman operator to
be used for the extraction and analysis of _invariant_ and _periodic_
structures in the state space [5]. Moreover, Mezić provided a Hilbert space
setting for spectral analysis of _disspative dynamical systems_ , and proved
that the spectrum of the Koopman operator on these spaces is the closure of
the product of the “on-attractor” and “off-attractor” spectra [18].

On the computational side, most existing numerical schemes motivated by the
Koopman operator are categorized as data-driven methods, as they use
spatiotemporal data to approximate a few of the leading Koopman eigenvalues,
eigenfunctions, and modes. In particular, the emerging computational method
dynamics mode decomposition (DMD) [24, 25, 30, 23, 14, 20, 1] as well as its
variant such as extended DMD (EDMD) [34] uses snapshots of a dynamical system
to extract temporal features as well as correlated spatial activity via matrix
decomposition techniques. DMD and EDMD produce results for any appropriately
formatted set of data, but connecting these outputs to the Koopman operator
requires additional knowledge about the nature of the underlying system in
that the system should be autonomous. Later, a modified EDMD [33] was proposed
to compensate for the effects of system actuation when it is used to explore
state space during the data collection, reestablishing the connection between
EDMD and the Koopman operator in this more general class of data sets. A
review of many of the DMD variants for approximating the Koopman operator can
be found in Brunton et al [4].

Our aim in this paper is to provide a numerical method based on the spectral-
collocation method (i.e., the pseudospectral method) to implement the Koopman-
operator approach to solving nonlinear ordinary differential equations (ODEs).
Unlike the data-driven methods, this approach is on the other end of the
“spectrum” of numerical methods, as it is based on the classical spectral
method [8, 29]. The main idea is to approximate eigenvalues, eigenfunctions,
and modes of the Koopman operator based on its discretized form. Specifically,
this method uses the differentiation matrix in spectral method to approximate
the generator of the Koopman operator, and then conducts eigendecomposition
numerically to obtain eigenvalues and eigenvectors that approximate Koopman
operator’s eigenvalues and eigenfunctions, respectively. Here, each element of
an eigenvector is the approximation of the associated eigenfunction evaluated
at a collocation point. The modes are approximated using the computed
eigenvalues, eigenvectors, and the initial state (or observable). This work
focuses on autonomous systems, and it would serve as a starting point for a
new framework of numerical methods for dynamical systems.

The paper is organized as follows. Background topics are introduced in section
2. Then, the adaptive spectral Koopman method is discussed in detail in
section 3. We present the experimental results in section 4, and the
discussion and conclusions follow in section 5.

##  2 Background

###  2.1 Koopman operator

Borrowing notions from [13], we consider an autonomous system described by the
ordinary differential equations

$\frac{\mathrm{d}\bm{x}}{\mathrm{d}t}=\operatorname{\mathbf{f}}(\bm{x}),$ (1)

where the state $\bm{x}=(x_{1},x_{2},\dotsc,x_{d})^{\top}$ belongs to a
$d$-dimensional smooth manifold $\mathscr{M}$, and the dynamics
$\operatorname{\mathbf{f}}:\mathscr{M}\rightarrow\mathscr{M}$ does not
explicitly depend on time $t$. Here, $\operatorname{\mathbf{f}}$ is a possibly
nonlinear vector-valued smooth function, of the same dimension as $\bm{x}$. In
many studies, we are concerned with the behavior of observables on the state
space. To this end, we define an observable to be a scalar function
$g:\mathscr{M}\rightarrow\mathbb{R}$, where $g$ is an element of some function
space $\mathcal{G}$ (e.g., $\mathcal{G}=L^{2}(\mathscr{M})$ as in [17]). The
flow map $\operatorname{\mathbf{F}}_{t}:\mathscr{M}\rightarrow\mathscr{M}$
induced by the dynamical system (1) depicts the evolution of the system as

$\bm{x}(t_{0}+t)=\operatorname{\mathbf{F}}_{t}(\bm{x}(t_{0}))=\bm{x}(t_{0})+\int_{t_{0}}^{t_{0}+t}\operatorname{\mathbf{f}}(\bm{x}(s))\,ds.$
(2)

Now we define the Koopman operator for continuous-time dynamical systems as
follows [18]:

######  Definition 2.1.

Consider a family of operators $\\{\mathcal{K}_{t}\\}_{t\geq 0}$ acting on the
space of observables so that

$\mathcal{K}_{t}g(\bm{x}_{0})=g(\mathbf{F}_{t}(\bm{x}_{0})),$

where $\bm{x}_{0}=\bm{x}(t_{0})$. We call the family of operators
$\mathcal{K}_{t}$ indexed by time t the Koopman operators of the continuous-
time system (1).

By definition, $\mathcal{K}_{t}$ is a _linear_ operator acting on the function
space $\mathcal{G}$ for each fixed $t$. Moreover, $\\{\mathcal{K}_{t}\\}$ form
a semi-group.

###  2.2 Infinitesimal generator

The Koopman spectral theory [17, 24] reveals properties that enable the
Koopman operator to convert nonlinear finite-dimensional dynamics into linear
infinite-dimensional dynamics. A key component in such spectral analysis is
the initesimal generator (or generator for brievity) of the Koopman operator.
Specifically, the generator of the Koopman operator $\mathcal{K}_{t}$, denoted
as $\mathcal{K}$, is given by

$\displaystyle\mathcal{K}g=\lim_{t\rightarrow 0}\frac{\mathcal{K}_{t}g-g}{t}.$
(3)

For any smooth function $g$, eq. 3 implies that

$\displaystyle\mathcal{K}g(\bm{x})=\frac{\mathrm{d}g(\bm{x})}{\mathrm{d}t}=\nabla
g(\bm{x})\cdot\frac{\mathrm{d}\bm{x}}{\mathrm{d}t}.$ (4)

Denoting $\operatorname{\varphi}$ an eigenfunction of $\mathcal{K}$ and
$\lambda$ the eigenvalue associated with $\operatorname{\varphi}$, we have

$\displaystyle\mathcal{K}\operatorname{\varphi}(\bm{x})=\lambda\operatorname{\varphi}(\bm{x}).$
(5)

Thus,

$\displaystyle\lambda\operatorname{\varphi}(\bm{x})=\mathcal{K}\operatorname{\varphi}(\bm{x})=\frac{\mathrm{d}\operatorname{\varphi}(\bm{x})}{\mathrm{d}t}.$
(6)

This implies that $\varphi(\bm{x}(t_{0}+t))=e^{\lambda
t}\varphi(\bm{x}(t_{0}))$, i.e.,

$\mathcal{K}_{t}\varphi(\bm{x}(t_{0}))=e^{\lambda t}\varphi(\bm{x}(t_{0})).$
(7)

Therefor, $\varphi$ is an eigenfunction of $\mathcal{K}_{t}$ associated with
eigenvalue $\lambda$. Of note, following the conventional notation, the
eigenpair for $\mathcal{K}_{t}$ is considered as $(\varphi,\lambda)$ instead
of $(\varphi,e^{\lambda t})$.

Now suppose $g$ exists in the function space spanned by all the eigenfunctions
$\operatorname{\varphi}_{j}$ (associated with eigenvalues $\lambda_{j}$) of
$\mathcal{K}$, i.e.,
$g(\bm{x})=\sum_{j}c_{j}\operatorname{\varphi}_{j}(\bm{x})$, then

$\displaystyle\mathcal{K}_{t}[g(\bm{x}(t_{0}))]=\mathcal{K}_{t}\left[\sum_{j}c_{j}\operatorname{\varphi}_{j}(\bm{x}(t_{0}))\right]=\sum_{j}c_{j}\mathcal{K}_{t}[\operatorname{\varphi}_{j}(\bm{x}(t_{0}))].$
(8)

Hence,

$g(\bm{x}(t_{0}+t))=\sum_{j}c_{j}\operatorname{\varphi}_{j}(\bm{x}(t_{0}))e^{\lambda_{j}t}.$
(9)

Similarly, if we choose a vector-valued observable
$\bm{g}:\mathscr{M}\rightarrow\operatorname{\mathbb{R}}^{d}$ with
$\bm{g}:=(g_{1}(\bm{x}),g_{2}(\bm{x}),\dotsc,g_{d}(\bm{x}))^{\top}$, the
system of observables becomes

$\displaystyle\frac{\mathrm{d}\bm{g}(\bm{x})}{\mathrm{d}t}=\mathcal{K}\bm{g}(\bm{x})=\begin{bmatrix}\mathcal{K}g_{1}(\bm{x})\\\
\mathcal{K}g_{2}(\bm{x})\\\ \vdots\\\ \mathcal{K}g_{d}(\bm{x})\\\
\end{bmatrix}=\sum_{j}\lambda_{j}\operatorname{\varphi}_{j}(\bm{x})\bm{c}_{j},$
(10)

where $\bm{c}_{j}\in\operatorname{\mathbb{C}}^{d}$ is called the $j$th Koopman
mode with $\bm{c}_{j}:=(c_{j}^{1},c_{j}^{2},\dotsc,c_{j}^{d})^{\top}$. In
general, there is no universal guide for choosing observables as this choice
is problem dependent. A good set of observables can lead to a system that is
significantly easier to solve. An example from [3, 16] is illustrated in
appendix A.

We finalize the introduction of the Koopam operator with the following simple
example. Consider the system $\frac{\mathrm{d}x}{\mathrm{d}t}=\mu x$ with
$x,\mu\in\mathbb{C}$ and $|\mu|\neq 1$. Then one can easily verify that
$\varphi_{n}(x):=x^{n}$ is an eigenfunction of the Koopman operator associated
with this dynamical system, and the corresponding eigenvalue is
$\lambda_{n}=\mu^{n}$ with $n\in\mathbb{N}^{+}$ (see Example 6 in [6]).
According to eq. 9, by setting $g(x)=x$ and let $x(0)=x_{0}$, we have

$x(t)=\sum_{j=1}^{\infty}c_{j}\operatorname{\varphi}_{j}(x_{0})e^{\lambda_{j}t}=\sum_{j=1}^{\infty}c_{j}x_{0}^{j}e^{\mu^{j}t}.$

We set $t=0$ to have $x_{0}=x(0)=\sum_{j=1}^{\infty}c_{j}x_{0}^{j}$, which
indicates $c_{1}=1$ and $c_{j}=0$ when $j\neq 1$. Therefore, we obtain the
solution of the ODE as $x(t)=x_{0}e^{\mu t}$.

##  3 Adaptive Spectral Koopman Method

In this section, we introduce the adaptive spectral Koopman (ASK) method,
which is a numerical method based on the Koopman operator and the spectral
method to solve ODE systems. Before describing details of this method, we
introduce the notations used in this algorithm. We consider solution of an ODE
system in time interval $[0,T]$ with $T>0$, i.e., we assume $t_{0}=0$ in eq.
2. The initial condition is $\bm{x}(t_{0})=\bm{x}_{0}$. We use $n$ to denote
the number of “check points” (see details in section 3.4). We denote the
radius of the neighborhood of $\bm{x}(t)$ by $r$, and use $\gamma$ as a
parameter to control the update of the neighborhood.

###  3.1 Finite-dimensional approximation

Based on the preliminaries introduced in section 2.2, the we aim to identify
the following truncated form of eq. 9

$g(\bm{x}(t))\approx
g_{N}(\bm{x}(t))=\sum_{j=0}^{N}c_{j}\operatorname{\varphi}_{j}(\bm{x}_{0})e^{\lambda_{j}t}.$
(11)

via approximating eigenpairs $\\{\varphi_{j},\lambda_{j}\\}_{j=0}^{N}$ and
computing modes $c_{j}$, where $N$ is a positive integer. For demonstration
purpose, we set the observable as the identity function (i.e.,
$g(\bm{x})=\bm{x}$) in this work, and the proposed method can be directly
implemented on other observables.

Because
$\frac{\mathrm{d}\bm{x}}{\mathrm{d}t}=\operatorname{\mathbf{f}}(\bm{x})$, we
obtain from eq. 4 and eq. 6:

$\displaystyle\mathcal{K}\operatorname{\varphi}$
$\displaystyle=\operatorname{\mathbf{f}}\cdot\nabla\operatorname{\varphi}=\left(f_{1}\frac{\partial\operatorname{\varphi}}{\partial
x_{1}}+f_{2}\frac{\partial\operatorname{\varphi}}{\partial
x_{2}}+...+f_{d}\frac{\partial\operatorname{\varphi}}{\partial x_{d}}\right)$
$\displaystyle=\left(f_{1}\frac{\partial}{\partial
x_{1}}+f_{2}\frac{\partial}{\partial x_{2}}+...+f_{d}\frac{\partial}{\partial
x_{d}}\right)(\operatorname{\varphi}).$

Thus,

$\displaystyle\mathcal{K}=f_{1}\frac{\partial}{\partial
x_{1}}+f_{2}\frac{\partial}{\partial x_{2}}+...+f_{d}\frac{\partial}{\partial
x_{d}}.$ (12)

Here we consider the case with $d\leq 3$, and adopt the approaches in the
spectral-collocation method. Specifically, we use Gauss-Lobatto points for the
interpolation of $\varphi$ and use differentiation matrices to approximate
(partial) derivatives (see e.g., [10, 26, 9]). Consequently, the first step of
our algorithm is to discretize $\mathcal{K}$.

  1. (1)

When $d=1$. Let $\\{\xi_{i}\\}_{i=0}^{N}$ be the Gauss-Lobatto points and the
Lagrange interpolation of $\varphi(\bm{x})$ is

$\varphi(x)\approx\sum_{j=0}^{N}\varphi(\xi_{i})P_{j}(x),$

where $P_{j}$ are Lagrange polynomials satisfying $P_{j}(\xi_{i})=\delta_{ij}$
and $\delta_{ij}$ is the Kronecker delta function. Let
$\bm{\varphi}=[\varphi(\xi_{0}),\varphi(\xi_{1}),\dotsc,\varphi(\xi_{N})]^{\top}$,
we have

$\displaystyle\mathcal{K}\bm{\varphi}=\operatorname{diag}(\operatorname{\mathbf{f}}(\xi_{0}),\operatorname{\mathbf{f}}(\xi_{1}),\dotsc,\operatorname{\mathbf{f}}(\xi_{N}))\mathbf{D}\bm{\varphi}:=\mathbf{K}\bm{\varphi},$
(13)

where $\mathbf{D}$ is the differentiation matrix associated with
$\\{\xi_{i}\\}_{i=0}^{N}$ and $\mathbf{K}$ is an $(N+1)\times(N+1)$ matrix.
Here, we abuse the notation to use $\mathcal{K}\bm{\varphi}$ to denote
evaluating $\mathcal{K}\varphi$ at $\\{\xi_{i}\\}_{i=0}^{N}$, and similar
notations are used in the following $d=2,3$ cases.

  2. (2)

When $d=2$. Let $\\{\xi_{i}\\}_{i=0}^{N}$ and $\\{\eta_{j}\\}_{j=0}^{N}$ be
the Gauss-Lobatto points of $x_{1}$ and $x_{2}$, respectively. The
eigenfunction $\operatorname{\varphi}$ is now a bivariate function, whose
Lagrange interpolation is based on its values on collocation points
$(\xi_{i},\eta_{j})$. Hence, we define a matrix $\mathbf{\Phi}$ as

$\displaystyle\mathbf{\Phi}=\begin{bmatrix}\operatorname{\varphi}(\xi_{0},\eta_{0})&\operatorname{\varphi}(\xi_{0},\eta_{1})&\dots&\operatorname{\varphi}(\xi_{0},\eta_{N})\\\
\operatorname{\varphi}(\xi_{1},\eta_{0})&\operatorname{\varphi}(\xi_{1},\eta_{1})&\dots&\operatorname{\varphi}(\xi_{1},\eta_{N})\\\
\vdots&\vdots&\ddots&\vdots\\\
\operatorname{\varphi}(\xi_{N},\eta_{0})&\operatorname{\varphi}(\xi_{N},\eta_{1})&\dots&\operatorname{\varphi}(\xi_{N},\eta_{N})\\\
\end{bmatrix}.$

We denote the differentiation matrices for $x_{1}$ and $x_{2}$ by
$\mathbf{D}_{1}$ and $\mathbf{D}_{2}$, respectively. Let $\mathbf{F}_{1}$ and
$\mathbf{F}_{2}$ be the matrices of $f_{1}$ and $f_{2}$ evaluated at
$(\xi_{i},\eta_{j})$. Then, the values of $\mathcal{K}\varphi$ at the
collocation points can be approximated as

$\displaystyle\mathcal{K}\mathbf{\Phi}=\mathbf{F}_{1}\odot\left(\mathbf{D}_{1}\mathbf{\Phi}\right)+\mathbf{F}_{2}\odot\left(\mathbf{\Phi}\mathbf{D}_{2}^{\top}\right),$

where $\odot$ denotes the Hadamard product. Then we vectorize $\mathbf{\Phi}$
(along columns) to obtain

$\displaystyle\mathcal{K}\operatorname{vec}(\mathbf{\Phi})$
$\displaystyle=\operatorname{vec}(\mathbf{F}_{1})\odot\Big{(}(\mathbf{I}\otimes\mathbf{D}_{1})\operatorname{vec}(\mathbf{\Phi})\Big{)}+\operatorname{vec}(\mathbf{F}_{2})\odot\Big{(}(\mathbf{D}_{2}\otimes\mathbf{I})\operatorname{vec}(\mathbf{\Phi})\Big{)}$
$\displaystyle=\Big{[}\operatorname{diag}\big{(}\operatorname{vec}(\mathbf{F}_{1})\big{)}(\mathbf{I}\otimes\mathbf{D}_{1})+\operatorname{diag}\big{(}\operatorname{vec}(\mathbf{F}_{2})\big{)}(\mathbf{D}_{2}\otimes\mathbf{I})\Big{]}\Big{(}\operatorname{vec}(\mathbf{\Phi})\Big{)}$
$\displaystyle:=\mathbf{K}\operatorname{vec}(\mathbf{\Phi}),$

where $\otimes$ denotes the Kronecker product, $\mathbf{I}$ is the identity
matrix, and $\mathbf{K}$ is an $(N+1)^{2}\times(N+1)^{2}$ matrix.

  3. (3)

When $d=3$. Let $\\{\xi_{i}\\}_{i=0}^{N}$, $\\{\eta_{j}\\}_{j=0}^{N}$, and
$\\{\zeta_{k}\\}_{k=0}^{N}$ be the Gauss-Lobatto points of $x_{1}$, $x_{2}$,
and $x_{3}$, respectively. The collocation points are then
$(\xi_{i},\eta_{j},\zeta_{k})$. In this case, $\operatorname{\varphi}$ is a
trivariate function. Hence, the values of $\mathcal{K}\varphi$ at the
collocation points are approximated by a tensor $\mathbf{\Phi}$ whose frontal
slices can be written as

$\displaystyle\bm{\Phi}(:,:,k)=\begin{bmatrix}\operatorname{\varphi}(\xi_{0},\eta_{0},\zeta_{k})&\operatorname{\varphi}(\xi_{0},\eta_{1},\zeta_{k})&\dots&\operatorname{\varphi}(\xi_{0},\eta_{N},\zeta_{k})\\\
\operatorname{\varphi}(\xi_{1},\eta_{0},\zeta_{k})&\operatorname{\varphi}(\xi_{1},\eta_{1},\zeta_{k})&\dots&\operatorname{\varphi}(\xi_{1},\eta_{N},\zeta_{k})\\\
\vdots&\vdots&\ddots&\vdots\\\
\operatorname{\varphi}(\xi_{N},\eta_{0},\zeta_{k})&\operatorname{\varphi}(\xi_{N},\eta_{1},\zeta_{k})&\dots&\operatorname{\varphi}(\xi_{N},\eta_{N},\zeta_{k})\\\
\end{bmatrix}.$

With the n-mode multiplication in tensor algebra, we arrive at a compact
representation of the approximation,

$\displaystyle\mathcal{K}\mathbf{\Phi}=\operatorname{\mathbf{F}}_{1}\odot\Big{(}\mathbf{\Phi}\times_{1}\mathbf{D}_{1}\Big{)}+\operatorname{\mathbf{F}}_{2}\odot\Big{(}\bm{\Phi}\times_{2}D_{2}\Big{)}+\operatorname{\mathbf{F}}_{3}\odot\Big{(}\mathbf{\Phi}\times_{3}\mathbf{D}_{3}\Big{)},$

where $\times_{p}$ denotes the mode-$p$ tensor-matrix multiplication. Here,
$\mathbf{D}_{1},\mathbf{D}_{2},\mathbf{D}_{3}$ are the differentiation
matrices, and
$\operatorname{\mathbf{F}}_{1},\operatorname{\mathbf{F}}_{2},\operatorname{\mathbf{F}}_{3}$
denote the tensors resulting from $f_{1},f_{2},f_{3}$ evaluated at
$(\xi_{i},\eta_{j},\zeta_{k})$. Following the same idea of vectorization, we
rewrite the tensor representation as

$\displaystyle\mathcal{K}\operatorname{vec}(\mathbf{\Phi})$
$\displaystyle=\operatorname{vec}(\operatorname{\mathbf{F}}_{1})\odot\Big{(}(\mathbf{I}\otimes\mathbf{I}\otimes\mathbf{D}_{1})\operatorname{vec}(\mathbf{\Phi})\Big{)}$
$\displaystyle+\operatorname{vec}(\operatorname{\mathbf{F}}_{2})\odot\Big{(}(\mathbf{I}\otimes\mathbf{D}_{2}\otimes\mathbf{I})\operatorname{vec}(\mathbf{\Phi})\Big{)}$
$\displaystyle+\operatorname{vec}(\operatorname{\mathbf{F}}_{3})\odot\Big{(}(\mathbf{D}_{3}\otimes\mathbf{I}\otimes\mathbf{I})\operatorname{vec}(\mathbf{\Phi})\Big{)}$
$\displaystyle=\Big{[}\operatorname{diag}\big{(}\operatorname{vec}(\operatorname{\mathbf{F}}_{1})\big{)}(\mathbf{I}\otimes\mathbf{I}\otimes\mathbf{D}_{1})$
$\displaystyle+\operatorname{diag}\big{(}\operatorname{vec}(\operatorname{\mathbf{F}}_{2})\big{)}(\mathbf{I}\otimes\mathbf{D}_{2}\otimes\mathbf{I})$
$\displaystyle+\operatorname{diag}\big{(}\operatorname{vec}(\operatorname{\mathbf{F}}_{3})\big{)}(\mathbf{D}_{3}\otimes\mathbf{I}\otimes\mathbf{I})\Big{]}\Big{(}\operatorname{vec}(\mathbf{\Phi})\Big{)}$
$\displaystyle:=\mathbf{K}\operatorname{vec}(\mathbf{\Phi}).$

where $\mathbf{K}$ is an $(N+1)^{3}\times(N+1)^{3}$ matrix.

In all these cases, the discretized generator $\mathcal{K}$ can be represented
as a matrix $\mathbf{K}$. We note that for $d=2$ and $d=3$, the total number
of eigenfuncitons used in eq. 11 is $(N+1)^{2}$ and $(N+1)^{3}$, respectively,
instead of $(N+1)$. For brevity, we still use $g_{N}$ to denote the
approximated observable for different $d$. The derivation of higher
dimensional systems amounts to further extensions of the three-dimensional
case by the Kronecker product.

###  3.2 Eigen-decomposition

Now the eigenvalue problem of the Koopman operator in eq. 5 is discretized as
the eigenvalue problem of matrix $\mathbf{K}$, i.e.,
$\mathbf{K}\bm{v}=\lambda\bm{v}$, where $\lambda\in\operatorname{\mathbb{C}}$
and $\bm{v}$ is a complex vector. The vector $\bm{v}$ is an approximation of
$\mathcal{K}$’s eigenfunction $\varphi$ evaluated at the collocation points
and $\lambda$ is the approximation of the associated eigenvalue of
$\mathcal{K}$. The matrix form of the eigenvalue problem is

$\mathbf{K}\mathbf{V}=\mathbf{V}\mathbf{\Lambda},$ (14)

where $\mathbf{V}$ consists of columns $\bm{v}_{j}$ and the diagonal elements
of $\mathbf{\Lambda}$ are $\lambda_{j}$. By construction, for $d=1$,
$(\bm{v}_{j})_{i}\approx\varphi_{j}(\xi_{i})$, and for $d=2$ or $3$,
$\bm{v}_{j}=\operatorname{vec}(\mathbf{\Phi}_{j})$, where $\mathbf{\Phi}_{j}$
consists of approximations of eigenfunction $\varphi_{j}$ at the collocation
points. Of note, the collocation points in multi-dimensional cases are
constructed by the tensor product of one-dimensional collocation points, but
we have not specified how to obtain such points, which will be explained in
section 3.3. Also, we emphasize that these collocation points are related to
$\bm{x}$ instead of $t$. In other words, we discretize $\varphi(\bm{x})$ in
space instead of discretizing $\bm{x}(t)$ in time, which is different from
conventional spectral methods for ODEs.

###  3.3 Constructing the solution

We first consider $d=1$. By the eigen-decomposition, we can access values of
eigenfunctions at the Gauss-Lobatto points
$\bm{\Xi}:=\\{\xi_{i}\\}_{i=0}^{N}$, where $\xi_{0}<\xi_{1}<\dotsc<\xi_{N}$.
Therefore, we can approximate $\varphi(\bm{x}_{0})$ when
$\xi_{0}\leq\bm{x}_{0}\leq\xi_{N}$. To avoid polynomial interpolation, we use
an even number for $N$ and set $\xi_{N/2}=\bm{x}_{0}$ in our algorithm. Based
on this setting, we consider a neighborhood of $\bm{x}_{0}$ with radius $r$,
i.e., $[\bm{x}_{0}-r,\bm{x}_{0}+r]$, where $r$ is tunable. Then we generate
Gauss-Lobatto points that satisfy
$\bm{x}_{0}-r=\xi_{0}<\xi_{2}<\dotsc<\xi_{N/2}=\bm{x}_{0}<\dotsc<\xi_{N}=\bm{x}_{0}+r$.
Thus, $g_{N}$ is constructed as

$g_{N}(\bm{x}(t))=\sum_{j=0}^{N}c_{j}\varphi_{j}(\bm{x}_{0})e^{\lambda_{j}t}=\sum_{j=0}^{N}c_{j}\varphi_{j}(\xi_{N/2})e^{\lambda_{j}t}\approx\sum_{j=0}^{N}c_{j}(\bm{v}_{j})_{N/2}e^{\lambda_{j}t},$
(15)

where $\bm{v}_{j}$ are eigenvectors of matrix $\mathbf{K}$ computed in section
3.1.

To compute Koopman modes $c_{j}$, we set $t=0$ in eq. 15, which yields

$g(\bm{x}_{0})\approx
g_{N}(\bm{x}_{0})=\sum_{j=0}^{N}c_{j}\varphi_{j}(\bm{x}_{0}),$

which holds for different initial state $\bm{x}_{0}$, e.g.,

$g(\xi_{i})\approx
g_{N}(\xi_{i})=\sum_{j=0}^{N}c_{j}\varphi_{j}(\xi_{i}),\quad i=0,\dotsc,N,$

where $\xi_{i}$ are the aforementioned Gauss-Lobatto points. Thus, we can
obtain $c_{j}$ by solving a linear system $\mathbf{V}\bm{c}=g(\bm{\Xi})$,
where $\mathbf{V}$ is defined in eq. 14,
$g(\bm{\Xi})=(g(\xi_{0}),\dotsc,g(\xi_{N}))^{\top}$ and
$\bm{c}=(c_{0},\dotsc,c_{N})^{\top}$. In this work, since we set
$g(\bm{x})=\bm{x}$, thus $g(\bm{\Xi})=(\xi_{0},\dotsc,\xi_{N})^{\top}$.

For $d=2$, we consider the neighborhood of
$\bm{x}_{0}=(x_{0}^{1},x_{0}^{2})^{\top}$ as
$[x_{0}^{1}-r,x_{0}^{1}+r]\times[x_{0}^{2}-r,x_{0}^{2}+r]$. Similarly, for
$d=3$, we set the neighborhood as
$[x_{0}^{1}-r,x_{0}^{1}+r]\times[x_{0}^{2}-r,x_{0}^{2}+r]\times[x_{0}^{3}-r,x_{0}^{3}+r]$,
where $\bm{x}_{0}=(x_{0}^{1},x_{0}^{2},x_{0}^{3})^{\top}$. We then generate
$(N+1)$ Gauss-Lobatto points in each direction and use the tensor product rule
to construct multi-dimensional collocation points. In practice, one can use
standard Gauss-Lobatto points in the spectral method such as Legendre-Gauss-
Lobatto and Chebyshev-Gauss-Lobatto points. Now the set of all collocation
points is $\mathbf{\Xi}=\\{(\xi_{i},\eta_{j})\\}_{i,j=0}^{N}$ for $d=2$ and
$\mathbf{\Xi}=\\{(\xi_{i},\eta_{j},\zeta_{k})\\}_{i,j=0}^{N}$ for $d=3$. Of
note, we use the isotropic set up here for demonstration purpose, i.e., we use
a fixed $r$ in each direction and use the same number of Gauss-Lobatto points
in each dimension. This is not necessarily the optimal choice, and one can use
different $r$ and different numbers of Gauss-Lobatto points in different
directions.

Next, since we vectorize matrix (or tensor) $\mathbf{\Phi}$ column by column
(or slice by slice) as shown in section 3.1, $\varphi_{j}(\bm{x}_{0})$ is
again approximated by the “middle” element of vector
$\operatorname{vec}(\mathbf{\Phi}_{j})$, which leads to

$\bm{g}_{N}(\bm{x}(t))\approx\begin{cases}\displaystyle\sum_{j=0}^{(N+1)^{2}-1}\bm{c}_{j}\>(\bm{v}_{j})_{[(N+1)^{2}-1]/2}\>e^{\lambda_{j}t},&d=2;\\\
\displaystyle\sum_{j=0}^{(N+1)^{3}-1}\bm{c}_{j}\>(\bm{v}_{j})_{[(N+1)^{3}-1]/2}\>e^{\lambda_{j}t},&d=3.\\\
\end{cases}$ (16)

Here, each element of the modes
$\bm{c}_{j}=(c_{j}^{1},\cdots,c_{j}^{d})^{\top}$ corresponds to a component of
$\bm{g}$, and it is computed in the same manner as in the $d=1$ case. For
example, for $d=2$, i.e.,
$\bm{g}(\bm{x})=(g_{1}(\bm{x}),g_{2}(\bm{x}))^{\top}$ (correspondingly,
$\bm{g}_{N}(\bm{x})=(g_{N}^{1}(\bm{x}),g_{N}^{2}(\bm{x}))^{\top}$), we have
$g_{N}^{1}(\bm{x}(t))\approx\sum_{j=0}^{(N+1)^{2}-1}c^{1}_{j}\operatorname{vec}(\bm{\Phi}_{j})_{[(N+1)^{2}-1]/2}e^{\lambda_{j}t}$.
Consider matrix $g_{1}(\mathbf{\Xi})$ whose elements are
$(g_{1}(\mathbf{\Xi}))_{ij}=g_{1}(\xi_{i},\eta_{j})$. The modes
$\bm{c}^{1}=(c_{0}^{1},\dotsc,c_{N}^{1})^{\top}$ are obtained by solving a
linear system $\mathbf{V}\bm{c}^{1}=\operatorname{vec}(g_{1}(\mathbf{\Xi}))$.
Similarly, we can compute the coefficient for $g_{N}^{2}(\bm{x})$. In
practice, we solve the linear system $\mathbf{V}\mathbf{C}=\bm{g}(\bm{\Xi})$,
where $\mathbf{C}=(\bm{c}^{1},\bm{c}^{2})$ and
$\bm{g}(\bm{\Xi})=(\operatorname{vec}(g_{1}(\bm{\Xi})),\operatorname{vec}(g_{l}(\bm{\Xi})))$.
The modes for $d=3$ are computed in the same manner. We present a pseudocode
in appendix B to illustrate how the solution is constructed.

###  3.4 Adaptivity

Since we use a finite-dimensional approximation of the Koopman operator and
use the Lagrange interpolation to approximate the eigenfunctions, the accuracy
of the solution may decay as time evolves, especially for highly nonlinear
systems. To further improve the accuracy, we propose an adaptive approach to
update $\mathbf{V},\mathbf{\Lambda}$ and $\bm{c}_{j}$. The main idea is to
identify the time when we need to repeat the procedure described in section
3.1– section 3.3. To this end, we set check points
$0<\tau_{1}<\tau_{2}<\dotsc<\tau_{n}<T$ to examine the “validity” of the
neighborhood of $\bm{x}(\tau_{k})$. Specifically, the component of
$\bm{x}(\tau_{k})=(x_{1}(\tau_{k}),\dotsc,x_{d}(\tau_{k}))^{\top}$ is
acceptable if $x_{i}(\tau_{k})\in R_{i}$ where

$\displaystyle R_{i}:=[L_{i}+\gamma r_{i},U_{i}-\gamma r_{i}].$ (17)

Here, $L_{i}$ and $U_{i}$ are the lower and upper bounds, $r_{i}$ is the
radius in the $i$th direction, and $\gamma\in(0,1]$ is a tunable parameter.
Recall that we use the isotropic setup in this work, thus, $r_{i}\equiv r$. In
the initial step, we set $L_{i}=x_{0}^{i}-r_{i}$ and $U_{i}=x_{0}^{i}+r_{i}$,
i.e., $\gamma=1$. We note that, in practice, we can fix $\gamma=1$ (or other
real number in $(0,1]$) and tune $r_{i}$ only. We keep both $\gamma$ and
$r_{i}$ here for future extension to anisotropic design and more advanced
adaptivity criterion.

If $x_{i}(\tau_{k})\in R_{i}$ for all $i$, then $R_{1}\times\dotsc\times
R_{d}$ is a valid neighborhood of $\bm{x}(\tau_{k})$. Otherwise, we update all
$L_{i},U_{i}$ and reconstruct $\varphi_{j},\lambda_{j},c_{j}$ to obtain
$\bm{x}(t)~{}(t>\tau_{k})$ as follows:

  1. 1.

Set $L_{i}=x_{i}(\tau_{k})-r_{i},U_{i}=x_{i}(\tau_{k})+r_{i}$, $1\leq i\leq
d$.

  2. 2.

Generate Gauss-Lobatto points and the differentiation matrix in each interval
$[L_{i}^{k},U_{i}^{k}]$. Repeat the procedure in section 3.1 to compute matrix
$\mathbf{K}$.

  3. 3.

Repeat the eigendecomposition in section 3.2 to update $\mathbf{V}$ and
$\mathbf{\Lambda}$ in eq. 14.

  4. 4.

Compute coefficient $\bm{c}_{j}$ as in section 3.3 with the updated
$\mathbf{V}$.

  5. 5.

Construct solution $\bm{x}(t)$ by replacing $e^{\lambda_{j}t}$ with
$e^{\lambda_{j}(t-\tau_{k})}$ in eq. 15 (or eq. 16 for $d=2,3$).

We note that the modification of constructing the solution in step 5 is
necessary because when an update is performed, we need to set $t_{0}=\tau_{k}$
and $\bm{x}_{0}=\bm{x}(t_{0})=\bm{x}(\tau_{k})$.

The parameter $\gamma$ decides how often we update the neighborhood and
reconstruct the solution. By construction, if we use a larger $\gamma$, then
we need to update the eigendecomposition more frequently. The extreme case
$\gamma=1$ enforces the update at every check point. In this work, we set
$\tau_{k+1}-\tau_{k}\equiv\Delta\tau$. Notably, since we discretize the
solution _in space_ instead of in time as in conventional ODE solvers, the
check points are different from time grids $0<t_{1}<t_{2}<\dotsc$ in those
solvers. If we set $k=0$, then no update is made, which indicates that the
solution $\bm{x}(t)$ only relies on the eigendecomposition based on
$\bm{x}_{0}$ (see the example pseudocode in Appendix appendix B).

###  3.5 Properties of the algorithm

In this work $\bm{x},\operatorname{\mathbf{f}},\bm{g}$ are real-valued
functions. Now we show that the solutions obtained by ASK are real numbers,
although $\mathbf{V},\mathbf{\Lambda},\bm{c}_{j}$ may contain complex values.
We start with reiterating a well-known conclusion:

######  Lemma 3.1.

If a real matrix has complex eigenvalues, then they always occur in complex
conjugate pairs. Furthermore, a complex conjugate pair of eigenvalues have a
complex conjugate pair of associated eigenvectors.

###### Proof.

Suppose the matrix $\mathbf{K}\in\operatorname{\mathbb{R}}^{n\times n}$ has an
eigenpair $\bm{v}$ and $\lambda$ such that $\mathbf{K}\bm{v}=\lambda\bm{v}$.
Let $\bar{\cdot}$ operator denote the complex conjugate. Taking the complex
conjugate of both sides of the equation, we have
$\bar{\mathbf{K}}\bar{\bm{v}}=\bar{\lambda}\bar{\bm{v}}$. However,
$\mathbf{K}=\bar{\mathbf{K}}$ since $\mathbf{K}$ has real entries. Thus,
$\mathbf{K}\bar{\bm{v}}=\bar{\lambda}\bar{\bm{v}}$. The claim follows. ∎

Next we present our theorem:

######  Theorem 3.2.

ASK yields real-valued solutions for dynamical systems with real-valued
$\bm{x},\operatorname{\mathbf{f}}$ and $\bm{g}$.

###### Proof.

We only need to consider the $d=1$ case since the solution for high-
dimensional cases are constructed in the same manner. Let $\bm{v}$ be a
eigenvector, then it is a column of matrix $\mathbf{V}$ in eq. 14. We only
need to consider the case when $\bm{v}$ is not a real-valued vector. According
to lemma 3.1 $\bar{\bm{v}}$ is also a column of $\mathbf{V}$. Let $\bm{u}$ be
a row of $\mathbf{V}^{-1}$ such that $\bm{u}\cdot\bm{v}=1$ and
$\bm{u}\cdot\tilde{\bm{v}}=0$, where $\tilde{\bm{v}}$ is any column of
$\mathbf{V}$ other than $\bm{v}$. It is clear that
$\bar{\bm{u}}\cdot\bar{\bm{v}}=1$ and $\bar{\bm{u}}\cdot\tilde{\bm{v}}=0$.
Therefore, $\bar{\bm{u}}$ is also a row of $\mathbf{V}^{-1}$. Next, as shown
in section 3.3, we compute the modes $\bm{c}$ as
$\bm{c}=\mathbf{V}^{-1}\bm{g}(\mathbf{\Xi})$. Let $c_{m}$ be the element of
$\bm{c}$ such that $c_{m}=\bm{u}\cdot\bm{g}(\mathbf{\Xi})$, then
$\bar{\bm{u}}\cdot\bm{g}(\mathbf{\Xi})=\overline{\bm{u}\cdot\bm{g}(\mathbf{\Xi})}=\bar{c}_{m}$
is also an element of $\bm{c}$.

In the numerical solution, it suffices to consider $c_{m}\nu e^{\lambda
t}+\bar{c}_{m}\bar{\nu}e^{\bar{\lambda}t}$, where
$\nu\in\operatorname{\mathbb{C}}$ denotes the middle element of the
eigenvector $\bm{v}$. For convenience, we dentoe
$\nu=A+Bi,\lambda=C+Di,c_{m}=E+Fi$. Here,
$A,B,C,D,E,F\in\operatorname{\mathbb{R}}$. Then,

$\displaystyle c_{m}\nu e^{\lambda t}+\bar{c}_{m}\bar{\nu}e^{\bar{\lambda}t}$
$\displaystyle=(E+Fi)(A+Bi)e^{(C+Di)t}+(E-Fi)(A-Bi)e^{(C-Di)t}$
$\displaystyle=(P+Qi)e^{(C+Di)t}+(P-Qi)e^{(C-Di)t}$
$\displaystyle=\big{(}Pe^{Ct}+Qe^{Ct}i\big{)}e^{Dti}+\big{(}Pe^{Ct}-Qe^{Ct}i\big{)}e^{-Dti}$
$\displaystyle=\big{(}Pe^{Ct}+Qe^{Ct}i\big{)}[\cos(Dt)+\sin(Dt)i]$
$\displaystyle\quad+\big{(}Pe^{Ct}-Qe^{Ct}i\big{)}[\cos(Dt)-\sin(Dt)i]$
$\displaystyle=2Pe^{Ct}\cos(Dt)+2\big{(}Qe^{Ct}i)\sin(Dt)i$
$\displaystyle=2Pe^{Ct}\cos(Dt)-2Qe^{Ct}\sin(Dt)\in\operatorname{\mathbb{R}},$

among which $P=AE-BF$ and $Q=AF+BE$.

∎

###  3.6 Algorithm summary

As a summary, section 3.1 to section 3.3 present a numerical scheme that
solves an autonomous ODE system using the eigendecomposition and a linear
system solver. Section 3.4 introduces a heuristic adaptivity criterion to
repeat the aforementioned procedure at appropriate time points to further
enhance the accuracy. We conclude the algorithm in algorithm 1.

Algorithm 1 Adaptive spectral Koopman method

0: $n,T,N,\bm{x}_{0},r,\gamma$

1: Set check points at $0=\tau_{0}<\tau_{1}<...<\tau_{n}<T$.

2: Let $L_{i}=x_{0}^{i}-r_{i},R_{i}=x_{0}^{i}+r_{i}$ and set neighborhood
$R_{i}$ as $R_{i}=[L_{i}+\gamma r_{i},U_{i}-\gamma r_{i}]$ for $i=1,2,...,d$,
where $r_{i}=r$.

3: Generate Gauss-Lobatto points and differentiation matrix $\mathbf{D}_{i}$
in $[L_{i},U_{i}]$ for $i=1,2,...,d$. Construct collocation points $\bm{\Xi}$
for $d>1$ using the tensor product rule. (For $d=1$, $\bm{\Xi}$ is the set of
the Gauss-Lobatto points.)

4: Construct matrix $\mathbf{K}$ using the formulas in section 3.1

5: Compute eigen-decomposition $\mathbf{K}\mathbf{V}=\mathbf{V}\Lambda$

6: Solve linear system $\mathbf{V}\mathbf{C}=\bm{g}(\bm{\Xi})$, where the
$l$th column of matrix $\bm{g}(\bm{\Xi})$ consists of the $l$th component of
all collocation points (see section 3.3).

7: for $k=1,2,3,...,n$ do

8: Let $\nu_{j}$ be the middle element of the $j$th column of $\mathbf{V}$.
Construct solution at time $\tau_{k}$ as
$\bm{x}(\tau_{k})=\displaystyle\sum_{j}\mathbf{C}(j,:)\nu_{j}e^{\lambda_{j}(\tau_{k}-\tau_{k-1})}$,
where $\mathbf{C}(j,:)$ is the $j$th row of $\mathbf{C}$.

9: if $(\bm{x}(\tau_{k}))_{i}\notin R_{i}$ for any $i$ then

10: Reset $L_{i}=x_{i}(\tau_{k})-r_{i},U_{i}=x_{i}(\tau_{k})+r_{i}$ and
$R_{i}=[L_{i}+\gamma r_{i},U_{i}-\gamma r_{i}]$.

11: Repeat steps $3-6$

12: end if

13: end for

14: return
$\bm{x}(T)=\displaystyle\sum_{j}\mathbf{C}(j,:)\nu_{j}e^{\lambda_{j}(T-\tau_{n})}$.

In the ASK scheme, the neighborhood for all components must be updated in the
adaptivity step. This is because we set the current state of each component to
be the midpoint of the corresponding neighborhood to avoid computing
interpolation. Also, following the standard practice in the spectral method,
we generate Gauss-Lobatto points $\xi_{i}$ and the associated differentiation
matrix $\mathbf{D}_{i}$ on $[-1,1]$ first, and then scale them to
$[L_{i},U_{i}]$ as $\frac{U_{i}-L_{i}}{2}(\xi_{i}+1\big{)}+L_{i}$ and
$\frac{2\mathbf{D}_{i}}{U_{i}-L_{i}}$ to improve the computational efficiency.
Moreover, we note again that the isotropic setup (i.e., using the same number
of Gauss-Lobatto points in each direction and fix $r_{i}\equiv r$) is not
necessarily the optimal choice and the adaptivity in different directions may
improve the efficiency of the algorithm. This is beyond the scope of this work
and will be included in the future study.

##  4 Numerical Results

In this section we present the performance of ASK on six nonlinear ODE systems
including $d=1,2,3$. In each example, we investigate the influence of number
of Gauss-Lobatto points $N$, number of check points $n$ and the radius $r$ on
the accuracy. The reference solution is generated by Verner’s ninth order
Runge-Kutta (RK9) method [31] if a close-form solution is not available. We
compare ASK with fourth order Runge-Kutta (RK4) method [7] since it is the
most widely used method in the Runge-Kutta family. In this work, we use
Chebyshev-Gauss-Lobatto points. We also tested Legendre-Gauss-Lobatto points
but there is no significant difference. Moreover, in section 4.7 we compare
the efficiency of ASK, RK4, and RK9 based on number of function evaluations
for systems with closed-form solutions. Herein, the function refers to
$\operatorname{\mathbf{f}}$ in the ODE. (All the codes can be downloaded at
https://github.com/Navarro33/Adaptive-Spectral-Koopman-Method)

###  4.1 Cosine model

The cosine model is a synthetic model invented for our demonstrative purposes.
The governing ODE is written as

$\frac{\mathrm{d}x}{\mathrm{d}t}=-0.5\cos^{2}(x).$

We set $x(0)=\frac{\pi}{4}$ and $T=20$ in this example. Despite the
nonlinearity, the system has a closed-form solution
$x(t)=\arctan(-0.5t+\tan(x_{0}))$. We aim to compute the solution at $T=20$.
The three experiments use the following parameters:

  1. (a)

test of $N$: $n=200,r=\frac{\pi}{20}$;

  2. (b)

test of $n$: $N=9,r=\frac{\pi}{20}$;

  3. (c)

test of $r$: $n=200,N=9$.

In all these tests, we set $\gamma=0.2$. Figure 1 summarizes these results in
plots (a), (b) and (c), respectively. The first test shows the exponential
convergence of ASK with respect to $N$, which is similar to the conclusions in
conventional spectral methods. Test (b) shows that the accuracy does not
change monotonically as $n$ varies given the parameter setting in this work.
On the other hand, using no more than $100$ check points is sufficient to
obtain good accuracy. The last test illustrates that the accuracy shows a “V
shape” with respect to the radius, i.e., $r$ can not be too large or too
small.

(a)

(b)

(c)

Figure 1: Cosine Model: (a) testing number of Gauss-Lobatto points $N$; (b)
testing number of checck points $n$; (c) testing radius $r$.

###  4.2 Lotka-Volterra model

The Lotka-Volterra equations model the interactive evolution of the population
of prey and predators [2]. Specifically, it is defined by

$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{d}t}=1.1x_{1}-0.4x_{1}x_{2},$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{d}t}=0.1x_{1}x_{2}-0.4x_{2}.$

We set $\bm{x}(0)=(10,5)^{\top}$ and $T=20$ in this example. The parameters
used in three different test are as follows:

  1. (a)

test of $N$: $n=200,r=1.5$;

  2. (b)

test of $n$: $N=5,r=1.5$;

  3. (c)

test of $r$: $n=200,N=5$.

In all the tests, we set $\gamma=0.5$. Note that for multi-dimensional systems
in the test on radius, all components share the same radius. Figure 2 presents
the results of these tests. Similar to the cosine model, the error decreases
exponentially with respect to $N$. The accuracy is quite stable with respect
to the number of check points in this case. Furthermore, fig. 2(c) shows that
the radius cannot not be too small as in the first example.

(a)

(b)

(c)

Figure 2: Lotka-Volterra Model: (a) testing number of Gauss-Lobatto points $N$
(total number of collocation points is $(N+1)^{2}$); (b) testing number of
checck points $n$; (c) testing radius $r$. $\medcircle,\medsquare$ denote
$x_{1},x_{2}$, respectively.

###  4.3 Simple pendulum

The simple pendulum is well studied in physics and mechanics. The movement of
the pendulum is described by a second order ordinary differential equation,

$\displaystyle\frac{\mathrm{d}^{2}\theta}{\mathrm{d}t^{2}}=-\frac{g}{L}\sin(\theta).$

Here, $\theta$ is the displacement angle and $L$ denotes the length of the
pendulum. The parameter $g$ is the gravity acceleration. This second order
equation can be converted to a two-dimensional first-order ODE system. To make
the notations consistent, we define $x_{1}:=\theta$ and
$x_{2}:=\frac{d\theta}{dt}$. Also, we set $L=g=9.8$ in our numerical
experiments. Correspondingly,

$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{d}t}$ $\displaystyle=x_{2},$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{d}t}$ $\displaystyle=-\sin
x_{1}.$

We set $\bm{x}(0)=(-\frac{\pi}{4},\frac{\pi}{6})^{\top}$ and $T=20$. The
parameters in three tests are as follows:

  1. (a)

test of $N$: $n=200,r=\frac{\pi}{8}$;

  2. (b)

test of $n$: $N=7,r=\frac{\pi}{8}$;

  3. (c)

test of $r$: $n=200,N=7$.

We set $\gamma=0.2$ in all these tests. The resutls are presented in fig. 3.
Again, we observe exponential convergence with respect to $N$ in fig. 3(a).
fig. 3(b) implies that more check points can improve the accuracy but the
difference is not very large. fig. 3(c) indicates there exists an “optimal”
$r$ as in the Lotka-Volterra example.

(a)

(b)

(c)

Figure 3: Simple Pendulum: (a) testing number of Gauss-Lobatto points $N$
(i.e., $(N+1)^{2}$ collocation points in total); (b) testing number of checck
points $n$; (c) testing radius $r$. $\medcircle,\medsquare$ denote
$x_{1},x_{2}$, respectively.

###  4.4 Limit cycle

The limit cycle is applied to model oscillatory systems in multiple research
fields [32]. Here, we follow the definition,

$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{d}t}=-x_{1}-x_{2}+\frac{x_{1}}{\sqrt{x_{1}^{2}+x_{2}^{2}}},$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{d}t}=x_{1}-x_{2}+\frac{x_{2}}{\sqrt{x_{1}^{2}+x_{2}^{2}}}.$

The closed-form solution is $x_{1}(t)=\cos(t+\arcsin(x_{1}(0)))$ and
$x_{2}(t)=\sin(t+\arcsin(x_{2}(0)))$. We set
$\bm{x}(0)=(\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2})^{\top}$ and $T=20$ in this
example. The parameters in the experiments are specified as follows:

  1. (a)

test of $N$: $n=200,r=\frac{\sqrt{2}}{6}$;

  2. (b)

test of $n$: $N=7,r=\frac{\sqrt{2}}{6}$;

  3. (c)

test of $r$: $n=200,N=7$.

We set $\gamma=0.2$ in all these tests. The results shown in fig. 4 reveal
similar patterns to the results of the simple pendulum, except that a very
small $r$ can still lead to accurate results.

(a)

(b)

(c)

Figure 4: Limit Cycle: (a) testing number of Gauss-Lobatto points $N$ (i.e.,
$(N+1)^{2}$ collocation points in total); (b) testing number of checck points
$n$; (c) testing radius $r$. $\medcircle,\medsquare$ denote $x_{1},x_{2}$,
respectively.

For this example, we also compared ASK with RK4 at various time within
$[0,T]$. Given the closed-form solution $x_{C}(t)$, we computed the errors by
$|x_{ASK}(t)-x_{C}(t)|$ and $|x_{RK4}(t)-x_{C}(t)|$. Here, RK4 employed
$M=200$ equidistant time points on $[0,T]$. As for ASK, we used
$N=9,r=\frac{\sqrt{2}}{8},\gamma=0.2$ and the check points are set to be the
same as the time points in RK4. With this set of parameters, ASK constantly
outperforms RK4 significantly, as illustrated in fig. 5. For both components
$x_{1}$ and $x_{2}$, the errors of ASK remain almost constant at the level of
$10^{-10}$. In comparison, the error of RK4 exhibited a periodic pattern,
rising slowly from $10^{-6}$ to $10^{-5}$. Moreover, fig. 6 illustrates the
evolution of the limit cycle model along time. The path decided by the two
components elevates spirally as time goes. If seen from above, the cross
section is an exact circle.

(a)

(b)

Figure 5: Errors of the limit cycle solution on $[0,T]$: $\medcircle$ denotes
ASK and $\medtriangleup$ denotes RK4. Figure 6: Limit cycle solution
trajectory: $\medcircle$ denotes ASK and $-$ dentoes the closed-form
solutions.

###  4.5 Kraichnan-Orszag model

The Kraichnan-Orszag model comes from the problem raised in [21]. This system
is nonlinear and three-dimensional, defined by

$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{d}t}=x_{2}x_{3},$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{d}t}=x_{1}x_{3},$
$\displaystyle\frac{\mathrm{d}x_{3}}{\mathrm{d}t}=-2x_{1}x_{2}.$

We set $\bm{x}(0)=(1,2,-3)^{\top}$ and $T=20$. In the three experiments, we
employed the following parameters:

  1. (a)

test of $N$: $n=400,r=0.1$;

  2. (b)

test of $n$: $N=3,r=0.1$;

  3. (c)

test of $r$: $n=400,N=3$.

Also, in all the tests, we set $\gamma=0.15$. The results are presented in
fig. 7(a). In particular, different from previous examples, fig. 7(b)
demonstrates that $n$ significantly influences the accuracy. This is because
the Kraichnan-Orszag model exhibited strong oscillations, so it requires more
frequent update of eigenpairs to guarantee high accuracy. Figure 7(c)
indicates that there is an “optimal” $r$ as in the Lokta-Volterra example.

(a)

(b)

(c)

Figure 7: Kraichnan-Orszag model: (a) testing number of Gauss-Lobatto points
$N$ (i.e., $(N+1)^{3}$ collocation points in total); (b) testing number of
checck points $n$; (c) testing radius $r$.
$\medcircle,\medsquare,\medtriangleup$ denote $x_{1},x_{2},x_{3}$,
respectively

###  4.6 Lorenz attractor

The Lorenz attractor was first introduced by Lorenz [15]. It is a highly
chaotic system that models the turbulence in dynamic flows. The governing
equations are as follows,

$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{d}t}=10(x_{2}-x_{1}),$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{d}t}=x_{1}(28-x_{3})-x_{2},$
$\displaystyle\frac{\mathrm{d}x_{3}}{\mathrm{d}t}=x_{1}x_{2}-3x_{3}.$

We set $\bm{x}(0)=(5,5,5)^{\top}$ and $T=20$ in this example. Parameters used
in the experiments are listed here:

  1. (a)

test of $N$: $n=500,r=1$;

  2. (b)

test of $n$: $N=5,r=4$;

  3. (c)

test of $r$: $n=500,N=5$.

In all the tests, we set $\gamma=0.5$. The results are summarized in fig. 8.
As the Loren attractor exhibits chaotic behaviour, it requires a larger number
of check points. Meanwhile, a relatively large radius favored the convergence
of the algorithm.This is probably because the eigenfunctions need to be
approximated in a larger neighborhood of the solution to include sufficient
information of the dynamics.

(a)

(b)

(c)

Figure 8: Loren attractor: (a) testing number of Gauss-Lobatto points $N$
(i.e., $(N+1)^{3}$ collocation points in total); (b) testing number of checck
points $n$; (c) testing radius $r$. $\medcircle,\medsquare,\medtriangleup$
denote $x_{1},x_{2},x_{3}$, respectively.

Next we compare ASK and RK4 in a larger time span $[0,20]$, where RK4 uses
$M=2000$ time steps, i.e., step size $\Delta t=0.01$. Since the Lorenz
attractor does not have closed-form solutions, RK9 is used to compute the
reference. To guarantee accuracy, RK9 used step size $\Delta t=0.001$, i.e.,
$M=20000$ time steps. On the other hand, ASK was implemented with
$n=2000,N=5,r=(1,1,1),\gamma=0.75$. For the aforementioned reason, we required
a large number of check points and small tolerance for the acceptable range.
Figure 9 reveals the accuracy of ASK in all three components. However, unlike
the limit cycle case, the error increases as time evolves. Although the error
increases to around $10^{-3}$ at $T=20$, ASK still yields an acceptable
accuracy for such a chaotic system. In comparison, RK4’s error ascends to a
level that makes it impractical. To obtain an insight of how the Lorenz system
evolves, we plot each of its component in fig. 10. Up to time $T=10$,
solutions given by ASK, RK4, and RK9 almost coincide. Nevertheless, RK4
deviates from the other two completely starting at $t=11$. The evolution
vibrates violently and does not possess periodicity, which imposes difficulty
on numerical solvers.

(a)

(b)

(c)

Figure 9: Error of Lorenz attractor solutions on $[0,T]$: $\medcircle$ denotes
ASK and $\medtriangleup$ dentoes RK4.

(a)

(b)

(c)

Figure 10: Lorenz Attractor Evolution: $\medcircle$ and $\medtriangleup$
denote ASK and RK4, respectively; $--$ donotes the reference solutions given
by RK9

The chaos can also be observed in a three-dimensional graph depicting the
trajectory, using the numerical solutions given by ASK. As in fig. 11, the
lemniscate shape demonstrates the complexity of the system.

Figure 11: Lorenz Attractor 3D Visualization

###  4.7 Computational complexity

By construction, the computational complexity of conventional ODE solvers that
approximate time integral is $\mathcal{O}(M)$ where $M$ is the number of time
steps. In other words, it is $M$ multiplied by a constant that represents the
number of operations in each step and it varies according to the accuracy of
the scheme. The computational complexity of ASK depends on the number of times
that eigenfunctions are constructed (and corresponding eigenvalues as well as
coefficients are computed). In this construction procedure, ASK needs to
perform the eigendecomposition and solve a linear system. For $d=1$ this is
not costly because empirically we set $4\leq N\leq 10$, and the size of matrix
in the eigendecomposition as well as the linear system is $N\times N$. But
when $d>1$, the complexity will increase exponentially with the dimension of
the current setting because we use the tensor product rule to construct the
collocation points and the matrix size is $(N+1)^{d}\times(N+1)^{d}$. Hence,
ASK can be less efficient for high-dimensional systems.

However, compared with the Runge-Kutta methods, ASK has advantages when it is
costly to evaluate the dynamics $\operatorname{\mathbf{f}}$ and when ASK only
needs a small number of updates. Here, we consider an evaluation of
$\operatorname{\mathbf{f}}$ as a function call, and compare the accuracy of
ASK, RK4 and RK9 against number of function calls. We use the cosine model
($d=1$) and the simple pendulum ($d=2$) to perform the comparison and show the
results in fig. 12. Here, the error in the simple pendulum case was computed
by $\sqrt{\frac{e_{1}^{2}+e_{2}^{2}}{2}}$, where $e_{1},e_{2}$ are the errors
in $x_{1},x_{2}$, respectively. The parameters in ASK are
$N=[4,6,8,10],n=50,r=\frac{\pi}{20}$ for the cosine model and
$N=[2,4,6,8],n=60,r=\frac{\sqrt{2}}{12}$ for the simple pendulum. ASK is
superior to the other two methods in the cosine model. For the simple
pendulum, RK9 is the most efficient method, and ASK outperforms RK4 when
number of function calls is beyond $2000$.

(a)

(b)

Figure 12: Test of Complexity: $\medcircle,\medsquare$, and $\medtriangleup$
denote ASK, RK4, and RK9, respectively.

##  5 Conclusion and Discussion

The ASK method uses the spectral-collocation (i.e., pseudo spectral) method in
the state space instead of in time to solve nonlinear autonomous dynamical
systems. It discretizes the generator of Koopman operator and employs the
eigendecomposition to obtain approximations of the eigenfunctions and
eigenvalues to construct solutions. Therefore, like the spectral method, ASK
is an expansion-based method to solve ODE systems, in which the basis
functions in the expansion are approximated eigenfunctions of the Koopman
operator. In each numerical example presented in this work, ASK exhibits
_exponential convergence_ as the conventional spectral method. Therefore, it
is suitable for the circumstances where high-accuracy solutions are desired
and $\operatorname{\mathbf{f}}$ is expensive to evaluate. Different from
existing ODE solvers that obtain solutions on mesh grids, ASK does not need a
time mesh and can evaluate the solution at any time. Hence, the resolution of
the time mesh which impacts the solutions of conventional ODE solvers like
Runge-Kutta methods does not influence ASK.

In the ASK algorithm, adaptively updating the approximation of eigenfunctions
in the neighborhood of the solution is necessary because it is challenging to
obtain very accuracy approximation of the eigenfunctions, eigenvalues and
Koopman modes using the initial condition only, especially for highly
nonlinear systems. When no information (e.g., range of states, regularity of
the eigenfunctions) of the system is available _a priori_ , the adaptivity
criterion serves as a updating step based on “ _posterior_ error estimates”.
Furthermore, tunable parameters $r$ and $\gamma$ affect the accuracy as they
are related to eigenfunction approximations and the adaptivity criterion.
Numerical analysis based on the spectrum theorem as well as the spectral
method is required to systematically understand the convergence and the impact
of all parameters on the performance of ASK, which will be included in our
future work.

Regarding the computational complexity, as indicated in section 4.7, ASK’s
efficiency decreases (compared with Runge-Kutta) as the system dimension
increases since the tensor product rule is applied to construct high-
dimensional collocation points. A possible way of improving the efficiency is
to employ the sparse grid methods to construct collocation points, which has
shown success in solving partial differential equations with the spectral
method [27, 28]. Also, applying an anisotropic setting, e.g., different number
of Gauss-Lobatto points, different radius, different $\gamma$ in each
direction, can potentially enhance the computational efficiency.

##  Appendix A An example of the obervable

As an example, we consider the following nonlinear dynamical system [3, 16]:

$\displaystyle\frac{\mathrm{d}x_{1}}{\mathrm{d}t}=\alpha x_{1},$
$\displaystyle\frac{\mathrm{d}x_{2}}{\mathrm{d}t}=\beta(x_{2}-x_{1}^{2}).$

Here, $\alpha$ and $\beta$ are the inherent parameters of the system. For such
a system, appropriate observables lead to a closed-form solution. In
particular, let $\bm{y}:=(x_{1},x_{2},x_{1}^{2})^{\top}$ be a three-
dimensional observable. Then, the system can be converted to the following
linear system,

$\displaystyle\frac{\mathrm{d}\bm{y}}{\mathrm{d}t}=\begin{bmatrix}\alpha&0&0\\\
0&\beta&-\beta\\\ 0&0&2\alpha\end{bmatrix}\bm{y}.$

For simplicity, assume $x_{1}(0)=x_{2}(0)=1$. Then, we have the closed-form
solution

$\displaystyle\bm{y}=\begin{bmatrix}1\\\ 0\\\ 0\end{bmatrix}e^{\alpha
t}+\frac{-2\alpha}{\beta-2\alpha}\begin{bmatrix}0\\\ 1\\\
0\end{bmatrix}e^{\beta t}+\begin{bmatrix}0\\\ \frac{\beta}{\beta-2\alpha}\\\
1\end{bmatrix}e^{2\alpha t}=\begin{bmatrix}e^{\alpha t}\\\
\frac{-2\alpha}{\beta-2\alpha}e^{\beta
t}+\frac{\beta}{\beta-2\alpha}e^{2\alpha t}\\\ e^{2\alpha t}\end{bmatrix}.$

Equivalently,

$\displaystyle x_{1}=e^{\alpha t},\qquad
x_{2}=\frac{-2\alpha}{\beta-2\alpha}e^{\beta
t}+\frac{\beta}{\beta-2\alpha}e^{2\alpha t}.$

##  Appendix B An example pseudocode

We demonstrate a pseudo code (in MATLAB) of solving
$\frac{\mathrm{d}x}{\mathrm{d}t}=\cos^{2}(x)$, which summarizes the steps
presented in section 3.1–section 3.3. The MATLAB code generating Chebyshev-
Gauss-Lobatto points and the associated differentiation matrix can be found in
[29].

        f = @(x) cos(x).^2;  % Function f
        x0 = pi/4;           % Initial condition
        r = 0.1;             % Radius of the neighborhood (tunable)
        N = 4;               % Number of collocation points (N+1 in total)
        T = 5;               % Final time
        % Generate collocation points and the differentiation matrix
        % on [x0-r, x0+r]
        [quad_pnt, diff_mat] = cheb(N, x0-r, x0+r);
        % Compute eigenpairs of the Koopman operator
        K = diag(f(quad_pnt))*diff_mat;
        [eig_vec, eig_val] = eig(K, ’vector’);
        % Compute coefficients (Koopman modes)
        coef = eig_vec\quad_pnt;
        % Construct solutions at time T
        sol = real(eig_vec(N/2,:).*coef’*exp(eig_val*T));

When the adaptive update in ASK is activated (see section 3.4), we only need
to repeat this pseudocode (as a subroutine) with an updated initial condition
$x_{0}$ and final time $T$.

## Acknowledgments

We thank Professor Yue Yu and Professor Hong Qian for fruitful discussions on
the spectral method and properties of the Koopman operator. Funding:  XY was
supported by the U.S. Department of Energy (DOE), Office of Science, Office of
Advanced Scientific Computing Research (ASCR) as part of Multifaceted
Mathematics for Rare, Extreme Events in Complex Energy and Environment Systems
(MACSER). BL was partially supported by Los Alamos National Laboratory. YM was
supported in part by the National Science Foundation Grants NSF-SCALE
MoDL(2134209) and NSF-CCF-2112665, and the Facebook research award. The work
of JNK was supported in part by the US National Science Foundation (NSF) AI
Institute for Dynamical Systems (dynamicsai.org), grant 2112085.

## References

  * [1] Travis Askham and J Nathan Kutz.  Variable projection methods for an optimized dynamic mode decomposition.  SIAM Journal on Applied Dynamical Systems, 17(1):380–416, 2018\.
  * [2] Immanuel M Bomze.  Lotka-Volterra equation and replicator dynamics: a two-dimensional classification.  Biological cybernetics, 48(3):201–211, 1983.
  * [3] Steven L Brunton, Bingni W Brunton, Joshua L Proctor, and J Nathan Kutz.  Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control.  PloS one, 11(2):e0150171, 2016.
  * [4] Steven L Brunton, Marko Budišić, Eurika Kaiser, and J Nathan Kutz.  Modern Koopman theory for dynamical systems.  arXiv preprint arXiv:2102.12086, 2021.
  * [5] Marko Budišić and Igor Mezić.  An approximate parametrization of the ergodic partition using time averaged observables.  In Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, pages 3162–3168. IEEE, 2009.
  * [6] Marko Budišić, Ryan Mohr, and Igor Mezić.  Applied Koopmanism.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(4):047510, 2012.
  * [7] John Charles Butcher and Gerhard Wanner.  Runge-Kutta methods: some historical notes.  Applied Numerical Mathematics, 22(1-3):113–151, 1996.
  * [8] Bengt Fornberg.  A practical guide to pseudospectral methods.  Number 1. Cambridge university press, 1998.
  * [9] Jan S Hesthaven, Sigal Gottlieb, and David Gottlieb.  Spectral methods for time-dependent problems, volume 21.  Cambridge University Press, 2007.
  * [10] George Karniadakis and Spencer Sherwin.  Spectral/hp element methods for computational fluid dynamics.  OUP Oxford, 2005.
  * [11] Bernard O Koopman.  Hamiltonian systems and transformation in Hilbert space.  Proceedings of the National Academy of Sciences of the United States of America, 17(5):315, 1931.
  * [12] Milan Korda, Mihai Putinar, and Igor Mezić.  Data-driven spectral analysis of the Koopman operator.  Applied and Computational Harmonic Analysis, 48(2):599–629, 2020\.
  * [13] J Nathan Kutz, Steven L Brunton, Bingni W Brunton, and Joshua L Proctor.  Dynamic mode decomposition: data-driven modeling of complex systems.  SIAM, 2016.
  * [14] J Nathan Kutz, Xing Fu, and Steven L Brunton.  Multiresolution dynamic mode decomposition.  SIAM Journal on Applied Dynamical Systems, 15(2):713–735, 2016\.
  * [15] Edward N Lorenz.  Deterministic nonperiodic flow.  Journal of atmospheric sciences, 20(2):130–141, 1963.
  * [16] Bethany Lusch, J Nathan Kutz, and Steven L Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature communications, 9(1):1–10, 2018.
  * [17] Igor Mezić.  Spectral properties of dynamical systems, model reduction and decompositions.  Nonlinear Dynamics, 41(1):309–325, 2005.
  * [18] Igor Mezić.  Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry.  Journal of Nonlinear Science, 30(5):2091–2145, 2020.
  * [19] Hiroya Nakao and Igor Mezić.  Spectral analysis of the Koopman operator for partial differential equations.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(11):113131, 2020.
  * [20] J Nathan Kutz, Joshua L Proctor, and Steven L Brunton.  Applied Koopman theory for partial differential equations and data-driven modeling of spatio-temporal systems.  Complexity, 2018, 2018.
  * [21] Steven A Orszag and LR Bissonnette.  Dynamical properties of truncated Wiener-Hermite expansions.  The Physics of Fluids, 10(12):2603–2613, 1967.
  * [22] Jacob Page and Rich R Kerswell.  Koopman analysis of Burgers equation.  Physical Review Fluids, 3(7):071901, 2018.
  * [23] Joshua L Proctor, Steven L Brunton, and J Nathan Kutz.  Dynamic mode decomposition with control.  SIAM Journal on Applied Dynamical Systems, 15(1):142–161, 2016\.
  * [24] Clarence W Rowley, Igor Mezić, Shervin Bagheri, Philipp Schlatter, and Dan S Henningson.  Spectral analysis of nonlinear flows.  Journal of fluid mechanics, 641:115–127, 2009.
  * [25] Peter J Schmid.  Dynamic mode decomposition of numerical and experimental data.  Journal of fluid mechanics, 656:5–28, 2010.
  * [26] Jie Shen and Tao Tang.  Spectral and high-order methods with applications.  Science Press Beijing, 2006.
  * [27] Jie Shen and Haijun Yu.  Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems.  SIAM Journal on Scientific Computing, 32(6):3228–3250, 2010.
  * [28] Jie Shen and Haijun Yu.  Efficient spectral sparse grid methods and applications to high-dimensional elliptic equations ii. unbounded domains.  SIAM Journal on Scientific Computing, 34(2):A1141–A1164, 2012.
  * [29] Lloyd N Trefethen.  Spectral methods in MATLAB.  SIAM, 2000.
  * [30] Jonathan H Tu, Clarence W Rowley, Dirk M Luchtenburg, Steven L Brunton, and J Nathan Kutz.  On dynamic mode decomposition: Theory and applications.  Journal of Computational Dynamics, 1(2):391–421, 2014.
  * [31] James Hamilton Verner.  Explicit Runge–Kutta methods with estimates of the local truncation error.  SIAM Journal on Numerical Analysis, 15(4):772–790, 1978.
  * [32] Mathukumalli Vidyasagar.  Nonlinear systems analysis.  SIAM, 2002.
  * [33] Matthew O Williams, Maziar S Hemati, Scott TM Dawson, Ioannis G Kevrekidis, and Clarence W Rowley.  Extending data-driven Koopman analysis to actuated systems.  IFAC-PapersOnLine, 49(18):704–709, 2016.
  * [34] Matthew O Williams, Ioannis G Kevrekidis, and Clarence W Rowley.  A data–driven approximation of the Koopman operator: Extending dynamic mode decomposition.  Journal of Nonlinear Science, 25(6):1307–1346, 2015.
  * [35] Dan Wilson and Jeff Moehlis.  Isostable reduction with applications to time-dependent partial differential equations.  Physical Review E, 94(1):012211, 2016.


# HofstadterTools: A Python package for analyzing the Hofstadter model

Bartholomew Andrews

††margin:  DOI: https://doi.org/DOI unavailablehttps://doi.orgDOI unavailable
Software • N/Ahttps://doi.orgReview • NO_REPOSITORYhttps://doi.orgRepository •
DOI unavailablehttps://doi.orgArchive Editor:
https://example.comhttps://doi.orgPending Editor
Reviewers: • https://github.com/Pending Reviewershttps://doi.org@Pending
Reviewers Submitted: N/A
Published: N/A License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/https://doi.orgCC BY 4.0).

## Summary

The Hofstadter model successfully describes the behavior of non-interacting
quantum particles hopping on a lattice coupled to a gauge field, and hence is
ubiquitous in many fields of research, including condensed matter, optical,
and atomic physics. Motivated by this, we introduce HofstadterTools
(https://hofstadter.toolshttps://doi.org), a Python package that can be used
to analyze the energy spectrum of a generalized Hofstadter model, with any
combination of hoppings on any regular Euclidean lattice. The package can be
applied to compute key properties of the band structure, such as quantum
geometry and topology, as well as plot Hofstadter butterflies and Wannier
diagrams that are colored according to their Chern numbers.

## Statement of need

The purpose of HofstadterTools is to consolidate the fragmented theory and
code relevant to the Hofstadter model into one well-documented Python package,
which can be used easily by non-specialists as a benchmark or springboard for
their own research projects. The Hofstadter model [1, 2, 3] is an iconic
tight-binding model in physics and famously yields a fractal energy spectrum
as a function of flux density, as shown in Figs. 1, 2, 3, and 4. Consequently,
it is often treated as an add-on to larger numerical packages, such as pyqula
[4], DiagHam [5], and TeNPy [6], or simply included as supplementary code
together with research articles [7]. However, the Hofstadter model’s
generalizability, interdisciplinary appeal, and recent experimental
realization, motivates us to create a dedicated package that can provide a
detailed analysis of its band structure, in the general case.

  1. 1.

Generalizability. The Hofstadter model was originally studied in the context
of electrons hopping in a periodic potential coupled to a perpendicular
magnetic field. However, the model transcends this framework and can be
extended in numerous directions. For example, the Peierls phases that arise in
the Hamiltonian due to the magnetic field [8] can also be generated using
artificial gauge fields [9] or Floquet modulation [10]. Moreover, the full
scope of the Hofstadter model is still being revealed, with papers on its
application to hyperbolic lattices [11], higher-dimensional crystals [12], and
synthesized materials [7], all published within the last couple of years.

  2. 2.

Interdisciplinary appeal. Owing to its generalizability, interest in the
Hofstadter model goes beyond its well-known connection to condensed matter
physics and the quantum Hall effect [13]. In mathematics, for example, the
difference relation arising in the solution of the Hofstadter model, known as
the Harper equation [1], is a special case of an "almost Mathieu operator",
which is one of the most studied types of ergodic Schrödinger operator [14,
15]. Moreover, in other branches of physics, the Hofstadter model has growing
relevance in a variety of subfields, including: cold atomic gases [16],
acoustic metamaterials [17], and photonics [18].

  3. 3.

Recent experimental realization. Although the Hofstadter model was introduced
last century [8, 1], it has only been experimentally realized within the last
decade. Signatures of the Hofstadter spectrum were first observed in moiré
materials [19] and optical flux lattices [20], and they have since been
reproduced in several other experimental platforms [16, 17, 18, 21]. Not only
does this spur recent theoretical interest, but it also increases the
likelihood of experimental groups entering the field, with the need for a
self-contained code repository that can be quickly applied to benchmark data
and related computations.

A prominent use-case of HofstadterTools is to facilitate the study of a rich
landscape of many-body problems. The Hofstadter model is an infinitely-
configurable topological flat-band model and hence, is a popular choice among
theorists studying strongly-correlated phenomena, such as the fractional
quantum Hall effect [22, 23] and superconductivity [24, 25]. Since there is a
relationship between the quantum geometry and topology of single-particle band
structures and the stability of exotic strongly-correlated states [26, 27, 28,
29, 30, 31], HofstadterTools may be used to guide theorists who are
researching quantum many-body systems. More broadly, we hope that
HofstadterTools will find many interdisciplinary applications, and we look
forward to expanding the package in these directions, with help from the
community.

Figure 1: Square Lattice. (a) Hofstadter butterfly and (b) Wannier diagram for
the Hofstadter model defined with nearest-neighbor hoppings on the square
lattice. (a) The energy $E$, and (b) the integrated density of states below
the gap $N(E)$, are plotted as a function of flux density
$n_{\phi}=BA_{\mathrm{min}}/\phi_{0}=p/499$, where $B$ is the perpendicular
field strength, $A_{\mathrm{min}}$ is the area of a minimal hopping plaquette,
$\phi_{0}$ is the flux quantum, and $p$ is an integer. The $r$-th gap is
colored with respect to $t=\sum_{i=0}^{r}C_{i}$, where $C_{i}$ is the Chern
number of band $i$. The size of the points in the Wannier diagram is
proportional to the size of the gaps. [12]

Figure 2: Triangular Lattice. (a) Hofstadter butterfly and (b) Wannier diagram
for the Hofstadter model defined with nearest-neighbor hoppings on the
triangular lattice. [32] Figure 3: Honeycomb Lattice. (a) Hofstadter butterfly
and (b) Wannier diagram for the Hofstadter model defined with nearest-neighbor
hoppings on the honeycomb lattice. [33]

Figure 4: Kagome Lattice. (a) Hofstadter butterfly and (b) Wannier diagram for
the Hofstadter model defined with nearest-neighbor hoppings on the kagome
lattice. [34]

## Acknowledgements

We thank Gunnar Möller, Titus Neupert, Rahul Roy, Alexey Soluyanov, Michael
Zaletel, Johannes Mitscherling, Daniel Parker, Stefan Divic, and Mathi Raja,
for useful discussions. This project was funded by the Swiss National Science
Foundation under Grant No.
https://data.snf.ch/grants/grant/203168https://doi.orgP500PT_203168, and
supported by the U.S. Department of Energy, Office of Science, Basic Energy
Sciences, under Early Career Award No. DE-SC0022716.

## References

  * [1] P. G. Harper.  Single band motion of conduction electrons in a uniform magnetic field.  Proceedings of the Physical Society. Section A, 68(10):874, Oct 1955\.
  * [2] M. Y. Azbel.  Energy spectrum of a conduction electron in a magnetic field.  Journal of Experimental and Theoretical Physics, 19(3):634–645, 1964.
  * [3] D. R. Hofstadter.  Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields.  Physical Review B, 14:2239–2249, Sep 1976.
  * [4] J. Lado.  pyqula.  Code repository at https://github.com/joselado/pyqula, 2021.
  * [5] N. Regnault.  DiagHam.  Code repository at http://www.nick-ux.org/diagham, 2001.
  * [6] J. Hauschild and F. Pollmann.  Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy).  SciPost Physics Lecture Notes, page 5, 2018.  Code available from https://github.com/tenpy/tenpy.
  * [7] D. Bodesheim, R. Biele, and G. Cuniberti.  Hierarchies of Hofstadter butterflies in 2D covalent organic frameworks.  npj 2D Materials and Applications, 7(1):16, Mar 2023.
  * [8] R. Peierls.  Zur Theorie des Diamagnetismus von Leitungselektronen.  Zeitschrift für Physik, 80(11):763–791, Nov 1933.
  * [9] N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman.  Light-induced gauge fields for ultracold atoms.  Reports on Progress in Physics, 77(12):126401, Nov 2014.
  * [10] A. Eckardt.  Colloquium: Atomic quantum gases in periodically driven optical lattices.  Reviews of Modern Physics, 89:011004, Mar 2017.
  * [11] A. Stegmaier, L. K. Upreti, R. Thomale, and I. Boettcher.  Universality of Hofstadter butterflies on hyperbolic lattices.  Physical Review Letters, 128:166402, Apr 2022.
  * [12] F. Di Colandrea, A. D’Errico, M. Maffei, H. M. Price, M. Lewenstein, L. Marrucci, F. Cardano, A. Dauphin, and P. Massignan.  Linking topological features of the Hofstadter model to optical diffraction figures.  New Journal of Physics, 24(1):013028, Jan 2022.
  * [13] J. E. Avron, D. Osadchy, and R. Seiler.  A Topological Look at the Quantum Hall Effect.  Physics Today, 56(8):38–42, Aug 2003.
  * [14] B. Simon.  Schrödinger operators in the twenty-first century, pages 283–288.  World Scientific, 2000.
  * [15] A. Avila and S. Jitomirskaya.  The Ten Martini Problem.  Annals of Mathematics, 170(1):303–342, 2009.
  * [16] N. R. Cooper, J. Dalibard, and I. B. Spielman.  Topological bands for ultracold atoms.  Reviews of Modern Physics, 91:015005, Mar 2019.
  * [17] X. Ni, K. Chen, M. Weiner, D. J. Apigo, C. Prodan, A. Alù, E. Prodan, and A. B. Khanikaev.  Observation of Hofstadter butterfly and topological edge states in reconfigurable quasi-periodic acoustic crystals.  Communications Physics, 2(1):55, Jun 2019.
  * [18] O. Zilberberg, S. Huang, J. Guglielmon, M. Wang, K. P. Chen, Y. E. Kraus, and M. C. Rechtsman.  Photonic topological boundary pumping as a probe of 4D quantum Hall physics.  Nature, 553(7686):59–62, Jan 2018.
  * [19] C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, T. Taniguchi, K. Watanabe, K. L. Shepard, J. Hone, and P. Kim.  Hofstadter’s butterfly and the fractal quantum Hall effect in moiré superlattices.  Nature, 497(7451):598–602, May 2013.
  * [20] M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch.  Realization of the Hofstadter Hamiltonian with ultracold atoms in optical lattices.  Physical Review Letters, 111:185301, Oct 2013.
  * [21] P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis.  Spectroscopic signatures of localization with interacting photons in superconducting qubits.  Science, 358(6367):1175–1179, 2017.
  * [22] B. Andrews and A. Soluyanov.  Fractional quantum Hall states for moiré superstructures in the Hofstadter regime.  Physical Review B, 101:235312, Jun 2020.
  * [23] B. Andrews, M. Mohan, and T. Neupert.  Abelian topological order of $\nu=2/5$ and $3/7$ fractional quantum Hall states in lattice models.  Physical Review B, 103:075132, Feb 2021.
  * [24] D. Shaffer, J. Wang, and L. H. Santos.  Theory of Hofstadter superconductors.  Physical Review B, 104:184501, Nov 2021.
  * [25] R. Sahay, S. Divic, D. E. Parker, T. Soejima, S. Anand, J. Hauschild, M. Aidelsburger, A. Vishwanath, S. Chatterjee, N. Y. Yao, and M. P. Zaletel.  Superconductivity in a topological lattice model with strong repulsion, 2023.  arXiv:2308.10935 [cond-mat.str-el].
  * [26] T. S. Jackson, G. Möller, and R. Roy.  Geometric stability of topological lattice phases.  Nature Communications, 6(1):8629, Nov 2015.
  * [27] B. Andrews, M. Raja, N. Mishra, M. Zaletel, and R. Roy.  Stability of fractional Chern insulators with a non-Landau level continuum limit, 2023.  arXiv:2310.05758 [cond-mat.str-el].
  * [28] P. J. Ledwith, A. Vishwanath, and D. E. Parker.  Vortexability: A unifying criterion for ideal fractional Chern insulators, 2022.  arXiv:2209.15023 [cond-mat.str-el].
  * [29] C. H. Lee, M. Claassen, and R. Thomale.  Band structure engineering of ideal fractional Chern insulators.  Physical Review B, 96:165150, Oct 2017.
  * [30] H. Tian, X. Gao, Y. Zhang, S. Che, T. Xu, P. Cheung, K. Watanabe, T. Taniguchi, M. Randeria, F. Zhang, C. N. Lau, and M. W. Bockrath.  Evidence for Dirac flat band superconductivity enabled by quantum geometry.  Nature, 614(7948):440–444, Feb 2023.
  * [31] J. Wang, J. Cano, A. J. Millis, Z. Liu, and B. Yang.  Exact Landau level description of geometry and interaction in a flatband.  Physical Review Letters, 127:246403, Dec 2021.
  * [32] J. E. Avron, O. Kenneth, and G. Yehoshua.  A study of the ambiguity in the solutions to the Diophantine equation for Chern numbers.  Journal of Physics A: Mathematical and Theoretical, 47(18):185202, Apr 2014.
  * [33] A. Agazzi, J.-P. Eckmann, and G. M. Graf.  The colored Hofstadter butterfly for the honeycomb lattice.  Journal of Statistical Physics, 156(3):417–426, Aug 2014.
  * [34] H. Jing-Min.  Light-induced Hofstadter’s butterfly spectrum of ultracold atoms on the two-dimensional kagomé lattice.  Chinese Physics Letters, 26(12):123701, Dec 2009.