Nonlinear Dyn (2023) 111:8439–8466
https://doi.org/10.1007/s11071-023-08247-7

ORIGINAL PAPER

Sorting-free Hill-based stability analysis of periodic
solutions through Koopman analysis

Fabia Bayer · Remco I. Leine

Received: 19 October 2022 / Accepted: 20 December 2022 / Published online: 6 February 2023
© The Author(s) 2023

Abstract In this paper, we aim to study nonlinear
time-periodic systems using the Koopman operator,
which provides a way to approximate the dynamics
of a nonlinear system by a linear time-invariant system
of higher order. We propose for the considered sys-
tem class a specific choice of Koopman basis functions
combining the Taylor and Fourier bases. This basis
allows to recover all equations necessary to perform
the harmonic balance method as well as the Hill anal-
ysis directly from the linear lifted dynamics. The key
idea of this paper is using this lifted dynamics to formu-
late a new method to obtain stability information from
the Hill matrix. The error-prone and computationally
intense task known by sorting, which means identify-
ing the best subset of approximate Floquet exponents
from all available candidates, is circumvented in the
proposedmethod. TheMathieu equation and an n-DOF
generalization are used to exemplify these findings.

Keywords Time-periodic systems · Harmonic
balance method · Koopman lift · Floquet multipliers ·
Monodromy matrix

F. Bayer (B) · R. I. Leine
Institute for Nonlinear Mechanics, University of Stuttgart,
Pfaffenwaldring 9, 70569 Stuttgart, Germany
e-mail: bayer@inm.uni-stuttgart.de

R. I. Leine
e-mail: leine@inm.uni-stuttgart.de

1 Introduction

The objective of this paper is to introduce a novel stabil-
ity method based on the Hill matrix, which differs from
the state-of-the-art methods in that a matrix projection
is applied before computing an eigenvalue problem.
The structure of this projection is obtained by consid-
ering the Hill matrix to be a result from the Koopman
lift for well-chosen basis functions.

The Koopman framework [1,2] has gained immense
popularity in recent years as a versatile tool for
various engineering applications, such as system iden-
tification [3], model order reduction [4] and feedback
control [5]. This is due to an auspicious promise: global
representation of a nonlinear system by a linear oper-
ator. To this end, in the Koopman framework, the
dynamical system is defined through the propagation
of functions on the state space, also called observables,
over time. Bernard Koopman first described a unitary
linear operator which evolves a class of measurable
functions along the flow of a conservative system [6],
which would later be named the Koopman operator.
After some relatively quiet years, the Koopman oper-
ator experienced a revival around the turn of the cen-
tury, when it was shown that its spectral characteristics
contain global properties for the underlying dynami-
cal system [7]. This sparked generalizations to non-
conservative systems [8,9] and at the same time, numer-
ical and data-driven methods like the Arnoldi method
[10] and extended dynamic mode decomposition [2]
emerged.

123

http://crossmark.crossref.org/dialog/?doi=10.1007/s11071-023-08247-7&domain=pdf
http://orcid.org/0000-0002-2038-3374
http://orcid.org/0000-0001-9859-7519


8440 F. Bayer, R. I. Leine

Classically, the Koopman framework is applied to
time-autonomous systems ẋ = f(x) and the approx-
imate linear dynamics obtained by the Koopman lift
then takes the form ż = Az. The incorporation of
a time-dependent input v(t) into the dynamics, i.e.,
ẋ = f(x, v(t)) or simply ẋ = f(x, t), generally poses
problems in the Koopman framework as the system
can only be approximated by a linear time-invariant
(LTI) system ż = Az + Bu(t) if products of state
and input are neglected. In this paper we focus on
non-autonomous systems for which the input is time-
periodic, i.e., v(t) = v(t + T ). In particular, we pro-
pose a specific choice of observable functions which
contains observables depending both on state and time,
opening the possibility to include products of state and
input.

The numerical computation of periodic solutions
in time-periodic non-autonomous systems is a task of
greatest interest in engineering application. Periodic
solutions are of prime importance for, e.g., nonlin-
ear vibration analysis in structural dynamics [11,12],
in acoustics [13] and thermo-acoustics [14], nonlinear
oscillation analysis in electronic circuits [15] as well as
heart rhythm analysis in cardiology [16]. Therefore, it
is important to find and characterize periodic solutions,
assess their stability properties and also track these
quantities along varying system parameters, a process
which is called continuation [17].

Naively, attractive periodic solutions can be found
by simply simulating numerically over a long time
interval, until transient effects are negligible. However,
since this method is very dependent on the initial con-
dition, can only find attractive solutions and is com-
putationally expensive, more sophisticated methods
have been developed. There are a multitude of periodic
solution solvers available, including finite differences
[18], shooting [19], multiple shooting [20], collocation
and generalized collocation [21] and harmonic
balancing [22]. We highlight two methods:

• The shooting method [19] uses the Newton method
and a numerical ODE solver to find initial condi-
tions x0 such that x(t0 + T ) = x(t0) = x0, i.e.,
whose trajectory over a period fulfills the periodic-
ity constraint. Since themonodromymatrix appears
in the update step of the Newton method, stabil-
ity information comes (almost) for free with this
method. This monodromy matrix also plays a role

in the tangent prediction in an arc-length continua-
tion method.

• The harmonic balance method (HBM) [22,23], in
contrast to the othermentionedmethods, is a purely
frequency-based method. The periodic solution is
parameterized globally by afinite set of trigonomet-
ric functions, and the coefficients are determined
using a residual in the frequency domain. The trans-
fer of nonlinear terms into the frequency domain
is non-trivial. In practice, this is achieved either
using an alternating frequency and time (AFT)
method [24] or by providing a recast of the sys-
tem in quadratic form [25]. One advantage of the
HBM is that it automatically provides a filtering
effect on the identified periodic solution.

While stability information about the periodic solu-
tion branches comes almost for free in the shooting
methods as the monodromy matrix is calculated as a
necessary continuation step, stability information and
hence information about the bifurcations that occur in
a system are not so trivial in the HBM. The Hill matrix
provides a frequency-based method to assess the sta-
bility of periodic solutions [26,27]. It can be found and
constructed easily in the numerical asymptotic method
(ANM) [25,28], being a continuation method based on
the HBM. Hence, computation of stability through its
eigenvalues,which approximate theFloquet exponents,
is a viable option.

However, two critical problems make the Hill
method often unattractive in practice. On the one hand,
fromanumerical viewpoint, computing the eigenvalues
of the large Hill matrix is computationally expensive
and potentially inaccurate [29]. On the other hand, for
correct assertion of stability, only a non-trivial subset
of these eigenvalues must be considered. This process
is known in the literature as sorting of Floquet expo-
nent candidates. The determination of this subset is still
an active area of research, with the approaches being
based on the imaginary parts [23,26,30] and poten-
tially in addition the real parts [31] of the eigenvalues,
or alternatively symmetry considerations of the eigen-
vectors [27,28].

In this work, we propose a different approach for
obtaining stability information from the Hill matrix,
which circumvents both issuesmentioned above.Using
the Koopman framework we motivate a novel dynam-
ical systems interpretation of the Hill matrix, which
allows to compute an approximation of themonodromy

123



Sorting-free Hill-based stability analysis 8441

matrix directly (i.e., without computing a large number
of eigenvalues and subsequent sorting). The proposed
method to find the monodromy matrix from the Hill
matrix involves the action of the matrix exponential of
the Hill matrix applied to a smaller sparse matrix, fol-
lowed by a projection to the n× n monodromy matrix.
Finally, the stability of the periodic solution candirectly
be assessed from the n eigenvalues of the monodromy
matrix, known to be the Floquet multipliers.

Parts of the research in this work were presented
in a preliminary form at the ENOC2020+2
Conference [32], in particular concerning the proposed
choice of the Koopman basis functions in Sect. 3 and
the connection to the harmonic balance equations and
the Hill matrix. The main (and original) contributions
of this work are the novel stability method of Sect. 4,
the considerations with respect to the matrix projection
and the formal proofs for the theorems in Sect. 3.2.

The paper is structured as follows. Section2 pro-
vides an overview of the notation and gives a theo-
retical background for the concepts that are central to
this work, in particular concerning a selection of top-
ics from Koopman theory as well as frequency-based
methods for periodic systems. Section3 introduces
the chosen basis for a Koopman lift on time-periodic
systems and states the three central theorems which
relate this Koopman lift to the classical frequency-
based methods. The proofs of these theorems can be
found in Appendix B. Section4 presents the novel sta-
bility method based on the findings from Sect. 3, and
the projection to the monodromy matrix as well as the
computational effort are discussed. These results are
illustrated in Sect. 5 using numerical investigations on
two exemplary dynamical systems. Finally, concluding
remarks are given in Sect. 6.

2 Theoretical background

In this section, the reader is provided with an overview
over the theory in the Koopman framework and the
Floquet theory that is necessary for the later parts of
the paper.

2.1 Notation and terminology

The frequency-based methods considered in this work
rely heavily on being represented as (generalized or

classical) Fourier series. Hence, a short overview over
Fourier series and the notation that is employed will be
given.

Scalar quantities will be represented by Greek or
Latin slanted lower case letters. This includes scalar-
valued functions as elements of a function space. Vec-
tors in Euclidean space will be represented by Latin
bold lowercase letters and matrices will be represented
by Latin bold uppercase letters. Tuples of functions
are represented by bold font and the distinction to
Euclidean space can be drawn from context.

The choice of index is connected to its meaning. The
index l ∈ N is used for indexing over the states of a
system, whereas the index k ∈ Z is used for index-
ing over frequency harmonics. While j may appear as
arbitrary auxiliary index, the letter i usually denotes
the imaginary number i2 = −1 and is only used for
indexing purposes if the indexing context is obvious.
The Kronecker delta is denoted by δ jl .

If an inner product 〈·, ·〉 with the usual inner prod-
uct properties (see, e.g., [33]) is defined on a vector
space F , elements f, g ∈ F of this space are orthogo-
nal if their inner product is zero. An orthonormal sys-
tem is a set

{
ξ j
}D
j=1 , D ∈ N ∪{∞}with 〈ξi , ξ j

〉 = δi j ,
and a maximal orthonormal system whose span is a
dense subset of the considered vector space is called an
orthonormal basis.

By slight abuse of notation, the inner product is
extended in this paper to tuples of functions g ∈ F l ,
h ∈ Fm element-wise via

〈g,h〉 :=
⎛

⎜
⎝

〈g1, h1〉 . . . 〈g1, hm〉
...

. . .
...

〈gl , h1〉 . . . 〈gl , hm〉

⎞

⎟
⎠ . (1)

It can be easily verified that, as a generalization of the
conjugate symmetry of the inner product, the relation

〈h, g〉 = 〈g,h〉∗ (2a)

holds, whereU∗ denotes the conjugate transpose of the
matrix U. Moreover, the sesquilinearity of the inner
product extends to matrices via

〈Ug,h〉 = U 〈g,h〉 (2b)

〈g,Vh〉 = 〈g,h〉V∗ (2c)

for constant complexmatricesU,V of appropriate size.

123



8442 F. Bayer, R. I. Leine

If a space F has an orthonormal basis
{
ξ j
}D
j=1

stacked into a tuple ξ , it is well-known [33] that ele-
ments g ∈ F admit a (generalized) Fourier series

g =
D∑

j=1

〈
g, ξ j
〉
ξ j = 〈g, ξ〉 ξ , (3)

where the matrix inner product notation was used in
the last term. The constant, scalar coefficient

〈
g, ξ j
〉
is

called the j-th Fourier coefficient of g. In particular,
in the space of trigonometric functions of period T ,
the functions uk : [0 T ) → C, t 	→ eikωt , k ∈ Z con-
stitute an orthonormal basis with respect to the inner
product

〈 f, g〉 =
∫ T

0
f (t)ḡ(t)dt. (4)

This is the classical Fourier series. For a finite-
dimensional subspace spanned by finitely many basis
functions {uk}Nu

k=−Nu
, a function g can thus be expressed

by

g = 〈g,u〉 u (5)

with the matrix-valued inner product notation intro-
duced earlier.

2.2 Koopman theory overview

A short introduction to the Koopman framework to set
the notation and help the reader understand the follow-
ing sections is given below. It is, however, not intended
to give a comprehensive understanding of the current
overall state of the art in the field that would be suitable
for a wider range of application. For such a more gen-
eral and in-depth treatment of the considered method-
ology, the authors rather recommend [10,34].

Consider a non-autonomous time-periodic finite-
dimensional dynamical system governed by

ẋ = f(x, t), (6)

where t ∈ R is the time, x(t) ∈ X ⊆ R
n is the state

trajectory starting at x(t0) = x0 and f : X × R → X
is a smooth vector field which is T -periodic in t . The
family ofmaps φt (x0, t0) = x(t) characterizes the flow

of the system and assigns to each (initial) configuration
(x0, t0) the resulting state at time t ≥ t0.

The Koopman framework [10] considers output
functions g(x, t), also called observables. Any Banach
spaces of functions over the complex or real numbers
are permitted in the general Koopman framework. In
this work, we consider in particular the space F of
complex-valued functions g : X × R → C which
are real analytic on X and T -periodic in the last argu-
ment t . Given any function g, it may be of interest how
its function values evolve along the trajectories of the
system. For instance, in the Lyapunov framework, it is
desired that function values of a Lyapunov candidate
decrease over time for any starting point. The operator
K t : F → F; g 	→ g ◦φt performs this shift along the
trajectory for arbitrary functions g from the considered
function space. The family of all these operators for
any t is called the Koopman semigroup of operators.
Indeed, the semigroup properties with respect to the
time parameter can be verified easily.

For suitable function spacesF , this Koopman semi-
group contains all information about the system with-
out explicitly knowing the vector field f or the flow φt .
In particular, ifF is chosen such that the identity func-
tion id is contained in the vector space, then the flowcan
be recovered easily by simply evaluating K t (id). As a
trivial counterexample, consider the one-dimensional
vector space of constant functions. Any constant func-
tion will not change its function value while being eval-
uated along an arbitrary trajectory of arguments. There-
fore, in this case, the Koopman operator semigroup is
well defined, albeit trivial. No information about the
underlying system is retained. This example shows that
an appropriate choice of function space is a crucial part
of the Koopman framework.

Under the aforementioned assumptions for the par-
ticular function space F and the vector field f , the
Koopman semigroup is continuous with respect to time
and there also exists the operator L : F → F with

g 	→ ġ = limt→0
K t g−g

t , mapping an observable g to
its total time derivative ġ along the flow with

ġ(x, t) = ∂g(x, t)
∂x

f(x, t) + ∂g(x, t)
∂t

. (7)

The operator L is called the infinitesimal Koopman
generator. Again, for suitable F , this representation
alone is a sufficient way to describe the behavior of the

123



Sorting-free Hill-based stability analysis 8443

Fig. 1 Schematic drawing of the infinitesimal Koopman genera-
tor L , its finite-dimensional approximation L N̂ and theKoopman

lift Â

dynamical system. In particular, if id ∈ F , the vector
field f is easily recovered.

In addition, the infinitesimal Koopman generator
and the Koopman semigroup of operators are linear
in the argument g, even if the governing differential
equation is nonlinear. This comes at the cost of deal-
ing with a mapping on an (infinite-dimensional) func-
tion space F instead of the (finite-dimensional) state
space X .

For practical reasons,wewill be forced to project the
dynamics on a finite-dimensional subspace FN̂ ⊂ F
spanned by N̂ linearly independent basis functions
{
ψ j
}N̂
j=1. Any projection �N̂ : F → FN̂ defines a

finite-dimensional approximation L̂ : FN̂ → FN̂ of
L on FN̂ by L N̂ := �N̂ L . The approximation pro-
cess and the subsequent approximation error are visu-
alized in Fig. 1. As the subspace FN̂ generally is not
closed w.r.t. L , the result of Lg must be projected back
onto FN̂ , introducing some approximation error. As
forF , the choice of the finite spaceFN̂ and the projec-
tion onto it is crucial.

To evaluate a system trajectory in the usual state-
space form, the initial condition is fixed first, the state
evolution is computed afterward and the output func-
tion is computed last. For the Koopman infinitesi-
mal generator (and its approximation on the finite-
dimensional function space), this is reversed. First, an
observable, or output, is fixed, the observable is prop-
agated along the flow and the initial condition is deter-
mined in the last step by evaluating (Lg)(x) for an ini-
tial condition x. To arrive back at a linear system rep-
resentation in the usual state-space form, this behavior
must be reversed again. This is achieved by evaluating
the action of the (approximate) Koopman infinitesimal
generator on the basis functions. Setting

ẑ j (0) := ψ j (x(0)) (8a)

˙̂z j (0) := (L N̂ ψ̂ j (x(0)) (8b)

and keeping this dynamics for increasing t , the linear
autonomous system

˙̂z = Âẑ (9a)

ẑ(0) = �̂(x(0)) :=
⎛

⎜
⎝

ψ1(x(0))
...

ψN̂ (x(0))

⎞

⎟
⎠ (9b)

results. This linear system is called the Koopman lift
and describes the dynamics represented in L N̂ . With
the lifted states

ẑ(t) ≈ �̂(x(t)) (9c)

it is a finite-dimensional linear approximation of the
original system dynamics. The Koopman lift matrix Â
can be derived from the original nonlinear dynamics
manually using a column vector � := (ψ1, . . . , ψN̂

)T

of the basis functions of FN̂ , computing d�
dt and iden-

tifying terms linear in elements of � after applying the
projection �N̂ from F onto FN̂ . If this projection is
orthonormal and � is an orthonormal system, then the
matrix entry Âi, j at i-th row and j-th column is given

by Âi, j =
〈
dψi
dt , ψ j

〉
, where 〈., .〉 denotes the corre-

sponding inner product.
Depending on the basis structure chosen, various

popular embedding techniques emerge as a Koopman
lift for specific system classes. For instance, if a mono-
mial basis is chosen for an autonomous system, the
Carleman linearization [35] results as Koopman lift.
For smooth, polynomial systems, this is often the first
choice of basis dictionary [36]. Alternatively, delay
coordinates are often employed in Koopman-based
applications [2]. Based on the Takens embedding theo-
rem, this can capture weakly nonlinear dynamics [37].
For periodic systems, the Fourier embedding [38] has
been known, although its properties have mainly been
analyzed in the frequency domain.

2.3 Harmonic balance method

Consider the non-autonomous time-periodic finite-
dimensional dynamical system (6) as above. Often one
is interested in finding a T -periodic solution, i.e., a
solution to the dynamical system (6) which fulfills

123



8444 F. Bayer, R. I. Leine

x(t + T ) = x(t) for all t ≥ 0. This constitutes a
boundary value problem (BVP) and there are meth-
ods to solve this type of BVP in the time domain and in
the frequency domain. Shooting, multiple shooting and
collocation methods all rely on an interplay between
time-integration (or finite differencing) of the ODE and
solving nonlinear functions for periodicity and conti-
nuity constraints [21].

In contrast, the HBM is a frequency-based method.
Under suitable smoothness assumptions, the periodic
solution has a convergent Fourier series. Hence the
periodic solution can be approximated by its Fourier
expansion up to order NHBM with unknown parame-
ters via

xp(t) =
NHBM∑

k=−NHBM

pkeikωt (10)

with ω = 2π
T , u(t) = (e−i NHBMωt , . . . , ei NHBMωt )

being a vector of Fourier base functions and{
p−NHBM , . . . ,pNHBM

}
gathering the corresponding

(unknown) coefficients. These coefficients pk are then
determined by substituting this ansatz into the system
equation (6). The comparison of coefficients for the
Fourier expansions of dxp

dt from the definition (10) and
f(xp, t) for every order up to NHBM transforms theBVP
into a system of n(2NHBM + 1) algebraic equations.
Existence and convergence of these HBM approxima-
tions has been shown [22]. While the left-hand side
of the equation as well as linear terms in f are easy
to handle, the frequency component of the nonlinear
terms can usually not be expressed in closed form. The
individual equations for each order are thus usually
determined and simultaneously solved using the fast
Fourier transform with an alternating frequency and
time (AFT) method [24]. The equations for each order
can also be isolated by projecting onto the correspond-
ing basis function from the collection in u through the
classical inner product (4). Hence, the HBM approxi-
mates a periodic solution by solving the n(2NHBM+1)
algebraic equations collected in

〈
dxp
dt

,u
〉

= 〈f(xp(·), ·),u
〉
. (11)

With this notation, the numerically cumbersome task
of calculating the Fourier coefficients of the nonlinear

components of f is hidden in the definition of the inner
product.

2.4 Floquet theory: stability of periodic solutions

When a periodic orbit xp is found (via HBMor by other
means), the next interesting question is that of its stabil-
ity properties; that is, whether trajectories that start suf-
ficiently close to the periodic orbit will approach it, stay
in a vicinity of it or tend away from it with increasing
time. To evaluate the stability properties, the dynamics
of a perturbation y = x − xp from the periodic solu-
tion is considered. Substitution of this definition into
the original system dynamics yields

ẏ = f(xp + y, t) − ẋp := J(t)y + O(‖y‖2), (12)

where J(t) = ∂f
∂x

∣∣∣
xp(t),t

is the Jacobian of the system

evaluated along the periodic solution. The approximate
linear time-varying (LTV) system

ẏ(t) = J(t)y(t) (13a)

y(0) = x(0) − xp(0) (13b)

has an equilibrium at zero, which corresponds to the
periodic orbit of the original system, and the stability
analysis of the periodic orbit reduces (in the hyperbolic
case) to the stability analysis of this equilibrium. This
will be the convention for the remainder of this paper,
unless stated otherwise.

The fundamental solutionmatrix�(t) is the solution
to the variational equation

�̇(t) = J(t)�(t) ; �(0) = I (14)

and any state can be obtained via y(t) = �(t)y0. In par-
ticular, the fundamental solution matrix �(T ) =: �T

evaluated after one period is called the monodromy
matrix of the system and its eigenvalues {λl}nl=1 are
called Floquet multipliers [39]. The Poincaré map
yk+1 = �T yk provides snapshots for the evolution of
the perturbation y, spaced at a time distance of T , i.e.,
yk = y(kT ). For the long-term behavior, it is sufficient
to consider the evolution of these snapshots. There-
fore, stability analysis of the periodic solution reduces
to stability analysis of the Poincaré map. Hence, if all
Floquet multipliers are of magnitude strictly less than

123



Sorting-free Hill-based stability analysis 8445

one, the equilibrium of the perturbed LTV system and
thus the periodic solution of the original system are
asymptotically stable; if at least one eigenvalue has a
magnitude strictly larger than one, they are unstable.
If there exist Floquet multipliers with magnitude equal
to one, but none with a magnitude larger than one, the
equilibrium is non-hyperbolic and further investigation
is necessary to give conclusive statements about stabil-
ity of the originally considered periodic solution.

Alternatively to the Floquet multipliers, the stabil-
ity properties of a time-periodic linear system can be
characterized by the Floquet exponents. In the lin-
earized perturbed system (13a), if the matrix �T is
diagonalizable, there exist n solutions yl(t) = pl(t)eαl t

which form a basis of the solution space, where each
function pl is T -periodic [39]. Hence, stability is char-
acterized by the real parts of the Floquet exponents
{αl}nl=1. If at least one Floquet exponent lies in the open
right half plane, i.e., if at least one real part is larger than
zero, the equilibrium is unstable. The Floquet multipli-
ers can be determined by substituting t = T in the
Floquet solution and it follows that

λl = eαl T , l = 1, . . . , n. (15)

In contrast to the Floquet multipliers, the Floquet
exponents are not uniquely defined. It is easy to
see that if the pair (pl(t), αl) generates a solution
yl(t), the same solution is generated by (p̃l(t), α̃l) =
(pl(t)e−ikωt , αl + ikω) with k ∈ Z. Hence, in total,
there are infinitely many valid Floquet exponents,
which can be categorized into n distinct groups. All
elements of one group have the same real part and
differ in the imaginary part by multiples of iω. As
stability is determined by the real part only, it is suf-
ficient for stability analysis to know any one element
from each of the n groups. All elements of one group
map to the same Floquet multiplier.

When a periodic orbit is determined using the purely
time-domain-based shooting method, the monodromy
matrix usually is a direct byproduct of the continua-
tionmethod [17]. In this case, the numerically obtained
monodromy matrix can be evaluated directly to obtain
the Floquet multipliers and their stability information.

When the HBM is computed in the standard way,
however, stability information about the identified limit
cycle is unclear without further investigation. The Hill
method [22,40] offers a frequency-domain-based way

to approximate the Floquet exponents of the linearized
perturbation equation.

The Floquet exponents are eigenvalues of the infi-
nite Hill matrix H∞ [30], which is constructed from
the Fourier coefficients of the periodic system matrix
J(t) =∑∞

k=−∞ Jkeiωkt and reads as

H∞ =

⎛

⎜⎜⎜⎜⎜
⎝

. . .
...

...
... . .

.

. . . J0 + iωI J−1 J−2 . . .

. . . J1 J0 J−1 . . .

. . . J2 J1 J0 − iωI . . .

. .
. ...

...
...

. . .

⎞

⎟⎟⎟⎟⎟
⎠

. (16)

Introducing a vector v∞ of appropriate (infinite) length,
the eigenproblem

H∞v∞ = α̃v∞ (17)

can be formulated. This infinite-dimensional prob-
lem has infinitely many discrete eigenvalues α̃, which
solve (17). They correspond identically to the Floquet
exponents α̃ for all k ∈ Z as introduced above and
can be sorted into n groups, where the entries of each
group differ by multiples of iω [30]. However, in prac-
tice, only the eigenvalues of a finite-dimensionalmatrix
approximation of H∞ can be computed numerically.
The matrix

H =
⎛

⎜
⎝

J0 + i NuωI . . . J−2Nu
...

. . .
...

J2Nu . . . J0 − i NuωI

⎞

⎟
⎠ (18)

of size n(2Nu + 1) × n(2Nu + 1) consists of the
n(2Nu + 1) most centered rows and columns of H∞
and approximates the original infinite-dimensional Hill
matrix. In the absence of truncation error, the eigenval-
ues of H would be a subset of the eigenvalues of H∞,
i.e., the Floquet exponents. Due to the inevitable error
which generally comes with truncation, however, this
does not quite hold. The N eigenvalues ofH, which do
not identically coincide with Floquet exponents, will
be called Floquet exponent candidates below.

The matrix H has a block Toeplitz structure except
for the middle diagonal, and, for sufficiently large Nu,
the bands near the diagonal dominate as the Fourier
coefficients of J tend to zero. Loosely speaking, some
eigenvalues affiliated most with the central rows of
H are less impacted by the truncation and provide a

123



8446 F. Bayer, R. I. Leine

better approximation to the Floquet exponents than
others [28]. Up until the change of the last century,
this property was neglected and stability of the peri-
odic solution was asserted based on the real parts of
all Floquet exponent candidates [40]. This naive Hill
method without any additional steps would often assert
instability for stable solutions due to spurious Floquet
exponent candidates without physical meaning, giving
it a reputation of being inaccurate [22,29].

However, the accuracy of the Hill method can be
improved significantly if only a subset of Floquet expo-
nent candidates is considered, instead of all of them.
Hence, the search for a selection criterion which deter-
mines the best approximation to the Floquet exponents
from the Floquet candidates has received much atten-
tion in the literature [23,26,28,30].

For sufficiently large Nu, it is proven that the can-
didates with minimal imaginary part in modulus con-
verge to the true Floquet exponents [26]. Since the
convergence may only occur for very large truncation
orders, an addition to this method was very recently
proposed [31]. This modified method first sorts the
Floquet exponent candidates based on their real parts,
before applying the imaginary part criterion to themost
highly populated groups. However, in this real-part-
based method, it is unclear how to proceed if two or
more true Floquet exponents have the same real parts,
and thus not enough groups are available.

An alternative criterion selects those candidates
whose eigenvectors are most symmetric [27,28], as
they should correspond most to the middle rows of
the matrix H. This symmetry is computed based on a
weighted mean. Even though there currently is no for-
mal convergence proof for this symmetry-based sorting
method, some numerical results indicate faster conver-
gence than with the aforementioned eigenvalue crite-
rion [12], while other results do not support this claim
[31].

For all these criteria, there are currently no meth-
ods to efficiently and accurately compute only those
eigenvalues which fit the criterion. Rather, all eigen-
pairs of the large matrix have to be computed first and
then most of them are discarded. As the cost of solv-
ing an eigenvalue problem of a N × N matrix is of the
order O(N 3), the computational cost of the approach
is usually dominated by determining the eigendecom-
position of a large matrix [29]. In addition, it is well
known that the accuracy of the computed eigenvalues
is not too high if all eigenvalues of a sparse matrix are

sought. Sparsity of H cannot be reasonably exploited
to decrease the computational cost of solving the com-
plete eigenproblem [41].

3 A Koopman dictionary for time-periodic systems

For smooth autonomous systems in the Koopman
framework, it is customary to choose as basis functions
the so-called Carleman basis [2,42], i.e., a finite set of
monomials ψβ(x) = xβ , where β ∈ N

n is a multi-
index and standard multi-index calculation rules (see
AppendixA) apply.As time-periodic functions are con-
sideredhere,wepropose in this paper to include as basis
functions combinations of monomial terms as well as
Fourier terms of the base frequency, i.e., basis functions
of the formψβ,k := xβeikωt , whereω = 2π

T . The func-
tions

{
ψβ,k |k ∈ Z,β ∈ N

n
0

}
are an orthonormal system

within the initially considered vector spaceF w.r.t. the
inner product

〈g, h〉 :=
∫ T

0

⎛

⎝ 1

(β!)2
∑

β∈Nn

∂βg

∂xβ

∣∣∣∣∣
0,t

∂β h̄

∂xβ

∣∣∣∣∣
0,t

⎞

⎠ dt.

(19)

This inner product contains the standard inner product
for Fourier series, and the derivatives serve as an inner
product for the monomials. The inner product proper-
ties can be readily verified.

Let Nz, Nu ∈ N be integers which describe the
assumed maximum polynomial and frequency order,

respectively. There is a set B = {β j }Nβ

j=1 collecting
all multi-indices with 1 ≤ ‖β‖ ≤ Nz. These are all
multi-indices that create monomials xβ of degree Nz

and less. By (A4), it holds that Nβ = (Nz+n
n

) − 1.
For the sake of brevity, define N = Nβ(2Nu + 1) and
N̂ = N + (2Nu + 1). The set

{
ψβ,k |β ∈ B ∪ {0} , |k| ≤ Nu

}
(20)

of orthogonal basis functions spans a specific finite-
dimensional subspace FN̂ ⊂ F . These basis functions
are the monomials up to degree Nz, multiplied to the
Fourier base functions up to frequency order Nu. With
k = 0 and ‖β‖ = 1, the set (20) includes the identity
function for each state.Moreover, since themulti-index
β = 0 is permitted in the set (20), the purely time-
dependent classical Fourier base functions eikωt are

123



Sorting-free Hill-based stability analysis 8447

represented. These basis functions are collected into
two vectors. The vector

u(t)T := (e−iωNut , . . . , e0, . . . , eiωNut ) (21)

collects all basis functions that are not dependent on
the state, but only on time, such that the evolution of u
is not a product of the system dynamics, but is known
a priori. All other state-dependent basis functions are
collected into the vector

�z(x, t) :=

⎛

⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜
⎝

xβ1e−iωNut

...

xβ1e0
...

xβ1eiωNut

xβ2e−iωNut

...

x
βNβ eiωNut

⎞

⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟
⎠

. (22)

Here, the basis functions are ordered by monomial
exponent first, and then all frequency base functions
for the samemonomial are grouped together in ascend-
ing order. It is notable that any other orderings are also
applicable and the matrices can be transformed into
each other by similarity transforms. The ordering (22)
has been chosen purely for convenience in later argu-
mentation.

If the full basis as introduced above is considered,
the state x itself is included in the basis. Hence, there
is a selector matrix Cz ∈ R

n×N containing select rows
of the identity matrix withCz�z(x, t) = x. There exist
other options to recover x from �z(x, t). For instance,
monomials of higher (uneven) order can be used
by considering the corresponding root. Also, the
matrix Cz can be allowed to be time-dependent. This
second option and its implications will be investigated
in more detail in Sect. 4.3. Both vectors �z and u are
collected into one vector � := (�T

z ,uT)T for conve-
nience.

As introduced inSect. 2.2, theKoopman lift describes
the evolution of the basis functions under the approxi-
mate infinitesimal Koopman generator, and it is deter-
mined by expressing the time derivative along the flow
of all basis functions in the finite-dimensional basis,
projecting them back again to the subspace. Using the
projection defined by the inner product (19), this gives

�̇ = 〈�̇,�
〉
� + r , (23)

where r ∈ F N̂ is the remainder that is orthogonal to the
finite-dimensional subspace and will be projected out.
Substituting (23) into the definition (9) for the Koop-
man lift and separating the two vectors of basis func-
tions yields

˙̂z :=
(
ż
u̇

)
= 〈�̇,�

〉
ẑ (24a)

=
(〈

�̇z,�z
〉 〈

�̇z,u
〉

〈u̇,�z〉 〈u̇,u〉
)(

z
u

)
(24b)

=:
(
A B
0 ∗
)(

z
u

)
, (24c)

where A ∈ C
N×N and B ∈ C

N×(2Nu+1) are con-
stant coefficient matrices. The lower rows of the large
matrix in (24b) determine the dynamics of u. However,
since u is purely time-dependent and its time evolution
is known a priori, the lower rows of (24b) are superflu-
ous and the original nonlinear system is approximated
by the LTI system

ż = Az + Bu (25a)

z(0) = �z(x0, 0) (25b)

for the very specific input u as in (21). The lifted state
vector z(t) is an approximation to the state-dependent
basis functions �z(x(t), t) evaluated along the flow of
the original system.

This is only an approximation, however, due to the
projection onto the finite-dimensional function space.
In particular, for t ≥ 0, the lifted state z will cease to
adhere to the constraints posed by �z, meaning that
there may not exist any vector x̃ ∈ R

n which fulfills
z(t) = �z(x̃, t).

Many results of this work tie in to classical results
that are obtained fromJacobian linearizationof the state
space. Therefore, the maximum monomial order in the
basis will often be restricted to Nz = 1. In this case, the
vector of basis functions will be denoted by �z,lin of
size n(2Nu + 1) and for the sake of brevity, these func-
tions will be called linear basis functions, even though
the influence of time is still decidedly nonlinear. To
make the ordering unique, these linear basis functions
are ordered in an ascending order, first by state and then

123



8448 F. Bayer, R. I. Leine

by frequency such that the resulting vector reads as

�z,lin(x, t) =

⎛

⎜⎜⎜⎜⎜⎜⎜
⎝

x1e−i Nuωt

...

x1ei Nuωt

x2e−i Nuωt

...

xnei Nuωt

⎞

⎟⎟⎟⎟⎟⎟⎟
⎠

. (26)

3.1 Koopman lift on the perturbed system

While a nonlinear dynamical systemmay have an arbi-
trary number of attractors of any kind (equilibrium,
periodic, quasi-periodic, chaotic) located anywhere in
the state space, the finite-dimensional linear Koopman
lift (25a) has a very limited range of possible attrac-
tors, i.e., at most one single globally attractive solution,
which is an equilibrium point. This means that the sys-
tem (25a) will most likely not be able to describe the
system dynamics of a nonlinear time-periodic system
globally. For time-periodic systems, periodic solutions
are expected. In this section, the Koopman lift of the
perturbed system around such a periodic solution will
be regarded to assess its local behavior.

As in the standard HBM, a periodic solution is
approximated by its Fourier expansion (10) up to order
NHBM with unknown Fourier coefficients. The per-
turbed system is then given by y(t) = x(t) − xp(t),
with dynamics

ẏ(t) := f̃(y(t), t) (27a)

as in (12). This allows theKoopman lift to be performed
on the perturbed system, i.e., on functions of the state y
evaluated along the flow f̃ . Now, the origin of the lifted
state does correspond to the origin equilibrium of the
perturbed dynamics, or a periodic solution of the orig-
inal dynamics. This also means that the Koopman lift
matrices A and B now depend on the Fourier coeffi-
cients p of the periodic solution around which we are
linearizing. The following sections shed light on the
properties of the matrices A(p) and B(p).

3.2 Koopman-based harmonic balance and Hill
equations

One of our key findings is that the B(p) matrix of the
Koopman lift contains information about the parame-
terization of the periodic solution that is also encoded
in the HBM. This is made more precise in the theorems
given in this section.

Let xp =∑k pke
ikωt be a real-valued periodic solu-

tion candidate of the time-periodic system (6). The har-
monic residual of this candidate is given by

r(t) = f(xp(t), t) − ẋp(t). (28)

Since xp as well as f are periodic, the residual is peri-
odic as well, and therefore it has a Fourier series

r(t) =
∑

k

rkeikωt . (29)

The HBM of order M returns solution candidates for
which rk = 0 for all |k| ≤ M . Since r is real-valued, it
also holds that r−k = r̄k , where rk = [rk,1, . . . , rk,n]T
is a vector with n complex-valued entries.

This definition of the residual allows to concisely
relate the HBM to the Koopman lift.

Theorem 1 Let ż = A(p)z + B(p)u be the lifted
dynamics of frequency order NHBM of system (6)
around an unknown periodic ansatz of the form (10).
The NHBM-th order HBM equations (11), i.e., rk =0,
|k| ≤ NHBM, are given byCzB(p) = 0, whereCz is the
constant selection matrix that fulfills y = Cz�z(y, t)
for all t.

The formal proof of this theorem is given in
Appendix B.1. From an intuitive point of view, this
condition is not surprising. For a periodic solution, the
lifted dynamics has an equilibrium at zero, meaning
that if y = 0 it should also hold that ẏ = 0. Since
the lifted dynamics approximates z(t) ≈ �z(y, t)
and this approximation holds identically at the initial
point (25b), it should hold that ẏ ≈ Czż is zero if
y = 0. This equation can be evaluated for an arbitrary
initial point on the periodic ansatz via a time-shift.With
y(0) = 0 and z(0) = �z(y(0),0) = 0, it turns out that
the only remaining summand in the approximate y-
dynamics is CzB(p).

If only basis functions that are linear in the state
(Nz = 1) are considered, the above result still holds.

123



Sorting-free Hill-based stability analysis 8449

Moreover, we can state Theorem 2 about all entries of
B and not just specific rows.

Theorem 2 Let ż = A(p)z + B(p)u be the lifted
dynamics of system (6)with linear basis functions�z,lin

of frequency order Nu that are sorted as in (22), eval-
uated for the perturbed system around an unknown
periodic ansatz of the form (10) up to frequency order
at least NHBM = 2Nu. Then, the matrix B(p) ∈
C
n(2Nu+1)×(2Nu+1) consists of n stacked Toeplitz matri-

ces. The l-th Toeplitz matrix Bl contains as entries
(ignoring duplicates) precisely the 4Nu + 1 residuals
rk,l(p), |k| ≤ 2Nu that follow from the HBM w.r.t the
l-th state.

If B(p) = 0, then all these residuals of the HBM
vanish. Conversely, if p solves the HBM equations
rk(p) = 0, |k| ≤ 2Nu, then it holds that B(p) = 0.

The formal proof of this theorem is given in
Appendix B.1.

In addition to theBmatrix, theAmatrix of theKoop-
man lift also holds frequency information about stabil-
ity of the periodic solution. This is summarized in the
following theorem.

Theorem 3 Let ż = A(p)z + B(p)u be the lifted
dynamics around a periodic solution of system (6) with
linear basis functions�z,lin of frequency order Nu that
are ordered as in (22). Then the Hill matrix H, trun-
cated to frequency order Nu, for the periodic solution
parameterized by p results from the matrix A(p) by
the similarity transformH = UA(p)UT, where U is an
orthogonal permutation matrix that satisfiesU�z,lin =
(yTei Nuωt , . . . , yTe−i Nuωt )T.

The formal proof of this theorem, again based on
explicit evaluation of the inner product in the Koopman
lift, is given in Appendix B.2.

With the three above theorems, qualitative insight
about the accuracy of the presented Koopman lift can
be gained. Locally (in the vicinity of a periodic solu-
tion), the lifted system contains the same dynamical
information that is encapsulated in the Hill matrix of
the same frequency order. Convergence results about
the HBM [22] and the Hill method [26] can thus be
related to the accuracy of the Koopman lift. Within
the applied Koopman community, such a link is quite
unusual. For many applications, convergence results
for the finite-dimensional Koopman lift do not exist at
all. The convergence results that do exist (e.g., [43,44])

state that, under special conditions, a given error toler-
ance can be reached using a large number of specific
basis functions,without giving an explicit upper bound.
Hence, for time-periodic systems of the form (6), the
proposed class of basis functions is a favorable choice
due to its added connection with the Hill matrix, which
indicates that the Koopman lift retains valuable stabil-
ity information.

4 Sorting-free stability method

The Koopman lift with Theorems 2 and 3 gives a
dynamical systems interpretation for the Hill matrix,
allowing for a novel stability method based on the Hill
matrix. This is demonstrated in the following section.

4.1 Approximating the monodromy matrix

Consider the Koopman lift as in Theorems 2 and 3,
i.e., consider the linear basis �z,lin evaluated for the
perturbed system around a periodic solution which is
determined up to a frequency order of 2Nu. From The-
orem 2, we know that B = 0 and from Theorem 3, that
A = UTHU. The linear dynamical system ż = Az+Bu
resulting from the Koopman lift therefore reduces to

ż(t) = Az(t) (30a)

y(t) ≈ zy(t) = C(t)z(t) (30b)

z(0) = �z,lin(y(0), 0) , (30c)

where C(t) ∈ C
n×N is a possibly time-dependent pro-

jection matrix that satisfies C(t)�z,lin(y, t) = y for all
t ∈ [0, T ). There is a nNu-parameter family of choices
for C if it is allowed to be time-dependent, which will
be investigated further in Sect. 4.3. However, the naive
choice would be to pick the entries in (26) that corre-
spond to frequency zero. In this case, the matrix C is
constant and given by

C = In×n ⊗ (0 . . . 0 1 0 . . . 0
)

, (31)

where In×n is the n × n identity matrix, ⊗ denotes the
Kronecker product and the second matrix in (31) is a
row vector ∈ R

2Nu+1 with zeros everywhere except for
the middle column.

Since for t = 0 all exponential terms are 1 andvanish
in (26), the vector �z,lin(y, 0) can be expressed as a

123



8450 F. Bayer, R. I. Leine

matrix product

�z,lin(y, 0) = Wy (32)

for all y. The matrixW ∈ R
N×n consists of (repeated)

rows of the identity matrix. More specifically, it can be
expressed as

W = In×n ⊗
⎛

⎜
⎝

1
...

1

⎞

⎟
⎠ , (33)

where the second matrix in (33) is a column vector
∈ R

2Nu+1 filledwith ones.As the system (30a)–(30c) is
a linear time-invariant (LTI) system (except for the pos-
sibly time-dependent matrix C), its closed form solu-
tion can be explicitly computed as

y(t) ≈ zy(t) = C(t)eAtz(0) (34a)

= C(t)eAtWy(0). (34b)

As a key finding of the current paper, the matrix
C(t)eAtW ∈ R

n×n is an approximation of the funda-
mental solution matrix �(t), which is the matrix that
satisfies y(t) = �(t)y(0). In particular, for t = T , the
monodromy matrix is approximated via

�T ≈ C(T )eATW. (35)

If the Hill matrixH is already known due to other com-
putations, e.g., as a by-product of a frequency-based
continuation method such as MANLAB [27], then the
similarity transform can be substituted into the matrix
exponential to yield

�T ≈ C(T )UTeHTUW =: C̃(T )eHT W̃ , (36)

where C̃, W̃ can also be computed directly via

W̃ =
⎛

⎜
⎝

In×n
...

In×n

⎞

⎟
⎠ (37a)

C̃ = (0 . . . 0 In×n 0 . . . 0
)

(37b)

by making use of the permutation properties of U (see
Thm. 3). The naive choice for C̃ was employed here
for demonstration purposes, although all other choices

are also applicable. Equations (35)–(37b) show that the
choice of basis function ordering, and thus the exact
numerical contents of the matrix, differ only formally
and can be easily transformed into each other. For the
sake of clarity, only the approximation (35)will be used
in the following sections, unless otherwise stated. All
results can, however, be transferred analogously to the
formulation (4.1).

4.2 Stability and computational effort

As (35) is an approximation of the monodromy matrix,
the Floquet multipliers can be approximated by its
eigenvalues. This constitutes a novel projection-based
stability method based on the Hill matrix. It is labeled
sorting-freebecause the approachdoes not use the same
steps as standard stability methods based on the Hill
matrix [28,30,31].

In all standard approaches, the complete set of eigen-
values of the Hill matrix are determined in the first
step. This is a computationally intense operation which
may lack accuracy [41]. The result of this first oper-
ation is N = n(2Nu + 1) Floquet exponent candi-
dates, of which only n are the best approximations to
the true Floquet exponents. There are various sorting
algorithms that vary in computational expense, which
aim to find these n best approximations as introduced
in Sect. 2.4.

In contrast to these approaches, we now propose
a sorting-free method. The Hill matrix is constructed
identically to the aforementioned approaches. How-
ever, once the Hill matrix is obtained, the sorting-
free projection-based stability method goes a different
route. The approximation to the monodromy matrix
is determined first via (35). In the most general case,
a matrix exponential is of the same complexity range
O(N 3) as an eigenvalue problem [41,45], however this
is not the case for the specific operations presented here.
We give two reasons:

1. The matrix exponential in (35) is always multiplied
by the matrix W, which is smaller and relatively
sparse. The matrix exponential can thus be com-
puted by multiple evaluations of the action of the
matrix exponential on a sparse vector. For this oper-
ation, efficient scaling and squaring approaches
which utilize these properties exist [46].

2. In many applications the constructed Hill matrix
will be relatively sparse. This is the case if the fre-

123



Sorting-free Hill-based stability analysis 8451

Fig. 2 Flowchart comparing the three general stability
approaches. For each of the methods, the most computationally
intense step is located in the second row

quencies of the dominant harmonics of the original
system are small in comparison to the frequency
order. This sparsity in H, or equivalently in A, can
be exploited in the scaling and squaring algorithm,
in addition to sparsity inW.

As a second step in the sorting-free approach, it remains
to solve an eigenvalue problem for the approximation
of�T . This matrix is only of the size n×n, i.e., in gen-
eralmuch smaller than the size n(2Nu+1)×n(2Nu+1)
of the Hill matrixH. The result are n Floquet multipli-
ers and no a posteriori sorting of candidates is neces-
sary. The reduction of candidates is performed implic-
itly by the projectionmatricesC andW. This reduction
through projection is essentially different from sorting,
as the resulting Floquetmultipliers from the projection-
based method do generally not coincide identically
with any Floquet multiplier candidates obtained as
eigenvalues of the Hill matrix, transformed from Flo-
quet exponents to Floquet multipliers via (15).

Alternatively to theHillmatrix approaches, themon-
odromy matrix can be determined via time-integration
of the variational equation (14) [29].Afterward, it again
remains to solve the eigenvalue problem on the n × n
monodromy matrix. The three general approaches are
visualized in Fig. 2. While the starting point for the
novel method is the Hill matrix as in the standard Hill

methods, the final steps of our method are instead iden-
tical to the time-integration method. However, the sec-
ond indicated step,whichneeds themost computational
effort in all three approaches, is different. The presented
sorting-free projectionmethod therefore has the advan-
tage of being aHill-basedmethod,which is favorable in
an HBM setting, and at the same time only requiring to
compute the smaller eigenvalue problem of the mon-
odromy matrix, similar to the time-integration-based
method.

4.3 Choice of projection matrix

Until now, it has been established that C(T )eATW
approximates the monodromy matrix if C�z,lin = y.
However, the choice of the matrix C to achieve this is
not unique. In addition to the naive choice (31), many
more matrices are admissible. Recall that the vector
�z,lin (26) is sorted such that it contains n consecu-
tive blocks of length 2Nu + 1, each block collecting
all functions of a single state. Therefore, the l-th row
of C, which singles out the l-th state, consists of n − 1
blocks of zeros, each of length 2Nu + 1, and one block
cl(t) ∈ R

1×(2Nu+1) that is possibly filled with nonzero
entries, yielding a non-square block-diagonal structure

C(t)=

⎛

⎜⎜⎜
⎝

c1(t) 0 . . . 0
0 c2(t) . . . 0
...

...
. . .

...

0 0 . . . cn(t)

⎞

⎟⎟⎟
⎠

∈Cn×n(2Nu+1) (38)

with cl(t) ∈ C
1×(2Nu+1); l = 1, . . . , n. With this

notation, the requirement C(t)�z,lin(y, t) = y for all
t ∈ [0, T ] can be simplified to

cl(t)u(t) = 1 ∀t ∈ [0, T ]; l ∈ 1, . . . , n (39)

by noting that yl can be pulled out on both sides of the
equation. As all entries in u are linearly independent,
the time dependency of cl is uniquely given via

cl(t) = ĉl

⎛

⎜
⎝

ei Nuωt . . . 0
...

. . .
...

0 . . . e−i Nuωt

⎞

⎟
⎠ =: ĉlV(t) (40)

to cancel out the time dependency in (39), such that
V(t)u(t) is constant. In this expression, ĉl ∈ C

1×(2Nu+1)

123



8452 F. Bayer, R. I. Leine

is a constant row vector and V(t) is given by V(t) =
diag
(
ei Nuωt , . . . , e−i Nuωt

)
. Due to (39) and y being

real, two additional conditions on ĉl are

Nu∑

k=−Nu

ĉl,k = 1 (41a)

ĉl,k = ¯̂cl,−k . (41b)

All choices for ĉl that satisfy these conditions are
admissible for the projection matrix. If a set

{
ĉl,k
}Nu
k=1

of Nu arbitrary independent complexpositive-frequency
coefficients is given, they can be easily extended to an
admissible choice ĉl using the conditions (41), which
implies that the admissible projection matrices form
a nNu-parameter family. For the sake of readability,
the constraints will be left explicit in the considera-
tions below. In Fig. 6–8 of Sect. 5.2, two choices for ĉ
are compared and it turns out that this choice indeed
strongly influences the approximation quality.

For numerical computations, handling the sparse
matrix C with only few unknown coefficients is cum-
bersome. It is easier to collect all unknown variables
ĉl , l = 1, . . . , n into a long row vector ĉall with

ĉall = (ĉ1 . . . ĉn
)
. (42)

If the true monodromy matrix �T would be known,
then (35) can be viewed as a fitting problem for the
unknown ĉall. The squared matrix residual

L true(ĉall) =
∥∥∥C(T )eATW − �T

∥∥∥
2

(43)

should be minimized. Setting Q = eATW and uti-
lizing the block diagonal structure (38) of C(T ), this
expression decouples row-wise into n (nonlinearly)
constrained least-squares problems

min
cl (T )

∥∥cl(T )T (Q)l -th block − (�T )l -th row

∥∥2

(44a)

subject to
Nu∑

k=−Nu

ĉl,k = 1 (44b)

ĉl,k = ¯̂cl,−k (44c)

for l = 1, . . . , n. In this equation, the matrix Q ∈
Cn(2Nu+1)×(2Nu+1) is divided into n square blocks and
due to the diagonal structure of C(T ), only the l-th
block of Q influences the l-th row of the least-squares
problem (43). A matrix C that was determined for
comparison purposes using this least-squares problem
under knowledge of the truemonodromymatrix will be
referred to as Ctrue below. This least-squares problem
has n equations for Nu independent unknowns. There-
fore, to enable an optimal solution, it is advisable that
the frequency order Nu is chosen to be at least n.

As the aim of this method is approximating the true
monodromy matrix, which is unknown in applications,
the least-squares problem (44) cannot be utilized in
practice.As a practicable condition, thematrixC(t) can
be chosen such that it optimally satisfies the variational
equation (14). The residual of (14) for the approximated
fundamental matrix is

R(t, ĉall) = (C(t)eAtW)· − J(t)C(t)eAtW (45)

and a cost function can be defined via

Lvar(ĉall) =
∫ T

0

∥∥R(t, ĉall)
∥∥2 dt. (46)

Similarly to the block-based approach for the least-
squares problem, this function with a matrix argument
can be transformed into a quadratic cost function

Lvar(ĉall) = ĉall

(∫ T

0
�(t)dt

)
ĉ∗
all (47)

with the vector ĉall as argument. This transformation of
‖R(t)‖ is detailed below.

By product rule, the total time derivative of the first
summand in (45) yields

(C(t)Q(t))· = C(t) (A + D(t))Q(t) =: C(t)L(t) ,

(48)

where D = In×n ⊗ diag(i Nuωt, . . . ,−i Nuωt) (cf.
(40)) is the matrix that satisfies Ċ = CD and Q̇ = AQ
follows from its definition. Below, the dependency on t
in the matrices is omitted for the sake of brevity. Sub-
stitution of Q, L into (45) yields

R = CL − JCQ. (49)

123



Sorting-free Hill-based stability analysis 8453

To exploit the diagonal structure (38) ofC, thematrices
Q,L ∈ C

n(2Nu+1)×(2Nu+1) are segmented into stacks
of column vectors of length 2Nu + 1 via

L =:
⎛

⎜
⎝

L11 . . . L1n
...

. . .
...

Ln1 . . . Lnn

⎞

⎟
⎠ , (50)

and Q analogously. For the first summand in (49), this
gives

CL =
⎛

⎜
⎝

c1L11 . . . c1L1n
...

. . .
...

cnLn1 . . . cnLnn

⎞

⎟
⎠ , (51)

where each of the indicated entries is a time-dependent
scalar inC. Similarly, the second summand can be sep-
arated into its entries to yield

JCQ =
⎛

⎜
⎝

∑n
l=1 cl J1lQl1 . . .

∑n
l=1 cl J1lQln

...
. . .

...∑n
l=1 cl JnlQl1 . . .

∑n
l=1 cl JnlQln

⎞

⎟
⎠ . (52)

The still time-dependent cl(t) is generated from the
constant coefficients ĉl via the diagonal matrix V(t)
as in (40). Collecting all unknown coefficients into one
large vector ĉall, the (i, j)-th scalar entry of thematrices
in (49) can be expressed by

(CL)i j = ĉall
(
0 . . . LT

i jV . . . 0
)T =: ĉallli j (53a)

(JCQ)i j = ĉall
(
Ji1QT

1 jV, . . . , JinQT
njV
)T

=: ĉallqi j . (53b)

In (53a), only the i-th block is nonzero.
If the Frobenius norm is used for the residual, all

absolute squared values of the entries of thematrix (49)
are summed, and this gives the expression

‖R‖2 =
n∑

i, j=1

ĉall
(
li j − qi j

) (
li j − qi j

)∗ ĉ∗
all . (54)

Finally, noting that ĉall is not time-dependent, it can be
pulled outside the integral to yield the quadratic cost
function (47) with

�(t) =
n∑

i, j=1

(
li j − qi j

) (
li j − qi j

)∗
. (55)

Therefore, the search for the best choice for the pro-
jection matrix Cvar can be formulated as a quadratic
program

min
ĉall

ĉall

(∫ T

0
�(t)dt

)
ĉ∗
all (56a)

subject to ĉallW = (1 . . . 1
)

, (56b)

where the equality constraint encodes the normaliza-
tion condition (41a). The condition (41b) is not explic-
itly stated in the quadratic program. This is because
minimizers of (56) always fulfill this condition. From
an arbitrary candidate ĉall, one can construct a symmet-
ric candidate ĉsym which fulfills (41b) and for whichR
has the same real part and zero imaginary part, meaning
that ‖R‖2 can not be larger for the symmetric candi-
date than for the non-symmetric one. The proof of this
is sketched in Appendix B.3.

As thematrix�(t) is of sizen(2Nu+1)×n(2Nu+1),
and therefore rather large, efficient numerical determi-
nation of the integral (56a) is crucial to retain a numer-
ically efficient stability method. Because the time-
dependency of � is introduced by terms of the form
eikωt , the integral can be reformulated into a (infinite)
sum of inner products with Fourier base functions if
the power series expression for the matrix exponential
is used. This suggests that the integral could be com-
puted efficiently using FFT methods. However, due to
the matrix exponential there is no limit in frequency of
these terms and aliasing effects must be considered.

Alternatively, the integral can be determined using
numerical quadrature schemes. This is accurate, but
computationally expensive. Finally, it is also an option
to require the integrand to be zero only at specific time
instants. This reduces the quadratic program to a linear
equation system. While this approach does not mini-
mize the quadratic program (56) in an integral sense,
the resulting approximate monodromymatrix seems to
be relatively accurate in application.

5 Numerical examples

The theoretical results of the previous sections will be
illustrated in this section using some numerical exam-
ples, which will allow us to demonstrate the numeri-
cal efficiency and accuracy of the proposed projection-
based Hill method. The Mathieu equation is utilized as
a simple linear time-periodic system of the form (14) to

123



8454 F. Bayer, R. I. Leine

explicitly illustrate theKoopman lift and the projection-
based stability approach. To demonstrate the computa-
tional advantage for larger numerical orders, the lin-
earization of the vertically excited n-pendulum is then
considered as a generalization of the Mathieu equation
with arbitrary degrees of freedom.

5.1 Mathieu equation

The Mathieu equation

ẍ + (a + 2b cos 2ωt)x = 0 (57)

is an example of a Hill differential equation [47], which
has become very well-known since it results from lin-
earization of a number of applications, among them
rolling of container ships [48] and a vertically excited
pendulum [49]. After bringing the system into first-
order-form

ẋ(t) = J(t)x(t) =
(

0 1
−a − 2b cos(2ωt) 0

)
x(t) (58)

with x = [x, ẋ]T, (58) is a linear periodic time-varying
homogeneous system of the form (12) with system
period T̃ = π

ω
, meaning that it is suitable to explore

the stability methods of Sect. 4. Necessarily, the sys-
tem is also T -periodic with T = 2T̃ = 2π

ω
. It is known

that the only parameter combinations of (58) which
admit non-trivial T -periodic solutions are located at
the stability boundaries [39,50]. The stability bound-
aries of the Mathieu equation can therefore be identi-
fied using the HBM or the shooting method. Below, the
base frequency ω (i.e., including the first subharmon-
ics of the system) has been chosen for investigation
using the frequency-based methods. This is due to two
reasons: First, as the stability boundaries may admit
T -periodic or T

2 -periodic solutions, this choice allows
to identify all stability boundaries using one unified
HBM approach without further distinction. Second, in
the projection-based Hill method, this choice of base
frequency allows to better showcase the influence of
the projection matrix.

The equivalence of the Koopman lift matrices, the
Hill matrix and the HBM are explicitly verified for the
smallest possible frequency order Nu = 1. First, the
Hill matrix is constructed in the standard way. The
Fourier decomposition of J(t) can be read directly
from (58). The nonzero Fourier coefficients with base
frequency ω are

J0 =
(

0 1
−a 0

)
(59)

J2 = J−2 =
(

0 0
−b 0

)
(60)

and all others vanish. By construction using (18), the
Hill matrix of frequency order 1 reads

H =

⎛

⎜⎜⎜⎜⎜⎜
⎝

iω 1 0 0 0 0
−a iω 0 0 −b 0
0 0 0 1 0 0
0 0 −a 0 0 0
0 0 0 0 −iω 1

−b 0 0 0 −a −iω

⎞

⎟⎟⎟⎟⎟⎟
⎠

. (61)

The block-Toeplitz-like structure with additional diag-
onal elements is visible in (61). Furthermore, it is obvi-
ous through the parameter b that terms of harmonic 2
influence the Hill matrix, even if the frequency order
was chosen to be Nu = 1 < 2.

Further, assume there exists a (potentially non-
trivial) periodic solution xp = ∑ j p jei jωt . Because
the original dynamics (58) is already linear, this also
holds for the perturbed dynamics

ẏ(t) = J(t)y(t) + J(t)xp(t) − ẋp(t). (62)

With (62) expressed explicitly in Fourier series form

ẏ1 = y2 +
∑

j

(p j,2 − i jωp j,1)e
i jωt (63a)

ẏ2 = −ay1 − by1e
−2iωt − by1e

2iωt

+
∑

j

[−ap j,1−b(p j+2,1+ p j−2,1)−i jωp j,2
]
ei jωt ,

(63b)

the explicit determination of the time derivatives of the
Koopmanbasis functions (26) for the perturbed dynam-
ics yields with the index shift j̃ = j + k

d

dt
[y1eikωt ] = ẏ1e

ikωt + ikωy1e
ikωt

= y2e
ikωt + ikωy1e

ikωt

+
∑

j̃

[
p
( j̃−k,2) − i( j̃−k)ωp

( j̃−k,1)
]
ei j̃ωt

(64a)

for the first state and

123



Sorting-free Hill-based stability analysis 8455

Fig. 3 Floquet multiplier
locations for ω = 1 and
b = 1.2, determined using a
Hill matrix of order 4

d

dt
[y2eikωt ] = − ay1e

ikωt + ikωy2e
ikωt

− by1e
i(k−2)ωt − by1e

i(k+2)ωt

+
∑

j̃

[
−ap

( j̃−k),1−iω( j̃−k)p
( j̃−k),2

]
ei j̃ωt

− b
∑

j̃

[
p
( j̃+2−k),1+ p

( j̃−2−k),1

]
ei j̃ωt

(64b)

for the second state. The coefficients for the basis
functions can be identified in (64) by inspection. For
the state-dependent basis functions with Nu = 1, this
yields

A =

⎛

⎜⎜⎜⎜⎜⎜⎜⎜
⎝

−iω 0 0 1 0 0
0 0 0 0 1 0
0 0 iω 0 0 1

−a 0 −b −iω 0 0
0 −a 0 0 0 0

−b 0 −a 0 0 iω

⎞

⎟⎟⎟⎟⎟⎟⎟⎟
⎠

. (65)

Theblocks ofAhavebeenvisually separated to indicate
the dependence on the individual states. The selection
matrix

U =

⎛

⎜⎜⎜⎜⎜⎜
⎝

0 0 1 0 0 0
0 0 0 0 0 1
0 1 0 0 0 0
0 0 0 0 1 0
1 0 0 0 0 0
0 0 0 1 0 0

⎞

⎟⎟⎟⎟⎟⎟
⎠

(66)

necessary to reorder the basis�z,lin as required by The-
orem 3 can be determined by inspection. With this, it

is indeed easy to verify that H = UAUT holds for
Eqs. (61), (65), (66).

The state-independent terms in (64) can be collected
into the B matrix. For the considered case Nu = 1,
only the values k, j ∈ {−1, 0, 1} are considered. As
expected from Theorem 2, this yields two stacked
Toeplitz matrices BT = [BT

1 ,BT
2 ] with

B1=
⎛

⎝
p0,2 p1,2−iωp1,1 p2,2−2iωp2,1

p−1,2+iωp−1,1 p0,2 p1,2−iωp1,1
p−2,2+2iωp−2,1 p−1,2+iωp−1,1 p0,2

⎞

⎠

(67)

and B2 omitted for the sake of brevity. The expected
Toeplitz structure and the HBM equations up to order
2Nu are clearly visible. It is also easy to see that B1

holds the HBM equations up to order 2 for the first
(time-independent) row of (58).

The accuracy of a stability traverse is investigated
in Fig. 3 for the Mathieu equation with ω = 1 and
b = 1.21. The parameter a takes the values −0.367
and−0.3673, respectively. Between these values, a sta-
bility change from unstable to stable occurs, where
the pair of Floquet multipliers meet at 1, and con-
tinue on the unit circle as a complex conjugated pair.
The Floquet multipliers labeled as true in Fig. 3 were
determined using the time-integration method with a
high relative and absolute tolerance of 10−14 each. In
contrast, all Hill-matrix-based methods rely on a Hill
matrix of relatively low frequency order Nu = 4, such
that the differences in performance are more visible.

From the figures, it is apparent that the projection-
basedmethodwith an optimized projectionmatrixCvar

(i.e., by solving (56)) can be able to track the true Flo-
quet multipliers more closely than any Floquet mul-
tiplier candidates obtained directly from the eigenval-

123



8456 F. Bayer, R. I. Leine

Fig. 4 Ince–Strutt
diagrams created with
various methods based on
the Hill matrix of order
Nu = 4. Color indicates
absolute value of largest
considered Floquet
multiplier candidate. Green
lines indicate true stability
boundaries

ues of the Hill matrix allow. Even more, both consid-
ered sorting methods choose another pair of eigenval-
ues than the closest one.

The stability regions of the Mathieu equation are
often visualized in a so-called Ince–Strutt diagram
[51]. Figure 4 showcases the properties of the different
approaches for drawing this stabilitymap using theHill
matrix of a relatively low frequency order Nu = 4. The
color indicates the absolute value of the largest Floquet
multiplier (in magnitude). Dark blue colors indicate a
maximum value of exactly 1, i.e., stable regions. White
color indicates a maximum value larger than 1, i.e.,
unstable regions.

For Hill equations, it is known that the product of
both Floquet multipliers must always equal to 1 [39],
i.e., the largest Floquet multiplier may never admit an
absolute value smaller than one. In the diagrams, red
and purple regions correspond to parameter combina-
tions where the largest identified Floquet multiplier has
an absolute value smaller than one, meaning that the
identified Floquet multipliers can certainly not be the
true ones. In every figure, the true stability boundaries
are given in green by the solution of an accurate shoot-
ing method.

If all eigenvalues of the Hill matrix are considered
for stability (naive approach, Fig. 4a), there are large
regions that are wrongly classified as unstable, while
none of the unstable regions are wrongly classified as
stable. This is expected since in the unstable case, the

best approximations will always be among the consid-
ered eigenvalues, but in the stable case the additional
candidates add additional possibilities of asserted insta-
bility. In contrast, the classical Hill method with sort-
ing procedures, which is the current state-of-the-art,
classifies most of the stable regions correctly. This is
visualized in Fig. 4b for the imaginary-part-based cri-
terion [23,30] and in Fig. 4c for the symmetry-based
criterion [27,28]. However, in the instability tongue
that is separated by non-trivial solutions of period T ,
red artifacts that are classified as stable are visible in
both approaches. These correspond to cases where the
sorting algorithm did not choose the correct Floquet
multipliers (of which one would be stable, i.e., real-
valued and< 1, and the other unstable, i.e., real-valued
and > 1), but rather chose two instances of the stable
class, both with real parts < 1. With the naive projec-
tion matrix, the novel projection-based approach pre-
servesmost of the stability regions,while not exhibiting
any stable artifacts. However, at the edges of the sta-
ble regions, in particular around a = 0.5 and b = 1,
slightly larger values of the magnitude are visible. The
stable regions are slightly underestimated. A very sim-
ilar behavior can also be observed with the optimized
projectionmatrixCvar. The examplewith the naive pro-
jection matrix shows that the presented method can be
very accurate, even in its simplest form which is easy
to implement since only the matrix equation (4.1) is
needed.

123



Sorting-free Hill-based stability analysis 8457

Fig. 5 The vertically excited multiple pendulum

5.2 Vertically excited multiple pendulum

TheMathieu equation considered in Sect. 5.1 can result
from linearization of a vertically excited mathematical
pendulum, also called the Kapitza pendulum [52]. As
a scalable generalization for arbitrary degrees of free-
dom, the linearized dynamics of a vertically excited
multiple pendulum is considered. A sketch of the con-
sidered mechanical system is given in Fig. 5. The pen-
dulum consists of np joints, each of mass m, with
viscous absolute damping d̂ , linked by np rods of

length l. The minimal coordinates θ = (θ1, . . . , θnp
)T

are the absolute angles of the individual joints. The sus-
pension point of the pendulum moves vertically with
y0(t) = ŷ0 cos(2ωt). Gravitation acts in the vertical
direction.

The equations of motion for the vertically excited
pendulum can be derived similar to [53]. We define the
auxiliary vectors

s(θ)T:=(np sin θ1, (np−1) sin θ2, . . . , sin θnp
)

(68a)

c(θ)T:=(np cos θ1, (np−1) cos θ2, . . . , cos θnp
)

(68b)

θ̇
2 :=

(
θ̇21 , . . . , θ̇2np

)T
(68c)

as well as matrices S(θ), C(θ) ∈ R
np×np with

Si j (θ) := [np + 1 − max(i, j)
]
sin(θi − θ j ) (68d)

Ci j (θ) := [np + 1 − max(i, j)
]
cos(θi − θ j ). (68e)

Dropping the arguments for the sake of brevity,
the potential energy V (θ, t) and the kinetic energy
T (θ , θ̇ , t) are given by

V = −npmgy0 − mgl
(
1 . . . 1

)
c (69a)

T = 1

2
m
[
np ẏ

2
0 + l2θ̇

TCθ̇ − 2l ẏ0sTθ̇
]
. (69b)

Using the Lagrange equations of the second kind (see,
e.g., [54, p. 76]), the equations of motion are given by

Cθ̈ + d̂

ml2
θ̇ + S θ̇

2 +
(

− ÿ0(t)

l
+ g

l

)
s(θ) = 0 (70)

after some algebra. Introducing the abbreviations

M := C(0) (71a)

D := diag(np, . . . , 1) (71b)

a := g/ l (71c)

2b := 4ŷ0ω
2/ l (71d)

d := d̂

ml2
(71e)

and linearizing around the origin, the first-order lin-
earized dynamics of the vertically excited pendulum
is

(
I 0
0 M

)(
θ̇

θ̈

)
=
(

0 I
−[a + 2b cos(2ωt)]D −dI

)(
θ

θ̇

)
,

(72)

which, after inversion of the left matrix, is of the
form ẏ = J(t)y that can be analyzed using the pre-
sented methods. The total number of system states is
n = 2np. The relation between the classical Mathieu
equation (58) and the linearized vertically excited mul-
tiple pendulum is visible from (72). The total num-
ber n of states as considered throughout this paper is
n = 2np.

To analyze the convergence behavior of our pro-
posed method, the accuracy of the Floquet multi-
pliers of the equilibrium for one specific parame-
ter set (a, b, d) is evaluated for various truncation
orders Nu and using the various approaches discussed
in this paper. As a basis for comparison, the “true”
Floquet multipliers are determined by integrating the
variational equation (14) using the MATLAB func-
tion ode45 with absolute and relative tolerances of

123



8458 F. Bayer, R. I. Leine

Fig. 6 Accuracy of Floquet multipliers against Floquet multi-
pliers by time-integration over frequency order for the linearized
6-pendulum with (a, b, d) = (5, 0.5, 0.2)

10−13. The total Floquet multiplier error is defined
as the square root of the sum of the squared absolute
differences between the true Floquet multipliers and
the obtained Floquet multipliers, while the latter are
ordered such that this error is minimal. More formally,
letP := {π : {1, . . . , n} → {1, . . . , n} |π bijective} be
the (finite) set of permutations, i.e., reorderings, of the
index set {1, . . . , n}. Then, the total Floquet multiplier
error is given by

εtotal = min
π∈P

√√√√
n∑

l=1

∣∣λl,true − λπ(l),cand
∣∣2. (73)

For the two standard Hill approaches, the eigenval-
ues and eigenvectors of the Hill matrix H were deter-
mined using the MATLAB procedure eig. The sort-
ing procedure based on the imaginary part as described
in [23,26,30] then singles out the n eigenvalues with
least imaginary part in modulus, while the symmetry-
based sorting procedure singles out the n eigenval-
ues whose eigenvectors have lowest weighted mean
according to [27]. The Floquet multipliers are then
determined from the Floquet exponents (i.e., the identi-
fied eigenvalues) using (15). The corresponding errors
are depicted in Fig. 6 in dashed blue and dotted red for
the imaginary part criterion and the symmetry criterion,
respectively.

The presented novel Koopman-based method is
evaluated for two choices of the projection matrix C
(see Sect. 4.3). The Floquet multipliers are determined
as the eigenvalues of the approximate monodromy

matrix (35). The matrix exponential in (35) is evalu-
ated directly with its action onto thematrixW using the
MATLAB function exmpv [46,55]. The dash-dotted
yellow line shows the accuracy of the Floquet mul-
tipliers obtained from the monodromy approximation
using the naive projection choice C0 (31), while the
solid purple line indicates the accuracy obtained by a
more informed projection matrix

C1 := In×n ⊗ (1,−1, 1,−1, . . . ,−1, 1
)

, (74)

which was obtained by inspection after running the
optimization (56) for a few test cases. From the fig-
ure, it is visible that all considered approaches eventu-
ally converge to the true Floquet multipliers. The final
error of order 10−12 can be attributed to inaccuracies
of the numerical integration. The imaginary-part-based
criterion, while being the only one with a rigorous con-
vergence proof for n → ∞, performs highly inaccu-
rate for smaller Nu. Moreover, the choice of projec-
tion matrix significantly influences the performance of
the novel projection-based approach. While the naive
choice of projection matrix does converge toward the
correct value, the optimization-based projectionmatrix
exhibits a significantly better convergence rate which
can compete with the symmetry-based criterion or out-
performs it.

The data set of Fig. 6 is again plotted in Fig. 7,
but against the computation time instead of the fre-
quency order, giving a work-precision diagram. The
strength of the projection-basedmethods for higher fre-
quency orders is apparent in this figure. While conver-
gence of the imaginary-part-basedmethodwith respect
to the frequency order Nu is better than that of the
new projection-based method for the naive projec-
tion choice C0, the accuracy of the projection-based
method is higher than that of the imaginary-part-based
method for any computation time. The reason for this
is that the eigenproblem becomes more expensive than
the matrix exponential for higher matrix dimensions,
i.e., for higher Nu. The same property can also be
observed regarding the symmetry-based approach and
the novel projection-based approach with optimized
projectionmatrixC1. The projection-based approach is
able to achieve maximum accuracy already at approx.
0.1s compared to approx. 0.25s for the symmetry-
based approach. This is a ratio of 2

5 , which is con-
siderably lower than the corresponding ratio 15

18 of

123



Sorting-free Hill-based stability analysis 8459

Fig. 7 Accuracy of Floquet multipliers over computation time
for the linearized 6-pendulumwith (a, b, d) = (5, 0.5, 0.2). The
data points correspond to values of Nu in the range between 1
and 40

frequency orders. If the system dimension increases,
this efficiency difference increases as well as the size
n(2Nu + 1) × n(2Nu + 1) of the Hill matrix is influ-
enced by a product of n and Nu. Therefore, the superior
efficiency of the projection-based approaches against
the standard approaches is even more prominent. Fig-
ure8 again compares computation time against accu-
racy for the various approaches, but for a pendulum
with 15 instead of 6 links. Since the number of states
is increased, the computation time for all the consid-
ered approaches is larger than in Fig. 7. However, the
increase in number of states impacts the standard Hill
methods which rely on solving the large eigenvalue
problem more than the projection-based method. This
is because of the key feature of the newmethod: a mod-
erate growth of computational cost as sparsity can be
exploited in the action of thematrix exponential. In fact,
in Fig. 8, the break-even point between the symmetry-
basedmethod and the projection-basedmethodwith the
optimized projection matrix already occurs at Nu = 6.
The imaginary-part-basedHillmethod fails to converge
within the depicted computation time interval.

6 Conclusions

In this paper, a relationship between the HBM, the Hill
method and the Koopman lift with specific basis func-
tions has been addressed. It has been shown that theHill
matrix is the system matrix of a linear time-invariant
dynamical system of higher order which approximates

Fig. 8 Accuracy of Floquet multipliers over computation time
for the linearized 15-pendulum with (a, b, d) = (5, 0.5, 0.2).
The data points correspond to values of Nu in the range between
1 and 40

the perturbed dynamics. This relationship has been
used to derive a novel stability method based on pro-
jecting the Hill matrix down to original system size,
instead of computing all its eigenvalues.

The resulting stability criterion in its naive form
is remarkably simple to implement. It only needs a
matrix exponential and two projection matrices, which
are the same for all systems and easy to construct
(see (31), (33)). It has been shown in the examples of
Sect. 5 that this simple method can compete with the
state-of-the-art algorithms in terms of Floquet multi-
plier accuracy over computational effort.

To further improve the accuracy of the proposed
method, an optimization-based approach for the choice
of the projection matrix through a quadratic program
has been presented in Sect. 4.3. While the determina-
tion of the matrix integral in (56) is more computa-
tionally expensive, it has been shown in Sect. 5 that
this optimized projection matrix allows for even better
accuracy.

There are many topics for further research expand-
ing from the results of this work: The approach itself
offers opportunities for additional development both in
its theoretical as well as for its computational proper-
ties. Further, there are many possible interesting appli-
cations of the approach to technical systems, and exten-
sions to wider classes of systems would be valuable.

With respect to theoretical properties, there is cur-
rently no rigorous convergence proof for the novel
approach for any choice of projection matrix. Such a

123



8460 F. Bayer, R. I. Leine

convergence proof would aid in an a priori determina-
tion which method to employ for the best results.

Considering computational properties, while the
presented approach does reduce computational effort
compared to the state-of-the-art approaches, some bot-
tlenecks remain. In particular, the efficient sparsity-
promoting determination of the matrix exponential
with projection from both sides as well as the integral
in the optimization problem (56) would benefit from
more efficient computation techniques. As the latter is
constructed by products of the Fourier coefficients of J
in A and additional Fourier terms, there could exist an
operation to obtain this integral directly from the (FFT)
Fourier series of J(t).

From an application point of view, going beyond the
simple and mainly academical examples presented in
this work, it would also be expedient to apply the novel
stabilitymethod to systems ofmore practical relevance.
To achieve this goal, it would be beneficial to integrate
the novel method within an established continuation
framework. The authors believe that the MANLAB
continuation framework [25,28] would be especially
suitable because the Hill matrix can be constructed in
this framework without much additional effort. Large
systems that were already analyzed extensively using
the MANLAB framework (e.g., [56]) could benefit
from the stability insights of our proposed method,
while simultaneously serving as practical examples of
the performance of the method.

Regarding possible extensions, while some relations
were established in Theorem 1 for polynomial basis
functions of the Koopman lift, most sections of this
work rely on basis functions that are linear in the state.
Some recent results [57,58] make statements about
return time and period for autonomous systems based
on the Carleman linearization. As our (full) basis pro-
vides a natural extension of the Carleman linearization,
these results could possibly be integrated into the pre-
sented framework.

Finally, it would be of practical value to extend the
method to a wider class of systems, for example, possi-
bly including delay differential equations (DDEs) with
delayed coordinates in the basis functions.

Author contributions Both authors contributed to the concep-
tualization, design and methodology. The formal analysis was
performed by Fabia Bayer. The original draft of the manuscript
was prepared by Fabia Bayer and both authors commented on
previous versions of the manuscript. Remco Leine was respon-

sible for funding acquisition and supervision. Both authors read
and approved the final manuscript.

Funding Open Access funding enabled and organized by Pro-
jekt DEAL. The authors declare that no funds, grants, or other
support was received during the preparation of this manuscript.

Data availability All numerical data in this work have been
generated from the considered system Eqs. (58), (72), with
the approaches described in this work and the cited works,
using MATLAB standard methods and the MATLAB function
expmv [55]. Therefore, it is possible to completely reproduce
the data from the information given in this work, except for the
computation times in the work-precision diagrams. These com-
putation times have been evaluated on a Windows PC, and more
detailed computation time data can be obtained from the authors
upon request.

Declarations

Conflict of interest The authors have no relevant financial or
non-financial interests to disclose.

Open Access This article is licensed under a Creative Com-
mons Attribution 4.0 International License, which permits use,
sharing, adaptation, distribution and reproduction in anymedium
or format, as long as you give appropriate credit to the original
author(s) and the source, provide a link to the Creative Com-
mons licence, and indicate if changes were made. The images or
other third partymaterial in this article are included in the article’s
Creative Commons licence, unless indicated otherwise in a credit
line to thematerial. If material is not included in the article’s Cre-
ative Commons licence and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need
to obtain permission directly from the copyright holder. To view
a copy of this licence, visit http://creativecommons.org/licenses/
by/4.0/.

Appendix A: Multi-index calculation rules

For convenience, the standard multi-index calculation
rules (see, for example, [59, p. 319]) are revisited. A
multi-index is a tuple β ∈ N

d
0 with nonnegative integer

entries. In particular, a multi-index β = (β1, . . . , βd)

has the 1-norm

‖β‖ =
d∑

i=1

βi (A1)

and a factorial

β! =
d∏

i=1

βi ! (A2)

123

http://creativecommons.org/licenses/by/4.0/
http://creativecommons.org/licenses/by/4.0/


Sorting-free Hill-based stability analysis 8461

given by the product of the factorials of all its com-
ponents. For a vector x ∈ C

d , exponentiation is given
by

xβ =
d∏

i=1

xβi
i (A3)

and for a function f : Rd → C, the β-th derivative is
given by

∂β f

∂xβ
= ∂‖β‖ f

∂xβ1
1 . . . ∂xβd

d

. (A4)

It is known in combinatorics that the number of multi-
indices β ∈ N

d
0 with ‖β‖ ≤ N is given by

(N+d
d

)
[60].

This is also the number of unique monomials of the
form (A3) with degree ≤ N , including x0.

Appendix B: Proofs for the theorems

In this section, the theorems of Sect. 3.2 are proven for-
mally. The slightly unusual matrix-valued inner prod-
uct notation (1)–(2c) is heavily utilized throughout the
proofs for the sake of brevity.

B.1 Theorems 1 and 2

Theorem 2 will be proven first, and Theorem 1 will
result as a consequence from that proof. Recall Theo-
rem 2:

Theorem 2 Let ż = A(p)z + B(p)u be the lifted
dynamics of system (6)with linear basis functions�z,lin

of frequency order Nu that are sorted as in (22), eval-
uated for the perturbed system around an unknown
periodic ansatz of the form (10) up to frequency order
at least NHBM = 2Nu. Then, the matrix B(p) ∈
C
n(2Nu+1)×(2Nu+1) consists of n stacked Toeplitz matri-

ces. The l-th Toeplitz matrix Bl contains as entries
(ignoring duplicates) precisely the 4Nu + 1 residuals
rk,l(p), |k| ≤ 2Nu that follow from the HBM w.r.t the
l-th state.

If B(p) = 0, then all these residuals of the HBM
vanish. Conversely, if p solves the HBM equations
rk(p) = 0, |k| ≤ 2Nu, then it holds that B(p) = 0.

Proof For the sake of brevity, bounds of sums are
omitted in this proof if they take the values of ±∞.
Recall that the HBM equations can be expressed as
the harmonic residual of the dynamics (29). Let (26)
be a linear basis for the Koopman lift as in Theo-
rem 2 with z(t) ≈ �z,lin(y(t), t). With the order-
ing (26), the first 2Nu + 1 entries of z approxi-
mate all basis functions connected to y1, the next
2Nu + 1 entries approximate all basis functions con-
nected to y2 et cetera. With this separation, divide the
matrix B(p) ∈ C

n(2Nu+1)×(2Nu+1) into n square matri-
ces {Bl}nl=1 such that Bl ∈ C

(2Nu+1)×(2Nu+1) influ-
ences the derivative of yl . Below, the argument is made
for one Bl and B can be stacked together in the end.

Consider now the time derivative along the system
flow of a Koopman basis function corresponding to the
l-th state and frequency j , evaluated on the perturbed
system:

ψ j,l(y, t) = yle
i jωt (B5)

= xle
i jωt −

∑

k

pk,le
i( j+k)ωt (B6)

ψ̇ j,l(y, t) = ẋle
i jωt + i jωxle

i jωt

−
∑

k

i( j + k)ωpk,le
i( j+k)ωt (B7)

= fl(y + xp, t)ei jωt + i jωyle
i jωt

−
∑

k

ikωpk,le
i( j+k)ωt . (B8)

The inner product
〈
ψ̇ j,l ,�z,lin

〉
is a row of A and〈

ψ̇ j,l ,u
〉
is a row of Bl . As the inner product for A

contains all summands that are dependent on y and
the inner product for B collects all parts that are not
dependent on y, we can set y = 0 in this proof about
properties of B. With this, it results

ψ̇ j,l = fl(xp(t), t)ei jωt −
(
∑

k

ikωpk,le
ikωt

)

ei jωt

(B9)

= rl(t)e
i jωt (B10)

=
∑

k

rk,le
i(k+ j)ωt (B11)

=
∑

k̃

rk̃− j,le
i k̃ωt . (B12)

123



8462 F. Bayer, R. I. Leine

The inner product with u is just a projection onto the
Fourier basis functions between ±Nu, and therefore,
the rowofBl that corresponds to frequency j is givenby

〈
ψ̇ j,l ,u

〉 = (r−Nu− j,l , . . . , rNu− j,l). (B13)

Comparing two successive rows, i.e., the j-th and the
( j + 1)-th row of Bl , it is obvious from (B13) that they
contain the same entries, but shifted to the right by one
element. Adequately, going from row j to row j+1, the
entry r−Nu−( j+1),l is added on the left, while the entry
rNu− j,l is pushed out on the right. The matrix Bl thus
has a Toeplitz structure, where all residuals between
r−2Nu,l (in the bottom left entry of Bl with j = Nu)
and r2Nu,l (in the top right entry with j = −Nu) occur,
yielding 4Nu + 1 distinct entries in total.

In particular, if the coefficients p describe the HBM
solution for order 2Nu, then all these residuals vanish
and it holds that B(p) = 0. ��

The arguments of the previous proof can be reused
for Theorem 1 with slightly modified assumptions.
Recall the theorem:

Theorem 1 Let ż = A(p)z + B(p)u be the lifted
dynamics of frequency order NHBM of system (6)
around an unknown periodic ansatz of the form (10).
The NHBM-th order HBM equations (11), i.e., rk =0,
|k| ≤ NHBM, are given byCzB(p) = 0, whereCz is the
constant selection matrix that fulfills y = Cz�z(y, t)
for all t.

Proof In contrast to the previously proven Theorem 2,
in this theorem the basis functions are not restricted to
be linear in the state. Due to sesquilinearity of the inner
product (2b), however, the selection matrix Cz can be
pulled into the expression for B(p) via

CzB(p) = Cz
〈
�̇z,u

〉
(B14)

= 〈Cz�̇z,u
〉

(B15)

=
〈⎛

⎜
⎝

ψ̇0,1
...

ψ̇0,n

⎞

⎟
⎠ ,u

〉

. (B16)

These expressions can again be evaluated row-wise
with (B13) and since j = 0 for all entries of (B14),
the condition CzB(p) = 0 is equivalent to rk = 0 for
all |k| ≤ NHBM, which are exactly the HBM equations.

��

B.2 Theorem 3

Similarly to the theoremsprovenpreviously, it is central
in the proof for Theorem 3 to evaluate the inner product
explicitly. The theorem is restated for convenience.

Theorem 3 Let ż = A(p)z + B(p)u be the lifted
dynamics around a periodic solution of system (6) with
linear basis functions�z,lin of frequency order Nu that
are ordered as in (22). Then, the Hill matrix H, trun-
cated to frequency order Nu, for the periodic solution
parameterized by p results from the matrix A(p) by
the similarity transformH = UA(p)UT, where U is an
orthogonal permutation matrix that satisfiesU�z,lin =
(yT ei Nuωt , . . . , yT e−i Nuωt )T.

Proof The matrixU is derived from the identity matrix
by reordering its rows, and it reorders the entries of
�z,lin such that the resulting vector has entries descend-
ing in frequency, where all states corresponding to the
same frequency are collected together. Denote �̃z :=
U�z,lin =: ((�̃z)

T−Nu
, . . . , (�̃z)

T
Nu

)T and separate it

into 2Nu+1 blocks of length n with (�̃z)k := ye−ikωt .
First, we show that the Koopman lift performed with
basis �̃z identically yields H as its system matrix Ã,
and afterward we demonstrate how performing the lift
with �z,lin instead of �̃z demands the additional sim-
ilarity transform.

As in the proof for Theorem 2, we start by deter-
mining the time derivative of the basis functions. For
the vector of time derivatives, by the chain rule it holds
that

˙̃
�z = ∂�̃z

∂y
f̃ + ∂�̃z

∂t
, (B17)

where f̃ is the nonlinear perturbed dynamics (27a).
Evaluating the individual summands of (B17) yields

∂�̃z

∂t

∣∣∣∣∣
y,t

=
⎛

⎜
⎝

i Nuω yei Nuωt

...

−i Nuω ye−i Nuωt

⎞

⎟
⎠ (B18)

∂�̃z

∂y

∣∣∣∣∣
y,t

=
⎛

⎜
⎝

Iei Nuωt

...

Ie−i Nuωt

⎞

⎟
⎠ (B19)

f̃(y, t) =
∞∑

j=−∞
J je

i jωty + O(‖y‖2) + b(t). (B20)

123



Sorting-free Hill-based stability analysis 8463

The terms of order O(‖y‖2) in (B20) can be dropped
because terms of higher polynomial order are not rep-
resented in the chosen basis �z,lin. The terms in b(t)
appear if the candidate ansatz is not a periodic solu-
tion. They are purely time-dependent and will remain
so after multiplication by (B19). Hence they are col-
lected in the Bmatrix. The column vector in (B18) and
the matrix in (B19) can be decomposed into 2Nu + 1
blocks with n rows and 1 or n columns as indicated
above, each block corresponding to a specific fre-
quency. These blocks are labeled in ascending order
by k, −Nu ≤ k ≤ Nu, such that the k-th block corre-
sponds to the term e−ikωt . Equations (B18)–(B20) are
substituted into the k-th block in (B17) to yield

( ˙̃
�z

)

k
= e−ikωt

⎛

⎝
∞∑

j=−∞
J je

i jωty

⎞

⎠− ikωye−ikωt

(B21)

=
∞∑

j=−∞
J je

i( j−k)ωty − ikωye−ikωt (B22)

= −ikω
(
�̃z

)

k
+

∞∑

j=−∞
J j

(
�̃z

)

k− j
, (B23)

where the purely time-dependent terms and the terms
of order O(‖y‖2) have been dropped for the sake of
legibility. To obtain the k-th block of the Koopman lift

matrix Ã, the inner product
〈( ˙̃

�z

)

k
, �̃z

〉
is considered.

Due to the sesquilinearity of the inner product, both
summands in (B23) can be computed separately. It is
easy to see that the first summand yields a sparsematrix
where only the diagonal in the k-th column block is
nonzero. With an index shift

∞∑

j=−∞
J j

(
�̃z

)

k− j
=

∞∑

j̃=−∞
Jk− j̃

(
�̃z

)

j̃
, (B24)

the inner product with the second summand of (B23)
yields the matrix

(
Jk−(−Nu) . . . Jk−Nu

)
. (B25)

Collecting all row blocks and both summands together,
the Koopman lift matrix à for the re-sorted basis is

given by

à =
〈 ˙̃
�z, �̃z

〉
=
⎛

⎜
⎝

J0 + i Nuω . . . J−2Nu
...

. . .
...

J2Nu . . . J0 − i Nuω

⎞

⎟
⎠ .

(B26)

Comparison of (B26) to the truncated Hill matrix (18)
shows the identity à = H for the Koopman lift with
the re-sorted basis.

Finally, the similarity transform obtained by the re-
ordering is considered. The permutation matrix U is
constant and can thus be pulled out of the time deriva-
tive.With the sesquilinearity of thematrix-valued inner
product (2b), (2c), it then follows

H = Ã =
〈 ˙̃
�z,�z

〉
(B27)

= 〈U�̇z,lin,U�z,lin
〉

(B28)

= U
〈
�̇z,lin,�z,lin

〉
U∗ (B29)

= UA(p)UT , (B30)

whereA(p) is the Koopman lift with the standard basis.
In the last step, it has been exploited that U is a real-
valued matrix and thus transpose and conjugate trans-
pose coincide. ��

B.3 Symmetry properties of the optimization problem

Below, it will be shown that there always exist solu-
tions to the optimization problem (56) fulfilling (41b),
even if this constraint is not enforced explicitly. This is
summarized in the following proposition.

Proposition 4 For an arbitrary collection ĉall ∈
C
n(2Nu+1) of coefficients which admit the residual R

in the variational equation (45), there exists another
collection ĉsym ∈ C

n(2Nu+1) constructed by

ĉsym,l,k = 1

2

(
ĉl,k + ¯̂cl,−k

)
, (B31)

which fulfills (41b). The residual Rsym of the varia-
tional equation (45) for ĉsym is real for arbitrary t with
Rsym = Re(R).

The proof for this proposition is only sketched here.
Some mathematical steps are only indicated but not
performed in detail for the sake of brevity.

123



8464 F. Bayer, R. I. Leine

Proof Below, vectors d ∈ C
(2Nu+1) which fulfill (41b)

will be called symmetric. The symmetric vectors form
a vector space over R (but not over C).

First, the structure of the Koopman lift matrix A is
revisited. Applying the similarity transform of Theo-
rem 3 to the Hill matrix explicitly yields for the matrix
A the structure

A =
⎛

⎜
⎝

A11 . . . A1n
...

. . .
...

An1 . . . Ann

⎞

⎟
⎠ (B32)

Al j =
⎛

⎜
⎝

Jl j,0 . . . Jl j,2Nu
...

. . .
...

Jl j,−2Nu . . . Jl j,0

⎞

⎟
⎠

+ δl jdiag(−iωNu, . . . , iωNu). (B33)

Hence, the l j-th block of the matrix A contains the
Fourier coefficients of the l j-th entry of the system
matrix J. In particular, since J(t) admits only real
values, Jl j,k = Jl j,−k and each block Al j fulfills a
variant of the symmetry property: The −k-th row is
the complex conjugate of the k-th row, with the order
of elements reversed. A little algebra shows that for
matrices with these properties and symmetric vectors
d ∈ C

(2Nu+1), the symmetry is retained though multi-
plication: Al jd is symmetric if d is symmetric.

If a matrix P = [dl j ] is constructed block-wise from
symmetric column vectors

{
dl j
}n
l, j=1, then the matrix

AP is again constructed from symmetric column vec-
tors since all block entries ofAP can be decoupled into
sums of productsAlidi j , which each retain the symme-
try as described above. In particular, since the matrix
W in (33) is constructed in the considered fashion, the
product AW again has block-wise symmetric entries.
Iterative application of this multiplicative invariance to

Q = eATW =
∞∑

k=0

T kAkW (B34)

L = (A + D)Q (B35)

shows that the column vectors Ql j and Ll j in (50) are
also symmetric.

Next, the scalar product bTd ∈ C of a symmetric
vector d and an arbitrary vector b of appropriate size
is considered. The operator flipb reverses the order of
entries of a vector b. Element-wise evaluation shows

that

(
flip b

)T d =
Nu∑

k=−Nu

b−kdk (B36a)

=
Nu∑

k=−Nu

bkdk (B36b)

= (bTd). (B36c)

Defining ĉl,sym = 1
2

(
ĉl + flip ĉl

)
and using the above

relation yields

ĉl,symQi j = Re(ĉlQi j ) , (B37)

and analogously forLi j . AsR is constructed from sum-
mands of this form via (49)–(52), the proposition fol-
lows. ��

A purely real matrix Rsym will necessarily have a
Frobenius norm smaller or equal than the norm of a
matrix with same real part and arbitrary imaginary part.
Hence, Proposition 4 allows to easily construct a min-
imizer that satisfies all constraints in case a different
minimizer is returned in the quadratic program (56).

References

1. Mauroy,A.,Mezić, I., Susuki,Y. (eds.): TheKoopmanOper-
ator in Systems and Control: Concepts, Methodologies and
Applications. Lecture Notes in Control and Information Sci-
ences, vol. 484. Springer, Cham (2020)

2. Brunton, S.L., Budišić, M., Kaiser, E., Kutz, J.N.: Modern
Koopman theory for dynamical systems. SIAM Rev. 64(2),
229–340 (2022). https://doi.org/10.1137/21m1401243

3. Williams, M.O., Kevrekidis, I.G., Rowley, C.W.: A
data-driven approximation of the Koopman operator:
extending dynamic mode decomposition. J. Nonlinear
Sci. 25(6), 1307–1346 (2015). https://doi.org/10.1007/
s00332-015-9258-5

4. Kono, Y., Susuki, Y., Hikihara, T.: Modeling of advec-
tive heat transfer in a practical building atrium via Koop-
man mode decomposition. In: Mauroy, A., Mezić, I.,
Susuki, Y. (eds.) The Koopman Operator in Systems and
Control: Concepts, Methodologies, and Applications, pp.
481–506. Springer, Cham (2020). https://doi.org/10.1007/
978-3-030-35713-9_18

5. Korda, M., Mezić, I.: Koopman model predictive control
of nonlinear dynamical systems. In: Mauroy, A., Mezić, I.,
Susuki, Y. (eds.) The Koopman Operator in Systems and
Control: Concepts, Methodologies and Applications, pp.
235–255. Springer, Cham (2020)

123

https://doi.org/10.1137/21m1401243
https://doi.org/10.1007/s00332-015-9258-5
https://doi.org/10.1007/s00332-015-9258-5
https://doi.org/10.1007/978-3-030-35713-9_18
https://doi.org/10.1007/978-3-030-35713-9_18


Sorting-free Hill-based stability analysis 8465

6. Koopman, B.O.: Hamiltonian systems and transformation
in Hilbert space. Proc. Natl. Acad. Sci. USA 17(16577368),
315–318 (1931). https://doi.org/10.1073/pnas.17.5.315

7. Mezić, I.: Spectral properties of dynamical systems, model
reduction and decompositions. Nonlinear Dyn. 41(1), 309–
325 (2005). https://doi.org/10.1007/s11071-005-2824-x

8. Mauroy, A., Mezić, I., Moehlis, J.: Isostables, isochrons,
and Koopman spectrum for the action-angle representation
of stable fixed point dynamics. PhysicaD261, 19–30 (2013).
https://doi.org/10.1016/j.physd.2013.06.004

9. Mohr, R., Mezić, I.: Construction of eigenfunctions for
scalar-type operators via Laplace averages with connections
to the Koopman operator (2014). https://doi.org/10.48550/
arXiv.1403.6559

10. Mauroy, A., Susuki, Y., Mezić, I.: Introduction to the Koop-
man operator in dynamical systems and control theory. In:
Mauroy, A., Susuki, Y., Mezić, I. (eds.) The KoopmanOper-
ator in Systems and Control: Concepts, Methodologies and
Applications, pp. 3–33. Springer, Cham (2020)

11. Nayfeh, A.H., Mook, D.T., Holmes, P.: Nonlinear oscilla-
tions. J. Appl. Mech. 47(3), 692–692 (1980). https://doi.
org/10.1115/1.3153771

12. Bentvelsen, B., Lazarus, A.: Modal and stability analysis of
structures in periodic elastic states: application to the Ziegler
column. Nonlinear Dyn. 91(2), 1349–1370 (2018). https://
doi.org/10.1007/s11071-017-3949-4

13. Karkar, S., Vergez, C., Cochelin, B.: Os