arXiv CSopen access

On the algebra of Koopman eigenfunctions and on some of their infinities

On the algebra of Koopman eigenfunctions and on some of their infinities

arXiv:2604.21825v1 [math.DS] 23 Apr 2026

Zahra Monfared1,2,3 , Saksham Malhotra3 , Sekiya Hajime3 , Ioannis Kevrekidis4 , Felix Dietrich3* 1

Interdisciplinary Center for Scientific Computing, University of Heidelberg, Germany. 2 Department of Mathematics and Computer Science, University of Heidelberg, Germany. 3 School of Computation, Information and Technology, Technical University of Munich, Germany & MDSI & MCML. 4 Departments of Chemical and Biomolecular Engineering and of Applied Mathematics and Statistics, Johns Hopkins University, Baltimore, USA. *Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected]; [email protected]; [email protected]; [email protected]; Abstract For continuous-time dynamical systems with reversible trajectories, the nowherevanishing eigenfunctions of the Koopman operator of the system form a multiplicative group. Here, we exploit this property to accelerate the systematic numerical computation of the eigenspaces of the operator. Given a small set of (so-called “principal”) eigenfunctions that are approximated conventionally, we can obtain a much larger set by constructing polynomials of the principal eigenfunctions. This enriches the set, and thus allows us to more accurately represent application-specific observables. Often, eigenfunctions exhibit localized singularities (e.g. in simple, one-dimensional problems with multiple steady states) or extended ones (e.g. in simple, two-dimensional problems possessing a limit cycle, or a separatrix); we discuss eigenfunction matching/continuation across such singularities. By handling eigenfunction singularities and enabling their continuation, our approach supports learning consistent global representations from locally sampled data. This is particularly relevant for multistable systems and applications with sparse or fragmented measurements.

1

Keywords: Nonlinear dynamical systems, Koopman operator, Koopman eigenfunctions, Singularities, Data-driven modeling.

1 Introduction Much research on Koopman operator approximation for dynamical systems focuses on single basins of attraction around fixed points [1–3], something that is exploited, at least locally, in the celebrated Hartman-Grobman theorem [4, 5]; for more recent results extending this to the entire basin of attraction, see [6, 7]. This is not surprising, because the spectrum of the Koopman operator around hyperbolic fixed points can be characterized using the relation to the local linearization around the fixed point [8]. However, even seemingly simple systems in one dimension with multiple steady states have complicated Koopman eigenfunction structure. Consider, for example, the following system defined through an ordinary differential equation with a polynomial vector field, ẋ = (x − a)(x − b)(x − c), x ∈ R. (1)

Figure 1 shows the three steady states of the system (Equation (1) for the case a = −1, b = 0, c = 3): It also shows the vector field itself, as well as three Koopman eigenfunctions, whose selection we now discuss. Let us concentrate at the interval

Fig. 1: Three Koopman eigenfunctions ϕ1 , ϕ2 , ϕ3 for a system with vector field ẋ = (x − a)(x − b)(x − c), with a = −1, b = 0, and c = 3. The color of the dashed lines at the steady states indicates where the correspondingly colored eigenfunction is zero.

between a and b. Using the linearization around steady state a, and solving the Koopman PDE ⟨∇ϕ, ẋ⟩ = λϕ, results in the blue eigenfunction ϕ1 ; we can pin, without loss of generality, the derivative of ϕ1 at x = a to an arbitrary value; this is just a scaling to select a single eigenfunction among all acceptable scalings. This eigenfunction helps us describe trajectories of the dynamical system between a and b. The same trajectories between a and b can also be described using the eigenfunction ϕ2 , similarly

2

constructed by focusing on the steady state b and its linearization eigenvalue. Cruλ /λ cially, the two eigenfunctions are simple transformations of each other: ϕ1 = ϕ2 1 2 . This is not surprising, as we will discuss below (Proposition 1): Powers of Koopman eigenfunctions are also Koopman eigenfunctions (a well-known fact, see [8]). The second important point to notice is that ϕ1 asymptotes to infinity at b, while ϕ2 asymptotes to infinity at a; so, the power relation between the two should only be valid in the open interval (a, b), where they are both finite (for details we refer to the concept of “open” [6] or, similarly, “primary” eigenfunctions [9]). Clearly, ϕ1 is also defined in (−∞, a), while ϕ2 is also defined in (b, c). Using the power relation λ /λ ϕ1 = ϕ2 1 2 we can now extend (in the spirit of analytic continuation) ϕ2 even to the left of a, and ϕ1 even to the right of b, up to c, where ϕ2 asymptotes to infinity. The same process allows us to realize that ϕ3 (computed based on the eigenvalue of the steady state c and the solution of the Koopman PDE in the interval (b, c)) can be transformed in that same interval to ϕ2 , using the ratio of λ2 and λ3 as the requisite power. And, finally, the same extension argument can be used to extend each of these eigenfunctions over the entire real line, with the possible exception of the singularities at the various steady states. In a sense, we can talk about a single, unique Koopman eigenfunction over the entire domain. The power transformation we discussed generates the other ones, modulo the (occasional) point-wise singularities. Both the items we discussed, the power transformation and the singularities, will naturally play a role in computing the eigenfunctions; that is the subject of the paper. In the first part of the paper we focus on the power relation between eigenfunctions (separating it from the singularities) and only discuss Proposition 1 and how it can be used to generate new eigenfunctions. The paper is organized as follows. In Section 2, we motivate how the group structure of the Koopman eigenfunctions can be used to systematize and speed up eigencomputations. In Section 3, we provide an overview of the Koopman operator framework and review various existing approximation algorithms. Section 4 focuses on the methodology for extending the computed set of eigenfunctions, with a detailed analysis presented in Section 4.1. Then, in Section 5, we illustrate our approach through several numerical examples, showcasing the construction of an expanded set of eigenfunctions. We also discuss how to approximate eigenfunctions with singularities, which is necessary for, e.g., systems with transients to multiple attractors. We further discuss how to consistently extend them across the singularities and how, when possible, to transform between them despite the singularities. Finally, in Section 6, we explore potential extensions and use cases of our work.

2 Computation of Koopman eigenfunctions Numerical computation of the spectral decomposition of the Koopman operator from observation data is an important yet challenging problem [10–13]. The main goal of Koopman operator approximation methods for dynamical systems with pure point spectrum is the accurate representation of a sufficiently large set of eigenfunctions. In the ideal case, the observables of interest lie in the span of these eigenfunctions.

3

Then, prediction of future values of these observables is trivial and only amounts to multiplication. The notion of “principal” eigenfunctions has been proposed [9, 14]: a minimal set of functions, such that the entire eigenspace of the operator can be constructed by systematically combining them (e.g., taking integer-valued powers). A prominent approximation algorithm, mainly for large-scale linear systems, is the Dynamic Mode Decomposition [15] (DMD). Its extension [10] (EDMD) approximates the spectral properties from the action of the operator on a larger set of functions. Many variants of these algorithms are available, see for example the recent work for ergodic systems [16]. Yet, most of these algorithms only employ the linearity of the operator in order to approximate a matrix representation of it: its own approximate Koopman matrix. This means that the chosen dictionary must be expressive enough to represent many Koopman eigenfunctions, so that the action of the operator on the dictionary does not deviate too much from the finite-dimensional dictionary subspace. In this paper, we propose to employ another property of the Koopman operator, which is not common to all linear operators: the group structure of its eigenfunctions [17]. Given any sufficiently well-approximated eigenfunction of the operator, (even negative) integer powers of the function are also eigenfunctions. This is true even if these powers cannot be represented in the original dictionary. It is therefore possible to derive eigenfunctions that cannot be well approximated using the initial Koopman matrix representation. As a motivating example (from [3]), consider a nonlinear ODE system of the form

(

x˙1 = ax1 ,

x˙2 = b x2 − x21 .

(2)

By an appropriate choice of observable functions, this system can be embedded in a higher-dimensional space where the dynamics become linear, although for a general nonlinear system it is not always possible to achieve linearization in finite dimensions. For instance, for y = (y1 , y2 , y3 ) = (x1 , x2 , x21 ), (2) can be represented as a 3D linear system that remains closed under the action of the Koopman operator. However, to better illustrate our point with this example, we expand it using more observable functions y = (y1 , y2 , y3 , y4 , y5 ) = (x1 , x2 , x21 , x1 x2 , x31 ) for which system (2) will be transformed into a 5D linear system:   y˙1 = ay1    
  y˙2 = b y2 − y3    . (3) y˙3 = 2ay3      y˙4 = (a + b)y4 − by5     y˙ = 3ay 5 5

4

The flow map F ∆t : M −→ M is given by   a0 0 0 0  0 b −b 0  0  
 ∆t y 0 = exp A ∆t y 0 , 0 0 2 a 0 0 F ∆t (y 0 ) = exp    0 0 0 a + b −b 0 0 0 0 3a

(4)

where y 0 = y (0). The matrix A represents the finite-dimensional Koopman generator associated with the chosen observable vector y . Indeed, the observables evolve according to the linear system

ẏ = A y,

(5)

The discrete-time Koopman operator over a sampling interval ∆t is then given by

K = eA∆t .

(6)

If w is a left eigenvector of A satisfying

wT A = λwT ,

(7)

ϕ(y ) = ⟨y, w⟩ = wT y

(8)

d ϕ(y ) = ⟨∇ϕ, ẏ⟩ = wT Ay = λϕ(y ). dt

(9)

then the function

satisfies

Therefore, ϕ(y ) is a Koopman eigenfunction associated with the eigenvalue λ. The left eigenvectors of the matrix A are

w1 = (1 0 0 0 0)T w2 = (0 0 1 0 0)T w3 = (0 0 0 0 1)T 2a − b w4 = (0 1 0 0)T b 2a − b w5 = (0 0 0 1)T , b and hence the eigenfunctions of the associated Koopman operator are

ϕ1 (y ) = ⟨y, w1 ⟩ = y1 , 5

(10)

ϕ2 (y ) = ⟨y, w2 ⟩ = y3 = y12 = ϕ21 (y ) ϕ3 (y ) = ⟨y, w3 ⟩ = y5 = y13 = ϕ31 (y ) ϕ4 (y ) = ⟨y, w4 ⟩ =

2a − b 2a − b y2 + y3 = y2 + y12 b b

ϕ5 (y ) = ⟨y, w5 ⟩ =

2a − b 2a − b y4 + y5 = y1 y2 + y13 . b b

(11)

Thus, the DMD algorithm creates a 5 × 5 Koopman matrix where the eigenfunction ϕ2 is the square of ϕ1 and ϕ3 is the cube of ϕ1 . However, the squares of the functions ϕ2 , ϕ3 , ϕ4 , ϕ5 are also eigenfunctions of the Koopman operator; yet, they cannot be derived from the given Koopman matrix (see Figure 2). Similarly, higher powers of the eigenfunctions cannot be derived from this matrix either.

Fig. 2: Our approach allows for the computation of eigenfunctions that cannot be derived directly from the given Koopman matrix we started with. For system (2), the eigenfunctions ϕ1 , ϕ2 , ϕ3 , ϕ4 , ϕ5 can be obtained from the Koopman matrix A. While the second and third powers of ϕ1 are also eigenfunctions that can be derived from this matrix, the powers of ϕ2 , ϕ3 , ϕ4 , and ϕ5 cannot be obtained directly from it. Therefore, our approach enables the derivation of additional eigenfunctions by utilizing higher powers (p ≥ 4) of ϕ1 and by considering powers p ≥ 2 of ϕ2 , ϕ3 , ϕ4 , and ϕ5 . Blue lines represent eigenfunctions obtained from the Koopman matrix, while orange lines correspond to integer powers of eigenfunctions. Solid lines indicate eigenfunctions derived directly from the Koopman matrix, whereas dashed lines denote those that are not.

6

3 Related work We now briefly discuss research related to the algebra of Koopman eigenfunctions, and their extension across singularities. The Koopman operator is often studied in the ergodic setting, which we introduce in Section A. The considerations on singular eigenfunctions concern the approach to attractors and not the attractors themselves (on which the dynamics can be ergodic), and so we require a more general setting. In general, a dynamical system is defined by a set M called the state space or the phase space and a map F : M → M called the flow (or evolution function). Typically, one considers the case where M is a measurable space, with a σ -algebra B, and F is B-measurable (for more information see [8]).

3.1 Koopman operator and its multiplicative property Let M be a B-measurable set and F = {f : M → C} a space of measurable, complexvalued functions.

Definition 1. (Koopman operator (discrete-time)); Consider a discrete-time dynamical system defined by the nonlinear (non-singular) map F : M → M s.t. x(n + 1) = F (x(n)), n ∈ N, x(1) ∈ M. The Koopman operator KF : F → F associated with the map F : M → M acts on observables g ∈ F and is defined through the composition [KF g ](x(n)) = (g ◦ F )(x(n)).

(12)

The composition on the left is a linear operation, so KF is a linear operator on F and can be spectrally decomposed.

Definition 2. (Koopman operator (continuous-time)); Consider a continuous-time dynamical system ẋ = F (x), x ∈ M, described by the one-parameter family of (nonsingular) flow maps F t : M → M, t ∈ R+ . The family of Koopman operators KFt t : F → F associated with the family of flow maps F t acts on observables g ∈ F and is defined as [KFt t g ](x) = (g ◦ F t )(x).

(13)

Definition 3. (Koopman eigenfunction and eigenvalue (discrete-time)); An eigenfunction of the Koopman operator KF associated with the map F : M → M is a nonzero observable ϕk ∈ F \ {0} such that KF ϕk = ϕk ◦ F = λk ϕk ,

(14)

where λk ∈ C is the corresponding eigenvalue.

Definition 4. (Koopman eigenfunction and eigenvalue (continuous-time)); An eigenfunction of the Koopman operator KFt t associated with the semigroup of flow maps

7

(F t )t≥0 is a nonzero observable ϕk ∈ F \ {0} such that

KFt t ϕk = ϕk ◦ F t = eλk t ϕk

∀t ≥ 0,

(15)

where λk ∈ C and eλk t is the corresponding eigenvalue.

If the semigroup of operators KFt t is strongly continuous, Koopman eigenfunctions and eigenvalues are characterized by the relation Lϕk = λk ϕk . In the case where the semigroup (F t )t≥0 arises from the flow of the system ẋ = F (x), they are obtained by solving the following equation [2]

F · ∇ϕk = λk ϕk .

(16)

A key property of the Koopman operator, central to this work, is its multiplicative structure: Proposition 1 Products and (integer) powers of eigenfunctions are also eigenfunctions of m1 m2 m1 m2 1 m2 KF , if the combination is in F: KF [ϕm k1 ϕk2 ] = λk1 λk2 [ϕk1 ϕk2 ], with k1 , k2 ∈ N and t m1 , m2 ∈ Z. The same property holds for KF t , corresponding to continuous-time systems. Proof Straightforward application of the definition of the Koopman operator.

□

Note that Proposition 1 requires that the powers of the eigenfunctions must be an element of the function space F . This means for negative m1 or m2 the associated eigenfunctions must not be zero on their domain. In the paper, whenever |m1 | > 1 or |m2 | > 1, we will call the corresponding functions monomial eigenfunctions to distinguish them from principal eigenfunctions. This characteristic imposes a lattice or group structure on the Koopman operator’s spectrum [2, 8, 18]. Building on this property, in [18] the authors introduced Multiplicative Dynamic Mode Decomposition (MultDMD) to incorporate the Koopman operator’s multiplicative structure into its finite-dimensional approximation. They developed a specialized dictionary of basis functions and an efficient optimization algorithm to maintain this structure, ensuring that the computed eigenfunctions include powers. Instead, in this work, we consider an arbitrary dictionary and generate a large set of eigenfunctions by taking their powers, either by multiplying them with themselves or with other eigenfunctions.

3.2 Extending eigenfunctions beyond infinity To create the powers of eigenfunctions in order to generate a larger set of them, we need a small set of initial eigenfunctions. However, approximating these initial eigenfunctions can be quite challenging in some cases. For instance, in certain systems, eigenfunctions may inherently asymptote to ± infinity on specific subsets of the entire space. This can occur, for example, for systems with multiple steady states or limit cycles [1, 19], where eigenfunctions may diverge at the boundaries between basins of

8

attraction [20, 21]. Therefore, it is crucial to develop effective techniques for extending (analytically continuing) eigenfunctions “beyond infinity”. We will explore how to extend eigenfunctions across singularities for different examples.

3.3 Existing algorithms for the approximation of the Koopman operator Dynamic Mode Decomposition (DMD), which was originally introduced by Schmid and Sesterhenn [15], can be seen as a variant of a standard Arnoldi method [22]. It can be used to approximate the Koopman operator for linear systems [23]. The approximations can also be useful for nonlinear systems, such as fluid flows [24], if a sufficiently rich set of nonlinear observations of the state of nonlinear system is provided. The latter was formalized as Extended Dynamic Mode Decomposition (EDMD) by Williams, Kevrekidis, and Rowley [10], through the choice of a proper truncated basis of the function space the operator acts on. Li et al. [11] use neural networks to construct a flexible, problem-dependent dictionary for this purpose. For measurepreserving dynamical systems, the mpEDMD algorithm as introduced by Colbrook [16] constrains the spectrum of the approximating matrix to lie on the unit circle.

3.4 General eigensolvers: Power method and QR algorithm We briefly outline the power method, based on Stoer and Bulirsch [25]. See Also Section B. Write the eigenpair of A ∈ Cn×n Pnas (λ1 , v 1 ), (λ2 , v 2 ), . . . (λn , v n ) and assume |λ1 | > |λ2 | ≥ · · · ≥ |λn |. Then for q = i=1 ci v i ∈ Cn , it is easy to see that

m n  n X X 1 λi Am vi m A q= = ci v i → c1 v 1 (m → ∞). ci |λ1 |m |λ1 |m |λ1 | i=1 i=1 Note that



λi |λ1 |

m

|λi | < 1. Therefore, the convergence of this → 0 if i ̸= 1 due to |λ 1|

2| algorithm for any given q is linear with respect to the ratio |λ |λ1 | . The second important algorithm is the QR algorithm. For the initial matrix A1 , define a matrix A2 as A2 = R1 Q1 where A1 = Q1 R1 is the QR decomposition computed by a standard algorithm. By definition, A2 = R1 Q1 = Q∗1 A1 Q1 , and as Q is unitary, A1 and A2 are unitary similar. Thus, they have the same eigenvalues. By repeating this procedure,

An+1 = Pn∗ A1 Pn , where Pn = Q1 Q2 . . . Qn . Since Pn is a product of unitary matrices, it is again a unitary matrix. Hence, An+1 and A1 are again unitary similar and have the same eigenvalues. An converges to an upper triangular matrix as n → ∞.

9

4 Mathematical framework 4.1 Numerical algorithms for Koopman operator approximation In this section, we discuss how the algebra of Koopman operator eigenfunctions can enable their computation, enhancing traditional eigensolvers.

4.1.1 Constructing eigenfunctions with integer exponents Once a single eigenfunction of the Koopman operator has been approximated, the multiplicative property from Proposition 1 can be used. If ϕ̄1 and ϕ̄2 are eigenfunction approximations of the Koopman operator KF with eigenvalues λ̄1 and λ̄2 , then

ϕ̄pq = ϕ̄p1 ϕ̄q2 , p, q ∈ Z

(17)

4 2 0 2 4

(x) 1

(x)

is also an eigenfunction approximation of KF with eigenvalue λ̄pq = λ̄p λ̄q . It is interesting at this point to consider negative powers p, q . Figure 3 shows how a simple eigenfunction ϕ(x) = x will lead to a singular eigenfunction if the power q = −1 is used. The singular transformation (·)q induced by taking this negative power also provides the simplest demonstration of how continuations across infinity may be rationalized.

2.5

0.0 x

2.5

2.5

0.0 x

2.5

Fig. 3: Left: eigenfunction ϕ(x) = x for a linear system ẋ = −x. Right: inverse of the same eigenfunction, with singularity at the steady state marked in orange.

In the following, we assume access to data within the basin of attraction of a hyperbolic fixed point of a (possibly nonlinear) dynamical system. In later sections, we will discuss examples with multiple basins, and discuss the related numerical pathologies. In the single basin case, assume we used EDMD with dictionary Ψ to construct the approximating Koopman matrix K ∈ Rd×d , with left eigenvector {wi }di=1 and eigenvalues {λ̄i }di=1 . Given i, j ∈ {1, . . . , d} and p, q ∈ N ∪ {0}, the extended set of eigenfunctions is given by

ϕ̄pq (x) = (wiT Ψ(x))p (wjT Ψ(x))q ,

10

(18)

with corresponding eigenvalues λ̄pq = λ̄pi λ̄qj . For simplicity, we will focus our analysis only on extending powers of known eigenfunctions and not those of combinations of known eigenfunctions (i.e., we set q = 0).

4.1.2 Constructing eigenfunctions with real-valued exponents Up to now, we considered integer exponents q, p in constructing new eigenfunctions. Yet, the procedure (also Proposition 1) applies to the case of real-valued exponents, and we will take advantage of this later below. The exponential of the complex number x is defined by (ex )y := ey log(e ) , x, y ∈ C, where the complex logarithm of non-zero complex number z , denoted as log(z ), is defined by elog(z) := z ∈ C. Note that if z is given by the polar form z = reiθ with r > 0 and θ ∈ R, then the complex logarithm is of the form log(z ) = ln(r) + i(θ + 2πk ), k ∈ Z,

where ln(·) is the natural logarithm, i.e., ln(·) := loge (·). Also remember that the principal value of the log(z ) is defined as the logarithm whose imaginary part lies in the (−π, π ], i.e. the principal value is ln(r)+ iθ′ such that θ′ = θ +2πk ∈ (−π, π ], k ∈ Z. In general, log(z ) is taken to denote the principal value. Now, assume there exists an eigenpair (λ, φ) with |λ| = 1, λ ̸= 1 + 0i. Then, one can represent it as λ = eci where c ∈ (−π, π ]. For any λ′ = edi ∈ S 1 with d ∈ (−π, π ], one can consider the power p := dc ∈ R so that d

λp = ep log(λ) = eipc = ei c c = eid = λ′ . Thus, by Proposition 1, (λ′ , φp ) is also an eigenpair.

4.1.3 Matching eigenfunctions across steady states Proposition 1 provides the mathematical basis for understanding the nature of computationally identified, “complicated” monomial eigenfunctions and how they arise from multiplicative combinations of principal ones. Using the logarithm helps reduce the multiplicative structure and identify the principal eigenfunctions that synthesize it. This is because for all x ∈ M where ϕk1 (x) ̸= 0 and ϕk2 (x) ̸= 0, m2 1 log |ϕm k1 (x)ϕk2 (x)| = m1 log |ϕk1 (x)| + m2 log |ϕk2 (x)|.

Now, consider a (finite) collection of known eigenfunctions {ϕk }, with their domain in a neighborhood B of a particular steady state of F , s.t. |ϕk (x)| > 0 ∀x ∈ B . Note that this is generally the case, because eigenfunctions are typically non-zero away from steady states (unless they are identically zero on a large portion of the state space). Then, we can systematically filter the collection by removing linear subspaces from it, only leaving the “principal spectrum” [14]. We will see below that computing logarithms will allow us to extend the domain of eigenfunctions far beyond the neighborhood B of the steady state around which they were originally computed; Remarkably, this extension can even “jump across” multiple singularities.

11

4.1.4 Algorithms to generate Koopman operator eigenfunctions Given a left eigenvalue λ and eigenvector ϕ of the Koopman operator matrix, we can define an algorithm to find monomial eigenfunctions ϕp (i.e., ϕp ≡ ϕp ) with prescribed trajectory error/upper bounds. For a discrete system, there is no time-integration error, so only the error δv due to eigenvector approximations plays a role. In this case, given ϵ > 0, to get EF G (ϕ̄p , λ̄p ) ≤ ϵ, we require

∥δv∥ ≤

ϵp . CF G (p, λ̄)

(19)

Using this relation, Algorithm 1 shows how to take powers of eigenpairs of the Koopman operator for a discrete time system.

4.1.5 Error metric for generated eigenfunctions in two dimensions Taking powers of eigenfunction approximations accumulates errors. Therefore, we need to define an error function for a monomial eigenfunction. To measure this error in examples with two spatial dimensions, we define a grid G on our domain Ω ⊆ R2

 G = (a + nh, b + mh) ∈ Ω | a, b ∈ R; n, m ∈ Z; h ∈ (0, 1) ,

(20)

and a norm ∥·∥G on the grid

s ∥ϕ∥G =

1 X (ϕ(x))2 . |G|

(21)

x∈G

For systems with higher-dimensional states x, we extend this notion accordingly, to higher-dimensional grids of equidistant nodes. We introduce the trajectory metric (EF G ) in Definition 5.

Definition 5. (Trajectory metric) Given an extended eigenfunction approximation ϕ̄p with an extended eigenvalue approximation λ̄p for the Koopman operator, the trajectory error of the eigenpair approximation (λ̄p , ϕ̄p ) is defined, for discrete systems, by (1/p) EF G (λ̄p , ϕ̄p ) = ϕ̄p (F (·)) − λ̄p ϕ̄p (·) G , (22) and for continuous systems by (1/p)

EF G (λ̄p , ϕ̄p ) = ϕ̄p (F ∆t (·)) − λ̄p ϕ̄p (·) G

.

(23)

For a continuous system, we also have an integration error in the flow calculation. Given ϵ > 0, to obtain EF G (ϕ̄p , λ̄p ) ≤ ϵ, we require

ϵG ≤

1 L



ϵp + (λ̄M )p

12

1/p

 − λ̄M .

(24)

Algorithm 1 Computing extended eigenpairs for discrete system. Given ∥δw∥, left eigenpair (λ̄, wc ) of the Koopman matrix K , Ψ as the dictionary basis and desired trajectory error bound ϵ p←1 while p ∈ N do Pp−1 p−1−i i CF G (p, λ̄) ← Ψ(F (·)) − λ̄Ψ(·) ∥Ψ(·)∥ λ̄i i=0 ∥Ψ(F (·))∥ G ϵp then if ∥δw∥ > CF G (p, λ̄) break end if ϕ̄p ← (wcT Ψ)p λ̄p = (λ̄)p p←p+1 end while Using this relation, Algorithm 2 outlines how to extend eigenpairs of the Koopman operator for a continuous time system.

Algorithm 2 Computing extended eigenpairs for a continuous time system. Given ϵG , left eigenpair (λ̄, wc ) of the Koopman matrix K and ψ as the dictionary basis and desired trajectory error bound ϵ p←1 while p ∈ N do 1 if ϵG > ((ϵp + (λ̄M )p )1/p − λ̄M ) then L break end if ϕ̄p ← (wcT ψ )p λ̄p = (λ̄)p p←p+1 end while

4.1.6 Constructing all admissible eigenvectors for continuous time systems Given the principal eigenfunctions, Algorithm 3 describes how to construct all monomial eigenfunctions and corresponding eigenvalues of a continuous time system that have a bounded error metric.

4.2 Theoretical analysis The two main sources of error in eigenfunction approximation are (a) the error in the eigenvector of the Koopman matrix K due to the eigensolver, and (b) the error in

13

Algorithm 3 (Iterative Koopman eigensolver) Algorithm for computing extending eigenvalues λ̄ip and extended eigenfunctions ϕ̄ip of a continuous time system, given integration error ϵG , desired trajectory error ϵ, the Koopman matrix K , Ψ as the dictionary basis and constant L i←0 A←K while i < n do λ̄i , vi ← power iteration complex(A) λ̄i , wi ← power iteration complex(AT ) p←1 while true do  1 if ϵG > (ϵp + (λ̄i M )p )1/p − λ̄i M then L break end if ϕ̄ip ← (wiT Ψ)p λ̄ip ← (λ̄i )p p←p+1 end while wi wi ← T wi vi A ← A − λ̄i vi wiT if imag(λ̄i ) > 10−6 then vi+1 ← v̄i wi+1 ← w̄i ¯ λ̄i+1 ← λ̄ i wi+1 wi+1 ← T wi+1 vi+1 T A ← A − λ̄i+1 vi+1 wi+1 i←i+1 end if i←i+1 end while the integration for the flow computation in continuous systems. Eigenvector error is present for both discrete and continuous system computations. We can estimate upper bounds for the trajectory error EF G with respect to these errors.

4.2.1 Integration error in continuous systems For a continuous time system, the trajectory error will depend on the error introduced while integrating the system to get the flow F ∆t . Let

F ∆t (x) − x∆t = ε(x),

14

(25)

where ϵ(x) ∈ Rn is the integration error at x, x∆t is the accurate flow at t = ∆t, and F t (x) is the computed flow. We consider the upper bound of the trajectory error with respect to the quantity ϵG = max ∥ε(x)∥ . (26) x∈G

Using the above, Proposition 2 gives the upper bound for eigenfunction approximation ϕ̄p and eigenvalue approximation λ̄p with respect to the integration error ϵG . Proposition 2 EF G (ϕ̄p , λ̄p ) ≤



λ̄∆t M + LϵG

p

− (λ̄∆t M )p

(1/p)

.

(27)

where M = maxx∈G ∥ψ(x)∥ and L is an upper bound on the spectral norm of the Jacobian of Ψ with respect to the l2 norm on the grid G, ∥JΨ (x)∥2 ≤ L. (28) Proof See Section C.1

□

Remark 1. To keep the trajectory error below ϵ for a power p, we will require that the integration error ϵG has the bound ϵG ≤

1 p [(ϵ + (λ̄∆t M )p )1/p − λ̄∆t M ]. L

(29)

Remark 2. To calculate L, we compute the maximum singular value of the matrix JΨ (x) over the grid G and take the maximum value σmax (x) =

q

λmax JΨ (x)T JΨ (x) ,

L = max σmax (x), x∈G

L ≥ ∥JΨ (x)∥2 ∀x ∈ G.

(30)

4.2.2 Eigenvector approximation error Let w be a left eigenvector of the Koopman matrix K ∈ Rd×d . Consider a left eigenvector wc ∈ Rd of K computed by an eigensolver. Let

wc = w + δw, where δw ∈ Rd is the error vector introduced due to the eigensolver. We can obtain an a posteriori upper bound on the trajectory error in (22) for the pth -power eigenfunction approximation

ϕ̄p (x) = (wTc Ψ(x))p .

(31)

We assume that the true left eigenvector is normalized so that ∥w∥ = 1. Assuming that λ̄ is the computed eigenvalue corresponding to wc and λ̄p = λ̄p , Proposition 3 gives the error bound for the trajectory error with respect to the eigenvector error ∥δw∥. 15

Proposition 3 EF G (ϕ̄p , λ̄p ) ≤ CF G (p, λ̄)(1/p) ∥δw∥(1/p) ,

(32)

where CF G (p, λ̄) =

Ψ(F (·)) − λ̄Ψ(·)

p−1 X i=0

∥Ψ(F (·))∥p−1−i ∥Ψ(·)∥i λ̄i

.

(33)

G

Proof See Section C.2

□

Remark 3. To keep the trajectory error below ϵ for a power p we will require that the error in computed eigenvector δw has norm such that ∥δw∥ ≤

ϵp . CF G (p, λ̄)

(34)

5 Computational experiments We now demonstrate the construction of additional (monomial) eigenfunctions in a series of computational experiments. In Section 5.1, we discuss how Koopman eigenfunctions of a linear system are constructed up to a certain accuracy. Note that even though the system matrix of the linear system is finite-dimensional and thus only has a finite number of eigenvectors, the Koopman operator of the system still acts on an infinite-dimensional space and has an infinite number of eigenfunctions. The benefit of analyzing it for a linear system is that all eigenfunctions are available analytically, so we can compute the approximation error of our numerical procedure to the ground truth. In order to validate the procedure for nonlinear systems, in Section 5.2, we construct an explicit, nonlinear example by transforming the state space of the linear system in Section 5.1 with a nonlinear diffeomorphism. This allows us to obtain analytic formulas for Koopman eigenfunctions even in this nonlinear setting. Section 5.3.1 discusses what happens in general for systems with separatrices. Section 5.4.1 shows the concept of isochrons for limit cycles. Section 5.5 shows how the concept of isochrons connects with/relates to the computation of eigenfunctions for systems with separatrices associated with saddle points.

5.1 Example - linear system in 2D, non ergodic Consider the continuous time linear system

ẋ = Ax,

(35)

where x ∈ R2 , A ∈ R2×2 with left eigenpairs (λ1 , w1 ) and (λ2 , w2 ). Let the system be sampled with fixed sampling interval ∆t. Then the Koopman eigenfunctions and eigenvalues of the system are given by

p
q λp = (eλ1 ∆t )p , ϕp (x) = ⟨w1 , x⟩ , λq = (eλ2 ∆t )q , ϕq (x) = ⟨w2 , x⟩ . 16

(36)

If λ̄ is a computed eigenvalue of the Koopman matrix K and λ is an eigenvalue of the Koopman generator of the continuous system, then λ̄ ≈ eλ∆t where ∆t is the temporal sampling interval for the system. We can compare the numerically constructed monomial eigenfunctions with the true eigenfunctions of the system. Let ϕ̄p be a constructed, monomial eigenfunction and ϕp be the true eigenfunction. Then on a grid G as defined in (20), the error given by EG is defined similarly to the discrete system:

G0 = {x ∈ G| |ϕp | < ε} ∪ {x ∈ G| |ϕ̄p | < ε}   ϕp|G/G0 cmode = mode ϕ̄p|G/G0 1/p

EG (ϕp , ϕ̄p ) = ϕp − cmode ϕ̄p G ,

(37) (38) (39)

where ε is some tolerance. The trajectory error for the computed eigenpair (λ̄p , ϕ̄p ) is given by 1/p EF G (λ̄p , ϕ̄p ) := ϕ̄p (F ∆t (·)) − λ̄p ϕ̄p (·) G . (40) Consider the case

 A=

 −0.9 0.1 . 0 −0.8

(41)

The matrix has left eigenpairs (−0.9, [1, −1]T ) and (−0.8, [0, 1]T ). The system is sampled with sampling interval ∆t. Therefore, the true Koopman eigenfunctions of the system are given by

x1 − x 2 √ ϕp (x) = 2 ϕq (x) = xq2 , 

p (42) (43)

with eigenvalues

λp = (e−0.9∆t )p

(44)

λq = (e−0.8∆t )q .

(45)

To approximate the eigenfunctions using DMD we collect 400 snapshot pairs (with the state variables as our Koopman observables), where the initial conditions are uniformly randomly distributed between [−2, 2] × [−2, 2], with temporal sampling interval ∆t = 0.2. Figure 4 shows the results of the DMD approximation. The out-ofsample prediction shows that the DMD approximation is able to predict trajectories accurately. We now create a grid G with a = −1, b = 1, n = 100, h = 0.01. On G, we calculate the flow approximated by the forward Euler method using step size h = 0.001. Then, we calculate ϵG using (26) and the true solution, x∆t = eA∆t x, and approximated

17

DMD model reconstruction (Identity state dictionary)

2.0

2.0

1.5

1.5

1.0

1.0

0.5

0.5

x2

x2

Training data used during fit

0.0

0.0

−0.5

−0.5

−1.0

−1.0

−1.5

−1.5

−2.0

−2.0 −2

−1

0

1

−2

2

x1

−1

out-of-sample prediction (x1 ) 2.00

0

1

2

x1

out-of-sample prediction (x2 ) 1.0

DMD true system

DMD true system

1.75 0.8 1.50 0.6

x2

x1

1.25 1.00

0.4

0.75 0.50

0.2 0.25 0.00

0.0 0

1

2

3

4

5

6

7

0

t

1

2

3

4

5

6

7

t

Fig. 4: DMD model for the continuous linear system (35). Top left shows the training data used for approximation, top right shows the reconstructed data, bottom shows the comparison between predicted trajectory and actual trajectory for a sample initial condition. solution using the forward Euler method,

x(0) = x, N =

∆t , x(i+1) = x(i) + hAx(i) , i = 0, . . . , N, T ∆t (x) = x(N ) . h

We then compute the trajectory error for increasing powers, p and q and their upper bound. For DMD, as Ψ = I , upper bound, L for ∥JΨ (x)∥2 ≤ L where L = 1. Figure 5 shows the trajectory error and upper bound with respect to the Euler integration error for extended eigenfunctions ϕ̄p and ϕ̄q . Assuming that flow F ∆t (x) = eA∆t x we can calculate the trajectory error with respect to eigenvector error by adding a random error vector δv with ∥δv∥ = 10−6 . 18

Figure 6 shows the trajectory error and the upper bound of the trajectory error due to this error in the eigenvector for extended eigenfunctions ϕ̄p and ϕ̄q .

0.0

0.0

−0.5

−0.5

−1.0

−1.0

−1.5

−1.5

log10

log10

Trajectory error and upper bound with respect to integration error G

−2.0 −2.5

−2.5

−3.0 −3.5

−2.0

−3.0

trajectory error: ET G (φ̄p , λ̄p ) upper bound: (λ̄M + LG )p − (λ̄M )p

1

2

3

4

5

6

7

8

9

1/p

−3.5

10

trajectory error: ET G (φ̄q , λ̄q ) upper bound: (λ̄M + LG )q − (λ̄M )q

1

2

3

4

p

5

6

7

8

9

1/q

10

q

Fig. 5: Integration error analysis for DMD eigenfunctions of continuous linear system (35). The trajectory error for computed eigenfunctions (blue) and upper bound (orange) with respect to the Euler integration error ϵG for powers p and q . Finally, we employ Algorithm 2, computing trajectory errors and upper bounds with respect to integration error to construct monomial eigenfunctions ϕ̄p and ϕ̄q with powers p and q such that the trajectory error stays below a desired upper bound ϵ = 0.2. Figure 7 shows the results of Algorithm 2 for a desired trajectory error upper bound ϵ = 0.1. As seen in the figure, the powers p and q suggested by the algorithm are close to the actual powers up to which the eigenfunctions can be extended.

5.2 Example - Nonlinear (transformation of linear) system Consider the linear continuous system

ẋ = Ax,

(46)



 −0.9 0.1 where A = . 0 −0.8 We transform using the diffeomorphism y = h(x) = log(ex + 1) (functions applied coordinate-wise) to obtain the new system

 1 − e−y1 0 A log(ey − 1). ẏ = 0 1 − e−y2 

(47)

Using Proposition 1, the explicit eigenfunctions of this nonlinear system are then given by (ϕ ◦ h−1 ) with eigenvalue λ, where ϕ is an eigenfunction of the linear system (46) with eigenvalue λ. 19

−1

−1

−2

−2

−3

−3

−4

−4

log10

log10

Trajectory error and upper bound with respect to eigenvector error ||δw|| = 10−6

−5

−5

−6

−6

−7

−7

trajectory error: ET G (φ̄p , λ̄p ) upper bound: CT G (p, λ̄)1/p ||δw||1/p

−8 1

2

3

4

5

6

7

8

9

trajectory error: ET G (φ̄q , λ̄q ) upper bound: CT G (q, λ̄)1/q ||δw||1/q

−8

10

1

2

3

4

5

6

p

7

8

9

10

q

Fig. 6: Eigenvector error analysis for DMD eigenfunctions of continuous linear system (35). The trajectory error for computed eigenfunctions (blue) and upper bound (orange) with respect to the eigenvector error ∥δw∥ = 10−6 for powers p and q .

0

−2

−2

−4

−4

log10

log10

Finding powers for extending eigenpairs given  = 0.2 with integration error G 0

−6

ET G (φ̄p , λ̄p )

−6

ET G (φ̄q , λ̄q )

−8

 G

1 p p 1/p − λ̄M ] L [( + (λ̄M ) )

1 q q 1/q − λ̄M ] L [( + (λ̄M ) )

 G

−8 1

2

3

4

5

6

7

8

9

10

1

2

3

4

5

p

6

7

8

9

10

q

Fig. 7: Results of Algorithm 2 applied to the DMD approximation of continuous linear system (35) with Euler integration error ϵG and desired trajectory error ϵ = 0.2. The value of p and q suggested by the algorithm– where the upper bound for ϵG (orange) crosses ϵG (red line) is close to actual value of p and q where the trajectory error (blue) crosses the required ϵ (black). As the eigenfunctions of the linear system are given by ϕi (x) = ⟨wi , x⟩ with eigenvalues λi , (i = 1, 2) where wi is left eigenvector of A with eigenvalue λi , the eigenfunctions of system (47) are given by

ϕnonlin (y ) = ϕi ◦ h−1 (y ) = ⟨wi , log(ey − 1)⟩ i with eigenvalues λi for i = 1, 2.

20

(48)

We sample the linear system using the exact solution with ∆t = 0.02, collecting 400 snapshot pairs with initial conditions uniformly randomly distributed between [−2, 2] × [−2, 2]. Then, we transform the sampled data using the diffeomorphism h. We use this transformed data to perform EDMD with a radial basis function (RBF) dictionary with 40 RBF Gaussian kernel functions with centers calculated using kmeans clustering of the transformed data. We take the grid G with a = 1, b = 2, n = 100, h = 0.01. Some of the explicit eigenfunctions on the grid G are shown in Figure 24 (Appendix E), and the computed EDMD eigenfunctions are shown in Figure 25 (Appendix E). Then we calculate ϵG by integrating over the grid and using the integration method RK45 with the explicit system (47) [26]. We calculate the first nine eigenpairs of the Koopman matrix. Then, we use the Koopman eigensolver algorithm defined in Algorithm 3 to extend the eigenfunctions of the system. The desired trajectory error is set to ϵ = 0.01. We use the algorithm to get up to p = 3 extended eigenfunctions for the first nine eigenfunctions with trajectory error less than ϵ. Figure 9 shows the spectrum and the powers up to which each eigenpair can be safely extended. Figure 26 (Appendix E) shows some of the extended eigenfunctions. The extended eigenfunctions do not match the explicit eigenfunctions in this case. This might be expected as the EDMD eigenfunctions can differ from the explicit eigenfunctions, as the transformed non-linear system can have an infinite number of independent eigenfunctions; we feel that given the size of the matrix and the accuracy of the calculation, such results are not unreasonable.

5.3 Bridging eigenfunctions across singularities To generate new eigenfunctions by repeatedly multiplying them with themselves or with other approximated eigenfunctions, we need to approximate a small initial set of them. However, as discussed in Sect. 3.2, in some cases, obtaining such initial sets of eigenfunctions is challenging, especially in cases where they involve infinities. In the following examples, we will explore how to extend eigenfunctions beyond infinity for systems with (possibly multiple) singularities.

5.3.1 Nonlinear system with two steady states Consider an ODE on R, with steady states at a = 2, b = 3, s.t.

ẋ = (x − a)(x − b).

(49)

The eigenvalues of the linearization around a and b are (a−b) and (b−a), respectively. Since scalar multiples cϕ(x) of Koopman eigenfunctions ϕ(x) are also eigenfunctions, we are allowed to consistently select a single representative member of this family by “pinning” the slope of an eigenfunction at some convenient reference value (e.g., link it to the slope of the eigenvector of the linearization at that steady state that has the same eigenvalue). The Koopman eigenfunctions of this system associated with the

21

EDMD model reconstruction (RBF dictionary) 2.2

2.0

2.0

1.8

1.8

1.6

1.6

x2

x2

Training data used during fit 2.2

1.4

1.4

1.2

1.2

1.0

1.0

0.8

0.8 0.6

0.6 0.75

1.00

1.25

1.50

1.75

2.00

0.75

1.00

1.25

x1 out-of-sample prediction (x1 )

1.75

2.00

out-of-sample prediction (x2 )

EDMD true system

2.50

EDMD true system

2.50

2.25

2.25

2.00

2.00

1.75

1.75

x2

x1

1.50

x1

1.50

1.50

1.25

1.25

1.00

1.00

0.75

0.75 0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.0

0.5

1.0

1.5

t

2.0

2.5

3.0

t

Fig. 8: EDMD model for the non-linear system (47). Top left shows the training data used for approximation, top right shows the reconstructed data, bottom shows the comparison between predicted trajectory and actual trajectory for one initial condition. steady states are obtained by solving ∇ϕki · ẋ = λki ϕki , i = 1, 2. This results in

ϕk1 (x) =

x−a x−b

λk1 /(b−a)

,

and

ϕk2 (x) =

x−b x−a

λk2 /(b−a)

.

(50)

When λk = k (b − a) with k ∈ Z, the eigenfunctions (50) take the form

 ϕk1 (x) =

x−a x−b

k1

 ,

and

ϕk2 (x) =

x−b x−a

k2 ,

k1 , k2 ∈ Z.

(51)

Then, almost by construction, real-valued Koopman eigenfunctions, associated with one of the steady states become zero at that steady state and approach infinity at the

22

Fig. 9: Spectrum computed using Algorithm 3 for the nonlinear system (47) and powers up to which the eigenfunctions can be extended for the first 9 eigenvalues. other. For instance, for k1 = k2 = 1, the Koopman eigenfunctions (51) are plotted in Fig. 10.

Fig. 10: The Koopman eigenfunctions (51) associated with the steady states x = 2 (left) and x = 3 (right) for k1 = k2 = 1. As shown, ϕk1 becomes zero at x = 2 and approaches infinity at x = 3, while ϕk2 becomes zero at x = 3 and approaches infinity at x = 2.

Eigenfunctions of this system can be systematically extended beyond the “next adjacent” steady state and even beyond the “next once removed” steady state. In the interval between a and b, both steady states lead to properly defined eigenfunction values; we show how this can be exploited to “cross” the singularity of the individual eigenfunctions. We first approximate the eigenfunctions associated with the linearization of each steady state in its neighborhood (so, either in a region [a − ϵ, a + ϵ], or [b − ϵ, b + ϵ], for an ϵ < b − a). If we compute them in a region that is large enough, i.e., we choose ϵ large, we get approximations in a region U ⊂ [a, b] for both sets. At this 23

point, it is appropriate to discuss some observations from computational experiments. Clearly, an eigencomputation of the EDMD matrix may converge to any power of the principal eigenvectors (discretized eigenfunctions). In principle, all eigenfunctions, including powers of eigenfunctions corresponding to each of the alternative nearby steady states, would be computable from data if we had an appropriate dictionary. What we practically observe is that data-driven computations with data collected only in the neighborhood of one steady state tend to numerically converge to the principal eigenfunctions “corresponding to” that steady state (not even to their integer powers). We never (in our experiments) saw the procedure converge to eigenfunctions associated with the other steady state. In fact, it is quite challenging to numerically approximate the functions ϕk1 and ϕk2 when they approach their farther away, more remote,“unassociated” steady state, given that they approach infinity quite rapidly there. We compute the logarithm of all functions in a region that excludes points close to the steady states, and test whether they are linearly dependent. The test consists of computing principal components of the logarithm of all eigenfunctions (evaluated on points in [1, 4]), stacked together in one large dataset L:



 | | | | | | L = log |ϕk1 | log |ϕ2k1 | · · · log |ϕ5k1 | log |ϕk2 | log |ϕ2k2 | · · · log |ϕ5k2 | . | | | | | |

10 3

0.1

10 8

0.0

PC 1

Singular value

Fig. 11 (left) shows that only a single direction is present in this dataset, i.e. all logarithms of the eigenfunctions are linearly dependent. The corresponding principal component of L is also shown (Fig. 11, right). Note that it is computed in log space, i.e. negative values here mean very small (but positive) values when considered in the original space. We cannot compute this with EDMD yet, as we show below. One of the goals of this paper is to discuss how such a global eigenfunction could be constructed numerically from partial EDMD results in the neighborhood of different steady states.

10 13

0.1 1 2 3 4 5 6 7 8 9 10 Index of singular value

1.0

1.5

2.0

2.5 3.0 Space

3.5

4.0

Fig. 11: Left: Energy (singular values) of the principal components of the logarithm dataset. Only one component is relevant, as expected—all the others have numerically zero energy. Right: the principal component associated to the largest singular value.

Approximation with EDMD. In practice, explicit forms of eigenfunctions are usually not available, and must be approximated from data. We thus now show how to approximate the eigenfunctions 24

Eigenfunction

numerically, with EDMD, using a radial basis function dictionary (Gaussian kernels), centered at points in the region [1, 4]. We already discussed that data in the neighborhood of each steady state tend to lead to computationally identified eigenfunctions “corresponding” to the linearization around that steady state. Computations combining data across the entire interval tend to provide very poor eigenfunction approximations due to the limitations of the dictionary. We therefore separately approximate eigenfunctions around x = a = 2 and x = b = 3, and then try to map them to each other in the intermediate region [a, b]. In practice, we cannot accurately approximate eigenfunctions close to their singularities. We thus only work in the regime where their values are comparatively small. In the regime where one eigenfunction approaches infinity, we can instead work with the eigenfunction that is its appropriate inverse power - which then approaches zero. Fig. 12 shows the approximations. Several spurious eigenfunctions are visible (shown are three in the last three columns), which are incorrectly identified.

0.75 0.50 0.25 0.00 0.25 0.50 0.75

Eigenfunction

0.50

KEF #1

1.8

2.0

2.2 2.4 Space

KEF #1

2.6

KEF #2

2.8 1.8 KEF #2

2.0

2.2 2.4 Space KEF #3

2.6

KEF #3

2.8 1.8

2.0

2.2 2.4 Space

KEF #4

2.6

KEF #4

2.8 1.8

KEF #5

2.0

2.2 2.4 Space KEF #6

2.6

KEF #5

2.8 1.8 KEF #7

2.0

2.2 2.4 Space

2.6

2.8

KEF #8

0.25 0.00 0.25 0.50 0.75 2.5

Space

3.0

2.5

Space

3.0

2.5

Space

3.0

2.5

Space

3.0

2.5

Space

3.0

2.5

Space

3.0

2.5

Space

3.0

2.5

Space

3.0

Fig. 12: Eigenfunctions approximated with EDMD (Gaussian radial basis function dictionary, kernel bandwidth ϵ = 0.05 and ϵ = 0.15 for top and bottom rows). Eigenfunctions associated to eigenvalues with nonzero imaginary part are excluded, and the eigenfunctions are normalized to have norm 1. The steady state is denoted by an orange vertical dashed line, indicating that the top row is computed around steady state x = a = 2, the bottom row around steady state x = b = 3. The last three eigenfunctions in each row are spurious and incorrectly oscillate around zero.

In the last step, we try to find a linear map from one set of the eigenfunctions (associated with the left steady state, for example) to the other set (here, to the right one) in logarithmic space. This can only practically be identified (again, due to dictionary limitations) in the intermediate region of the interval between the steady states; we choose the domain [2.25, 2.75]. Fig. 13 (top) shows the logarithm of the absolute value of several non-spurious eigenfunctions in this region. We then solve the following linear systems with a least squares method (incl. Tikhonov regularization)

25

for the coefficients c = (c1 , c2 ): log |ϕk1 |c1 = log |ϕk2 |,

log |ϕk2 |c2 = log |ϕk1 |.

The coefficients c1 and c2 allow us to map from one set of (logarithms of) eigenfunctions to the other. The result is shown in Fig. 13. 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5

2 3

2.3

2.4

Original Approximation

4 5 6

2.3

2.4

2.5 Space

2.6

2.7

2.5

Log Abs Eigenfunctions (right s.s.)

Log Abs Eigenfunctions (left s.s.)

6.0

2.6

2.7

2

Original Approximation

3 4 5 6

2.3

2.4

2.5 Space

2.6

2.7

Fig. 13: Top: Logarithms of the absolute value of selected, non-spurious eigenfunctions in the intermediate region [2.25, 2.75] between the steady states. The eigenfunctions associated to a = 2 are orange, the ones associated to b = 3 are blue. Bottom: A linear map from the one set of eigenfunctions (associated to left/right steady state) to the other (right/left) provides a good approximation, meaning we can accurately construct the value of eigenfunctions “on the other side of the steady state” (not shown).

5.4 Singularities in two dimensions and the use of isochrons We saw that in one dimension, Koopman eigenfunction singularities arise naturally as we approach adjacent steady states (which also naturally constitute the boundary of basins of attraction). We now proceed to explore singularities of Koopman eigenfunctions for continuous time dynamical systems in two dimensions. Here, the singularities 26

will arise along entire one-dimensional curves: (a) separatrices, consisting of the stable manifolds of saddle steady states, or (b) limit cycles, surrounding a steady state in the neighborhood of which the Koopman eigenfunction computation is initiated. Such a limit cycle is also the boundary of the basin of attraction of the steady state (e.g., forward in time if the steady state inside is stable, and the limit cycle unstable; or backward in time in the opposite case). In both cases, the important notion is that of an isochron, first defined in terms of limit cycles [1], but straightforward to also extend in the case of saddle-associated separatrices.

Definition 6. (Dominant Koopman eigenvalue and eigenfunction) Consider a dynamical system ẋ = F (x), where x ∈ M, and F ∈ C 2 in the neighborhood of an asymptotically stable equilibrium x∗ . That is, the eigenvalues λj of the Jacobian matrix J (x∗ ) satisfy Re(λj ) < 0 for all j . In this case, the eigenvalues λj are directly related to the eigenvalues exp(tλj ) of the associated Koopman operator KFt t . The dominant eigenvalue λ1 is defined as the eigenvalue satisfying Re(λ1 ) > Re(λj )

for all λ1 ̸= λj .

(52)

We assume the existence of such an eigenvalue, which is guaranteed in monotone systems [27, Chapter 5.5]. The eigenfunction associated with λ1 is referred to as the dominant eigenfunction and is denoted by φλ1 . Under these assumptions, the dominant eigenfunction φλ1 can be computed using the Laplace average: 1 T →∞ T

fλ∗1 (x) = lim

Z T 0

f ◦ F t (x) e−λ1 t dt,

(53)

where f ∈ C 1 satisfies f (x∗ ) = 0 and (∇f (x∗ ))T v1 ̸= 0, with v1 being the right eigenvector of the system Jacobian J (x∗ ) corresponding to λ1 . The function fλ∗1 represents the dominant eigenfunction φλ1 (x) up to a scalar multiple [2]. For a system with a stable limit cycle Γ, the dominant Koopman eigenfunction associated with the linearization around the limit cycle is linked with the leading nontrivial Floquet exponent λ1 , which governs the slowest transverse decay of trajectories towards the cycle. In an action-angle coordinate system, the function fλ∗1 provides a global action coordinate that captures the attractivity of the limit cycle. For a dynamical system in d dimensions with an exponentially stable limit cycle Γ of period ω ∈ R+ , there exist Koopman eigenfunctions [2, Chapter 15] associated with eigenvalues λi , i = 1, . . . , d, such that λ1 = ω , and λi are the (d − 1) Floquet exponents.

Definition 7. (Isostables, [2, Chapter 15]) The isostables of the system are level sets of the Koopman eigenfunctions associated to the Floquet exponents. Formally, for any value r ∈ R and any point x in the state space M, an isostable associated to one of

27

the eigenvalues λi is the set

Ir = {x ∈ M : |φλi (x)| = erλi },

(54)

where φλi (x) is the principal Koopman eigenfunction associated with the eigenvalue λi . For an equilibrium, isostables describe trajectories that decay at the same exponential rate towards x∗ . For a limit cycle, isostables represent sets of points that converge to the cycle with the same phase-independent transient behavior. Isostables provide a global partition of the state space and are closely linked to the stable and unstable manifolds of the system invariant sets, offering valuable insights into the transient dynamics of the system. They can be computed throughout the entire basin of attraction of an attractor using methods such as Laplace averages [1].

Definition 8. (Isochrons) The isochrons of the limit cycle are the level sets of the Koopman eigenfunction φiω , formally defined as: Iθ = {x ∈ M | ∠φiω (x) = θ},

θ ∈ [0, 2π ).

(55)

Each isochron represents a set of initial conditions that asymptotically converge to the same phase on the limit cycle. The evolution of these sets follows the periodic relation (F t )k (Iθ ) = Iθ+ωt (mod 2π) . In particular, after a full period, the isochrons satisfy (F t )k (Iθ ) = Iθ , where k = 2π ω [2]. Isochrons provide a global phase coordinate system for the dynamics, ensuring that all points on a given isochron approach the same asymptotic phase on the attractor.

5.4.1 Isochrons and isostables of limit cycles Consider the Van der Pol system with the parameter µ = 0.3. For this parameter value, the system exhibits a stable limit cycle, shown in red in Figure 14, surrounding an unstable (source, focus) steady state. This figure also illustrates the isochrons (left) and isostables (right) for this Van der Pol system. Isochrons and isostables are (separately) Koopman eigenfunctions for the system; here, they are both anchored on the linearization around the stable limit cycle and their eigenvalues embody the attractivity of the limit cycle as well as its period. Clearly, the isostable component goes to zero on the limit cycle and to infinity at the steady state in the middle (see Figure 17). A negative power of the isostable would then go to infinity on the limit cycle. If we also computed the isochrons and isostables associated with the unstable spiral steady state inside the limit cycle, we would encounter a situation similar to the one in the interval between different steady states (different invariant objects) in one dimension. There, we used the ratios of respective eigenvalues to transform between eigenfunctions. Doing this in the context of two dimensional systems with limit cycles is the context of current research (see [28] for possible structures for “inner” and “outer” isochrons).

28

Inside the limit cycle, trajectories converge to it and cannot escape or cross it. Therefore, it is important to examine the region interior, and the one exterior of the limit cycle, along with how these two areas are connected. One effective method for approximating them in the entire domain is to use Laplace averaging [1] of the observable g (x1 , x2 ) = sin(x1 + x2 ).

3

2.7

3

2

1.8

2

1

0.9

1

0

0.0

0

1

0.9

1

2

1.8

2

2.7

3

3 3

2

1

0

1

2

3

3

2

1

0

1

2

3

Fig. 14: Isochrons (left) and isostables (right) associated with the limit cycle of the Van der Pol system (µ = 0.3), computed with Laplace averaging of the observable g (x1 , x2 ) = sin(x1 + x2 ).

5.4.2 Transformation of trajectories across the limit cycle We now describe a procedure to transform trajectories from within the limit cycle invertibly into trajectories outside of the limit cycle. The key idea is to represent the trajectories in terms of eigenfunctions (isochrons and isostables), and to construct the transformation by transforming between eigenfunctions “inside” and eigenfunctions “outside” of the limit cycle. The procedure works as follows. 1. Compute two complex-valued eigenfunctions ϕ1 , ϕ2 : R2 → C (cf. Figure 15, leftmost and rightmost plot), where ∠ϕ1 and ∠ϕ2 are isochrons associated with the limit cycle and eigenfunctions associated with the steady state, respectively, and |ϕ1 |, |ϕ2 | are the corresponding isostables. The computation of ϕ1 is performed in the neighborhood of the limit cycle, the computation of ϕ2 in the neighborhood of the steady state. 2. Extend the functions ϕ1 , ϕ2 as much as possible (accurately) within the limit cycle, i.e., ϕ1 so that its domain approaches the steady state, and ϕ2 so that its domain approaches the limit cycle. This is possible by solving the Koopman PDE with the method of characteristics. Note that |ϕ1 | approaches infinity towards the steady state (cf. Figure 14, right), while |ϕ2 | is zero at the steady state and approaches infinity towards the limit cycle.

29

0.54 0.50 0.46 0.42 0.38 0.34 0.30 0.26 0.22 0.18

3. Define an invertible transformation function Ti : C → C (subscript i for “inside”) between the values of ϕ1 and the values of ϕ2 within the limit cycle. This is possible, because both ϕ1 and ϕ2 separately can be used to represent the open 2D domain within the limit cycle (excluding it and the steady state). 4. As isostables (and isochrons) of the limit cycle are defined on either side of it, and both sets have levels from zero to infinity, it is possible to identify isostables inside with isostables outside (cf. Figure 15, rightmost plot, black circles). Using this identification, it is possible to create a separate transformation To (subscript o for “outside”) that (a) maps isostables “inside” the limit cycle to isostables “outside” the limit cycle invertibly and (b) does so across the isochrons with identical values inside and outside. 5. For any trajectory within the limit cycle, it is possible to (a) transform it from Cartesian coordinate representation into an isochron/isostable representation (i.e., range of ϕ2 ), then (b) map this representation to the corresponding one in the range of ϕ1 using Ti , and finally (c) map this representation to isochrons/isostables outside of the limit cycle using To . limit cycle isochrons isostables

limit cycle x coordinate y coordinate invertible transform (around limit cycle)

invertible transform (inside)

Steady state eigenfunction

Limit cycle eigenfunction

Cartesian coordinates

Fig. 15: Isochrons and isostables (combined as a single complex-valued Koopman eigenfunction) can be used to construct invertible transforms inside the limit cycle, by transforming to the Cartesian coordinates (x, y ). The transformation for the steady state breaks down close to the limit cycle, while the transformation for the limit cycle breaks down close to the steady state.

We now provide explicit formulations for the eigenfunctions and the transformations in a specific example.

Example 1. Consider a planar system (in polar coordinates) given by ( ṙ = r(µ − r2 ) θ̇ = ω + α (r −

√

µ),

30

µ, ω > 0, α ∈ R

.

(56)

At r = 0, the origin is an unstable equilibrium, and the circle r = cycle (Fig. 16).

√

µ is a stable limit

Fig. 16: Top: From left to right, the phase portrait of system (56), for µ = ω = α = 1, with an unstable steady state (green point) and a stable limit cycle (red circle), along with trajectories (light blue). Isostables and isochrons of the limit cycle (top, for C = 1) and the steady state (bottom row) are shown in four separate plots. Bottom: The isostables and isochrons of the steady state (for C = 1). The way this example is constructed means the linearization of the limit cycle and the steady state are identical, for general dynamical systems with limit cycles this is not the case.

I. Koopman eigenfunctions : The separable eigenfunctions of the form

ϕλ,k (r, θ) = hλ,k (r) eikθ ,

k ∈ Z,

(57)

with eigenvalue λ ∈ C, can be obtained by solving the Koopman PDE using the method of characteristics:

ϕλ,k (r, θ) = C



r2  |µ − r2 |

λ−ikω 2µ

 √µ + r −ikα/√µ r

eikθ .

(58)

We consider C > 0 for the sake of simplicity. Using linearization around the limit cycle √ r = µ yields a radial Floquet exponent −2µ and angular frequency ω . Thus, for 31

Fig. 17: Isostables (Koopman eigenfunctions) for a limit cycle system, plotted in three and two dimensions. Two representative trajectories are shown, with the z-coordinate equal to the respective isostable value (blue inside, orange outside of limit cycle). Left: An isostable associated to the steady state. It approaches (negative) infinity on the limit cycle, and is zero on the steady state (also see Fig. 16, bottom left). Outside of the limit cycle, it is extended using the limit cycle isostable. Right: An isostable associated to the limit cycle. It is zero on the cycle, and approaches infinity towards the steady state (also see Fig. 16, top, center). k = 1 and λ1 = −2µ + iω , the eigenfunction associated with the limit cycle becomes ϕλ1 (r, θ) = C

√ √ |µ − r2 |  µ + r −iα/ µ iθ · e , r2 r

C > 0.

(59)

Therefore, the isostables of the limit cycle are

Ir,LC = |ϕλ1 (r, θ)| = C

32

|µ − r2 | , r2

(60)

and the corresponding isochrons are given by

 √µ + r  α (mod 2π ). Iθ,LC = ∠ϕλ1 (r, θ) = arg(ϕλ1 (r, θ)) = θ − √ ln µ r

(61)

Both isostables and isochrons of the limit cycle are defined on either side of it. For √ isostables, we have Ir,LC = 0 on the limit cycle r = µ. We consider two branches √ corresponding to the regions inside and outside the limit cycle (0 < r < µ and √ r > µ):

• Interior branch: ϕλ1 (0,√µ)×[0,2π) where |µ − r2 | = µ − r2 > 0. On this punctured open disk

√ bijective ϕλ1 (0,√µ)×[0,2π) : (0, µ) × [0, 2π ) −−−−−−→ (0, ∞) × [0, 2π ),
and Ir,LC (r) = C rµ2 − 1 ∈ (0, ∞). • Exterior branch: ϕλ1 (√µ,∞)×[0,2π) where |µ − r2 | = r2 − µ > 0. On this exterior region √ bijective ϕλ1 (√µ,∞)×[0,2π) : ( µ, ∞) × [0, 2π ) −−−−−−→ (0, C ) × [0, 2π ),
and Ir,LC (r) = C 1 − rµ2 ∈ (0, C ). Moreover

• If α = 0, then the isochrons are ∠ϕλ1 (r, θ) = θ , which correspond to straight radial lines. • If α ̸= 0, then the isochrons are curved lines. Similarly, linearization around the unstable steady state r = 0 gives a radial expo√ √ nent µ and angular frequency ω − α µ. Hence, for k = 1 and λ1 = µ + i(ω − α µ), the eigenfunction associated with the steady state is

 √µ + r −iα/√µ · p eiθ , µ − r2 µ − r2

ϕλ2 (r, θ) = C p

r

0≤r<

√

µ,

C > 0.

(62)

µ,

(63)

The isostables of the steady state can be expressed as

r

Ir,SS = |ϕλ2 (r, θ)| = C p

µ − r2

∈ [0, ∞),

0≤r<

√

and the associated isochrons are

 √µ + r  α (mod 2π ), Iθ,SS = ∠ϕλ2 (r, θ) = arg(ϕλ2 (r, θ)) = θ − √ ln p µ µ − r2

33

0≤r< (64)

√

µ,

√ which are straight radial lines for α = 0. Furthermore, ϕλ2 is defined for 0 ≤ r < µ √ (inside the limit cycle) and satisfies ϕλ2 (0, θ) = 0. Its restriction to (0, µ) × [0, 2π ) is a bijection onto C \ {0}, while the full map is surjective onto C but not injective at 0 (all θ map to 0). From eq. (60) and eq. (63), we can see that the isostables |ϕλ1 (r, θ)| of the limit cycle are zero at the limit cycle and approach infinity toward the steady state inside (and toward infinity outside). In contrast, the isostables |ϕλ2 (r, θ)| of the steady state are zero at the steady state and approach infinity towards the limit cycle. For more details, see Fig. 16.

II. Invertible transformation Ti :

√ For fixed parameters µ, C > 0, and α ∈ R, the maps ϕλ1 (0,√µ)×[0,2π) : (0, µ) × √ [0, 2π ) −→ C \ {0} and ϕλ2 (0,√µ)×[0,2π) : (0, µ) × [0, 2π ) −→ C \ {0} are bijec√ tive on (0, µ) × [0, 2π ) inside the limit cycle on the punctured open disk where the steady state and the limit cycle are excluded. The transformation Ti can be built by composition as  −1 : C \ {0} −→ C \ {0} Ti := ϕλ2 (0,√µ)×[0,2π) ◦ ϕλ1 (0,√µ)×[0,2π) s     C3 |z| α Ti (z ) = ln exp i arg z + 2√ , µ C |z|

(65)

which is invertible with the inverse

 −1 : C \ {0} −→ C \ {0} Ti−1 = ϕλ1 (0,√µ)×[0,2π) ◦ ϕλ2 (0,√µ)×[0,2π)     C3 √α ln C Ti−1 (v ) = exp i arg v − , µ |v| |v|2

(66)

which implies that we can transform the values of ϕ1 and the values of ϕ2 within the limit cycle.

III. Invertible transformation To : We construct a transformation To that maps the values of interior isostables of the limit cycle to their exterior counterparts while preserving the associated isochron values. For the isostables of the limit cycle, we have

 µ  √  C 2 − 1 ∈ (0, ∞), for 0 < r < µ (inside) r Ir,LC = .    C 1 − µ ∈ (0, C ), for r > √µ (outside) r2

34

(67)

If we define a monotone bijection

T̃o : (0, ∞)×[0, 2π ) −→ (0, C ) × [0, 2π )   Cr (r, θ) 7−→ , θ , 1+r

(68)

then the map

To =

ϕλ1 |(√µ,∞)×[0,2π)

−1

√ √ ◦ T̃o ◦ ϕλ1 (0,√µ)×[0,2π) : (0, µ) × [0, 2π ) −→ ( µ, ∞) × [0, 2π ),

defines a bijection between the interior and exterior regions of the limit cycle while preserving the Koopman phase along the isochrons by construction.

IV. Mapping a trajectory inside the limit cycle to exterior isostables/isochrons: Transforming trajectories within the limit cycle to trajectories outside of it is possible if both isostables and isochrons are available. If only isochrons of limit cycle and steady state are available, then CIFT-bifurcations [28] can prevent this type of mapping. In the former case, any trajectory x(t) = (r(t), θ(t)) inside the limit cycle, with √ 0 < r(t) < µ and t ∈ I , can be mapped into an isochron/isostable representation of the steady state using the eigenfunction ϕλ2 . Next, this representation can be transferred to the isochron/isostable representation associated with the limit cycle by applying Ti−1 . Finally, we can map this representation to isostables/isochrons outside
−1 of the limit cycle using ϕλ1 |(√µ,∞)×[0,2π) ◦ T̃o . Thus the overall transformation is (r, θ) 7−→

ϕλ1 |(√µ,∞)×[0,2π)

−1

◦ T̃o ◦ Ti−1 ◦ ϕλ2 (0,√µ)×[0,2π) (r, θ).

(69)

5.5 Isochrons for systems with saddle points Figures 20–23 show Koopman eigenfunctions of the bi-stable system (Equation (70)) whose two basins of attraction are separated by the stable manifold of the saddle. Starting with each of the stable steady states, and constructing the Koopman eigenfunctions based on their linearization, we can clearly see these eigenfunctions going to infinity at the separatrix. Two important questions arise: is there a way of extending each one of them across the separatrix to the “other” basin? Is there Is there (and if yes, can we construct) an analogy with the power transformation in one dimension? To explore this, we first performed a simple numerical experiment in a linear system that only possesses a single saddle (not multiple basins of attractions, and accordingly no separatrices). The results are shown in Figure 18. We visualize that the level sets of the two principal eigenfunctions of the system intersect the stable/unstable manifolds of the saddle exactly in the same way that isochrons and isostables would for a limit cycle. Here, the limit cycle corresponds to the unstable manifold of the saddle, the level sets of the eigenfunction associated to its dynamics (Figure 18, left) corresponds to the isochrons, and the level sets of the eigenfunction associated to the stable dynamics (Figure 18, left) corresponds to the isostables converging to the limit cycle. Points 35

on these level sets will, backward in time, approach the unstable eigendirection of the saddle and will asymptotically approach backward in time the point at which the level sets cross the unstable eigendirection. In that sense, these level sets are isochrons of the unstable eigendirection. They will help us bridge eigenfunctions defined on alternate sides of it. For more complicated systems with limit cycles (cf. section 5.4.2) or separatrices (cf. section 5.5.2), the same concept can, in principle, be applied.

5.5.1 Isochrons for a nonlinear systems with a single saddle point Consider a nonlinear 2D system with a saddle at the origin (red dot, with trajectories as gray curves in Figure 18). For this saddle, we can construct one eigenfunction associated with the unstable direction/manifold and another associated with the stable direction/manifold, both of which precisely correspond to the saddle’s eigenvalues by construction. To construct this example, we start with a linear system associated with the eigenvalues (λ1 , λ2 ) = (−1, 1.5) of a saddle point and then transform the state space to obtain a nonlinear system. For the original linear system, the eigenfunctions are analytically available. Here, this linear system is defined by

ẋ1 = λ1 x1 , ẋ2 = λ2 x2 . Hence, the eigenfunctions are simply ϕ1 (x) = x1 , ϕ2 (x) = x2 , and associated to the eigenvalues λ1 , λ2 of the saddle. To calculate these eigenfunctions for the saddle in the nonlinear system, a nonlinear transformation (diffeomorphism) is used to transform the state space. Here, we choose x 7→ exp( π1 Rx) − (1, 1)T , with R a matrix that rotates x by 60 degrees and the exponential function applied coordinate-wise. This also transforms the eigenfunctions (through the same function). Figure 18 shows these transformed Koopman eigenfunctions associated with the stable and unstable manifolds of the saddle point of the nonlinear system. If we take any points on one contour line associated with the unstable manifold (left plot), they will diverge in the same manner towards infinity in the unstable direction. The right plot illustrates the parameterization of the stable direction. All points on one of the contour lines will converge backwards in time to the same point on the stable manifold. This is the analog to isochrons being associated with limiting points on the limit cycle, and to isostables being associated with convergence to the cycle.

36

2.5

stable unstable saddle

2.0 1.0

y2

y2

1.5

stable unstable saddle

0.5 0.0 0.5 0

1 y1

2

0

1 y1

2

Fig. 18: Koopman eigenfunctions (contour lines) of a system with a saddle (red dot) at the origin, and example trajectories as thin, gray curves in the background. The eigenfunctions associate level sets of the state space with points on the unstable (left) and stable (right) manifold of the saddle.

5.5.2 Separatrices for multistable systems Consider now a 2D nonlinear system (c.f. Figure 19, left) given by

( ẏ1 = −(y1 − 1/4)(y1 + 1/4) y1 , ẏ2 = −y2

.

(70)

This system has a single saddle point at P0 = (0, 0) with eigenvalues λ0,1 = 1/16 and λ0,2 = −1, and two stable nodes at P1 = (−1/4, 0) and P2 = (1/4, 0), both having eigenvalues λ1,1 = λ2,1 = −1/8 and λ1,2 = λ2,2 = −1. To make the discussion more visually compelling, we transform this system using an analytic, non-linear diffeomorphism



x1 x2



 y1 + y14 + 2y12 y2 + y22 . = x = h(y ) = 2 y12 + y2 

(71)

The inverse of h is also analytic, and given by

    y1 x1 /2 − (x2 /2)2 −1 . = y = h (x) = y2 −(x1 /2)2 + x2 /2 + 1/4(x1 x22 ) − (x2 /2)4

(72)

The transformed system has a saddle at P̃0 = (0, 0), and two stable nodes at P̃1 = (0.5078125, 1/8) and P̃2 = (−0.4921875, 1/8). Nearly all initial conditions except those on the saddle point’s stable manifold - will converge towards one of the stable equilibria. The stable manifold of the saddle point is represented by a green curve, and it constitutes the boundary separating the basins of attraction of the two stable nodes. Meanwhile, the two sides of the unstable manifold of the saddle are shown in red and asymptote each to one of the two stable nodes (Fig. 19). The system in (x1 , x2 ) coordinates is constructed through a diffeomorphism h from a simpler system 37

0.5

0.5

0.0

0.0

0.5 1.0

1.0

y1 y2

0.5 x2

1.0

x2

y2

1.0

0.5 1.0

0.5

0.0 y1

0.5

1.0

1.0

0.0 0.5

1.0

0.5

0.0 x1

0.5

1.0

1.0

1.0

0.5

0.0 x1

0.5

1.0

Fig. 19: Phase portrait of the 2D nonlinear system (70) with two stable nodes and one saddle. The stable and unstable manifolds of the saddle point are depicted as green and red curves respectively. The left plot shows the system in original coordinates y , the center plot shows the nonlinear transformation h−1 from x to y , and the right plot shows the transformed system in x coordinates.

Fig. 20: Koopman eigenfunction ϕ41 of the system (70) in the coordinates y = (y1 , y2 ) for λ1 = 1/16 (top left); and λ3 = −1/8 (top right). The transformed eigenfunction ϕ̃41 is represented in the coordinates x = (x1 , x2 ) for λ1 = 1/16 (bottom left); and λ3 = −1/8 (bottom right). in (y1 , y2 ) coordinates, whose Koopman eigenfunctions can be constructed explicitly (and thus also transformed explicitly). Although the equilibrium points P0 , P1 and P2 give rise to six local eigenvalues, there are only three distinct eigenvalues λ1 = 1/16, λ2 = −1, λ3 = −1/8, since some of them are equal. For system (70) in coordinates y = (y1 , y2 ), Koopman eigenfunctions ϕi , associated with the distinct eigenvalues λi , i = 1, 2, 3, can be constructed separately for each coordinate, so that

ϕk1 (y ) = (y1 − 1/4)(y1 + 1/4) y1−2

38

k λi /λ3

, i = 1, 3,

k ∈ Z,

(73)

Fig. 21: The transformed eigenfunction ϕ̃41 in 3D space with coordinates (x1 , x2 , ϕ̃41 ) for λ1 = 1/16 (left); and λ3 = −1/8 (right). Note that the functions extend to positive infinity, but the plot shows a flat plateau at ϕ = 50 (where our plotting saturates) instead.

1 ,m2 Fig. 22: Combinations of the transformed eigenfunctions ϕ̃m in the coordinates 12 x = (x1 , x2 ) for m1 = −8/9, m2 = 8/9 and λ3 = −1/8 (left); and m1 = m2 = −2 and λ1 = 1/16 (right).

1 ,m2 Fig. 23: Combinations of the transformed eigenfunctions ϕ̃m in the coordinates in 12 m1 ,m2 3D space with coordinates (x1 , x2 , ϕ̃12 ) for m1 = −8/9, m2 = 8/9 and λ3 = −1/8 (left); and m1 = m2 = −2 and λ1 = 1/16 (right).

39

and

ϕk2 (y ) = y2

−k λ2

k

= y2 ,

k ∈ Z.

(74)

m2 1 ,m2 1 Of course, by Proposition 1, all combinations ϕm = ϕm with m1 , m2 ∈ R 12 1 · ϕ2 are also eigenfunctions.

The eigenfunctions of the transformed system in x = (x1 , x2 ) are given by

ϕ̃k1 (x) = ϕ1 ◦ h−1 (x) =

(75)

x1 /2 − (x2 /2)2 − 1/4 x1 /2 − (x2 /2)2 + 1/4 x1 /2 − (x2 /2)2

i = 1, 3,

−2 k λi /λ3

k ∈ Z,

and k

ϕ̃k2 (x) = ϕ2 ◦ h−1 (x) = − (x1 /2)2 + x2 /2 + 1/4(x1 x22 ) − (x2 /2)4 ,

k ∈ Z. (76)

m2 1 ,m2 1 Similarly, all combinations ϕ̃m = ϕ̃m for m1 , m2 ∈ R are again eigenfunc12 1 · ϕ̃2 tions for the transformed system. See the eigenfunctions plotted in Figures 20–23 and Figure 32 in Appendix G.

6 Discussion We introduce an approach to numerically construct a large family of eigenfunctions of the Koopman operator, given a small number of eigenvectors from a Koopman matrix, approximated using extended dynamic mode decomposition. Because integer powers amplify numerical errors exponentially fast, we also introduce error metrics to track how many new eigenfunctions we can credibly (i.e. accurately enough) construct. We demonstrate the approach in several computational experiments. The construction of a larger set of eigenfunctions has several benefits. Since integer powers are typically not linearly related to each other, the newly constructed eigenfunctions help to span a larger subspace of functions over the data. This allows us to approximate general observables of dynamical systems more accurately. Another important observation concerns the construction of eigenfunctions across singularities. We demonstrate how such eigenfunctions arise and how they approach infinity as the input approaches the singularity. We then show how numerical approximations of the functions on both sides of singularities allow us to construct numerical approximations of the entire function, across the singularity. If we used neural networks with singularities to approximate such eigenfunctions, we would need tools such as rational activation functions (see [29] for a recent pre-print; [30] discuss rational functions, but without singularities). The representation of trajectories in terms of eigenfunctions poses challenging new questions for systems with bifurcations related to those eigenfunctions (e.g., so-called “cubic isochron foliation tangency” bifurcations, CIFT, see [28]). Such qualitative 40

,

changes of the system are only visible through tangential crossing of isochrons, but not apparent in the system’s dynamic behavior. The representation of the state in terms of isochron eigenfunctions becomes singular at the bifurcation point, and thus may provide a path to analyze this behavior.

Declarations • Funding: The authors are grateful to Prof. Alex Townsend for several comments on the manuscript. F.D. was partially funded by the German Research Foundation/DFG, project 468830823. Z.M. is grateful to the Bundesministerium für Forschung, Technologie und Raumfahrt (BMFTR, Federal Ministry of Research, Technology and Space) for funding through project OIDLITDSM, No. 01IS24061. The work of IGK was partially supported by the US Department of Energy and the US National Science Foundation. • Conflict of interest/Competing interests: The authors declare that they have no Conflict of interest. • Ethics approval and consent to participate: Not applicable. • Consent for publication: Not applicable. • Data availability: No datasets were generated during the current study. • Materials availability: Not applicable. • Code availability: All codes developed in this study will be made publicly available upon publication. • Author contribution: F.D.; Y.G.K. devised the main ideas, Z.M.; F.D.; Y.G.K. devised the experiments and wrote the paper; S.M. and S.H. conducted experiments and computed the numerical error estimates.

References [1] Mauroy, A., Mezić, I., Moehlis, J.: Isostables, isochrons, and Koopman spectrum for the action–angle representation of stable fixed point dynamics. Physica D: Nonlinear Phenomena 261, 19–30 (2013) https://doi.org/10.1016/j.physd.2013. 06.004 [2] Mauroy, A., Susuki, Y., Mezic, I.: Koopman Operator in Systems and Control vol. 484. Springer, ??? (2020)

41

[3] Kaiser, E., Kutz, J.N., Brunton, S.L.: Data-driven discovery of Koopman eigenfunctions for control. Machine Learning: Science and Technology 2(3), 035023 (2021) [4] Hartman, P.: A lemma in the theory of structural stability of differential equations. Proceedings of the American Mathematical Society 11(4), 610–620 (1960) https://doi.org/10.1090/S0002-9939-1960-0121542-7 [5] Grobman, D.M.: Homeomorphism of Systems of Differential Equations. Doklady Akademii Nauk SSSR (128), 880–881 (1959) [6] Mezić, I.: Spectrum of the Koopman Operator, Spectral Expansions in Functional Spaces, and State-Space Geometry. Journal of Nonlinear Science 30(5), 2091– 2145 (2019) https://doi.org/10.1007/s00332-019-09598-5 [7] Kvalheim, M.D., Sontag, E.D.: Global linearization of asymptotically stable systems without hyperbolicity. Systems & Control Letters 203, 106163 (2025) https://doi.org/10.1016/j.sysconle.2025.106163 [8] Budišić, M., Mohr, R.M., Mezić, I.: Applied Koopmanism. Chaos: An Interdisciplinary Journal of Nonlinear Science 22, 047510 (2012) https://doi.org/10.1063/ 1.4772195 [9] Bollt, E.M.: Geometric considerations of a good dictionary for Koopman analysis of dynamical systems: Cardinality, “primary eigenfunction,” and efficient representation. Communications in Nonlinear Science and Numerical Simulation 100, 105833 (2021) https://doi.org/10.1016/j.cnsns.2021.105833 [10] Williams, M.O., Kevrekidis, I.G., Rowley, C.W.: A data-driven approximation of the koopman operator: Extending dynamic mode decomposition. Journal of Nonlinear Science 25(6), 1307–1346 (2015) https://doi.org/10.1007/ s00332-015-9258-5 [11] Li, Q., Dietrich, F., Bollt, E.M., Kevrekidis, I.G.: Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the Koopman operator. Chaos: An Interdisciplinary Journal of Nonlinear Science 27(10), 103111 (2017) https://doi.org/10.1063/1.4993854 [12] Schmid, P.J.: Dynamic Mode Decomposition and Its Variants. Annual Review of Fluid Mechanics 54(1), 225–254 (2022) https://doi.org/10.1146/ annurev-fluid-030121-015835 [13] Colbrook, M.J., Townsend, A.: Rigorous data-driven computation of spectral properties of Koopman operators for dynamical systems. Communications on Pure and Applied Mathematics 77(1), 221–283 (2024) https://doi.org/10.1002/ cpa.22125

42

[14] Mezic, I.: Koopman Operator Spectrum and Data Analysis. arXiv (1702.07597) (2017) arxiv:1702.07597v1 [15] Schmid, P.J.: Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28 (2010) https://doi.org/10.1017/ S0022112010001217 [16] Colbrook, M.J.: The mp-EDMD algorithm for data-driven computations of measure-preserving dynamical systems. SIAM Journal on Numerical Analysis 61(3), 1585–1608 (2023) [17] Mohr, R., Mezić, I.: Construction of eigenfunctions for scalar-type operators via Laplace averages with connections to the Koopman operator. arXiv (2014) arXiv:1403.6559v2 [18] Lee, M., Park, J.: An optimized dynamic mode decomposition model robust to multiplicative noise. SIAM Journal on Applied Dynamical Systems 22(1), 235– 268 (2023) [19] Bollt, E.M., Li, Q., Dietrich, F., Kevrekidis, I.: On matching, and even rectifying, dynamical systems through Koopman operator eigenfunctions. SIAM Journal on Applied Dynamical Systems 17(2), 1925–1960 (2018) [20] Bakker, C., Ramachandran, T., Rosenthal, W.S.: Learning bounded Koopman observables: Results on stability, continuity, and controllability. arXiv preprint arXiv:2004.14921 (2020) [21] Dietrich, F., Thiem, T.N., Kevrekidis, I.G.: On the Koopman operator of algorithms. SIAM Journal on Applied Dynamical Systems 19(2), 860–885 (2020) [22] Rowley, C.W., Mezić, I., Bagheri, S., Schlatter, P., Henningson, D.S.: Spectral analysis of nonlinear flows. Journal of Fluid Mechanics 641, 115–127 (2009) https: //doi.org/10.1017/S0022112009992059 [23] Mezić, I.: Spectral Properties of Dynamical Systems, Model Reduction and Decompositions. Nonlinear Dynamics 41(1), 309–325 (2005) https://doi.org/10. 1007/s11071-005-2824-x [24] Mezić, I.: Analysis of Fluid Flows via Spectral Properties of the Koopman Operator. Annual Review of Fluid Mechanics 45(1), 357–378 (2013) https://doi.org/ 10.1146/annurev-fluid-011212-140652 [25] Stoer, J., Bulirsch, R.: Introduction to Numerical Analysis, 3rd edn. Springer, ??? (2010). https://doi.org/10.1007/978-0-387-21738-3 [26] Dormand, J.R., Prince, P.J.: A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics 6(1), 19–26 (1980) https:

43

//doi.org/10.1016/0771-050X(80)90013-3 [27] Hirsch, M.W., Smith, H.: Chapter 4 Monotone Dynamical Systems. In: Handbook of Differential Equations: Ordinary Differential Equations vol. 2, pp. 239–357. Elsevier, ??? (2006). https://doi.org/10.1016/S1874-5725(05)80006-9 [28] Langfield, P., Krauskopf, B., Osinga, H.M.: Forward-Time and Backward-Time Isochrons and their interactions. SIAM Journal on Applied Dynamical Systems 14(3), 1418–1453 (2015) https://doi.org/10.1137/15m1010191 [29] Derevianko, N., Kevrekidis, I.G., Dietrich, F.: Neural network-based singularity detection and applications. arXiv (arXiv:2509.10110) (2025) https://doi.org/10. 48550/arXiv.2509.10110 2509.10110 [30] Boullé, N., Nakatsukasa, Y., Townsend, A.: Rational neural networks. Advances in neural information processing systems 33, 14243–14253 (2020) [31] Cornfeld, I.P., Fomin, S.V., Sinai, Y.G.: Ergodic Theory, 1st edn. Springer, ??? (2012). https://doi.org/10.1007/978-1-4615-6927-5

44

Appendix A Ergodic dynamical systems Ergodic systems provide a setting in which the Koopman operator is well-studied. For these dynamical systems, the evolution function T is invertible and measurepreserving, i.e., there exists an invariant measure µ such that for any S ∈ B, µ(S ) = µ(F −1 (S )), where F −1 (S ) is understood as the pre-image of the set S ∈ B, i.e., F −1 (S ) := {x ∈ M | F (x) ∈ S}. Ergodic theory then describes statistical properties of such dynamical systems and their long-term behavior. The theory is developed in the context of measure-preserving transformations in measure theory. Ideas related to the Koopman operator are based on dynamical systems theory (work from Budišić, Mohr, and Mezić [8]), specifically, ergodic theory (work from Cornfeld, Fomin, and Sinai [31]). Ergodic theory has applications in physics, probability theory, and information theory. We now briefly introduce ergodic theory following Cornfeld, Fomin, and Sinai [31]. A measure-preserving function g is called invariant with respect to the evolution function F if for all x ∈ M, g (F (x)) = g (x) = g (F −1 (x)). If this is true almost everywhere instead of for all x ∈ M, namely, if it is true for x ∈ M \ {B ∈ B | µ(B ) = 0}, then it is said to be invariant mod 0. Analogously, a set A ∈ B is called invariant, or invariant mod 0, if the indicator function χA (x) (with χA (x) = 1 if x ∈ A and χA (x) = 0 otherwise) is an invariant function, or an invariant mod 0 function, respectively. A dynamical system (M, B, µ, F ) is called ergodic if the measure µ(A) of any invariant set A equals 0 or 1. There are several equivalent conditions of ergodicity, and we use the following definition.

Definition 9. (Equivalent definition of ergodic system) A dynamical system (M, B, µ, F ) is called ergodic if any function f ∈ L2 (M, B, µ) that is invariant with respect to the Koopman operator KF is constant almost everywhere.

B Algorithms B.1 Eigensolvers Algorithm 4 Deflation based algorithm for real non-Hermitian matrix A with real eigenvalues i ← 0. while i < n do λi , v i ← power iteration(A). λi , wi ← power iteration(AT ). wi wi ← T . wi v i A ← A − λi v i wTi . i ← i + 1. end while

45

Algorithm 5 (power iteration complex) Arnoldi-based power iteration algorithm for real non-Hermitian matrix A with complex dominant eigenvalues λmax , λ̄max and eigenvectors v, v̄ , given tolerance and maximum iterations N x ← power iteration(A, max iterations = 500). λo ← ∞. i ← 0. while i < N do V [:, 0] ← x. h, V ← modified Gram-Schmidt(A, V, degree = 2). λ1 , λ2 ← eigenvalue2D(h). k ← arg max1,2 {|λ1 |, |λ2 |}. λmax ← λk w ← eigvector2D(h, λmax ). v ← V [: 0 : 2]w. v . v← ∥v∥ if |λmax − λo | < tolerance then break. else λo ← λmax . end if i ← i + 1. end while

C Proofs of Propositions C.1 Proof of Proposition 2 Proof Using the definition of Koopman operator, ϕ(F ∆t (x)) = ϕ(x∆t + ε(x)) ≈ ϕ(x∆t ) + ε(x)T ∇ϕ(x∆t )

= λ∆t ϕ(x) + ε(x)T ∇ϕ(x∆t ).

Using the bound on the spectral norm on the Jacobian of the dictionary basis Ψ, we get ∥∇ϕ(x)∥ = ∥JΨ (x)w∥ ≤ ∥JΨ (x)∥2 ∥w∥ ≤ L, ∀x ∈ G,

(77)

where ∥.∥2 is the spectral norm. Then using the above equations, the Cauchy-Schwarz inequality, the Binomial theorem, and ϕp ≡ ϕp , we get 2

2

EF G (ϕ̄p , λ̄p )2p = ϕ̄p (F ∆t (·)) − λ̄p ϕ̄p (·) = ϕ̄(F ∆t (·))p − (λ̄∆t ϕ̄(·))p G G  2 1 X p ∆t p ∆t = (ϕ̄(F (x))) − (λ̄ ϕ̄(x))) |G| x∈G  2 1 X = (λ̄∆t ϕ̄(x) + ε(x)T ∇ϕ̄(x∆t ))p − (λ̄∆t ϕ̄(x)))p |G| x∈G

46

! 2 p (λ̄∆t ϕ̄(x))p−k (ε(x)T ∇ϕ̄(x∆t ))k − (λ̄∆t ϕ̄(x)))p k x∈G k=0 !  p 2 1 X X p (λ̄∆t ϕ̄(x))p−k (ε(x)T ∇ϕ̄(x∆t ))k = |G| k x∈G k=1 !  2 p 1 X X p = (λ̄∆t wT Ψ(x))p−k (ε(x)T ∇ϕ̄(x∆t ))k |G| k x∈G k=1 !   p k 2 1 X X p ≤ (λ̄∆t )p−k ∥w∥p−k ∥Ψ(x)∥p−k ∥ε(x)∥k ∇ϕ̄(x∆t ) |G| k x∈G k=1 ! 2  p 1 X X p ≤ (λ̄∆t )p−k M p−k ∥ε(x)∥k Lk |G| k x∈G k=1 !  2 p 1 X X p (λ̄∆t )p−k M p−k ∥ε(x)∥k Lk − (λ̄∆t )p M p ≤ |G| k x∈G k=0 2 
p 1 X λ̄∆t M + L ∥ε(x)∥ − (λ̄∆t )p M p ≤ |G| 1 X = |G|

X p

x∈G

λ̄∆t M + L ∥ε(·)∥

p

2

− (λ̄∆t M )p G  2
p λ̄∆t M + LϵG − (λ̄∆t M )p . ≤

=

□

C.2 Proof of Proposition 3 Proof Using the identity (ap − bp )2 = (a − b)2 inequality, we get 2

p−1−i i
2 b , and the Cauchy–Schwarz i=0 a

Pp−1

EF G (ϕ̄p , λ̄p )2p = ϕ̄p (F (·)) − λ̄p ϕ̄p (·) G
p
p 2 = wcT Ψ(F (·) − λ̄wcT Ψ(·) G  2 1 X T p T = (wc Ψ(F (x))) − (wc λ̄Ψ(x))p |G| x∈G

=

1 X wcT Ψ(F (x)) − wcT λ̄Ψ(x))2 |G| x∈G

 p−1 X

(wcT Ψ(F (x)))p−1−i (wcT λ̄Ψ(x))i

i=0

2

 p−1 2 X T 1 X = (δwT (Ψ(F (x)) − λ̄Ψ(x))2 (wc Ψ(F (x)))p−1−i (wcT Ψ(x))i λ̄i |G| i=0

x∈G

≤

1 X 2 ∥δw∥2 Ψ(F (x)) − λ̄Ψ(x) |G| x∈G

 p−1 X i=0

∥wc ∥p−1−i ∥Ψ(F (x))∥p−1−i ∥wc ∥i ∥Ψ(x)∥i λ̄i

 p−1 2 X 1 X 2 = ∥δw∥2 Ψ(F (x)) − λ̄Ψ(x) ∥Ψ(F (x))∥p−1−i ∥Ψ(x)∥i λ̄i |G| i=0

x∈G

47

2

2

= ∥δw∥

Ψ(F (·)) − λ̄Ψ(·)

= CT G (p, λ̄)2 ∥δw∥2 .

 p−1 X i=0

p−1−i

∥Ψ(F (·))∥

i

∥Ψ(·)∥ λ̄

i

 2

G

□

D Analysis of real powers of complex numbers Recall that the real-exponential of a complex number is defined for α ∈ R and z = reiθ ∈ C, r, θ ∈ R as z α := eα log(z) , where the complex logarithm log(z ) is defined as log(z ) = ln(r) + i(θ + 2πk ), k ∈ Z. Note that for r = 1, log(z ) = i(θ + 2πk ) ≡ θ (mod 2π ) but log(ex ) ̸= x in general. Furthermore, consider eikαx = eixπ+2nπ , n ∈ Z. For (eikx )α , take θ = arg eikx = kx + 2nπ ∈ [−π, π ], and return eiθα . However, for α ∈ R \ Z, ′

eiαxπ = (e(ixπ+2nπ) )α = eiαxπ+2αnπ = eiαxπ+2n π , n′ = αn ∈ R \ Z. Thus, the angle is not unique (mod 2π ). The following is a small numerical experiment for calculating eikα , where i is the imaginary unit, k ∈ Z, and α ∈ R \ Z. As an example, consider k = 4 and α = 0.5. Define the following three functions in Python:

• e1(x) = eikαx • e2(x) = (eiαx )k • e3(x) = (eikx )α These functions should represent the same function, however, the function, which does real power later, shows a different function as is shown in the following image. Furthermore, consider for f (x) = e3(x) = (eikx )α . Since eikx is 2π k periodic function ikx 2π (with respect to x), f (x) is also k periodic. By definition, f (x) = eα log(e ) = eiα(kx+2nπ) with n ∈ Z, kx + 2nπ ∈ [π, π ). Thus, consider the value of f (x) in one period, for m ∈ Z, 2π 2π m≤x< (m + 1) k k ∴ 2πm ≤ kx < 2π (m + 1)

∴ 2π (m + n) ≤ kx + 2πn < 2π (m + n + 1) ∴ 2π (m + n)α ≤ (kx + 2πn)α < 2π (m + n + 1)α

∴ e2π(m+n)α ≤ e(kx+2πn)α = f (x) < e2π(m+n+1)α . Here, for small α ≈ 0, both e2π(m+n)α and e2π(m+n+1)α are close to 1, and thus, the values f (x) are close to 1. Hence, for small α, f (x) always takes value close to 1. 48

λ0,2 = 0.98410 .0.98222 = 0.9646

λ1,1 = 0.98411 .0.98221 = 0.9666

λ0,4 = 0.98410 .0.98224 = 0.9305

x2

x2

x2

E Additional figures for the non-linear system (47) in Sect. 5.2 2.0

1.5

2.0

1.0

1.5

2.0

1.0

1.0 1.0

1.5

2.0

1.0

1.5

2.0

x1 λ1,3 = 0.98411 .0.98223 = 0.9324

x1 λ2,2 = 0.9841 .0.98222 = 0.9343

x2

2.0

x1 λ2,0 = 0.98412 .0.98220 = 0.9685

x2

1.5

x2

1.0

1.5

2.0

1.5

2.0

1.0

1.5

2

2.0

1.0

1.0 1.0

1.5

2.0

1.0

1.5

2.0

x1 λ4,0 = 0.98414 .0.98220 = 0.9380

x1 λ5,0 = 0.98415 .0.98220 = 0.9231

x2

2.0

x1 λ3,1 = 0.98413 .0.98221 = 0.9361

x2

1.5

x2

1.0

1.5

2.0

1.5

2.0

1.0

1.5

2.0

1.0

1.0 1.0

1.5

2.0

1.0

1.5

2.0

x1 λ3,2 = 0.98413 .0.98222 = 0.9194

x1 λ2,3 = 0.98412 .0.98223 = 0.9176

x2

2.0

x1 λ4,1 = 0.98414 .0.98221 = 0.9213

x2

1.5

1.5

x2

1.0

2.0

1.5

2.0

1.0

1.5

0.0

1.0 1.5

2.0

1.0

1.5

2.0

x1 λ1,4 = 0.98411 .0.98224 = 0.9158

x1 λ2,1 = 0.98412 .0.98221 = 0.9512

x2

1.0

x1 λ3,0 = 0.98413 .0.98220 = 0.9531

x2

2.0

1.5

x2

1.5

2.0

1.5

2.0

1.0

1.5

2.0

1.0

1.0 1.0

1.5

2.0

1.0

1.5

2.0

x1 λ1,2 = 0.98411 .0.98222 = 0.9493

x1 λ0,3 = 0.98410 .0.98223 = 0.9474

x2

2.0

x1 λ0,5 = 0.98410 .0.98225 = 0.9139

x2

1.5

2.0

1.5

2.0

1.0

1.5

1.0 2.0

1.5

2.0

x1 λ0,1 = 0.98410 .0.98221 = 0.9822

2.0

1.5

1.0

−1.0

2.0

1.5

1.0 1.0

x1 λ1,0 = 0.98411 .0.98220 = 0.9841

x2

1.5

x2

1.0

−0.5

1.5

x2

1.0

0.5

2.0

1.0 1.0

1.0

1.0

1.5

2.0

x1

2.0

1.5

1.0 1.0

1.5

x1

2.0

1.0

1.5

2.0

x1

Fig. 24: Explicit eigenfunctions for the non-linear system (47) computed using ϕ◦h−1 .

49

x2

λ = 0.99993 + 0.00000j

λ = 0.98198 + 0.00032j

λ = 0.96132 + 0.00000j

2.0

λ = 0.98198 − 0.00032j

λ = 0.97169 + 0.00000j

2.0

2.0

2.0

2.0

1.5

1.5

1.5

1.5

1.5

1.0

1.0 1.0

1.5 x1

2.0

x2

2.0

1.0 1.0

1.5 x1

2.0

λ = 0.95606 + 0.03288j

2.0

2.0

1.5

1.5

1.0

1.0 1.0

x2

1.5 x1

λ = 0.95909 − 0.01175j

λ = 0.95909 + 0.01175j

1.5 x1

2.0

1.5 x1

2.0

1.5

1.5 x1

2.0

1.0

1.5 x1

2.0

2.0

2.0

2.0

1.5

1.5

1.5

1.0 1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

λ = 0.93423 + 0.05314j

2.0

λ = 0.94083 − 0.02258j 2.0

2.0

λ = 0.93423 − 0.05314j

1.5

1.5

1.5

1.5

1.0 1.0

1.5 x1

λ = 0.94369 + 0.00375j

2.0

λ = 0.94083 + 0.02258j

1.0 1.0

λ = 0.95606 − 0.03288j

1.0 1.0

λ = 0.94369 − 0.00375j

1.0 2.0

1.0 1.0

1.5 x1

2.0

2.0

1.0 1.0

1.5 x1

2.0

1.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

λ = 0.92944 − 0.02835j

λ = 0.92899 + 0.07828j

2.0

2.0

λ = 0.92899 − 0.07828j

λ = 0.92532 + 0.03684j

2.0

2.0

2.0

1.5

1.5

1.5

1.5

1.5

λ = 0.92944 + 0.02835j

x2

1.0 1.0

0.5 1.0

1.0

x2

1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

λ = 0.92532 − 0.03684j

λ = 0.92396 + 0.00000j

λ = 0.91014 + 0.07986j

2.0

2.0

λ = 0.91014 − 0.07986j

λ = 0.90992 + 0.04540j

2.0

2.0

2.0

1.5

1.5

1.5

1.5

1.5

1.0

1.0 1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

λ = 0.90992 − 0.04540j

λ = 0.90005 + 0.09108j

2.0

1.5

1.0 1.0

1.5 x1

2.0

0.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

λ = 0.89664 + 0.06120j

2.0

λ = 0.90005 − 0.09108j 2.0

2.0

λ = 0.89664 − 0.06120j

1.5

1.5

1.5

1.5

2.0

x2

−0.5

1.0

1.0 1.0

1.5 x1

2.0

x2

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

λ = 0.89173 + 0.10864j

2.0

λ = 0.89495 − 0.12371j 2.0

1.5

1.5

λ = 0.89495 + 0.12371j

1.0

1.0 1.0

x2

1.0 1.0

1.5 x1

2.0

1.5 x1

1.5 x1

2.0

1.0

1.5 x1

2.0

λ = 0.88893 + 0.09273j

2.0

λ = 0.89173 − 0.10864j 2.0

2.0

1.5

1.5

1.5

1.0 1.0

1.0 1.0

2.0

1.0 1.0

1.5 x1

2.0

−1.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

λ = 0.88893 − 0.09273j

λ = 0.88685 + 0.00460j 2.0

λ = 0.88685 − 0.00460j

λ = 0.88525 + 0.14987j

2.0

2.0

2.0

λ = 0.88525 − 0.14987j

1.5

1.5

1.5

1.5

1.5

1.0

1.0 1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

2.0

1.0 1.0

1.5 x1

2.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

x2

λ = 0.00000 + 0.00000j λ = −0.00000 + 0.00000j 2.0

2.0

1.5

1.5

1.0

1.0 1.0

1.5 x1

2.0

1.0

1.5 x1

2.0

Fig. 25: Eigenfunctions for the non-linear system (47) approximated using EDMD.

50

λ11 = 0.99993 + 0.00000j

λ21 = 0.99987 + 0.00000j

λ31 = 0.99980 + 0.00000j

1.5 1.0

2.0 x2

2.0 x2

x2

2.0

1.5 1.0

1.5 1.0

1.0 1.5 2.0

1.0 1.5 2.0

1.0 1.5 2.0

x1

x1

λ12 = 0.98198 + 0.00032j

λ22 = 0.96429 + 0.00062j

x1 λ32 = 0.94691 + 0.00092j

1.5

1.5

1.5

1.0 1.5 2.0

x1

x1

x1

λ13 = 0.98198 − 0.00032j

λ23 = 0.96429 − 0.00062j

λ33 = 0.94691 − 0.00092j x2

1.0

1.0 1.5 2.0

x2

1.0

x2

1.0

2.0 x2

2.0 x2

x2

2.0

2.0 1.5

2.0

2.0

1.0 1.5 2.0

1.0 1.5 2.0

x1

x1

λ14 = 0.97169 − 0.00000j

λ24 = 0.94418 − 0.00000j

x1 λ34 = 0.91745 − 0.00000j x2

1.0

1.0 1.5 2.0

x2

1.0

1.5

x2

1.0

1.5

1.0 1.5 2.0

2.0 1.5

2.0

1.0

1.5

2.0

1.0

1.5 0.5

1.0

x1

x1

x1

λ15 = 0.96132 − 0.00000j

λ25 = 0.92413 − 0.00000j

λ35 = 0.88839 − 0.00000j x2

1.0 1.5 2.0

x2

1.0 1.5 2.0

x2

1.0 1.5 2.0

2.0 1.5

2.0

1.0

1.5

2.0

1.0

1.5

0.0

1.0

1.0 1.5 2.0

1.0 1.5 2.0

1.0 1.5 2.0

x1 λ16 = 0.95909 + 0.01175j

x1 λ26 = 0.91971 + 0.02253j

x1 λ36 = 0.88181 + 0.03241j

1.5

1.5 1.0 1.5 2.0

1.0 1.5 2.0

x1 λ17 = 0.95909 − 0.01175j

x1 λ27 = 0.91971 − 0.02253j

x1 λ37 = 0.88181 − 0.03241j x2

1.0

1.0 1.5 2.0

x2

1.0

−0.5

1.5

x2

1.0

2.0 x2

2.0 x2

x2

2.0

2.0 1.5

2.0

1.0

1.5

2.0

1.0

1.0 1.0 1.5 2.0

x1

x1

x1

λ18 = 0.95606 + 0.03288j

λ28 = 0.91297 + 0.06287j

λ38 = 0.87079 + 0.09013j

1.5 1.0 1.5 2.0

1.0 1.5 2.0

x1

x1

λ19 = 0.95606 − 0.03288j

λ29 = 0.91297 − 0.06287j

x1 λ39 = 0.87079 − 0.09013j x2

1.0

1.0 1.5 2.0

x2

1.0

1.5

x2

1.0

2.0 x2

x2

x2

1.0 1.5 2.0

2.0

1.5

2.0 1.5 1.0

−1.0

1.5

1.0 1.5 2.0

2.0

1.0

2.0 1.5 1.0

2.0 1.5 1.0

1.0 1.5 2.0

1.0 1.5 2.0

1.0 1.5 2.0

x1

x1

x1

Fig. 26: Extended eigenfunctions for the first 9 eigenvalues computed using Algorithm 3 for the non-linear system (47).

51

F Separatrices for multistable systems Consider the unforced Duffing equation

ẍ = −δ ẋ − x(β + αx2 ),

(78)

with δ = 0.5, β = −1, and α = 0.1. For this setting, there exist two stable spirals at (±3.1623, 0) and a saddle at the origin. Thus, nearly all initial conditions, excluding those on the stable manifold of the saddle point, will be attracted to one of the spiral equilibria. The stable manifold of the saddle point is plotted as a blue curve, which separates the basins of attraction of the two stable spirals. The unstable manifold of the saddle is depicted in orange, along which the two stable spirals lie (see Fig. 27).

Fig. 27: Phase portrait of the unforced Duffing system (y = ẋ), with two stable spirals and one saddle. The stable and unstable manifolds of the saddle point are depicted as blue and orange curves respectively.

Fig. 28: Koopman spectrum computed using EDMD for the unforced Duffing system (78). Orange dots indicate real eigenvalues for which Im(λ) = 0.

52

Fig. 29: Pseudocolor plots of the first twenty eigenfunctions for the unforced Duffing system (78) approximated using EDMD.

We utilize EDMD to approximate the Koopman eigenfunctions associated with the attractors at (±3.1623, 0). Following [10], we use a dataset consisting of 3 × 103 trajectories, each with 11 samples taken at a sampling interval of ∆t = 0.25. The data 4 points are represented as X, Y ∈ R6×10 , with initial conditions uniformly distributed over x, ẋ = y ∈ [−6, 6]. We used a dictionary that included 500 radial basis functions (RBFs). The RBF centers were selected using k-means clustering on the entire data set. Fig. 28 and Fig. 29 illustrate the Koopman spectrum and the first twenty eigenfunctions for the unforced Duffing system(78). For the Koopman eigenfunction corresponding to the stable manifold, the magnitude and phase of the eigenfunction parameterize the basin of attraction for each stable spiral. But, the analytical eigenfunctions go to infinity at the boundary between the basins [1].

53

Fig. 30: Sampled points S on the unstable manifold of the saddle points with (x, y ) ∈ W u (0, 0) ∩ [−2, 2] × [−1.33, 1.3].

Fig. 31: Evaluated eigenfunctions ϕ associated with real eigenvalues (Im(λ) = 0) except for λ = 1 over the set S . 54

Theoretically, the eigenfunctions associated with real eigenvalues also approach infinity along the unstable manifold at the saddle point. To show this, first we accurately sample the unstable manifold W u (0, 0) around the saddle with 100 points S = {(xi , yi ) : (xi , yi ) ∈ W u (0, 0) ∩ [−2, 2] × [−1.33, 1.3], i = 1, . . . , 100} as shown in Fig. 30. Then, we evaluate the previously computed Koopman eigenfunctions (corresponding to the attractors) over the set S . These evaluated eigenfunctions associated with real eigenvalues ( Im(λ) = 0) except for λ = 1, are plotted in Fig. 31.

G Additional figures for the bi-stable system (70) in Sect. 5.5.2

Fig. 32: Four Koopman eigenfunctions for the system (70), plotted over the separatrix (y1 coordinate) or the connection between the two steady states through the saddle (y2 coordinate). From the left, the first function is zero on the connection and does not diverge. The second function diverges at the separatrix and is zero on the two steady states. The third function is the inverse of the second. The fourth function is the inverse of the product of the first and the second.

55