# A spectral characterization of nonlinear normal modes

Giuseppe Cirillo, Alexandre Mauroy, Ludovic Renson, Gaëtan Kerschen, Rodolphe Sepulchre

► To cite this version:

Giuseppe Cirillo, Alexandre Mauroy, Ludovic Renson, Gaëtan Kerschen, Rodolphe Sepulchre. A spectral characterization of nonlinear normal modes. Journal of Sound and Vibration, 2016, 377, pp.284-301. 10.1016/j.jsv.2016.05.016. hal-01336227

HAL Id: hal-01336227

https://hal.science/hal-01336227

Submitted on 31 May 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A spectral characterization of nonlinear normal modes

G.I. Cirillo<sup>a,*</sup>, A. Mauroy<sup>a</sup>, L. Renson<sup>b</sup>, G. Kerschen<sup>c</sup>, R. Sepulchre<sup>d</sup>

<sup>a</sup>Systems and Modeling research group, Department of Electrical Engineering and Computer Science, University of Liège, Liège, Belgium.

<sup>b</sup>Department of Engineering Mathematics, University of Bristol, Bristol, UK.

<sup>c</sup>Space Structures and Systems Lab (S3L), Department of Aerospace and Mechanical Engineering, University of Liège, Liège, Belgium.

<sup>d</sup>Control Group, Department of Engineering, University of Cambridge, Cambridge, UK.

## Abstract

This paper explores the relationship that exists between nonlinear normal modes (NNMs) defined as invariant manifolds in phase space and the spectral expansion of the Koopman operator. Specifically, we demonstrate that NNMs correspond to zero level sets of specific eigenfunctions of the Koopman operator. Thanks to this direct connection, a new, global parametrization of the invariant manifolds is established. Unlike the classical parametrization using a pair of state-space variables, this parametrization remains valid whenever the invariant manifold undergoes folding, which extends the computation of NNMs to regimes of greater energy. The proposed ideas are illustrated using a two-degree-of-freedom system with cubic nonlinearity.

## 1. Introduction

Nonlinear normal modes (NNMs) of vibration can be regarded as the extension of the concept of linear normal modes to nonlinear systems. They have proved useful in a number of applications going from model reduction [1] and the study of localization phenomena [2] to system identification [3, 4]. NNMs of undamped systems were first defined as a vibration of unison of the system [5]. The NNM definition was later generalized to encompass modal interactions [6]. Inspired by the center manifold technique, Shaw and Pierre extended the NNM concept to damped systems [7]. They defined NNMs as two-dimensional invariant manifolds in phase space which are tangent to the planes formed by the linear modes at the equilibrium point. Such manifolds are invariant under the flow (i.e. trajectories initialized in the manifold do not leave the manifold), a property which extends the invariance of linear modes to nonlinear systems.

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: gcirillo@ulg.ac.be (G.I. Cirillo), a.mauroy@ulg.ac.be (A. Mauroy), l.renson@bristol.ac.uk (L. Renson), g.kerschen@ulg.ac.be (G. Kerschen), r.sepulchre@eng.cam.ac.uk (R. Sepulchre)</footnote>
Operator-theoretic methods have a long and rich tradition in dynamical systems. In particular, a (nonlinear) system can be studied through the spectral properties of the associated (linear) Koopman operator. Koopman operators were first defined for Hamiltonian systems in [8]. In the seminal work of Mezić [9] the spectral properties of the operator were connected to geometrical features of ergodic (i.e. conservative) systems. More recently, this spectral approach was extended to the case of dissipative systems in the context of model reduction [10] and stability analysis [11]. The Koopman operator framework is also conducive to data analysis (see e.g. Dynamic Mode Decomposition [12]) and was successfully used to capture coherent structures (called Koopman modes) or detect instabilities in various fields (e.g. fluid flows [13], power grids [14], energy-efficient buildings [15]).

Our aim is to explore the connections that exist between NNMs and the Koopman operator framework. In this direction, a fundamental result reported by Mezić in [9] is that normal modes of linear oscillations have their natural analogs, termed the Koopman modes, in the context of nonlinear dynamics. Specifically, it was shown that, for a hyperbolic fixed point, the dynamics can be linearized in the entire basin of attraction [16] and that the coordinates for the linearization can be obtained through the eigenfunctions of the Koopman operator [17]. This framework therefore seems very appropriate to extend the concept of linear normal modes to nonlinear systems. Indeed, we will demonstrate that there is a direct connection between NNMs and the spectral expansion of the Koopman operator. More precisely, the invariant manifolds are related to the zero level sets of some eigenfunctions of the operator, in direct analogy with the linear case.

Using this new characterization, we will show how issues related to the parametrization of the invariant manifold can be overcome. In the original work of Shaw and Pierre, a single pair of variables (either in physical or modal coordinates) was used to obtain a parametrization of the manifold, which can possibly hold only locally. For instance, the parametrization may fail when the invariant manifold presents a complex geometry with foldings, as shown in [18, 19]. In contrast, interpreting the concept of NNMs in the Koopman operator framework allows us to obtain a parametrization that is valid globally.

The paper is organized as follows. Sections 2 and 3 review the theory of NNMs and of the Koopman operator, respectively. In Section 4 the connection between the two frameworks is investigated, first in the linear case and then in the nonlinear case. We also derive a partial differential equation describing the geometry of the invariant manifold and an approximation of the manifold using Taylor series. The theoretical developments are illustrated with a two-degree-of-freedom mechanical system with cubic nonlinearity in Section 5. The conclusions of the present study are summarized in Section 6.

## 2. Introduction to NNMs

Consider the first-order equations of motion of an N-degree-of-freedom nonlinear mechanical system

$$ \dot{x}_i = y_i, \quad \dot{y}_i = f(\mathbf{x}, \mathbf{y}) \quad i = 1, \dots, N \quad (1) $$

where $x_i$ represents a generalized coordinate (displacement or rotation), $y_i$ is the corresponding velocity, $\mathbf{x} = [x_1, \dots, x_N]$ and $\mathbf{y} = [y_1, \dots, y_N]$. From this point on, we assume that
system (1) is damped, $(\mathbf{x}, \mathbf{y}) = (0, 0)$ is an equilibrium point, and the system linearized around $(0, 0)$ has $n/2$ pairs of complex-conjugate eigenvalues with non-zero real parts (hyperbolic equilibrium point).

In the 1990s, Shaw and Pierre defined a NNM as a two-dimensional manifold in phase space [7] that is invariant under the flow, i.e. trajectories with initial conditions in the manifold remain in it for all time. As such, the invariance property of linear normal modes was extended to nonlinear systems: a motion that starts on a mode stays on this mode, while the other modes remain quiescent. Inspired by the center manifold technique [20], Shaw and Pierre mathematically described the two-dimensional invariant manifold using a pair of master coordinates, the other state-space coordinates being functionally related to the master pair. Choosing arbitrarily $(x_i, y_i)$ as master coordinates, the remaining coordinates follow the constraint equations

$$x_j = X_j(x_i, y_i), \quad y_j = Y_j(x_i, y_i), \quad j = 1, \dots, N; \quad j \neq i. \qquad (2)$$

In several other works (see, for instance, [21]), a linear change of coordinates was used to express Eq. (2) in modal space and describe the invariant manifold as a function of a linear modal coordinate and its corresponding modal velocity.

A set of equations for the functions $X_j$ and $Y_j$ can be derived by substituting the time derivative of Eq. (2) into the equations of motion (1). This leads to a set of $2N - 2$ coupled partial differential equations (PDEs)

$$\begin{aligned} Y_j(x_i, y_i) &= \frac{\partial X_j(x_i, y_i)}{\partial x_i} y_i + \frac{\partial X_j(x_i, y_i)}{\partial y_i} f_i, \\ f_j &= \frac{\partial Y_j(x_i, y_i)}{\partial x_i} y_i + \frac{\partial Y_j(x_i, y_i)}{\partial y_i} f_i, \end{aligned} \qquad (3)$$

with $f_j = f_j(x_i, \mathbf{X}, y_i, \mathbf{Y})$, $\mathbf{X} = \{X_j(x_i, y_i)\}$, $\mathbf{Y} = \{Y_j(x_i, y_i)\}$, and $j = 1, \dots, N; \ j \neq i$. The PDEs (3) do not depend on time and can be solved for $\mathbf{X}$ and $\mathbf{Y}$. Finally, after a substitution of the solutions in the ordinary differential equations governing the master coordinates, the dynamics on the NNM is reduced to a single-degree-of-freedom (SDOF) nonlinear oscillator dynamics

$$\dot{x}_i = y_i, \quad \dot{y}_i = f_i(x_i, \mathbf{X}, y_i, \mathbf{Y}). \qquad (4)$$

The PDEs (3) admit several solutions that correspond to the extension of the underlying linear normal modes. The constraint relations $(\mathbf{X}, \mathbf{Y})$ describe the geometry of a single NNM in phase space. For a linear system, the mode corresponds to a plane whose dependence in $x_i$ and $y_i$ is function of the inertia, stiffness and damping properties of the system. For a nonlinear system, the mode is no longer flat but curved due to nonlinear distortions. At the equilibrium point, the NNM is tangent to the plane formed by the linear mode of the underlying linear system. Figure 1 illustrates the in-phase NNM of a two-degree-of-freedom system with cubic nonlinearity. The surface shown in the figure corresponds to the position of the second DOF $x_2$ as a function of the master coordinates $(x_1, y_1)$. The solutions of (3) can be approximated in the form of a polynomial series expansion in $x_i$ and $y_i$ [7], whose
coefficients are the solutions of a set of algebraic equations. More recently, several numerical methods were proposed in the literature and applied to various conservative [21, 22] and nonconservative [19] systems, see [23] for a review.

<figure><img src="image_2.png" /><figcaption>Figure 1: In-phase NNM of a two-degree-of-freedom system with cubic nonlinearity. The manifold was computed using the Taylor series expansion presented in Section 4.4.</figcaption></figure>

Either derived using state space or modal coordinates, Eq. (2-4) are based on the fundamental underlying assumption that the invariant manifold can be described as a function of the chosen pair of variables. However, the invariant manifold may present a complex geometry with several foldings that appear for a number of reasons including internal resonances, localization, and multiple fixed points [23, 24, 25]. Such a manifold cannot generally be described using a simple parametrization that is linear in the state space coordinates. For instance, the invariant manifold shown in Figure 1 exhibits foldings, which prevent a parametrization in terms of the master coordinates $(x_1, y_1)$. In this case, a proper parametrization cannot be obtained as a solution of the PDEs (3).

In Ref. [24], the parametrization issue was addressed using a larger set of master coordinates. When considering $n_m$ pairs, the slave coordinates can be described as

$$x_j = X_j(\mathbf{x}_m, \mathbf{y}_m), \quad y_j = Y_j(\mathbf{x}_m, \mathbf{y}_m), \qquad (5)$$

where $(\mathbf{x}_m, \mathbf{y}_m)$ are the vectors of master coordinates. After solving the corresponding manifold-governing PDEs, the dynamics is reduced to $m$ coupled (nonlinear) oscillators. Although this method is elegant and proved effective on a nonlinear beam [24, 25], it does
not completely solve the intrinsic parametrization issue and may therefore fail in some regions of the phase space.

In this paper the parametrization issue is addressed using the relation that exists between NNMs and Koopman operator. Through a nonlinear change of coordinates, the dynamics on the NNM will be expressed as a single-degree-of-freedom linear oscillator. This novel approach bears strong resemblance with the normal form method [1, 26], whose main idea is to simplify the vector field using a series of coordinate changes. Unlike the normal form method which has only local validity, the Koopman operator framework provides a change of coordinates that is valid in the entire basin of attraction of a hyperbolic fixed point [16]. Therefore, it is possible to think of methods based on Koopman operator as a generalization of methods based on normal form, with the notable difference that the system of coordinates obtained in this case is valid in the whole basin of attraction of a fixed point. By doing so the parametrization issue is overcome, and coupled PDEs describing the NNM in the entire basin of attraction of the equilibrium point can be obtained.

## 3. Koopman Operator

### 3.1. Motivation and definition

<figure><img src="image_3.png" /><figcaption>Figure 2: 2DOF system with cubic nonlinearity.</figcaption></figure>

The system depicted in Figure 2 is considered in order to illustrate the Koopman operator definition and properties. The system dynamics can be expressed through a system of ordinary differential equations (ODEs)

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}), \quad \mathbf{x} = (x_1, x_2, x_3, x_4)^T; \qquad (6) $$

where $x_1$ and $x_2$ are, respectively, the positions of the first and second mass, and $x_3$ and $x_4$ the corresponding velocities. The general solution, or flow, of (6) is a function, $\phi(\mathbf{x}_0, t)$, that assigns a point in state space to an initial condition $\mathbf{x}_0$ and a time $t$.

It is often interesting to analyze the dynamics of a system in terms of specific quantities such as, for instance, the relative displacement between the first and second mass $g_1(\mathbf{x}) = x_2 - x_1$, the potential energy $g_2(\mathbf{x}) = \frac{1}{2}(k_a(x_1^2 + x_2^2) + k_b(x_1 - x_2)^2) + \frac{1}{4}k_{nl}x_1^4$ or the kinetic energy $g_3(\mathbf{x}) = \frac{1}{2}(m_1x_3^2 + m_2x_4^2)$. All these quantities are functions from the state space to $\mathbb{R}$ (or more generally $\mathbb{C}$) and are called *observables*.
Since state-space variables change with time, so do these observables. Given initial conditions $x_0$, the time evolution of an observable is obtained through its composition with the flow. For instance, the potential energy varies according to $V(x_0, t) = g_2(\phi(x_0, t))$, while the kinetic energy is given by $K(x_0, t) = g_3(\phi(x_0, t))$. At every time instant $t$, $K(x_0, t)$ and $V(x_0, t)$ can still be considered as observables (i.e. functions from the state space to $\mathbb{C}$) such that, generally speaking, the time evolution of an observable results in a family of observables given by its composition with the flow.

The Koopman operator theory concentrates on how arbitrary observables change with time, instead of how trajectories are organized in state space. Given a dynamical system

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}), \quad \mathbf{x} \in \mathbb{R}^n, \qquad (7) $$

and a vector space of observables $\mathcal{F}$, the Koopman operator (or more rigorously, semigroup of operators) $U^t$ is defined as the composition of the flow with $g(\mathbf{x}) \in \mathcal{F}$:

$$ U^t g(\mathbf{x}) = g(\phi(\mathbf{x}, t)). \qquad (8) $$

Equation (8) considers the case of a scalar observable $g$. The definition can be extended to a vector-valued function $g(\mathbf{x})$ by applying $U^t$ to each component $g_i(\mathbf{x})$ of the vector, provided that each vector component $g_i(\mathbf{x})$ is in $\mathcal{F}$.

The fundamental property of Koopman operator is linearity. Recalling the example in Figure 2, the evolution of the total energy of the system, $g_4(\mathbf{x}) = g_2(\mathbf{x}) + g_3(\mathbf{x})$, is given by the sum of the evolutions of $g_2$ and $g_3$, i.e. $E(x_0, t) = K(x_0, t) + V(x_0, t)$. More generally, applying $U^t$ to the linear combination of two observables $g_\alpha(\mathbf{x})$ and $g_\beta(\mathbf{x})$ that belong to $\mathcal{F}$ yields

$$ g_{\gamma}(\mathbf{x}) = a g_{\alpha}(\mathbf{x}) + b g_{\beta}(\mathbf{x}) \qquad (9) $$

$$ \begin{aligned} U^t g_{\gamma}(\mathbf{x}) &= g_{\gamma}(\phi(\mathbf{x}, t)) = a g_{\alpha}(\phi(\mathbf{x}, t)) + b g_{\beta}(\phi(\mathbf{x}, t)) = \\ &= a U^t g_{\alpha}(\mathbf{x}) + b U^t g_{\beta}(\mathbf{x}). \end{aligned} \qquad (10) $$

The linearity property of the Koopman operator is general and does not depend on the dynamical system. This property is very attractive, especially in the context of vibration analysis, where eigenvalues and eigenvectors form the basis of vibration modes. The Koopman operator allows such tools to be used with nonlinear systems, but the price to pay is that infinite dimensional vector spaces have to be considered.

### 3.2. Spectral properties

From this point on the space of observables $\mathcal{F}$ is fixed. Although the only requirement is invariance with respect to the action of $U^t$ (given $g \in \mathcal{F}$, $U^t g \in \mathcal{F}$ for every $t$), different choices of $\mathcal{F}$ result in different spectral properties. In this work, $\mathcal{F}$ is taken as the space of $C^1$ functions that are analytic in a neighborhood of the origin. Moreover, only systems with a stable hyperbolic fixed point (assumed to be the origin) are considered, and our investigations are restricted to the basin of attraction of that point.
Based on the above assumptions, this section presents the Koopman operator spectral properties that are essential to this paper. A more comprehensive introduction can be found in Ref. [17].

### 3.2.1. Eigenfunctions

Since $U^t$ is linear, it is possible to consider its eigenfunctions: those observables $s(\mathbf{x}) \in \mathcal{F}$ whose evolution in time is given by

$$U^t s(\mathbf{x}) = s(\phi(t, \mathbf{x})) = \exp(\lambda t) s(\mathbf{x}) \quad \forall t, \qquad (11)$$

in which $\lambda \in \mathbb{C}$ is the corresponding eigenvalue.

As an example consider a simple one-dimensional system

$$\dot{x} = -\eta x, \quad \eta > 0, \qquad (12)$$

with corresponding flow $\phi(x, t) = x \exp(-\eta t)$. Applying the Koopman operator to the observable $s_1(x) = x$ yields

$$U^t s_1(x) = s_1(\phi(x, t)) = x \exp(-\eta t) = \exp(-\eta t) s_1(x), \qquad (13)$$

thus, the identity function $s_1(x) = x$ is an eigenfunction with eigenvalue $-\eta$.

Now consider a general linear system defined by a matrix $\mathbf{A}$. The eigenvalues of $\mathbf{A}$ are also part of the spectrum of the Koopman operator. Moreover, this property extends to nonlinear systems: if $\mathbf{A}$ is the Jacobian matrix evaluated at the origin of the vector field defined in (7), that is

$$\mathbf{A} = \frac{\partial \mathbf{f}}{\partial \mathbf{x}} \Big|_{\mathbf{x}=0}, \qquad (14)$$

then its eigenvalues $\lambda_k$ are part of the spectrum of the Koopman operator [11]. In fact, if $\mathbf{A}$ admits $n$ linearly independent eigenvectors, its eigenvalues are associated with $n$ independent eigenfunctions $s_k(\mathbf{x})$ that are defined in the entire basin of attraction of the equilibrium.

However, the spectrum of $U^t$ is much richer than the spectrum of $\mathbf{A}$: it always contains an infinity of eigenvalues. For instance, in the case of system (12), the observable $s_2(x) = x^2$ is also an eigenfunction:

$$U^t s_2(x) = (\phi(x, t))^2 = (x \exp(-\eta t))^2 = x^2 \exp(-2\eta t) = \exp(-2\eta t) s_2(x). \qquad (15)$$

A similar reasoning can be applied to each integer power $n$ of $x$, so that $-n\eta$ and $s_n(x) = s_1^n(x) = x^n$ are eigenvalues and eigenfunctions of system (12), respectively.

It turns out that this construction can be applied to general systems. Let $s_1(\mathbf{x})$ and $s_2(\mathbf{x})$ be two eigenfunctions:

$$\begin{aligned} U^t s_1(\mathbf{x}) &= \exp(\lambda_1 t) s_1(\mathbf{x}) \\ U^t s_2(\mathbf{x}) &= \exp(\lambda_2 t) s_2(\mathbf{x}). \end{aligned} \qquad (16)$$
Applying the Koopman operator to their product $s_{1,2}(\mathbf{x}) = s_1(\mathbf{x})s_2(\mathbf{x})$ yields

$$ 
\begin{aligned}
U^t s_{1,2}(\mathbf{x}) &= s_{1,2}(\phi(\mathbf{x}, t)) = s_1(\phi(\mathbf{x}, t))s_2(\phi(\mathbf{x}, t)) = \\
&= \exp(\lambda_1 t)s_1(\mathbf{x})\exp(\lambda_2 t)s_2(\mathbf{x}) = \exp((\lambda_1 + \lambda_2)t)s_1(\mathbf{x})s_2(\mathbf{x}).
\end{aligned} 
\quad (17) $$

Thus, the sum of two eigenvalues is also an eigenvalue, and equivalently the product of two eigenfunctions is an eigenfunction. Using (17) it is possible to construct an entire system of eigenfunctions. Considering the basis eigenfunctions $s_k$ associated with the eigenvalues $\lambda_k$ of the Jacobian matrix, we can obtain an infinity of eigenfunctions and eigenvalues:

$$ s_{k_1, \dots, k_n}(\mathbf{x}) = s_1^{k_1}(\mathbf{x}) \dots s_n^{k_n}(\mathbf{x}), \quad \lambda_{k_1, \dots, k_n} = k_1 \lambda_1 + \dots + k_n \lambda_n. \quad (18) $$

In the following a particular property of the eigenfunctions will be central to establish a connection to NNMs defined in Section 2. Let $s(\mathbf{x})$ be an eigenfunction and $\lambda$ the corresponding eigenvalue; then, its zero level set $\{\mathbf{x} \in \mathbb{R}^n | s(\mathbf{x}) = 0\}$ is invariant for the dynamics (7) [11]. Indeed, given an initial condition $\mathbf{x}_0$ satisfying

$$ s(\mathbf{x}_0) = 0, \quad (19) $$

the corresponding trajectory satisfies

$$ s(\phi(t, \mathbf{x}_0)) = \exp(\lambda t)s(\mathbf{x}_0) = 0. \quad (20) $$

If the eigenvalue associated with $s(\mathbf{x})$ is real, then so is $s(\mathbf{x})$, and (19) corresponds to one equality. If $\lambda$ has a nonzero imaginary part, $s(\mathbf{x})$ cannot be purely real, so that (19) is a complex equation yielding two equalities $\mathrm{Re}(s(\mathbf{x})) = 0$ and $\mathrm{Im}(s(\mathbf{x})) = 0$. In the latter case, the complex conjugate $\bar{\lambda}$ is also an eigenvalue with eigenfunction $\bar{s}_k(\mathbf{x})$, so that (19) corresponds to the zero level set of two (complex conjugate) eigenfunctions.

### 3.2.2. Spectral expansion

One of the main tools in the analysis of linear systems is the expansion of the dynamics in terms of linear modes. For a general linear system

$$ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x}, \quad \mathbf{x} \in \mathbb{R}^n, \quad (21) $$

with $\mathbf{A}$ diagonalizable, the flow can be written as

$$ \phi(\mathbf{x}, t) = \sum_{k=1}^{n} s_k(\mathbf{x}) \mathbf{v}_k \exp(\lambda_k t), \quad (22) $$

where $\mathbf{v}_k$ is a (possibly complex) eigenvector of $\mathbf{A}$, $\lambda_k$ is the corresponding eigenvalue, and $s_k(\mathbf{x})$ is the projection of $\mathbf{x}$ onto the eigenspace generated by $\mathbf{v}_k$. Note the use of the same notation for the projections of $\mathbf{x}$ and the Koopman operator eigenfunctions; these two concepts will be linked in Section 4.1.
The Koopman operator theory allows the decomposition (22) to be generalized to nonlinear systems [17]. Let $g$ be an observable that belongs to the span of the eigenfunctions (18). The observable can be expressed as a (possibly infinite) sum of eigenfunctions

$$g(\mathbf{x}) = \sum_{\substack{(k_1, \dots, k_n) \in \mathbb{N}^n \\ k_1 + \dots + k_n > 0}} C_{k_1, \dots, k_n} s_1^{k_1}(\mathbf{x}) \dots s_n^{k_n}(\mathbf{x}). \quad (23)$$

Applying $U^t$ to (23), the time evolution of the observable is given by

$$U^t g(\mathbf{x}) = \sum_{\substack{(k_1, \dots, k_n) \in \mathbb{N}^n \\ k_1 + \dots + k_n > 0}} C_{k_1, \dots, k_n} s_1^{k_1}(\mathbf{x}) \dots s_n^{k_n}(\mathbf{x}) \exp((k_1 \lambda_1 + \dots + k_n \lambda_n)t). \quad (24)$$

The coefficients $C_{k_1, \dots, k_n}$ are the so-called Koopman modes. Each of them is the projection of the observable $g$ onto the eigenspace generated by the corresponding eigenfunction $s_1^{k_1} \dots s_n^{k_n}$. If a vector-valued observable $\mathbf{g}: \mathbb{R}^n \to \mathbb{C}^d$ is considered, the Koopman modes in (23–24) become vector-valued modes $\mathbf{C}_{k_1, \dots, k_n}$ and the scalar eigenfunctions remain unchanged.

As an example, consider for system (12) the observable $g_1(x) = x + x^2$. In the previous section it was shown that the functions $s_1^k(x) = x^k$ are eigenfunctions. Since $g_1(x)$ is a finite sum of eigenfunctions, we directly find the corresponding Koopman modes: $C_1 = 1$, $C_2 = 1$, $C_k = 0$, $k \ge 3$. A less trivial example of spectral decomposition is given by the function $g_2(x) = \exp(x) - 1$. In this case, the Taylor series

$$g_2(x) = \exp(x) - 1 = \sum_{k=1}^{\infty} \frac{1}{k!} x^k = \sum_{k=1}^{\infty} \frac{1}{k!} s_1^k(x), \quad (25)$$

shows that the Koopman modes are $C_k = \frac{1}{k!}$. In both cases (24) can be used to obtain the time evolution of the observables.

It is clear that (24) bears strong resemblance with the linear mode expansion (22). To obtain the same decomposition of the flow as in (22), it is necessary to find an observable whose evolution corresponds to the flow. This is achieved using the identity function $\mathbf{id}(\mathbf{x}) = \mathbf{x}$. In a neighborhood of the origin, $\mathbf{id}(\mathbf{x})$ can be decomposed as sum of eigenfunctions [10]

$$\mathbf{id}(\mathbf{x}) = \mathbf{x} = \sum_{\substack{(k_1, \dots, k_n) \in \mathbb{N}^n \\ k_1 + \dots + k_n > 0}} \mathbf{v}_{k_1, \dots, k_n} s_1^{k_1}(\mathbf{x}) \dots s_n^{k_n}(\mathbf{x}). \quad (26)$$

Using (26) together with (24), a decomposition of the flow in terms of eigenfunctions is obtained

$$U^t \mathbf{id}(\mathbf{x}) = \mathbf{id}(\phi(\mathbf{x}, t)) = \phi(t, \mathbf{x}) = \sum_{\substack{(k_1, \dots, k_n) \in \mathbb{N}^n \\ k_1 + \dots + k_n > 0}} \mathbf{v}_{k_1, \dots, k_n} s_1^{k_1}(\mathbf{x}) \dots s_n^{k_n}(\mathbf{x}) \exp((k_1 \lambda_1 + \dots + k_n \lambda_n)t). \quad (27)$$

The decomposition (27) generalizes (22) to nonlinear systems. This infinite series can be truncated to a finite number of terms to approximate the system trajectories. The Koopman modes $\mathbf{v}_{k_1, \dots, k_n}$ associated with the identity function can be computed by solving a sequence of linear systems; the details of the calculation are reported in appendix A. For linear systems, Eq. (27) coincides with the linear normal mode decomposition, as shown in Section 4.1.
# 4. Koopman Operator and Vibration Modes

## 4.1. The linear case

Consider the linear system

$$ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x}, \quad \mathbf{x} \in \mathbb{R}^n, \qquad (28) $$

with A diagonalizable, and the corresponding expansion in linear modes given by

$$ \phi(\mathbf{x}, t) = \sum_{k=1}^{n} s_k(\mathbf{x})\mathbf{v}_k \exp(\lambda_k t). \qquad (29) $$

The terms $s_k(\mathbf{x})$ are the projections onto the corresponding eigenspace; they can be written in the form

$$ s_k(\mathbf{x}) = \mathbf{w}_k^* \mathbf{x}, \qquad (30) $$

where * denotes the conjugate transpose and $\mathbf{w}_k$ are the eigenvectors of $\mathbf{A}^*$, normalized so that

$$ \mathbf{w}_k^* \mathbf{v}_h = \begin{cases} 1, & k = h \\ 0, & k \neq h \end{cases}. \qquad (31) $$

Applying the Koopman operator to $s_k(\mathbf{x})$, using the previous property, yields

$$ U^t s_k(\mathbf{x}) = s_k(\phi(\mathbf{x}, t)) = \mathbf{w}_k^* \sum_{i=1}^{n} s_i(\mathbf{x}) \mathbf{v}_i \exp(\lambda_i t) = s_k(\mathbf{x}) \exp(\lambda_k t). \qquad (32) $$

Thus, the projections $s_k(\mathbf{x})$ are eigenfunctions of $U^t$ (which explains the notation used for them). It follows that the expansions of the flow through linear modes and Koopman operator eigenfunctions are identical in the linear case. Also, the general decomposition (27) appears as the natural extension of linear normal mode decomposition to nonlinear systems.

The relation between vibration modes and spectral expansion is well-known in the linear case: given a pair of complex conjugate eigenvalues $\lambda_1$ and $\lambda_2$, the motion on the linear mode is obtained when (29) reduces to

$$ s_1(\mathbf{x})\mathbf{v}_1 \exp(\lambda_1 t) + s_2(\mathbf{x})\mathbf{v}_2 \exp(\lambda_2 t), \qquad (33) $$

i.e. when the following conditions on the eigenfunctions are verified:

$$ s_k(\mathbf{x}) = 0, \quad k = 3, \dots, n. \qquad (34) $$

Thus, the plane defining the linear mode is determined by Eq. (34). Given a pair of complex conjugate eigenvalues, the corresponding linear mode is the intersection of the zero level set of all the other eigenfunctions.

In the next two sections (34) is generalized to nonlinear systems, and it is shown that this generalization is equivalent to define NNMs in terms of invariant manifolds.
## 4.2. Nonlinear coordinate transformation

The properties of Koopman operator in the general case of a nonlinear system with hyperbolic fixed point can be studied through linearization techniques [16, 11]. The results of [16] are particularly useful, as they allow a local smooth change of coordinates (or diffeomorphism) $h(x)$ to be extended to the entire basin of attraction of the fixed point. Assuming that the eigenvalues are nonresonant (i.e. relations of the type $\sum_i n_i \lambda_i = 0$ with $n_i \in \mathbb{N}$ and $\sum_i n_i > 0$ are not possible), there exists in a neighborhood of the origin an analytic diffeomorphism of the form [27]

$$ \mathbf{z} = \mathbf{x} + \mathbf{Z}(\mathbf{x}), \quad \mathbf{Z}(\mathbf{0}) = \mathbf{0}, \quad \left. \frac{\partial \mathbf{Z}}{\partial \mathbf{x}} \right|_{\mathbf{x}=\mathbf{0}} = \mathbf{0}, \qquad (35) $$

such that

$$ \dot{\mathbf{z}} = \mathbf{A}\mathbf{z}. \qquad (36) $$

Note that normal form methods aim at approximating the inverse transformation $\mathbf{x} = \mathbf{z} + \mathbf{X}(\mathbf{z})$ through its Taylor series.

According to the results and methods given in [16], there exists a linearizing global $C^1$ diffeomorphism $h(x)$, which is analytic when restricted to a neighborhood of the origin and characterizes the eigenfunctions of the Koopman operator. In particular, the eigenfunction corresponding to the eigenvalue $\lambda_k$ of $\mathbf{A}$ is given by

$$ s_k(\mathbf{x}) = \mathbf{w}_k^* \mathbf{h}(\mathbf{x}), \qquad (37) $$

where $\mathbf{w}_k$ is a left eigenvector of $\mathbf{A}$ (here and in the following $\mathbf{A}$ is assumed diagonalizable). Moreover, it follows from (35) that

$$ \left. \frac{\partial \mathbf{h}}{\partial \mathbf{x}} \right|_{\mathbf{x}=\mathbf{0}} = \mathbf{I}, \qquad (38) $$

so that

$$ \left. \frac{\partial s_k}{\partial \mathbf{x}} \right|_{\mathbf{x}=\mathbf{0}} = \mathbf{w}_k^*, \qquad (39) $$

i.e. the gradient of the eigenfunction at the origin is a left eigenvector of $\mathbf{A}$. Note that the eigenfunctions $s_k(\mathbf{x}) = \mathbf{w}_k^* \mathbf{h}(\mathbf{x})$ are related to the coordinates $\mathbf{z}$. More precisely, they correspond to coordinates $\xi$ given by

$$ \begin{gathered} \xi = \mathbf{W}\mathbf{z}, \\ \mathbf{W} = \begin{pmatrix} \mathbf{w}_1^* \\ \vdots \\ \mathbf{w}_n^* \end{pmatrix}, \end{gathered} \qquad (40) $$

i.e. they are linearizing coordinates, in which the dynamics is diagonal.
### 4.3. The nonlinear case

This section generalizes the connection between the eigenfunctions of Koopman operator and linear normal modes of vibration to NNMs as defined by Shaw and Pierre [7]. Recall from Section 2 that NNM manifolds are characterized by the following properties:

1. the manifold is invariant;

2. it passes through a stable fixed point (the origin);

3. the plane tangent to it at the origin is a linear mode of the linearized dynamics.

Following the discussion of Section 4.1, consider the surface defined by the intersection of $n-2$ zero level sets of eigenfunctions

$$s_k(\mathbf{x}) = 0, \quad \forall k \neq k_1, k_2, \quad \lambda_{k_1} = \overline{\lambda_{k_2}}. \qquad (41)$$

Relying on the properties of the previous sections, we can show that the manifold defined by (41) verifies the properties 1-3:

1. As shown in Section 3.2.1 the level sets $s_k(\mathbf{x}) = 0$ are all invariant, therefore the intersection of any number of them is invariant as well.

2. From (35) and (37) it follows that

$$s_k(\mathbf{0}) = \mathbf{w}_k^*\mathbf{h}(\mathbf{0}) = 0, \qquad (42)$$

thus any zero level set contains the origin, and so does their intersection.

3. The tangent space at the origin of (41) is defined by the linear equations

$$\frac{\partial s_k}{\partial \mathbf{x}}(\mathbf{0})\mathbf{x} = 0, \quad \forall k \neq k_1, k_2 \qquad (43)$$

which, using (39), reduces to

$$\mathbf{w}_k^*\mathbf{x} = 0, \quad \forall k \neq k_1, k_2, \qquad (44)$$

and corresponds to the linear mode associated with $\lambda_{k_1}$ and $\lambda_{k_2}$ as shown in Section 4.1.

Therefore, the generalization of linear modes to nonlinear systems in the Koopman operator framework corresponds to the definition of NNMs proposed by Shaw and Pierre.

### 4.4. NNM approximation using Taylor series

Thanks to the connection established in the previous section, it is possible to use methods from the Koopman operator framework to compute NNMs. In this section a Taylor series approximation based on the spectral expansion (27) is proposed. The underlying idea is to use the values of the eigenfunctions as coordinates. This has the advantage that the dynamics in these coordinates is linear as illustrated in Section 4.2; furthermore, equations (41) have a particularly simple form in this coordinate system.
Let $\xi_k = s_k(\mathbf{x})$ so that (26) becomes

$$ \mathbf{id}(\mathbf{x}) = \mathbf{x} = \sum_{\substack{(k_1, \dots, k_n) \in \mathbb{N}^n \\ k_1 + \dots + k_n > 0}} \mathbf{v}_{k_1, \dots, k_n} \xi_1^{k_1} \dots \xi_n^{k_n}, \qquad (45) $$

and (41)

$$ \xi_k = s_k(\mathbf{x}) = 0, \quad \forall k \neq k_1, k_2, \quad \lambda_{k_1} = \overline{\lambda_{k_2}}. \qquad (46) $$

Suppose $k_1 = 1$ and $k_2 = 2$. Equation (45) expresses the transformation between the coordinate system defined by the Koopman operator eigenfunctions $\xi$ and the original coordinates $\mathbf{x}$, while (46) defines the NNM in the $\xi$ coordinates. Combining the two, we can express the points on the NNM as a function of the variables $\xi_1$ and $\xi_2$:

$$ \mathbf{x} = \sum_{\substack{(k_1, k_2) \in \mathbb{N}^2 \\ k_1 + k_2 > 0}} \mathbf{v}_{k_1, k_2} \xi_1^{k_1} \xi_2^{k_2} := \mathbf{\Psi}(\xi_1, \xi_2), \qquad (47) $$

where the notation $\mathbf{v}_{k_1,k_2} = \mathbf{v}_{k_1,k_2,0,\dots,0}$ is used. Equation (47) assigns to any pair of values of the eigenfunctions $s_1(\mathbf{x})$ and $s_2(\mathbf{x})$ a point on the manifold, and each point on the manifold corresponds to a value $\xi_1 = \bar{\xi}_2 \in \mathbb{C}$. In other words, the manifold is parametrized using all the admissible values for the eigenfunctions. Moreover, it follows from (27) that the time evolution of a trajectory on the manifold is given by

$$ \mathbf{x}(t) = \sum_{\substack{(k_1, k_2) \in \mathbb{N}^2 \\ k_1 + k_2 > 0}} \mathbf{v}_{k_1, k_2} \xi_1^{k_1} \xi_2^{k_2} \exp((k_1 \lambda_1 + k_2 \lambda_2)t) = \mathbf{\Psi}(\xi_1(t), \xi_2(t)). \qquad (48) $$

Equation (47) can be truncated to a fixed order to approximate the invariant manifold in a neighborhood of the origin, and correspondingly (48) to approximate the NNM motion. The vectors $\mathbf{v}_{k_1,k_2}$ are the Koopman modes of the identity function and can be computed as discussed in Appendix A.

To conclude this section some intrinsic limitations of approximating NNMs with the above Taylor series are discussed. First of all, the eigenfunctions might be analytic only in a neighborhood of the origin, so that the amplitude to which NNMs can be computed using (47) might also be limited (see also [11] for a similar issue in approximating eigenfunctions with Taylor series). Even in cases in which the series converges in the whole basin of attraction, the computation of the Koopman modes is more demanding as the order increases, and becomes unfeasible for very high orders.

Due to this the method proposed in this section is not meant to be a ready-to-use algorithm for the computation of NNMs of arbitrary systems, but rather a simple method to compute NNMs for limited amplitudes in small-dimensional systems. Its main purpose is to show that the theory of Koopman operator can lead to computational methods that use the coordinates given by the eigenfunctions. More generic algorithms could be derived using the results of the following section, in which a system of PDEs resembling (3) is derived.
### 4.5. Analogy with Shaw and Pierre's manifold governing PDE

In Section 2 a system of PDEs (3) is derived starting from Shaw and Pierre definition of NNMs. A similar system can be derived for $\Psi(\xi_1, \xi_2)$ defined in (47). To this end, note that the dynamics of $\xi_k$ is by definition

$$ \dot{\xi}_k = \lambda_k \xi_k. \qquad (49) $$

It follows that the time derivative of $\Psi(\xi_1(t), \xi_2(t))$ is given by

$$ \dot{\mathbf{x}}(t) = \frac{\partial \Psi(\xi_1(t), \xi_2(t))}{\partial \xi_1} \lambda_1 \xi_1(t) + \frac{\partial \Psi(\xi_1(t), \xi_2(t))}{\partial \xi_2} \lambda_2 \xi_2(t). \qquad (50) $$

Substituting (50) into the equations of motion (7), and evaluating at $t = 0$ yields

$$ \mathbf{f}(\Psi) = \frac{\partial \Psi}{\partial \xi_1} \lambda_1 \xi_1 + \frac{\partial \Psi}{\partial \xi_2} \lambda_2 \xi_2. \qquad (51) $$

Equation (51) is a set of $n$ PDEs very similar in principle to the one derived in Section 2. The approaches followed here and in Section 2 can be both interpreted as a way to solve

$$ s_k(\mathbf{x}) = 0, \quad i = 3, \dots, N. \qquad (52) $$

When using state-space variables as in (3), Eq. (52) is inverted, but the inversion might be valid only locally. In contrast, a different set of coordinates is used in (51) to parametrize the solution, which avoids issues arising when NNMs have a complex geometry (e.g. folding). The two systems of PDEs (3) and (52) have a similar structure, but (3) corresponds to $n-2$ equations while (52) consists of $n$ equalities. Moreover, the second system admits more solutions, due to the fact that the eigenfunctions are defined modulo a scaling factor. Indeed, if $\Psi(\xi_1, \xi_2)$ is a solution then $\hat{\Psi}(\xi_1, \xi_2) = \Psi(c\xi_1, c\xi_2)$ is a solution for any complex $c \neq 0$.

To obtain $\Psi(\xi_1, \xi_2)$ real-valued $\xi_1 = \bar{\xi}_2$ must be assumed in (51). Equivalently, a system of equations can be derived using as variables $u = \text{Re}(\xi_1)$ and $v = \text{Im}(\xi_1)$, whose dynamics is

$$ \dot{u} = \sigma_1 u - \omega_1 v, \qquad (53) $$

$$ \dot{v} = \omega_1 u + \sigma_1 v, \qquad (54) $$

with $\lambda_1 = \sigma_1 + i\omega_1$. This leads to

$$ \mathbf{f}(\psi) = \nabla \psi \begin{pmatrix} \sigma_1 & -\omega_1 \\ \omega_1 & \sigma_1 \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix}, \qquad (55) $$

with $\psi(u, v) = \Psi(u + iv, u - iv)$. Clearly, the above procedure can be followed for any pair of complex conjugate eigenfunctions, yielding a system of PDEs (51) or (55) for each NNM. These systems differ only for the value of $\lambda$ (or $\sigma$ and $\omega$).
### 4.6. Uniqueness of NNMs

The section is concluded by a remark about conditions 1-3 of Section 4.3 defining the NNMs. Generally these conditions do not determine a unique manifold in the dissipative case. For the sake of illustration, consider a two-dimensional linear system

$$ 
\begin{aligned}
\dot{x}_1 &= \sigma_1 x_1, \\
\dot{x}_2 &= \sigma_2 x_2,
\end{aligned}
\qquad (56) $$

with $\sigma_2 < \sigma_1 < 0$. The phase portrait, illustrated in Figure 3, shows that almost all trajectories approach the linear eigenspace $x_2 = 0$, leading to a continuum of tangent manifolds. Analogous phenomena are possible in higher dimensions. For completeness, Appendix B presents a family of invariant manifolds tangent to a linear normal mode, for a four-dimensional linear system. However, assuming nonresonant eigenvalues, the coefficients of a polynomial approximation are uniquely determined at any order, so that only one of the manifolds is analytic. This observation applies to Shaw and Pierre's approach as well (note that the existence of locally analytic manifolds is guaranteed by the existence of a local analytic change of coordinates, see Section 4.2). The issue of non-uniqueness is discussed in detail in [28].

We note that also in the context of the Koopman operator there is the possibility to obtain multiple manifolds for only one eigenvalue. Specifically, our construction of the eigenfunctions relies on the analyticity of the transformation $\mathbf{h}$ of Section 4.2, which is uniquely defined by its Taylor series. However, if this smoothness condition is dropped it is possible to find different transformations verifying (35) and (36), obtaining in this way multiple eigenfunctions for the same eigenvalue.

<figure><img src="image_4.png" /><figcaption>Figure 3: Phase portrait of a typical two-dimensional linear system</figcaption></figure>
## 5. A 2 DOF example

The Taylor series expansion derived in Section 4.4 and based on the Koopman operator framework is used to compute the NNMs of the system shown in Figure 2. The equations of motion are

$$ 
\begin{aligned}
m_1 \ddot{x}_1 + k_a x_1 + c \dot{x}_1 + k_{nl} x_1^3 + k_b (x_1 - x_2) + c(\dot{x}_1 - \dot{x}_2) &= 0, \\
m_2 \ddot{x}_2 + k_a x_2 + c \dot{x}_2 + k_b (x_2 - x_1) + c(\dot{x}_2 - \dot{x}_1) &= 0,
\end{aligned} 
\quad (57) $$

where $m_1 = m_2 = m = 1$ Kg, $c = 0.05$ Ns/m, $k_a = 1$ N/m and $k_{nl} = 0.5$ N/m³. With $y_1 = \dot{x}_1$, $y_2 = \dot{x}_2$, $\eta_1 = (k_a + k_b)/m$, $\eta_2 = k_b/m$, $\beta = c/m$, $\alpha = k_{nl}/m$, (57) is recast in first-order form

$$ 
\begin{pmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{y}_1 \\ \dot{y}_2 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\eta_1 & \eta_2 & -2\beta & \beta \\ \eta_2 & -\eta_1 & \beta & -2\beta \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ y_1 \\ y_2 \end{pmatrix} + \begin{pmatrix} 0 \\ 0 \\ -\alpha x_1^3 \\ 0 \end{pmatrix}. \quad (58) $$

The Jacobian matrix

$$ 
\mathbf{A} = \begin{pmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -\eta_1 & \eta_2 & -2\beta & \beta \\ \eta_2 & -\eta_1 & \beta & -2\beta \end{pmatrix} \quad (59) $$

has two pairs of complex conjugate eigenvalues $\lambda_1 = \bar{\lambda}_2$, $\lambda_3 = \bar{\lambda}_4$, which verify a 3 : 1 resonance condition $\lambda_1 = 3\lambda_3$ for $k_b = 4$. Exact linear resonance is eliminated by varying $k_b$ but folding is likely to occur in the neighborhood of those parameter values. Table 1 shows how the eigenvalues vary for different values of $k_b$.

<table><caption>Table 1: Eigenvalues, natural frequencies and damping ratios for different values of $k_b$.</caption><thead><tr><th>k<sub>b</sub></th><th>λ<sub>1</sub></th><th>f<sub>n,1</sub></th><th>ζ<sub>1</sub></th><th>λ<sub>2</sub></th><th>f<sub>n,2</sub></th><th>ζ<sub>2</sub></th></tr></thead><tbody><tr><td>4.7</td><td>-0.075+i3.224</td><td>3.224</td><td>0.0233</td><td>-0.025+i0.999</td><td>1</td><td>0.025</td></tr><tr><td>4.3</td><td>-0.075+i3.079</td><td>3.098</td><td>0.0242</td><td>-0.025+i0.999</td><td>1</td><td>0.025</td></tr><tr><td>4.1</td><td>-0.075+i3.032</td><td>3.033</td><td>0.0247</td><td>-0.025+i0.999</td><td>1</td><td>0.025</td></tr></tbody></table>

To determine the coefficients describing the NNM in Eq. (47), the trajectory on the manifold given by (48) is injected into the equations of motion. Evaluating at time $t = 0$ yields

$$ \sum_{\substack{(k_1, k_2) \in \mathbb{N}^2 \\ k_1 + k_2 > 0}} \mathbf{v}_{k_1, k_2} \xi_1^{k_1} \xi_2^{k_2} (k_1 \lambda_1 + k_2 \lambda_2) = $$
$$A \sum_{\substack{(k_1,k_2) \in \mathbb{N}^2 \\ k_1+k_2>0}} \mathbf{v}_{k_1,k_2} \xi_1^{k_1} \xi_2^{k_2} + \left( -\alpha \left( \left( \sum_{\substack{(k_1,k_2) \in \mathbb{N}^2 \\ k_1+k_2>0}}} \mathbf{v}_{k_1,k_2} \xi_1^{k_1} \xi_2^{k_2} \right)_1 \right)^3 \right), \quad (60)$$

where $(\mathbf{v})_i$ stands for the ith component of the vector $\mathbf{v}$. Matching the terms corresponding to the first-order powers results in two systems

$$(\lambda_1 \mathbf{I} - \mathbf{A}) \mathbf{v}_{10} = 0, \quad (\lambda_2 \mathbf{I} - \mathbf{A}) \mathbf{v}_{01} = 0, \quad (61)$$

so that the first-order modes $\mathbf{v}_{10}$ and $\mathbf{v}_{01}$ are the eigenvectors of $A$ corresponding to $\lambda_1$ and $\lambda_2$. If there is no resonance between the eigenvalues of $A$, the successive systems admit a unique solution, and are obtained by matching like-power terms (see the appendix for the complete derivation). Doing so for the second-order terms yields linear equations of the form

$$((k_1\lambda_1 + k_2\lambda_2)\mathbf{I} - \mathbf{A})\mathbf{v}_{k_1k_2} = 0, \quad k_1 + k_2 = 2, \quad (62)$$

where the right hand side is zero due to the absence of second-order terms in the vector field. We can proceed iteratively, solving for any mode $\mathbf{v}_{k_1,k_2}$ a linear system of the form

$$(r\mathbf{I} - \mathbf{A})\mathbf{v}_{k_1k_2} = \mathbf{b}, \quad r = k_1\lambda_1 + k_2\lambda_2 \quad (63)$$

where $\mathbf{b}$ is obtained as a sum of lower-order modes (see appendix A for the details). A change of coordinates using the Jordan basis is used to improve numerical stability.

The above procedure is used to compute the high-frequency (i.e. the out-of-phase) mode of the 2 DOF system for $k_b = 4.3$. The Taylor series used to generate the figures of this section were truncated at the 50th order. The maximum amplitude at which the manifolds are represented was chosen so that trajectories obtained with (48) agree with those computed through numerical integration. To have a quantitative measure of the difference between two trajectories the normalized mean square error (NMSE) was used:

$$\text{NMSE} = \frac{100}{\sigma_x^2} \sum_{k=0}^{K} \|\hat{\mathbf{x}}(k) - \mathbf{x}(k)\|^2, \quad (64)$$

where $\hat{\mathbf{x}}$ is the time series obtained with (48), $\mathbf{x}$ is the reference series obtained through numerical integration, and $\sigma_x^2$ is the variance of the latter. The NMSEs obtained to validate the manifolds were always lower than 1%, a value usually taken as threshold for a good match between two time series.

The result is compared with the finite element (FE) method developed in [19]. This numerical method is based on an explicit parametrization of the invariant surface using a pair of master coordinates and it is known to be accurate provided that a good mesh is employed. For the out-of-phase mode, the selected master coordinates correspond to the
<figure><img src="image_5.png" /><figcaption>Figure 4: Out-of-phase NNM; ■ FE method, □ Taylor series; (a) $x_1, y_1, x_2$; (b) $x_1, y_1, y_2$</figcaption></figure>

displacement and velocity of the first mass. As can be seen in Figure 4 the results from the two methods agree very well, but the FE method allows higher levels of amplitude to be reached, while the Taylor series validity is more limited.

The lower-frequency (in-phase) NNM was computed for the different values of $k_b$ reported in Table 1. The results are shown in Figure 5 and compared with the NNMs obtained using the FE algorithm. As $k_b$ is decreased from 4.7 (Figures 5a and 5b) to 4.3 (Figures 5c and 5d) and 4.1 (Figures 5e and 5f), the manifold geometry changes and starts to fold, which prevents the use of master coordinates for the parametrization. This is also confirmed by Figure 6, where the manifold obtained for $k_b = 4.1$ is shown with a viewpoint orthogonal to the $x_2y_2$ plane. The FE algorithm reaches the folding but cannot go beyond, since it cannot converge in region where the master coordinates parametrization is not valid. In contrast the folding does not pose any problem with the parametrization obtained through the Koopman eigenfunctions.

In Figure 7 one of the trajectories considered for the validation is compared with the corresponding numerical integration, in the case $k_b = 4.1$. As expected according to the values of the NMSE, the two trajectories cannot be distinguished by visual inspection.

Equation (48) also highlights a structure of the evolution that is reminiscent of that obtained using Fourier series with periodic trajectories. The main difference is that the number of terms in (48) is not the same at all orders, but increases with the order. Figure 8 shows the terms corresponding to the first three non-zero orders, for $k_b = 4.7$ and $k_b = 4.1$, and in Figure 9 their sum is compared with the trajectories obtained using numerical integration. When only terms up to order five are used the agreement between the two trajectories is not as good as in Figure 7, nonetheless, it is remarkable that the main aspects of the dynamics are retained by this approximation, at least in the case considered herein. It is also interesting to note that in the case $k_b = 4.7$ the trajectory presents a lower harmonic content, and using the same order a better approximation is obtained.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 5: In-phase NNM for different values of $k_b$; (a)-(b) $k_b = 4.7$, (c)-(d) $k_b = 4.5$, (e)-(f) $k_b = 4.1$; ■ FE method, ● Taylor series</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 6: Different views of the invariant manifold shown in Figure 5 ($k_b = 4.1$); (a) $x_1$ (Figure 5e) top view, (b) $x_2$ (Figure 5f) bottom view. Note that the manifold undergoes folding, which prevents a parametrization with state-space variables.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7: Comparison of time series;— Numerical integration,--- Taylor approximation; initial condition: $x_1(0) = -0.23$, $x_2(0) = -0.76$, $y_1(0) = -0.4$, $y_2(0) = -0.86$; (a) $x_1$, (b) $x_2$, (c) $y_1$, (d) $y_2$</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8: Decomposition of trajectories in different terms; —1st order term, ---3rd order term, ···5th order term; (a) $k_b = 4.7$ ($x_1(0) = -0.57$, $x_2(0) = -0.55$, $y_1(0) = -0.44$, $y_2(0) = -0.66$); (b) $k_b = 4.1$ ($x_1(0) = -0.23$, $x_2(0) = -0.76$, $y_1(0) = -0.4$, $y_2(0) = -0.86$)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9: Comparison of 5th order approximation and numerical integration; — Numerical integration, ··· 5th order approximation; (a) $k_b = 4.7$ ($x_1(0) = -0.57$, $x_2(0) = -0.55$, $y_1(0) = -0.44$, $y_2(0) = -0.66$), (b) $k_b = 4.1$ ($x_1(0) = -0.23$, $x_2(0) = -0.76$, $y_1(0) = -0.4$, $y_2(0) = -0.86$)</figcaption></figure>
## 6. Conclusion

This research investigates the intimate connection that exists between NNMs of damped systems and the Koopman operator associated with their dynamics. In particular, we have shown that there is a correspondence between the invariant manifolds defined by Shaw and Pierre and the zero level sets of particular eigenfunctions of the Koopman operator, and that this correspondence is a natural generalization of the linear case.

Through the Koopman operator eigenfunctions a new characterization of NNMs is possible, and a novel approach to their parametrization is proposed. The main advantage of this new approach is its global validity, obtained as a consequence of the linearizing properties provided by the eigenfunctions. This parametrization overcomes a known limitation in the use of master coordinates.

However, the Taylor series expansion used herein is local in nature and restricted to low-dimensional systems. In order to target arbitrarily high levels of amplitude, as well as high-dimensional systems, more efficient methods should be developed for solving the system of PDEs (51) or (55). As mentioned in Section 4.5 there are differences between these systems and the one obtained in the work of Shaw and Pierre. Nonetheless, all these systems share a common structure, which might allow the extension of numerical methods developed using master coordinates to this new framework. In addition, the Koopman operator approach is associated with recent data analysis techniques that are well-suited to high-dimensional systems [29]. This could potentially open the door to develop data-driven methods for computing NNMs of large systems.

The preliminary results presented here motivate a deeper investigation of NNMs for damped systems. For example it is well known that in conservative systems NNMs can bifurcate, and exceed the number of degrees of freedom of the system, but how these phenomena survive the presence of damping is still an open question.

## Acknowledgment

This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office.

This work was performed while A. Mauroy held a return grant from the Belgian Science Policy (BELSPO). The author L. Renson is a Marie-Curie COFUND Postdoctoral Fellow of the University of Liège, co-funded by the European Union.

## A. Computation of the Koopman Modes of the Identity

In this appendix we derive a method for computing the Koopman modes $v_{k_1, \dots, k_n}$ of the identity (see (26)). In particular, we show that these modes are the solutions of linear systems.

We fix some notation. For $\alpha = (\alpha_1, \dots, \alpha_n) \in \mathbb{N}^n$ let

$$ \alpha! = \prod_{j=1}^{n} \alpha_j!, \qquad (65) $$
$$|\alpha| = \sum_{j=1}^{n} \alpha_j, \qquad (66)$$

$$\partial^{\alpha} \mathbf{f}_l = \frac{\partial^{|\alpha|} \mathbf{f}_l}{\partial x_1^{\alpha_1} \dots \partial x_n^{\alpha_n}}, \qquad (67)$$

$$x^{\alpha} = \prod_{j=1}^{n} x_j^{\alpha_j}. \qquad (68)$$

Using this notation, we can write

$$s_{k_1, \dots, k_n}(\mathbf{x}) = s^k(\mathbf{x}), \quad k = (k_1, \dots, k_n), \qquad (69)$$

with

$$\mathbf{s}(\mathbf{x}) = \begin{pmatrix} s_1(\mathbf{x}) \\ \vdots \\ s_n(\mathbf{x}) \end{pmatrix}. \qquad (70)$$

In addition, (26) can be rewritten as

$$\mathbf{x} = \sum_{h \ge 1} \sum_{\substack{k \in \mathbb{N}^n \\ |k|=h}} \mathbf{v}_k s^k(\mathbf{x}), \qquad (71)$$

and

$$\dot{\phi}(t, \mathbf{x})|_{t=0} = \mathbf{f}(\mathbf{x}) = \sum_{h \ge 1} \sum_{\substack{k \in \mathbb{N}^n \\ |k|=h}} \mathbf{v}_k s^k(\mathbf{x}) k \lambda, \qquad (72)$$

where $k\lambda = \sum_{j=1}^n k_j \lambda_j$ and $\lambda_j$ are the eigenvalues of the Jacobian matrix

$$\mathbf{A} = \left. \frac{\partial \mathbf{f}}{\partial \mathbf{x}} \right|_{\mathbf{x}=\mathbf{0}}. \qquad (73)$$

Similarly, the Taylor series of $\mathbf{f}(\mathbf{x})$ can be written as

$$\mathbf{f}_l(\mathbf{x}) = \sum_{h \ge 1} \sum_{\substack{k \in \mathbb{N}^n \\ |k|=h}} \frac{1}{k!} \partial^k \mathbf{f}_l(\mathbf{0}) \mathbf{x}^k, \qquad (74)$$

and injecting (71) in (74) leads to

$$\mathbf{f}_l(\mathbf{x}) = \sum_{h_1 \ge 1} \sum_{\substack{k_1 \in \mathbb{N}^n \\ |k_1|=h_1}} \frac{1}{k_1!} \partial^{k_1} \mathbf{f}_l(\mathbf{0}) \prod_{j=1}^{n} \left( \sum_{h_2 > 1} \sum_{\substack{k_2 \in \mathbb{N}^n \\ |k_2|=h_2}} (\mathbf{v}_{k_2})_j s^{k_2}(\mathbf{x}) \right)^{k_1 j} . \qquad (75)$$

By comparison of (72) and (75), for any $k \in \mathbb{N}^n$ we obtain

$$(\mathbf{v}_k)_l k \lambda = \sum_{h=1}^{|k|} \sum_{|j|=h} \frac{1}{j!} \partial^j \mathbf{f}_l(\mathbf{0}) \sum_{\substack{(k_1 \dots k_j)| \\ \sum_i k_i = k}} \prod_{i=1}^n \prod_{j=1}^{s} (\mathbf{v}_{k_m})_i. \qquad (76)$$
These equations are linear in the terms with maximum value of $|k|$, and can be solved if all the $v_k'$ with $|k'| < |k|$ are known.

The matrix defining the linear system related to $v_k$ has the form $k\lambda I - A$. The first-order modes are the eigenvectors of $A$ and it is necessary to choose a scaling factor, which correspondingly fixes the value of the eigenfunctions (defined up to a multiplicative constant). If there are no resonances between the eigenvalues, the linear systems related to modes of higher order are uniquely determined, since they are defined by a nonsingular matrix. We can also use the known form of the matrices defining the linear systems to simplify the computation. In the case of a diagonalizable Jacobian matrix, the systems are diagonal when expressed in coordinates defined by a Jordan basis.

To obtain a manifold using (47) not all the modes are necessary, but just a subset of them: the $v_k$, $k = (k_1, \dots, k_n)$, for which only two indices $k_{j1}$ and $k_{j2}$ are nonzero. Since to obtain the linear system corresponding to a given mode $v_k$ only the modes $v_{k'}$, $k' = (k'_1, \dots, k'_n)$, with $k_l \ge k'_l$ are needed, the computation can be restricted to the modes that appear in (47).

As concluding remarks we note that approximating a solution of (51) with a Taylor series yields the same linear systems. Also, the procedure presented here corresponds in essence to a normal form approximation. In that context, starting from a dynamical system with diagonal linear part, a transformation can be obtained so that the resulting equations of motion are the simplest possible. The only algebraic condition needed to obtain a transformation that linearizes the dynamic is the nonresonant eigenvalues. The corresponding linearizing coordinates are the eigenfunctions of the Koopman operator.

## B. A family of invariant manifolds for a linear system

Consider a four-dimensional linear system

$$ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x}, \qquad (77) $$

with

$$ \mathbf{A} = \begin{pmatrix} -\sigma_1 & \omega_1 & 0 & 0 \\ -\omega_1 & -\sigma_1 & 0 & 0 \\ 0 & 0 & -\sigma_2 & \omega_2 \\ 0 & 0 & -\omega_2 & -\sigma_2 \end{pmatrix}. \qquad (78) $$

Note that any four-dimensional linear system with two couples of complex conjugate eigenvalues can be transformed into this form. It will be assumed that $\sigma_1 > 0$ and $\sigma_2 > 0$ (stability).

For any value of $c_1$ and $c_2$ consider the manifold $\mathbf{x}(u, v) = \mathbf{f}(u, v; c_1, c_2)$,

$$ f_1(u, v) = u $$

$$ f_2(u, v) = v $$

$$ f_3(u, v) = (u^2 + v^2)^{\frac{\sigma_2}{2\sigma_1}} \left( \cos \left( \frac{\log(u^2 + v^2)\omega_2}{2\sigma_1} \right) c_1 - \sin \left( \frac{\log(u^2 + v^2)\omega_2}{2\sigma_1} \right) c_2 \right) \qquad (79) $$

$$ f_4(u, v) = (u^2 + v^2)^{\frac{\sigma_2}{2\sigma_1}} \left( \sin \left( \frac{\log(u^2 + v^2)\omega_2}{2\sigma_1} \right) c_1 + \cos \left( \frac{\log(u^2 + v^2)\omega_2}{2\sigma_1} \right) c_2 \right). $$
The derivatives of $\mathbf{f}$ are

$$f_{u,1} = 1$$

$$f_{u,2} = 0$$

$$f_{u,3} = \frac{u(u^2 + v^2)^{\frac{\sigma_2}{2\sigma_1} - 1} \left( \cos\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_1\sigma_2 - c_2\omega_2) - \sin\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_2\sigma_2 + c_1\omega_2) \right)}{\sigma_1}$$

$$f_{u,4} = \frac{u(u^2 + v^2)^{\frac{\sigma_2}{2\sigma_1} - 1} \left( \cos\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_2\sigma_2 + c_1\omega_2) + \sin\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_1\sigma_2 - c_2\omega_2) \right)}{\sigma_1}$$

$$f_{v,1} = 0$$

$$f_{v,2} = 1$$

$$f_{v,3} = \frac{v(u^2 + v^2)^{\frac{\sigma_2}{2\sigma_1} - 1} \left( \cos\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_1\sigma_2 - c_2\omega_2) - \sin\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_2\sigma_2 + c_1\omega_2) \right)}{\sigma_1}$$

$$f_{v,4} = \frac{v(u^2 + v^2)^{\frac{\sigma_2}{2\sigma_1} - 1} \left( \cos\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_2\sigma_2 + c_1\omega_2) + \sin\left(\frac{\log(u^2+v^2)\omega_2}{2\sigma_1}\right) (c_1\sigma_2 - c_2\omega_2) \right)}{\sigma_1} \quad (80)$$

These are continuous for $\sigma_2 > \sigma_1$, and more regular if the distance between the two real parts of the eigenvalues is increased further.

Assuming sufficient regularity, $\mathbf{f}$ defines a manifold, with tangent plane at a point $(u, v)$ spanned by the vectors

$$\mathbf{f}_u = \frac{\partial \mathbf{f}}{\partial u}, \quad \mathbf{f}_v = \frac{\partial \mathbf{f}}{\partial v}. \qquad (81)$$

Since it holds

$$\mathbf{A}\mathbf{f}(u, v) = (-\sigma_1 u + \omega_1 v)\mathbf{f}_u + (-\omega_1 u - \sigma_1 v)\mathbf{f}_v, \qquad (82)$$

all the manifolds are invariant ((82)) is the analogous of the equation derived by Shaw and Pierre equation for the system (77)). Moreover, if $\sigma_2 > \sigma_1$

$$\mathbf{f}_u(0, 0) = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \quad \mathbf{f}_v(0, 0) = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \qquad (83)$$

so that at the origin the tangent plane is $x_3 = 0, x_4 = 0$, which corresponds to the linear normal mode, but only for $c_1 = c_2 = 0$ the whole manifold is the linear normal mode.

Note that apart from the points verifying $x_{1,0} = x_{2,0} = 0$ (i.e. points on the fast linear mode), for any other point $(x_{1,0}, x_{2,0}, x_{3,0}, x_{4,0})$, there exists a manifold which passes through
it, this is obtained using the parameters

$$c_1 = x_{3,0} \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)^{\frac{\sigma_2}{2\sigma_1}} \cos \left( \frac{\omega_2 \log \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)}{2\sigma_1} \right) - x_{4,0} \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)^{\frac{\sigma_2}{2\sigma_1}} \sin \left( \frac{\omega_2 \log \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)}{2\sigma_1} \right), \\ c_2 = x_{3,0} \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)^{\frac{\sigma_2}{2\sigma_1}} \sin \left( \frac{\omega_2 \log \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)}{2\sigma_1} \right) + x_{4,0} \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)^{\frac{\sigma_2}{2\sigma_1}} \cos \left( \frac{\omega_2 \log \left( \frac{1}{x_{1,0}^2 + x_{2,0}^2} \right)}{2\sigma_1} \right). \tag{84}$$

As a final remark we note that even if each couple of coefficients $c_1, c_2$ defines a different manifold, the derivatives of all the functions defined by (79) have the same values at the origin (up to the order at which they are defined): they all share the same polynomial approximation. Interestingly the same situation is encountered in the context of center manifolds.

## References

[1] C. Touzé, M. Amabili, Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures, Journal of Sound and Vibration 298 (2006) 958–981.

[2] A. F. Vakakis, L. I. Manevitch, Y. V. Mikhlin, V. N. Pilipchuk, A. A. Zevin, Normal modes and localization in nonlinear systems, Springer, 1996.

[3] M. Peeters, G. Kerschen, J. Golinval, Modal testing of nonlinear vibrating structures based on nonlinear normal modes: Experimental demonstration, Mechanical Systems and Signal Processing 25 (2011) 1227–1247.

[4] L. Renson, A. Gonzalez-Buelga, D. A. W. Barton, S. A. Neild, Robust identification of backbone curves using control-based continuation, Journal of Sound and Vibration 367 (2016) 145–158.

[5] R. M. Rosenberg, Normal modes of nonlinear dual-mode systems, Journal of Applied Mechanics 27 (1960) 263–268.

[6] Y. S. Lee, G. Kerschen, A. F. Vakakis, P. Panagopoulos, L. Bergman, D. M. McFarland, Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment, Physica D: Nonlinear Phenomena 204 (2005) 41–69.

[7] S. W. Shaw, C. Pierre, Normal modes for non-linear vibratory systems, Journal of Sound and Vibration 164 (1993) 85–124.

[8] B. O. Koopman, Hamiltonian systems and transformation in hilbert space, Proceedings of the National Academy of Sciences of the United States of America 17 (1931) 315.

[9] I. Mezić, Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynamics 41 (2005) 309–325.

[10] A. Mauroy, I. Mezić, J. Moehlis, Isostables, isochrons, and Koopman spectrum for the action–angle representation of stable fixed point dynamics, Physica D: Nonlinear Phenomena 261 (2013) 19–30.

[11] A. Mauroy, I. Mezić, Global stability analysis using the eigenfunctions of the Koopman operator, Automatic Control, IEEE Transactions on PP (2016).

[12] J. H. Tu, D. M. Rowley, Clarence W. andLuchtenburg, S. L. Brunton, J. N. Kutz, On dynamic mode decomposition: Theory and applications, Journal of Computational Dynamics 1 (2014) 391–421.

[13] C. W. Rowley, I. Mezić, S. Bagheri, P. Schlatter, D. S. Henningson, Spectral analysis of nonlinear flows, Journal of Fluid Mechanics 641 (2009) 115–127.
[14] Y. Susuki, I. Mezić, Nonlinear Koopman modes and power system stability assessment without models, IEEE Transactions on Power Systems 29 (2014) 899-907.

[15] B. Eisenhower, T. Maile, M. Fischer, I. Mezić, Decomposing building system data for model validation and analysis using the Koopman operator, in: Proceedings of the National IBPSAUSA Conference, New York, USA, 2010.

[16] Y. Lan, I. Mezić, Linearization in the large of nonlinear systems and Koopman operator spectrum, Physica D: Nonlinear Phenomena 242 (2013) 42-53.

[17] I. Mezić, Analysis of fluid flows via spectral properties of the Koopman operator, Annual Review of Fluid Mechanics 45 (2013) 357-378.

[18] D. Jiang, C. Pierre, S. Shaw, The construction of non-linear normal modes for systems with internal resonance, International Journal of Non-Linear Mechanics 40 (2005) 729 - 746.

[19] L. Renson, G. Deliége, G. Kerschen, An effective finite-element-based method for the computation of nonlinear normal modes of nonconservative systems, Meccanica 49 (2014) 1901-1916.

[20] J. Carr, Applications of Centre Manifold Theory, volume 35 of *Applied Mathematical Sciences*, Springer-Verlag, 1981.

[21] E. Pescheck, C. Pierre, S. W. Shaw, A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds, Journal of Sound and Vibration 249 (2002) 971-993.

[22] F. Blanc, C. Touzé, J. F. Mercier, K. Ege, A. S. Bonnet Ben-Dhia, On the numerical computation of nonlinear normal modes for reduced-order modelling of conservative vibratory systems, Mechanical Systems and Signal Processing 36 (2013) 520-539.

[23] L. Renson, G. Kerschen, B. Cochelin, Numerical computation of nonlinear normal modes in mechanical engineering, Journal of Sound and Vibration 364 (2016) 177-206.

[24] N. Boivin, C. Pierre, S. W. Shaw, Non-linear modal analysis of structural systems featuring internal resonances, Journal of Sound and Vibration 182 (1995) 6.

[25] D. Jiang, C. Pierre, S. W. Shaw, The construction of non-linear normal modes for systems with internal resonance, International Journal of Non-Linear Mechanics 40 (2005) 729-746.

[26] S. A. Neild, A. Cammarano, D. J. Wagg, Nonlinear modal decomposition using normal form transformations, in: Topics in Nonlinear Dynamics, Volume 1, Springer, 2013, pp. 179-187.

[27] V. I. Arnol'd, Geometrical methods in the theory of ordinary differential equations, volume 250, Springer New York, 1988.

[28] G. Haller, Nonlinear normal modes and spectral submanifolds: Existence, uniqueness and use in model reduction, Submitted to Nonlinear Dynamics, arXiv preprint arXiv:1602.00560 (2016).

[29] M. O. Williams, C. W. Rowley, I. G. Kevrekidis, A kernel approach to data-driven Koopman spectral analysis, arXiv preprint arXiv:1411.2260 (2015).