# Nonautonomous Spectral Submanifolds for Model Reduction of Nonlinear
Mechanical Systems under Parametric Resonance

Thomas Thurnher  Institute for Mechanical Systems, ETH Zürich, Leonhardstrasse
21, 8092 Zürich, Switzerland  George Haller  Institute for Mechanical Systems,
ETH Zürich, Leonhardstrasse 21, 8092 Zürich, Switzerland  Shobhit Jain
<EMAIL_ADDRESS>Delft Institute of Applied Mathematics, TU Delft,
Mekelweg 4, 2628 CD Delft, The Netherlands

###### Abstract

We use the recent theory of Spectral Submanifolds (SSM) for model reduction of
nonlinear mechanical systems subject to parametric excitations. Specifically,
we develop expressions for higher-order nonautonomous terms in the
parameterization of SSMs and their reduced dynamics. We provide these results
both for general first-order as well as second-order mechanical systems under
periodic and quasiperiodic excitation using a multi-index based approach,
thereby optimizing memory requirements and the computational procedure. We
further provide theoretical results that simplify the SSM parametrization for
general second-order dynamical systems. More practically, we show how the
reduced dynamics on the SSM can be used to extract the resonance tongues and
the forced response around the principal resonances in parametrically excited
systems. In the case of two-dimensional SSMs, we formulate explicit
expressions for computing the steady-state response as the zero-level set of a
two-dimensional function for systems that are subject to external as well as
parametric excitation. This allows us to parallelize the computation of the
forced response over the range of excitation frequencies. We demonstrate our
results on several examples of varying complexity, including finite-element
type examples of mechanical systems. Furthermore, we provide an open-source
implementation of all these results in the software package SSMTool.

> As the complexity of engineering systems, such as MEMS devices, is growing,
> accurate reduced-order models (ROMs) are in high demand. Not only do such
> models promise quick routines for system design and analysis, but they can
> also deliver specifically desired information that might be inaccessible
> from a full system simulation. Reduction of the full system to Spectral
> Submanifolds (SSM) offers a rigorous mathematical procedure to derive such
> ROMs even for strongly nonlinear systems. In this paper, we develop a
> methodology for computing nonautonomous approximations to these SSMs up to
> arbitrary order, generalizing, improving and expanding available results in
> the literature. Using the exact ROM provided by the SSM-reduced dynamics, we
> analyze instabilities that arise due to parametric resonance and obtain the
> resulting steady states for a number of examples. Our results are also
> implemented in an open-source software library.

##  I Introduction

Parametric excitation is a widely observed phenomenon in mechanical systems
that arises from time-dependence of parameters in a dynamical system. This
time dependence is either due to external modulation or inherent nonlinear
couplings Champneys2013 . Perhaps the simplest example is the periodic
modulation of the length of a pendulum that may even destabilize its stable
fixed point. This phenomenon is known as parametric resonance and occurs if
the modulation frequency and the eigenfrequency of the pendulum assume a
rational relationship.

Due to its destabilizing effect, parametric resonance is often avoided in the
operational regime of engineering systems Sanches2012 ; Coleman1956 . However,
it is also possible to take advantage of the resulting instability, which
leads to exponential growth of energy and oscillation amplitudes. In rotor
dynamics, for instance, parametric resonance may be used to virtually increase
damping Ecker2011 . Nonlinear energy harvesters show increased performance
when driven in the regime of parametric resonance compared to systems using
conventional resonance Tehrani2015 . Micro-Electromechanical systems (MEMS)
also frequently make use of parametric resonance. Examples include low-noise
parametric amplifiers Rhoads2006_MEMS , mass-sensing devices Li2014 , signal
amplifiers Welte2013 , gyroscopes Polunin2017 ; Sharma2012 , picoprojectors,
and microscanners Frangi2017 . Parametric resonance in combination with
nonlinearites can also enable operation in a wider frequency range Rhoads2009
.

Among various types of parametric resonances, the principal parametric
resonance, where the external excitation frequency $\Omega$ and a system’s
modal frequency $\omega_{i}$ obey a relation of the form $\Omega=2\omega_{i}$,
is of pronounced practical interest. Indeed, due to the presence of damping,
higher-order parametric resonances become effectively inaccessible in most
engineering structures Verhulst2002 ; Frangi2017 ; Verhulst2009 ;
Champneys2013 . Nonetheless, high-order parametric resonances have been
observed in specially engineered systems Jia2016 .

### Forced response curves

Forced response curves (FRC) visualize the relationship between a system’s
response amplitude and its excitation frequency. FRCs provide insights into
the nonlinear characteristics of a system when compared with the linear
frequency response, especially near resonance. Various techniques have been
used to analyze parametrically excited systems. For instance, the Poincaré-
Lindstedt method assumes the dependence of the period of the solution on a
small excitation amplitude and expresses the periodic solution as an
asymptotic series for weakly excited systems Verhulst1996 . The method of
multiple scales introduces different temporal variables for different
timescales. These are then used in perturbation expansions and can be applied
to parametrically excited systems as well Warminski2005 . Other common methods
found in the parametric excitation literature are averaging Rhoads2006_MEMS ;
Rhoads2006_Microbeam ; Li2020 ; Polunin2017 ; Zaghari2016 and the analytical
construction of approximate solutions, for instance, using a harmonic series
Szabelski1995 ; Szabelski1997 ; Sorokin2015VibrationModulation ;
Aghamohammadi2022DynamicExcitations or direct numerical simulation
Warminski2003 .

For a systematic computation of FRCs, numerical continuation packages, such as
AUTO AUTO , coco Dankowicz2013 , or the harmonic balance-based NLvibKrack2019
, may be used. These methods are based on collocation or spectral
discretization and typically suffer from the curse of dimensionality. Hence,
their direct applications to high-dimensional finite element models are
generally unfeasible Jain2021 .

### Stability Diagrams

Establishing the stability type of trivial fixed points is the first important
step in the study of parametrically excited systems. The regions of
instability and stability of these fixed points are plotted in an excitation
amplitude vs. frequency diagram, often referred to as Strutt or Ince Strutt
diagram. In these stability diagrams, the regions of instability are often
called Floquet tongues that extend in the direction of decreasing forcing
amplitude as the system’s damping decreases. The simplest technique for
constructing stability diagrams is numerical time integration over a grid of
amplitude and frequency parameters Welte2013 . More sophisticated approaches
employ Floquet theory and range from the method of infinte determinants in the
conservative case Hansen1985 to spectral methods for obtaining the eigenvalues
of the monodromy matrix Lindh1970 ; Hansen1985 ; Welte2013 . In the extended
phase space, a bifurcation of the trivial periodic orbits occurs along the
stability boundary, which enables its computation using continuation
techniques Champneys2013 .

### Model Reduction

The methods discussed above are well suited for computing FRCs and stability
diagrams in low-dimensional examples under parametric excitation. For the
analysis of realistic high-dimensional systems, however, reduced-order models
(ROM) are necessary. While projection onto modal subspaces is a powerful ROM
technique for linear systems, in nonlinear systems, those subspaces lose their
invariance. The dynamics and the trajectories projected on them thus no longer
correspond to those of the full dynamical system.

In constrast, nonlinear normal modes (NNM) Shaw1991 are invariant manifolds
that perturb from linear modal subspaces and serve as the nonlinear analogues
of those subspaces. Using projections to NNMs, Warminski et al. Warminski2012
analyzed parametrically and externally excited mechanical systems. Sinha and
coworkers Sinha2005 ; SinhaManifolds employed Lyapunov-Floquet transformations
to parametrically forced systems to obtain systems with autonomous linear part
before computing NNMs. However, this method is feasible only for low-
dimensional systems because of its reliance on the initial transformation.
Furthermore, the existence and uniqueness of such NNMs is a priori unclear
ssmexist .

The theory of Spectral Submanifolds (SSMs) ssmexist has clarified the
mathematical preliminaries required for rigorous model reduction of nonlinear
mechanical systems. Specifically, an SSM is the unique smoothest invariant
manifold that perturbs from a modal subspace under the addition of
nonlinearities and small-amplitude parametric or external forcing. The
existence and uniqueness of an SSM are guaranteed under appropriate non-
resonance conditions ssmexist . SSMs associated with slow eigenspaces (slow
SSMs) attract full system trajectories. Hence, the reduced dynamics on an SSM
serves as an exact reduced-order model for the full system.

In the recent years, significant advance has been made towards the automated
computation of SSMs. In a first implementation, Ponsioen et al. SSMTool
proposed a routine for automated computation of two-dimensional SSMs in modal
coordinates. More recent efforts have focused on multidimensional SSM
computation in physical coordinates, using only the eigenvectors associated
with the master modal subspace Jain2021 . This alleviated the computational
burden imposed by a modal transformation and enabled SSM computations for
high-dimensional finite element models.

In particular, for systems without internal resonances, the reduced dynamics
on two-dimensional SSMs have been used to compute the forced response of
periodically forced nonlinear systems without any numerical simulation, simply
by finding the zero-level set of a function isolatedresponsesten . Coupled
with the automated SSM computation procedure Jain2021 , this makes the
extraction of FRCs very fast for high-dimensional finite element problems as
well. If internal resonances are present, SSMs of higher dimension need to be
computed over the linear modes involved in the resonance MingwuPart1 ;
MingwuPart2 . More recently, SSMs have also been constructed directly from
individual trajectory data Cenedese2022Data-drivenSystems , used in control
theoryMahlkneckt2022 ; Alora2022Data-DrivenRobots and for systems with
algebraic constraintsLi2022ModelSubmanifolds . All these developments,
however, used only a leading-order approximation for the nonautonomous terms
in the SSM and its reduced dynamics.

To capture parametric effects, a high-order approximation of non-autonomous
SSMs and their reduced dynamics is necessary. Ponsioen et al. Ponsioen2020
have already presented such higher-order computations of the non-autonomous
SSM, albeit in modal coordinates. For efficient memory management and for
avoiding redundant computations, they employ a multi-index notation in their
numerical algorithm. More recently, Opreni et al. Opreni2022 have developed
expressions for the computation of nonautonomous SSMs in physical coordinates.
Their expressions, however, are limited to systems with proportional damping
and geometric nonlinearities up to cubic order. They use these expressions to
analyze several mechanical systems subject to external periodic excitation
(but no parametric forcing).

### Our contributions

Here we develop an automated algorithm for computing arbitrarily high-order
approximations for non-autonomous SSMs and their reduced dynamics for
mechanical systems in physical coordinates with general damping and
geometric/material nonlinearities of arbitrary polynomial order. Furthermore,
we develop separate expressions for general first-order dynamical systems in
the Appendix. While our applications focus on periodically excited systems,
our expressions and their implementations are developed for the general
quasiperiodic case. In contrast to the prior approaches Jain2021 ;
Vizzaccaro2022 ; Opreni2022 that use a tensor-based notation, we use a
multiindex notation, which minimizes the memory requirements and the number of
equations that need to be solved at each polynomial degree. Indeed, we argue
that for two-dimensional SSMs, the overhead of equations resulting from the
use of tensors increases exponentially with the order of the approximation.

In prior work Vizzaccaro2022 ; Opreni2022 , a special structure in the
invariance equations associated with second-order mechanical systems
containing up to cubic-order geometric nonlinearity and proportional damping
was used to reduce the number of equations in SSM computation. We now provide
a theoretical foundation and physical intuition for this structure, linked to
SSM theory and computations for general second-order dynamical systems.

For the case of two-dimensional periodic SSMs, we develop the following
results. For the first time, we demonstrate how SSM theory can be used to
obtain stability diagrams for the principal parametric resonance of finite
element examples directly from the ROM. Additionally, we analyze the forced
response arising from these instabilities using the reduced dynamics on the
SSM. We show how the FRCs of mechanical systems subject to external as well as
parametric periodic excitation can be obtained using an analytic
transformation of the reduced dynamics on the SSM, without the need for any
numerical continuation. This result allows for parallelization of the FRC
computation over any desired frequency range. Finally, we provide an open-
source implementation of our computational routines and numerical examples in
the MATLAB-based software package SSMTool 2.4 SSMTool2 .

The remainder of the paper is structured as follows. We discuss the general
class of mechanical systems we treat in Section II. Section III provides the
basic results from SSM theory relevant for our computations. In Section IV, we
present the autonomous and nonautonomous invariance equations associated to
the SSM, and develop specific resonance conditions which determine the normal
form of the SSM-reduced dynamics. We provide an automated routine for the
extraction of stability diagrams in Section V along with numerical examples.
Furthermore, we develop explicit expressions that enable analytical extraction
of FRCs from two-dimensional, nonautonomous SSMs along with examples in
Section VI. Finally, we draw conclusions in Section VII.

##  II Setup

We consider a quasiperiodic forced nonlinear mechanical system of dimension
$n$.

$\displaystyle\bm{M\ddot{y}+C}$
$\displaystyle\bm{\dot{y}+Ky+\bm{f}(y,\dot{y})}=\epsilon\bm{g}(\bm{\Omega}t,\bm{y},\bm{\dot{y}}),\
\ 0\leq\epsilon\ll 1,$ (1)

with $\bm{y}\in\mathbb{R}^{n}$ a generalized displacement; and
$\bm{M},\bm{C},\bm{K}\in\mathbb{R}^{n\times n}$ denoting the mass, damping and
stiffness matrices, respectively. The function $\bm{f}(\bm{y},\bm{\dot{y}})$
is a nonlinear function that is $r$ times continuously differentiable with
$r\geq 1$. The mechanical system is subject to periodic parametric excitation
via $\bm{f}_{\text{param}}(\bm{\Omega}t,\bm{y},\bm{\dot{y}})$ and possibly to
external periodic forcing $\bm{f}_{\text{ext}}(\bm{\Omega}t)$ with the
frequency vector $\bm{\Omega}\in\mathbb{R}^{K}$ for some $K\geq 0$. Here
$\bm{f}_{\text{param}}$ is a general function of $\bm{y}$ and $\bm{\dot{y}}$
and is taken to be $C^{r}$ in these variables. The two types of excitation are
included in the general excitation function

$\bm{g}(\bm{\Omega}t,\bm{y},\bm{\dot{y}})=\bm{f}_{\text{ext}}(\bm{\Omega}t)+\bm{f}_{\text{param}}(\bm{\Omega}t,\bm{y},\bm{\dot{y}}),$
(2)

where $\bm{g}$ will be autonomous for $K=0$, periodic in $t$ for $K=1$, and
quasiperiodic in $t$ for $K>1$ with $K$ rationally incommensurate frequencies.

The second-order system (1) can be transformed to the first-order form

$\displaystyle\bm{B\dot{z}=Az+F(z)}+\epsilon\bm{G}(\bm{\Omega}t,\bm{z}),$ (3)

where

$\displaystyle\bm{z}=\begin{bmatrix}\bm{y}\\\
\dot{\bm{y}}\end{bmatrix},\bm{A}=\begin{bmatrix}-\bm{K}&\bm{0}\\\
\bm{0}&\bm{M}\end{bmatrix},\bm{B}=\begin{bmatrix}\bm{C}&\bm{M}\\\
\bm{M}&\bm{0}\end{bmatrix},$ (4)
$\displaystyle\bm{F(z)}=\begin{bmatrix}-\bm{f(y,\dot{y})}\\\
\bm{0}\end{bmatrix},\bm{G}(\bm{\Omega}t,\bm{z})=\begin{bmatrix}\bm{g}(\bm{\Omega}t,\bm{y},\bm{\dot{y}})\\\
\bm{0}\end{bmatrix},$ (5)

with $\bm{z}\in\mathbb{R}^{N},N=2n$. This choice for the first-order form is
not uniqueJain2021 but all such forms can be used for the computation of SSMs.
The autonomous linear part of system (3),

$\displaystyle\bm{B\dot{z}=Az},$ (6)

has a fixed point at $\bm{z}=\bm{0}$, which we assume to be hyperbolic and
stable. We focus on an even subset of $M<<N$ eigenvalues and eigenvectors of
system (6), the choice of which is discussed below. Note that the task of
finding specific subsets of eigenvalues and eigenvectors can be carried out
efficiently with iterative algorithms, as presented, for instance, by Golub
and Van Loan Golub2013MatrixComputations .

The reduced left and right eigenvalue problems are given by

$\displaystyle(\bm{A}-\lambda_{j}\bm{B})\bm{v}_{j}=\bm{0},\qquad\bm{u}^{*}_{j}(\bm{A}-\lambda_{j}\bm{B})=\bm{0}$
(7)
$\displaystyle(\bm{A}-\overline{\lambda_{j}}\bm{B})\overline{\bm{v}_{j}}=\bm{0},\qquad\overline{\bm{u}^{*}_{j}}(\bm{A}-\overline{\lambda_{j}}\bm{B})=\bm{0}$
(8)

for $j=1,...,M/2$. We arrange these eigenvectors in the increasing order of
magnitudes of the real parts of the associated eigenvalues, i.e.,

$\displaystyle\text{Re}(\lambda_{M/2})\leq\text{Re}(\lambda_{M/2-1})\leq...\leq\text{Re}(\lambda_{1})<0.$
(9)

The choice of $\bm{A}$ and $\bm{B}$ in eq. (4) results in displacement and
velocity variables that are are inherently related, such that

$\displaystyle\bm{v}_{j}=\begin{bmatrix}\bm{\phi}_{j}\\\
\bm{\phi}_{j}\lambda_{j}\end{bmatrix},\
\bm{u}_{j}=\begin{bmatrix}\bm{\theta}_{j}\\\
\bm{\theta}_{j}\overline{\lambda_{j}}\end{bmatrix},$ (10)
$\displaystyle\overline{\bm{v}_{j}}=\begin{bmatrix}\overline{\bm{\phi}_{j}}\\\
\overline{\bm{\phi}_{j}\lambda_{j}}\end{bmatrix},\
\overline{\bm{u}_{j}}=\begin{bmatrix}\overline{\bm{\theta}_{j}}\\\
\overline{\bm{\theta}_{j}}\lambda_{j}\end{bmatrix},$ (11)

where the vectors $\bm{\phi}_{j},\bm{\theta}_{j}\in\mathbb{C}^{n}$ are the
right and left eigenvectors associated with the quadratic eigenvalue problem

$\displaystyle(\lambda_{j}^{2}\bm{M}+\lambda_{j}\bm{C}+\bm{K})\bm{\phi}_{j}=\bm{0},\quad\bm{\theta}_{j}^{*}(\lambda_{j}^{2}\bm{M}+\lambda_{j}\bm{C}+\bm{K})=\bm{0}$
(12)
$\displaystyle(\overline{\lambda_{j}^{2}}\bm{M}+\overline{\lambda_{j}}\bm{C}+\bm{K})\overline{\bm{\phi}_{j}}=\bm{0},\quad\overline{\bm{\theta}_{j}^{*}}(\overline{\lambda_{j}^{2}}\bm{M}+\overline{\lambda_{j}}\bm{C}+\bm{K})=\bm{0}$
(13)

for $j=1,...,M/2$.

Note that the above procedure is valid for any general damping matrix $\bm{C}$
and reduces to the computation of conservative eigenmodes when $\bm{C}$ is
simultaneously diagonalizable with the mass and stiffness matrices, e.g., in
the case of proportional damping.

##  III Non-autonomous invariant manifolds and model reduction

We now introduce the methodology for computing the SSMs on which the ROM of
the parametrically excited system (1) is constructed. The complex-conjugate
pairs of right eigenvectors span spectral the subspaces
$E_{j}=\text{span}\\{\bm{v}_{j},\overline{\bm{v}_{j}}\\}$, which are invariant
under the linear flow of system (6). To analyze the linear system response via
a ROM, we select a $M$-dimensional spectral subspace ${\cal E}=E_{1}\oplus
E_{2}\oplus...\oplus E_{{M/2}}$, the master spectral subspace. The modes
$E_{j}$ are chosen such that the corresponding spectrum

$\displaystyle\text{Spec}(\mathcal{E})=\\{\lambda_{1},\overline{\lambda_{1}}...,\lambda_{M/2},\overline{\lambda_{M/2}}\\}$
(14)

contains pairs of eigenvalues that potentially assume near-resonance
relationships with either the forcing frequencies or among each other. In the
absence of any (near) resonances, the subspace $\mathcal{E}$ typically
contains the slowest modes of the system. The corresponding eigenvectors are
normalised such that

$\displaystyle(\bm{u}_{i})^{*}\bm{B}\bm{v}_{j}=\delta_{ij}$ (15)

which implies
$\bm{\theta}_{i}^{*}\bm{C}\bm{\phi}_{j}+\lambda_{i}\bm{\theta}_{i}^{*}\bm{M}\bm{\phi}_{j}+\lambda_{j}\bm{\theta}_{j}^{*}\bm{M}\bm{\phi}_{i}=\delta_{ij}$.
Note that this normalization choice is different from the commonly adopted
mass normalization for proportionally damped systems, i.e.,
$\bm{\phi}_{i}^{*}\bm{M}\bm{\phi}_{j}=\delta_{ij}$.

Under the addition of nonlinearities to system (6), the subspace $\mathcal{E}$
does not remain invariant. Model reduction onto $\mathcal{E}$ then looses its
justification as the trajectories projected onto $\mathcal{E}$ no longer
correspond to the trajectories of the full dynamical system. A mathematically
rigorous alternative is to construct an invariant manifold that acts as a
nonlinear continuation of $\mathcal{E}$. Spectral Submanifolds (SSM) ssmexist
are such continuations with their existence conditioned on a set of resonance
conditions Ponsioen2020 . Specifically, an autonomous SSM,
$\mathcal{W}(\mathcal{E})$, is the unique smoothest invariant manifold that
perturbs smoothly from the spectral subspace $\mathcal{E}$ under the addition
of the nonlinear autonomous terms $\bm{f(x,\dot{x})}$ to the linearized
dynamical system (6). The dynamics on such an SSM are called its  reduced
dynamics and serve as an exact ROM of the full dynamical system (1) for
$\varepsilon=0$.

Under the addition of quasiperiodic excitation (i.e., for $\varepsilon>0$,
small), this autonomous manifold starts to oscillate and deform, becoming a
non-autonomous SSM, $\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$.
Instead of being anchored to the trivial fixed point $\bm{z=0}$, the
quasiperiodic SSM is now attached to $\mathcal{O}(\varepsilon)~{}C^{r}$-close
invariant torus $\bm{\gamma}_{\epsilon}$ into which this fixed point
perturbsGuckenheimer1983 . The tangency of the autonomous SSM to $\mathcal{E}$
implies that the nonautonomous SSM perturbs from the spectral subbundle
$\bm{\gamma}_{\epsilon}\times\mathcal{E}$. Now, we will use

$\bm{\phi}=\bm{\Omega}t\in{\mathbb{T}}^{K}$

to denote the phase of the quasiperiodic excitation term. For computing
$\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$ and its reduced dynamics,
the following results are useful ssmexist .

  1. 1.

$\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$ can be described by a time
dependent parametrization

$\displaystyle\bm{W}_{\epsilon}(\bm{p},\bm{\phi}):\mathcal{U}=U\times\mathbb{T}^{K}\mapsto\mathbb{R}^{N},\quad
U\in\mathbb{C}^{M},$ (16)

from an open set $\mathcal{U}$ onto the phase space of the full system (3).

  2. 2.

The reduced vector field
$\bm{R}_{\epsilon}(\bm{p},\bm{\phi}):\mathcal{U}\times\mathbb{T}^{K}\mapsto\mathcal{U}$
on $\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$ satisfies the
invariance equation

$\displaystyle\bm{B}((\partial_{\bm{p}}{\bm{W}}_{\epsilon})\bm{R}_{\epsilon}+(\partial_{\bm{\phi}}\bm{W}_{\epsilon})\bm{\Omega})=\bm{A}\bm{W}_{\epsilon}+\bm{F}\circ\bm{W}_{\epsilon}+\epsilon\bm{G}(\bm{\phi},\bm{W}_{\epsilon})$
(17)

and generates the reduced dynamics on
$\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$ via the ODE

$\displaystyle\dot{\bm{p}}=\bm{R}_{\epsilon}(\bm{p},\bm{\phi}),\quad\dot{\bm{\phi}}=\bm{\Omega}.$
(18)

If for a given dynamical system we wish to obtain a mathematically exact ROM,
it is therefore our task to try and construct representations of the invariant
manifold parametrization $\bm{W}_{\epsilon}$ and the reduced vector field
$\bm{R}_{\epsilon}$. Once this has been achieved, the low-dimensional reduced
dynamics can be analyzed in detail. Due to the invariance of the SSM
$\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$, any periodic orbit,
invariant torus, and bifurcation that is observed in eq. (18) also exists in
the full system and can be mapped back using the SSM parametrization
$\bm{W}_{\epsilon}$ to the full system (1).

##  IV SSM Computation

In this section, we address how the SSM parametrization $\bm{W}_{\epsilon}$
and the reduced vector field $\bm{R}_{\epsilon}$ can be efficiently computed
up to arbitrary orders of accuracy for general mechanical systems. We will
merely require knowledge of the eigenvectors and eigenvalues of the underlying
master modal subspace over which the SSM is constructed. All computations will
be carried out in physical coordinates, without the need of a modal
transformation. This will enable the treatment of large dynamical systems
arising, e.g., from finite-element models.

To eliminate the redundancy that would be introduced by a tensor-based
notation for the expansions, we compute the SSM parametrization and the
reduced dynamics with the use of multi-index notation (see Appendix A.1 and
Ref. Ponsioen2020, ). This optimizes the memory requirements for storing and
manipulating the expansion coefficients and reduces the number of computations
to be performed. We will present results for both the autonomous and
nonautonomous cases. The general treatment for first-order ODEs that do not
originate from an underlying second-order mechanical system is given in
Appendix A.3. The complete algorithm and the implementation of the examples
treated here are openly available in the MATLAB-based software package SSMTool
2.4 SSMTool2 .

In summary, our treatment here extends the available results Ponsioen2020 ;
Jain2021 ; Vizzaccaro2022 ; Opreni2022 on automated computation of SSMs.
Specifically, we combine the advantages of a multi-index approach Ponsioen2020
, yet use physical coordinates Jain2021 . We compute the nonautonomous SSM up
to arbitrary orders of approximation and exploit a special structure inherent
in phase space for second-order systems Vizzaccaro2022 ; Opreni2022 . We also
provide the theoretical foundation for this special structure (see Lemmas 2,
4), allowing extension of its advantages to any second-order system. At the
same time, we consider more general systems, including direct parametric
excitation, and geometric and system nonlinearities of arbitrary polynomial
order including nonlinear damping.

###  IV.1 SSM parametrization

For a second-order system, the phase space variables
$\bm{z}=(\bm{y},\dot{\bm{y}})^{T}$ are inherently connected as the velocity
components are simply the time derivative of the displacement variables.
Naturally, this relationship carries over to the SSM parametrization
${\bm{W}}_{\epsilon}(\bm{p},\bm{\phi})$, which maps to the phase space. This
allows us to express the SSM parametrization as

$\displaystyle{\bm{W}}_{\epsilon}(\bm{p},\bm{\phi})=\begin{bmatrix}{\bm{w}}_{\epsilon}(\bm{p},\bm{\phi})\\\
\dot{{\bm{w}}}_{\epsilon}(\bm{p},\bm{\phi})\end{bmatrix},$ (19)

where ${\bm{w}}_{\epsilon}$ and $\dot{{\bm{w}}}_{\epsilon}$ denote the
parametrizations for the displacement and the velocity variables. The smooth
dependence of the SSM $\mathcal{W}_{\bm{\gamma}_{\epsilon}}(\mathcal{E})$ on
$\epsilon$ allows us to further expand these parametrizations in $\epsilon$ as

$\displaystyle{\bm{W}}_{\epsilon}(\bm{p},\bm{\phi})$
$\displaystyle=\bm{W}(\bm{p})+\epsilon\bm{X}(\bm{p},\bm{\phi})+\cdots$ (20)
$\displaystyle=\begin{bmatrix}\bm{w}(\bm{p})\\\
\dot{\bm{w}}(\bm{p})\end{bmatrix}+\epsilon\begin{bmatrix}\bm{x}(\bm{p},\bm{\phi})\\\
\dot{\bm{x}}(\bm{p},\bm{\phi})\end{bmatrix}+\dots.$ (21)

Similarly, we expand the reduced dynamics parametrization as

$\displaystyle{\bm{R}}_{\epsilon}\bm{(}p,\bm{\phi})=\bm{R}(\bm{p})+\epsilon\bm{S}(\bm{p},\bm{\phi})+O(\epsilon^{2}).$
(22)

We will show that, in order to compute the SSM for second-order systems, it
suffices to solve the invariance equation for the displacement parametrization
${\bm{w}}_{\epsilon}$ only. Subsequently, the velocity parametrization
$\dot{{\bm{w}}}_{\epsilon}$ can be analytically reconstructed by
differentiating ${\bm{w}}_{\epsilon}$ with respect to time.

###  IV.2 Autonomous SSM ($\varepsilon=0$)

In order to compute autonomous SSMs, the autonomous invariance equation

$\displaystyle\bm{B}(\text{D}\bm{W}(\bm{p}))\bm{R}(\bm{p})=\bm{A}\bm{W}(\bm{p})+\bm{F}(\bm{W}(\bm{p}))$
(23)

must be solved. This equation is obtained by substituting the expansions (21)
and (22) into the invariance equation (17) and collecting terms at
$\mathcal{O}(\epsilon^{0})$. As we have noted, expansions of the SSM
parametrization and the reduced dynamics will be written using multi-indices.
A multi-index $\bm{m}$ is a vector $\bm{m}\in\mathbb{N}^{M}$ for which
addition, subtraction and other operations are defined elementwise. With the
help of multi-indices, multivariate monomials can be written as
$\bm{p}^{\bm{m}}=p_{1}^{m_{1}}\dots p_{M}^{m_{M}}$.

The autonomous parametrizations of the SSM are expanded using these multi-
indices

$\displaystyle\bm{w}(\bm{p})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{w}_{\bm{m}}\bm{p}^{\bm{m}},$
(24) $\displaystyle\dot{\bm{w}}(\bm{p})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\dot{\bm{w}}_{\bm{m}}\bm{p}^{\bm{m}}$
(25)

where $\bm{w}_{\bm{m}},\bm{\dot{w}}_{\bm{m}}\in\mathbb{C}^{n}$. As we
elaborate in Appendix A.1, when compared to tensor-based approachesJain2021 ;
Vizzaccaro2022 ; Opreni2022 , this notation results in a drastic reduction in
the number of systems of equations that need to be solved, depending on the
SSM dimension and the order of approximation. For instance, in the case of
two-dimensional SSMs, we observe an exponentially lower number of equations
for multi-index notation relative to the tensor-notation as the order of
approximation increases.

The autonomous reduced vector field $\bm{R}\bm{(p)}$ is also expanded using
multi-index notation as

$\displaystyle\bm{R}(\bm{p})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{R}_{\bm{m}}\bm{p}^{\bm{m}},\
\bm{R}_{\bm{m}}=\begin{bmatrix}R^{1}_{\bm{m}}\\\ \vdots\\\
R^{M}_{\bm{m}}\end{bmatrix}$ (26)

Substituting the expansions (24),(25) and (26) into the autonomous invariance
equation (23) results in linear system of equations that can be solved
recursively at each order for the unknown coefficients
$\bm{W}_{\bm{m}},\bm{R}_{\bm{m}}$. These equations are decoupled for distinct
multi-indices at each order. Furthermore, the splitting (24) of the
parametrization yields a single $n$-dimensional system of linear equations
that can be solved for the displacement coefficients $\bm{w}_{\bm{m}}$
according to the following statement.

######  Lemma 1.

The autonomous invariance equation (23) for a multi-index
$\bm{$}canbewrittenas\begin{aligned}
\underbrace{\bigg{(}\bm{K}+\Lambda_{\bm{m}}\bm{C}+\Lambda_{\bm{m}}^{2}\bm{M}\bigg{)}}_{:=\bm{L}_{\bm{m}}}\bm{w}_{\bm{m}}=\bm{D}_{\bm{m}}\bm{R}_{\bm{m}}+\bm{C_{m}},\end{aligned}where$C_m∈C^n,D_m∈C^n×M$aredefinedinAppendix~{}\ref{app_sec:SSMComp_2ndorder},and$Λ_m
= Λ⋅m$with$Λ = (λ_1, …, ¯λ_M/2) .

###### Proof.

See Appendix A.2. ∎

In system (1), the matrix $\bm{D_{m}}$ is computed solely based on the
eigenvalues and eigenvectors associated to the master subspace, whereas
$\bm{C_{m}}$ contains contributions from the lower-order terms in the
expansion of the SSM and its reduced dynamics. Based on Lemma 1, the system
(1) can be solved recursively up to any arbitrary order for the autonomous
coefficients $\bm{w}_{\bm{m}},\bm{R}_{\bm{m}}$ of the SSM parametrization and
the reduced dynamics. At any order $m$, we have $z_{m}=\binom{m+M-1}{M-1}$
distinct multi-indices in the expansions (24), (26). Thus, at each order, we
need to solve $z_{m}$ systems of $n$-dimensional linear equations. Indeed, for
real dynamical systems, symmetries allow for further reduction in the number
of equations, as we elaborate in Appendix A.3.

######  Remark.

In case of near-resonances of the form

$\displaystyle\lambda_{i}\approx\bm{\Lambda}_{\bm{m}},$ (27)

the coefficient matrix $\bm{L}_{\bm{m}}$ will be nearly singular Jain2021 . If
the resonance condition (27) is fulfilled for any eigenvalue
$\lambda_{i}\in\text{Spec}(\mathcal{E})$ then a nonzero choice for
$\bm{R}_{\bm{m}}$ is necessary to solve the linear system (1), where different
choices lead to different styles of parametrization Haro2016 ; Jain2021 .

To understand the resonance condition (27), consider a two-dimensional SSM
tangent to a modal subspace spanned by lightly damped (complex-conjugate)
modes with associated eigenvalues $\lambda$ and $\overline{\lambda}$
satisfying

$\displaystyle\lambda\approx l\lambda+(l-1)\overline{\lambda}$ (28)

for all small enough $l\in\mathbb{N}^{+}$. This relationship holds for the
pair $(\lambda,\overline{\lambda})$ if its real part is small, i.e.,
underlying physical system is lightly damped. In this case, for any multi-
index of the form $\bm{m}=\begin{bmatrix}l,&l-1\end{bmatrix}$, the coefficient
matrix $\bm{L}_{\bm{m}}$ in system (1) turn out to be nearly singular.

Finally, once we have computed the coefficients $\bm{w}_{\bm{m}}$ for the
displacement parametrization $\bm{w(p)}$, the coefficients for the velocity
parametrization can be directly obtained using the following statement.

######  Lemma 2.

The coefficients of the velocity parametrization corresponding to a multi-
index $\bm{m}$ can be obtained directly from the displacement parametrization
as

$\displaystyle\dot{\bm{w}}_{\bm{m}}=[\text{D}\bm{w}(\bm{p})\bm{R}(\bm{p})]_{\bm{m}}$
(29)

###### Proof.

See Appendix A.2. ∎

###  IV.3 Non-autonomous SSM ($\epsilon>0$)

In the non-autonomous setting, inserting the expansions (21),(22) into the
invariance equation (17), and collecting terms of $\mathcal{O}(\epsilon)$
leads to the non-autonomous invariance equation

$\displaystyle\bm{B}\bigg{(}(\text{D}\bm{W}(\bm{p}))\bm{S}(\bm{p},\bm{\phi})+[\partial_{\bm{p}}\bm{X}(\bm{p},\bm{\phi})]\bm{R}(\bm{p})+[\partial_{\bm{\phi}}\bm{X}(\bm{p},\bm{\phi})]\bm{\Omega}\bigg{)}=$
$\displaystyle\bm{A}\bm{X}(\bm{p},\bm{\phi})+\big{[}\text{D}\bm{F}(\bm{W}(\bm{p}))\big{]}\bm{X}(\bm{p},\bm{\phi})+\bm{G}(\bm{W}(\bm{p},\bm{\phi}),\bm{\phi}),$
(30)

which needs to be solved for the non-autonomous terms $\bm{X}$ and $\bm{S}$ in
the parametrizations for the SSM and its reduced dynamics. As in the
autonomous case, we now expand the nonlinear functions in Taylor series using
multi-indices. The coefficients of these series, however, depend
quasiperiodically on time. Hence, we futher expand these coefficients in
temporal Fourier series to obtain

$\displaystyle\bm{x}(\bm{p},\bm{\phi})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\sum_{\bm{\kappa}\in\mathbb{Z}^{k}}\bm{x}_{\bm{m},\bm{\kappa}}e^{i\langle\bm{\kappa},\bm{\phi}\rangle}\bm{p}^{\bm{m}},$
(31) $\displaystyle\dot{\bm{x}}(\bm{p},\bm{\phi})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\sum_{\bm{\kappa}\in\mathbb{Z}^{k}}\dot{\bm{x}}_{\bm{m},\bm{\kappa}}e^{i\langle\bm{\kappa},\bm{\phi}\rangle}\bm{p}^{\bm{m}}.$
(32)

The reduced dynamics are similarly expanded in a Taylor-Fourier Series as

$\displaystyle\bm{S}(\bm{p},\bm{\phi})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\sum_{\bm{\kappa}\in\mathbb{Z}^{k}}\bm{S}_{\bm{m},\bm{\kappa}}e^{i\langle\bm{\kappa},\bm{\phi}\rangle}\bm{p}^{\bm{m}},$
$\displaystyle\ \bm{S}_{\bm{m},\bm{\kappa}}$
$\displaystyle=\begin{bmatrix}S^{1}_{\bm{m},\bm{\kappa}}\\\ \vdots\\\
S^{M}_{\bm{m},\bm{\kappa}}\end{bmatrix}.$ (33)

The procedure of determining these coefficients recursively is similar to that
of the autonomous case. The contributions to distinct harmonics and to
different multi-indices can be considered independently as the corresponding
equations decouple. For second-order systems, the computation of these
coefficients can again be carried out by first obtaining the displacement
parametrization $\bm{x}(\bm{p},\bm{\phi})$ via a reduced invariance equation
of dimension $n$, followed by computing the velocity parametrization
$\dot{\bm{x}}(\bm{p},\bm{\phi})$ directly from the time derivative of the
displacement terms, as we detail below.

######  Lemma 3.

The non-autonomous invariance equation in terms of a multi-index $\bm{m}$ and
harmonic $\bm{\kappa}$ is given by

$\displaystyle\underbrace{\big{(}\bm{K}+\Lambda_{\bm{m},\bm{\kappa}}\bm{C}+\Lambda_{\bm{m},\bm{\kappa}}^{2}\bm{M}\big{)}}_{\bm{L}_{\bm{m},\bm{\kappa}}}\bm{x}_{\bm{m},\bm{\kappa}}=\bm{D}_{\bm{m},\bm{\kappa}}\bm{S}_{\bm{m},\bm{\kappa}}+\bm{C}_{\bm{m},\bm{\kappa}},$
(34)

$\bm{C}_{\bm{m},\bm{\kappa}}\in\mathbb{C}^{n}$,
$\bm{D}_{\bm{m},\bm{\kappa}}\in\mathbb{C}^{n\times M}$ are defined in Appendix
A.2, and
$\Lambda_{\bm{m},\bm{\kappa}}:=\bm{\Lambda}\cdot\bm{m}+i\bm{\Omega}\cdot\bm{\kappa}$.

###### Proof.

See Appendix A.2. ∎

In system (34), the matrix $\bm{D}_{\bm{m},\bm{\kappa}}$ depends on the
eigenvalues and eigenvectors of the master subspace as well as on the forcing
frequency $\bm{\Omega}$. The vector $\bm{C}_{\bm{m},\bm{\kappa}}$ contains
contributions from lower-order terms in the parametrization of the SSM and its
reduced dynamics. The invariance system (34) contains decoupled equations for
each pair of multi-indices and harmonics $(\bm{m},\bm{\kappa})$. Thus, for
$z_{m}$ distinct multi-indices and $z_{K}$ distinct harmonics, we must solve
$z_{m}z_{K}$ systems of $n-$dimensional linear equations.

######  Remark.

In case of near-resonances of the form

$\displaystyle\lambda_{i}\approx\bm{\Lambda}\cdot\bm{m}+i\bm{\Omega}\cdot\bm{\kappa}$
(35)

the coefficient matrix $\bm{L}_{\bm{m},\bm{\kappa}}$ will be nearly singular
Jain2021 . If the resonance condition (27) is fulfilled for any eigenvalue
$\lambda_{i}\in\text{Spec}(\mathcal{E})$ then a nonzero choice for
$\bm{S}_{\bm{m},\bm{\kappa}}$ is necessary to solve the linear system (1),
where different choices lead to different styles of parametrization Haro2016 ;
Jain2021 .

As an example of the resonance condition (35), consider a two-dimensional
periodic SSM. At the leading order, the resonance condition then reduces to
$\lambda\approx i\Omega$, which is a simple resonance created by the external
periodic excitation. Higher-order terms lead to more complicated resonances
which are able to capture the effects of parametric resonance. For instance,
the a principal parametric resonance of the form $\Omega\approx
2\text{Im}(\lambda)$ results in

$\displaystyle\lambda\approx\overline{\lambda}+i\Omega$ (36)

in terms of relationship (35). Finally, once we have computed the coefficients
$\bm{x}_{\bm{m},\kappa}$ for the displacement parametrization
$\bm{x(p,\phi)}$, the coefficients for the nonautonomous velocity
parametrization can be obtained directly using the following statement (cf.
Lemma 2).

######  Lemma 4.

The coefficients of the nonautonomous velocity parametrization
$\dot{\bm{x}}_{\bm{m},\bm{\kappa}}$ for a multi-index $\bm{m}$ and harmonic
$\bm{\kappa}$ can be computed directly from the displacement parametrization
as

$\displaystyle\dot{\bm{x}}_{\bm{m},\bm{\kappa}}=[\text{D}\bm{w}(\bm{p})\cdot\bm{S}+\partial_{\bm{p}}\bm{x}\cdot\bm{R}+\partial_{\bm{\phi}}\bm{x}\cdot\bm{\Omega}]_{\bm{m},\bm{\kappa}}$
(37)

###### Proof.

See Appendix A.2 ∎

Prior works Vizzaccaro2022 ; Opreni2022 identified a special structure in the
invariance equations of second-order mechanical systems with proportional
damping and up to cubic-order geometric nonlinearities. This special structure
led to reduced computational costs by halving the number of invariance
equations relative to a first-order setting. We have shown that such a
reduction in the number of equations arises from the underlying structure (19)
of the phase space variables in any second-order system. Indeed, according to
Lemmas 1-4, these computational advantages apply to general second-order
systems with nonlinearities of arbitrary polynomial order, including nonlinear
damping.

##  V Stability Diagrams and SSM

###  V.1 Mathematical Foundation

In this section, we consider mechanical systems subject to periodic parametric
excitation, which means that $\bm{f}^{ext}=\bm{0}$ in eq. (2). For such
systems, $\bm{z=0}$ remains a fixed point for system (3) even for
$\epsilon>0$. This fixed point may, however, be destabilized by resonances
between the eigenfrequencies of the linear system and the parametric
excitation frequency.

In order to study the stability of this fixed point for $\epsilon>0$, the
dynamical system (3) can be extended to an autonomous system of variables
$(\bm{z},\tau)\in\mathbb{R}^{N}\times S^{1}$ such that the trivial fixed point
can be interpreted as the periodic orbit $(\bm{z},\tau)=(\bm{0},t\ \text{mod}\
2\pi)$Champneys2013 . Changes in the stability of this periodic orbit result
in bifurcations with respect to parameters. Tools such as the po toolbox of
the continuation package coco Dankowicz2013 can detect and continue families
of such bifurcations using Floquet theory. We refer to Dankowicz and Schilder
Dankowicz2013 ) for a detailed account. Numerical continuation for finding an
initial bifurcation that characterizes the stability boundary may be carried
out by varying either the forcing frequency or the amplitude, as shown in
Figure 1.

Figure 1: Sketch of a stability diagram at the principal resonance. The thin
dashed line corresponds to the boundary of stable and unstable regions in the
conservative limit. When damping is added, the tongue starts to lift up,
bounded by the thick black line. The numerical continuation in the variables
$\epsilon$ (blue arrow) and $\Omega$ (red arrow), respectively, leads to the
detection of a bifurcation indicated by the blue and red dots. Choosing these
bifurcations as a starting point, the family of bifurcations on the boundary
can then be found via further continuation.

We employ SSM computation in order to efficiently obtain the _resonance
tongue_ of the principal resonance, which marks the stability boundary in
Figure 1. The SSM is constructed over the spectral subspace corresponding to
the resonant eigenfrequency $\omega_{i}$. Consequently, the $M$-dimensional
ROM

$\displaystyle\dot{\bm{p}}$ $\displaystyle=\bm{R}_{\epsilon}(\bm{p},\Omega t)$
(38)
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{R}_{\bm{m}}\bm{p}^{\bm{m}}+\epsilon\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{M}\\\
m\geq
1\end{subarray}}\sum_{\kappa\in\mathbb{Z}}\bm{S}_{\bm{m},\kappa}e^{i\kappa\Omega
t}\bm{p}^{\bm{m}}$ (39)

needs to be analyzed for bifurcations using the po toolbox of
cocoDankowicz2013 . Details on the corresponding continuation problem are
recounted in Appenidx C. As the forcing frequency changes, the ROM has to be
updated, which highlights the need for a fast computational algorithm to
calculate the SSM. The reduced dynamics are of low dimension, so continuation
is a suitable method, as it is fast on small systems and provides key insights
into the bifurcation behavior. The proposed approach thus combines the
advantages of reduced-order modeling and numerical continuation. This
procedure is now automated and available in SSMTool 2.4, which we use to
analyze the following examples. Unless stated otherwise, the computations are
performed using Matlab 2021a on a Macbook Air with a 1,4 GHz Intel Core i5
processor and 4GB of total RAM.

######  Remark.

Other resonances with nonmaster modes might also occur close to this principal
resonance and lead to resonance tongues that interact with the principal
resonance tongue. Naturally, these other resonances are not detected by the
ROM as it is ignorant of the spectrum of the other modes. Such near-resonance
relationships with the master mode are called internal resonances. For such
internally resonant systems, the existence and uniqueness of the underlying
SSM are guaranteed if these resonant modes are also included in the master
subspace, leading to a higher-dimensional SSM. We refer to Li et al.
MingwuPart1 ; MingwuPart2 for the details of SSM theory and computation for
internally resonant systems.

###  V.2 Examples

####  V.2.1 Mathieu equation

The most frequently discussed system in the parametric excitation literature
is the Mathieu equation

$\displaystyle m\ddot{y}+c\dot{y}+\omega_{0}^{2}(1+\epsilon\cos{\Omega
t})y=0,\quad 0\ll\epsilon\ll 1,$ (40)

whose stability is well understood. Various studies have appeared for the
damped and undamped casesVerhulst2009 , some with additional nonlinearities to
system Rhoads2006_MEMS and others in higher dimensions Hansen1985 . Here, we
extend the system (40) to a set of damped oscillators coupled via springs with
parametrically varying stiffness. Adding stiffness nonlinearities then results
in the mechanical system displayed in FIG. 2.

Figure 2:  Schematic model of two oscillators described by two coupled Mathieu
equations. The coupling springs have time-varying stiffness coefficients
$k(t)=k(1+\epsilon\cos(\Omega t))$. The induced parametric excitation may
destabilise the trivial fixed point if chosen in subharmonic resonance with an
eigenfrequency of the linearized dynamical system.

The resulting generalised Mathieu equation is given by

$\displaystyle\bm{M}\bm{\ddot{y}}+\bm{C}\bm{\dot{y}}+\bigg{(}\bm{K}+\epsilon\bm{Q}\cos(\Omega
t)\bigg{)}\bm{y}+\bm{f}(\bm{y})=\bm{0}.$ (41)

In the case considered here, these matrices and the nonlinearity are
explicitly given as

$\displaystyle\small\bm{M}$ $\displaystyle=\left[\begin{array}[]{cc}m_{1}&0\\\
0&m_{2}\end{array}\right],~{}\bm{C}=\left[\begin{array}[]{cc}2c&-c\\\
-c&2c\end{array}\right],~{}\bm{K}=\left[\begin{array}[]{cc}2k&-k\\\
-k&2k\end{array}\right],$ (48) $\displaystyle\bm{Q}$
$\displaystyle=\left[\begin{array}[]{cc}k&-k\\\
-k&k\end{array}\right],~{}\bm{f}({y_{1},y_{2}})=\kappa\left[\begin{array}[]{cc}-y_{1}^{3}-(y_{1}-y_{2})^{3}\\\
-y_{2}^{3}+(y_{1}-y_{2})^{3}\end{array}\right].$ (53)

We choose non-dimensional parameters $m_{i}=1,k=1,c=0.05$ and $\kappa=0.1$,
for which we obtain the eigenvalues of the $\epsilon=0$ system in the form
$\lambda_{1},\overline{\lambda_{1}}\approx-0.0250\pm 0.9997i$ and
$\lambda_{2},\overline{\lambda_{2}}\approx-0.0500\pm 1.7304i$. We seek to
analyze the nonlinear response under the principal parametric resonance
$\Omega\approx 2\omega_{2}$ by constructing the non-autonomous SSM
corresponding to the spectral subspace
$\mathcal{E}=\text{span}(\bm{v}_{2},\overline{\bm{v}_{2}})$. Using numerical
continuation in the reduced dynamics, we will detect the corresponding
resonance tongue.

First, we choose an initial parameter set $(\epsilon_{i},\Omega_{i})$, then
either of the two parameters is released, and the trivial
$2\pi/\Omega_{i}$-periodic response is continued in search of period-doubling
bifurcations. The detected bifurcation point is then continued in the
parameter space, yielding the stability diagram in FIG. 3a. The entire
procedure is automated, constantly updating the coefficients of the
nonautonomous reduced dynamics when the excitation frequency gets changed.

(a)

(b)

Figure 3: Stability diagram and forced response of the two coupled Mathieu
oscillators shown in FIG. 2. (a) Stability diagram of the second mode obtained
using an SSM-based ROM, which locally approximates the resonance tongue in a
region around the principal resonance where $\Omega\approx 2\omega_{2}$. (b)
The $||z_{1}||_{\infty}=0$ plane corresponds to the trivial response of the
full system. On this plane, the resonance tongue is displayed in gray shades.
From its boundary, orbits with period $4\pi/\Omega$ emerge that show hardening
behavior due to the nonlinearity.

To obtain a more complete picture of how these regions of stability influence
the system’s dynamics, we consider the non-trivial forced response of system
(41) which emanates from the instability. The orbits with period
$T=\frac{4\pi}{\Omega}\approx\frac{2\pi}{\omega_{0}}$ bifurcate out of the
trivial response at the stability boundary. The amplitude of these response
curves increases smoothly with increasing excitation amplitude $\epsilon$. The
reduced dynamics on the SSM serves as a ROM to explore periodic orbit families
via numerical continuation. The resulting stability diagram, as well as the
forced response, can be seen in FIG. 3b. The resonance tongue and the forced
response are accurately reproduced using the ROM, as we verify using the po-
toolbox of coco on the full dynamical system. The parameters used for
continuation of the full and reduced models are noted in Appendix C.

The computation time for the stability diagram using the full system was 22
seconds. The same computation using the ROM, including the continual update of
the SSM parametrization, takes 1 minute and 39 seconds. Similarly, for the
FRCs displayed in FIG. 3b the full system analysis took 1 minute and 57
seconds, while using the ROM amounted to 3 minutes and 43 seconds. Therefore,
because the full dynamical system is low-dimensional, the advantage of using a
ROM is outweighed by the computational cost of constructing the reduced
dynamics on the time-dependent SSM for each forcing frequency.

####  V.2.2 Bernoulli beam

(a)

(b)

Figure 4: (a) Model of a Bernoulli beam of length $L$ and height
$h$isolatedresponsesten . The beam is aligned with the $x$ axis, cantilevered
on one side, and parametrically excited on the other, as indicated by the red
arrow. A nonlinear (cubic) damper and a cubic spring are attached to the free
end of the beam. (b) Schematic depiction of a possible experimental setup for
parametrically exciting the beam. Two coils are used as tunable magnets via
the induced current $I$. Depending on its magnitude and direction, the
effective attraction and repulsion between the beam and the adjacent magnets
can be controlled. This results in a time varying stiffness Chen2001 .

We now consider the nonlinear Bernoulli beam shown in FIG. 4 and treated in
Ref. isolatedresponsesten, . It is cantilevered on one side and subject to a
cubic spring, a cubic damper, and parametric excitation at its free end. A
similar example has been experimentally studied by Chen and YehChen2001 , who
implemented parametric excitation through a pair of coils fixed to the beam.
Nearby magnets mounted on a support were then used to exert a force on these
coils when a current is applied.

We discretize the beam into 5 finite elements, resulting in $n=10$ degrees of
freedom. Damping is assumed to be proportional and of the form
$\bm{C}=\sigma(\alpha\bm{M}+\beta\bm{K})$. Here, the cubic nonlinearity
$\bm{f}(\bm{y},\dot{\bm{y}})$ includes the nonlinear damping coefficient
$\gamma$ and the nonlinear stiffness $\kappa$. The parameters and beam
properties are set according to TABLE LABEL:tab:BB_model, and the parametric
excitation amplitude is set to 1 Newton. This system is driven in the regime
of parametric resonance by targeting the principal resonance of the mode
associated with the leading pair of eigenvalues

$\displaystyle\lambda_{1},\overline{\lambda_{1}}\approx-0.1238\pm 6.9995i$
(54)

The prinicipal resonance is thus obtained for forcing frequency values
$\Omega\approx 13.999\ \text{rad}/s$. We use the reduced dynamics on the
corresponding SSM to construct the stability diagram of the principal
resonance of the selected mode. The resulting resonance tongue is shown in
FIG. 5a for various values of the damping coefficient $\sigma$. The forced
response obtained from SSM theory for a set of parametric excitation
amplitudes with SSM theory is shown in FIG.5b. We verified these reduced
results against the full system using coco, for which the relevant simulation
parameters can be found in Appendix C. The computation times for the stability
diagram and forced responses are reported in TABLE LABEL:tab:CompTimeBB_SD and
TABLE LABEL:tab:CompTimeBB_FRC, which show that the computational savings
arising from the use of SSM-based ROM in this higher-dimensional example are
significant, unlike in our first example.

Table 1: Physical parameters used for the model of the Bernoulli beam. Symbol | Meaning | Value
---|---|---
$L$ | Length of beam | 2.7 (m)
$h$ | Height of beam | 10 (mm)
$b$ | Width of beam | 10 (mm)
$E$ | Young’s modulus |  45 $\times 10^{6}$ (kPa)
$I$ | Area moment of inertia |  833.3 (mm4)
$\alpha$ | Structural damping parameter |  1.25 $\times$ $10^{-4}$ (s-1)
$\beta$ | Structural damping parameter |  2.5 $\times$ $10^{-4}$ (s)
$\kappa$ | Cubic stiffness coefficient |  50 (N/m3)
$\gamma$ | Cubic damping coefficient |  0.01 (Ns/m3)
Table 2: Computation time for for the stability diagrams in FIG. 5a with $n=10$ DOFs in the format hh:mm:ss. Damping ($\sigma$) | SSM-based ROM | Full system
---|---|---
10 | 00:01:10 | 00:02:07
20 | 00:00:35 | 00:01:42
30 | 00:00:30 | 00:01:31
40 | 00:00:37 | 00:01:21
Table 3: Computation time for the FRCs in FIG. 5b with $n=20$ DOFs in format hh:mm:ss. Damping ($\sigma$) | SSM-based ROM | Full system
---|---|---
10 | 00:14:58 | 06:12:31

(a)

(b)

Figure 5:  (a) Stability diagram for the principal resonance of the first mode
of the Bernoulli beam for 5 discrete elements with $n=10$ degrees of freedom.
The solid line denotes the border of the resonance tongues as obtained using
the dynamics of the full system and has been obtained using po. The reduced
model of these dynamics provided by a cubic truncation of the SSM accurately
locates the position and shape of the tongue for different values of the
damping parameter $\sigma$. (b) Forced response bifurcating out of the trivial
response at the boundary of the stability diagram for $n=20$ DOFs. The
amplitude of the response curve increases with growing parametric excitation
amplitude $\epsilon$.

####  V.2.3 Prismatic Beam

Next we analyze the dynamics of a hinged-clamped prismatic beam. Extending the
initial treatment of external excitation by Nayfeh et al.
Nayfeh1974NonlinearElements , Li et al MingwuPart1 have recently applied SSM-
theory to the case of a 1:3 internal resonance between the first two bending
modes of this beam model. Their implementation serves as the basis of this
example but we choose the nondimensionalized length in terms of the
characteristic length here as $l=1.7$ to avoid internal resonances.
Introducing nondimensional parameters and coordinates, Nayfeh et
al.Nayfeh1974NonlinearElements derive nonlinear partial equations for the
transverse displacement $w(y,t)$ of the beam. Under axial loading at the
hinged tip, these equations must be modified to Barari2011NonlinearBeams

$\displaystyle\frac{\partial^{4}w}{\partial
y^{4}}+\frac{\partial^{2}w}{\partial
t^{2}}=\epsilon\left(H\frac{\partial^{2}w}{\partial
y^{2}}-p_{a}(t)\frac{\partial^{2}w}{\partial y^{2}}\right),$ (55)

where $p_{a}(t)$ describes the axial loading, $H$ accounts for the axial
stretching forces and depends on $w$, and $c$ is the distributed damping
parameter. The axial excitation effectively leads to a parametric excitation
on the transverse displacements of the beam, as is evident in eq. (55). The
beam is clamped at $x=l$ and hinged at $y=0$, resulting in the boundary
conditions

$\displaystyle w(0)=w(l)=0,\quad w^{\prime}(l)=0,\quad w^{\prime\prime}(0)=0.$
(56)

After writing the displacement $w(y,t)$ as a modal expansion, we obtain a
system of second-order ODEs for the evolution of the modal amplitudes
$z_{j}(t)$.

(a)

(b)

Figure 6:  (a) Stability diagram for the principal resonance of the first mode
of the prismatic beam for a truncation with $n=10$ modal coordinates. The full
system response has been obtained with the use of collocation, specifically
po. The reduced model of these full dynamics given by the dynamics on the SSM
can accurately predict the stability behaviour of the full system. (b)
Nontrivial period $T=4\pi/\Omega$ forced response for a number of excitation
amplitudes for a fixed damping coefficient $c=200$. The ROM on the SSM
accurately reproduces the full system response.

$\displaystyle\ddot{z}_{j}+2\epsilon
c\dot{z}_{j}+\omega_{j}^{2}z_{j}=\sum_{i,k,s}\epsilon\alpha_{jiks}z_{i}z_{k}z_{s}+\epsilon\sum_{i}p_{a}(t)z_{i}a_{ji},$
(57)

with $j=1,\dots,n$. A detailed derivation of these equations using Galerkin
projection along with the definition of the coefficients
$\alpha_{jiks},a_{ji}$ are given in Appendix D. For this particular system,
the nonlinear internal forces appear premultiplied by $\epsilon$. Note,
however, that this not a requirement for SSM-theory to be applicable. The
shape of the SSMs depends smoothly on the magnitude of the nonlinearities,
they thus also allow for the analysis of systems that include large
nonlinearities. We take the axial forcing to be harmonic and of the form

$\displaystyle p_{a}(t)=\mu\cos(\Omega t).$ (58)

Furthermore, we set $\epsilon=1\times 10^{-4}$. The modal expansion is
truncated at $n=10$ modes which gives rise to a phase space dimensionality of
$N=20$.

We wish to analyze the stability of the beam when axially excited at the
principal resonance of its fundamental frequency. We construct the SSM over
the first pair of modes with eigenvalues

$\displaystyle\lambda_{1},\overline{\lambda_{1}}\approx-0.0200\pm 5.3361i.$
(59)

The instability of the trivial fixed point induced by the parametric resonance
occurs near the excitation frequency values $\Omega\approx 10.6672\ [rad/s]$.
FIG. 6 shows the corresponding stability diagram and the steady-state response
computed using the SSM-based ROM. The full system’s stability is again
verified using the po-toolbox of coco. The parameters chosen for numerical
continuation of the full and reduced system are documented in Appendix C.

TABLE LABEL:tab:CompTimePB shows the computation times for obtaining the FRCs
as well as the stability diagram for a series of damping parameters. Again,
the use of an SSM-based ROM brings a significant computational advantage over
the full model.

Table 4: Computation time for the prismatic beam example in hh:mm:ss. FIG. | DOFs (n) |  Damping ($c$) | SSM-based ROM | Full system
---|---|---|---|---
6a | 10 | 100 | 00:01:11 | 00:05:41
|  | 200 | 00:00:42 | 00:04:39
|  | 300 | 00:00:38 | 00:04:19
|  | 400 | 00:00:36 | 00:04:02
6b | 10 | 200 | 02:21:09 | 21:03:03

##  VI Forced response and Spectral submanifolds

###  VI.1 Mathematical Foundation

For externally excited systems, SSMs have proven highly effective in capturing
the forced response accurately at a low computational costSSMTool ;
isolatedresponsesten ; Ponsioen2020 ; Jain2021 . Next, we extend these results
to obtain FRCs of systems subject to both external and parametric excitation.
Our algorithm is completely parallelizable and can be carried out for each
forcing frequency independently. We consider periodic external and parametric
excitation and assume that their harmonic $\kappa_{0}$ is in 1:1 resonance
with the master mode, i.e.

$\displaystyle\mathrm{i}\kappa_{0}\Omega\approx\lambda_{j}$ (60)

for some $j$. We further assume that no (near) internal resonances are present
in the spectrum of the dynamical system. The periodic forced response and its
stability can then be explicitly computed from the SSM-based ROM according to
the following statement.

######  Lemma 5.

Assume that the near-resonance condition (60) holds for a mechanical system
subject to periodic external and parametric excitation. We then obtain the
following results for the dynamics on the two-dimensional SSM that serves as
the non-autonomous, nonlinear continuation of the $j^{\mathrm{th}}$ linear
mode:

  1. (i)

Polar dynamics: The reduced dynamics in polar coordinates $(\rho,\theta)$ is
given by

$\displaystyle\begin{bmatrix}\dot{\rho}\\\ \rho\dot{\psi}\end{bmatrix}$
$\displaystyle=\bm{r}(\rho,\psi,\Omega),$ (61) $\displaystyle\dot{\phi}$
$\displaystyle=\Omega,$ (62)

where

$\displaystyle\psi=$ $\displaystyle\ \theta-\kappa_{0}\phi,$ (63)
$\displaystyle\bm{r}(\rho,\psi,\Omega)=$
$\displaystyle\begin{bmatrix}a(\rho)\\\
b(\rho)\end{bmatrix}+\epsilon\bigg{(}\bm{Q}(\psi)\begin{bmatrix}\operatorname{Re}(S^{1}_{\bm{0},\kappa_{0}})\\\
\operatorname{Im}(S^{1}_{\bm{0},\kappa_{0}})\end{bmatrix}$
$\displaystyle+\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{2}\\\
\kappa\in\mathbb{Z}^{*}\end{subarray}}\rho^{|\bm{m}|}\bm{Q}((\kappa/\kappa_{0})\psi)\begin{bmatrix}\operatorname{Re}(S^{1}_{\bm{m},\kappa})\\\
\operatorname{Im}(S^{1}_{\bm{m},\kappa})\end{bmatrix}\bigg{)},$ (64)
$\displaystyle a(\rho)=$ $\displaystyle\
\operatorname{Re}\bigg{(}\sum_{\bm{m}}R^{1}_{\bm{m}}\rho^{|\bm{m}|}\bigg{)},$
(65) $\displaystyle b(\rho,\Omega)=$ $\displaystyle\
\operatorname{Im}\bigg{(}\sum_{\bm{m}}R^{1}_{\bm{m}}\rho^{|\bm{m}|}\bigg{)}-\kappa_{0}\rho\Omega,$
(66) $\displaystyle\bm{Q(\bullet)}:=$ $\displaystyle\
\begin{bmatrix}\cos\bm{(\bullet)}&\sin\bm{(\bullet)}\\\
-\sin\bm{(\bullet)}&\cos\bm{(\bullet)}\end{bmatrix}$ (67)

and the sum over $\kappa$ and $\bm{m}$ includes the harmonics and spatial
multi-indices associated to nonzero reduced dynamics coefficients in a normal-
form style of parametrization.

  2. (ii)

Forced response curve: The fixed points of the SSM-reduced dynamics lie on the
forced response curve that satisfies

$\bm{r}(\rho,\psi,\Omega)=\bm{0}.$ (68)

  3. (iii)

Stability of the forced response: For any hyperbolic fixed point $(\rho,\psi)$
of (61), the eigenvalues of the Jacobian

$\displaystyle J=\begin{bmatrix}\partial_{\rho}a(\rho)&0\\\
\partial_{\rho}(b(\rho)/\rho)&0\end{bmatrix}+\epsilon\bigg{(}\begin{bmatrix}0&d_{\bm{0},\kappa_{0}}(\psi)\\\
-\frac{d_{\bm{0},\kappa_{0}}(\psi)}{\rho^{2}}&-\frac{c_{\bm{0},\kappa_{0}}(\psi)}{\rho}\end{bmatrix}$
$\displaystyle+\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{2}\\\
\kappa\in\mathbb{Z}^{*}\end{subarray}}\begin{bmatrix}c_{\bm{m},\kappa}(\psi)|\bm{m}|\rho^{|\bm{m}|-1}&(\kappa/\kappa_{0})d_{\bm{m},\kappa}(\psi)\rho^{|\bm{m}|}\\\
d_{\bm{m},\kappa}(\psi)(|\bm{m}|-1)\rho^{|\bm{m}|-2}&-(\kappa/\kappa_{0})c_{\bm{m},\kappa}(\psi)\rho^{|\bm{m}|-1}\end{bmatrix}\bigg{)}$
(69)

determine the stability of the fixed point, where

$\displaystyle c_{\bm{m},\kappa}(\psi)$
$\displaystyle:=\operatorname{Re}(S^{1}_{\bm{m},\kappa})\cos((\kappa/\kappa_{0})\psi)$
$\displaystyle\ \
+\operatorname{Im}(S^{1}_{\bm{m},\kappa})\sin((\kappa/\kappa_{0})\psi),$ (70)
$\displaystyle d_{\bm{m},\kappa}(\psi)$
$\displaystyle:=-\operatorname{Re}(S^{1}_{\bm{m},\kappa})\sin((\kappa/\kappa_{0})\psi)$
$\displaystyle\ \
+\operatorname{Im}(S^{1}_{\bm{m},\kappa})\cos((\kappa/\kappa_{0})\psi).$ (71)

###### Proof.

See Appendix B ∎

Fixed points of the reduced dynamics in polar coordinates correspond to
periodic orbits of the full system as the SSM is a time-periodic invariant
manifold, i.e., $\bm{z(t)}=\bm{W}_{\epsilon}(\bm{p}(t),\phi)$. As statement
(ii) of Lemma 5 only requires the computation of zero level sets, the forced
response can be computed independently for each forcing frequency
$\\{\Omega_{1},...\Omega_{z}\\}$ in the interval of interest. Another
consequence is that isolated regions of forced response, the so-called isolas,
can be computed at the same time. This is in stark contrast to continuation-
based methods, which will only find such branches of response from initial
conditions that are specified close to isolas.

We make use Lemma 5 to analyze the forced response of dynamical systems that
are subject to direct and parametric excitation simultaneously. In our
examples, we will consider systems of the form

$\displaystyle\bm{M}\bm{\ddot{y}}+\bm{C}\bm{\dot{y}}+\bigg{(}\bm{K}+\epsilon\mu\bm{Q}\cos(2\Omega
t)\bigg{)}\bm{y}+\bm{f}(\bm{y,\dot{y}})=\epsilon\bm{g}(\bm{y,\dot{y}},\Omega
t),$ (72)

with the function $\bm{g}$ including nonlinear parametric excitation and
external excitation. Solving the invariance equation (17) up to
$\mathcal{O}(\epsilon)$ provides a parametrization of the nonautonomous SSM.
The computation of these parametrizations and the extraction of the FRCs will
be performed with SSMTool 2.4. For a given dynamical system of first or second
order, the algorithm automatically solves the invariance equations (23) and
(30) up to any desired order and computes the steady-states using Lemma 5.

###  VI.2 Examples

####  VI.2.1 Self-excited parametric oscillator

Self-excitation of mechanical systems is caused by negative linear damping
that leads to a destabilization of the trivial fixed point of the linearized
system. In contrast, nonlinear positive damping stabilizes the system for
large oscillation amplitudes and forms a limit cycle, resulting in behavior
described by a Hopf-type normal form Lu2017 . Such self-excitations may arise,
for instance, from dry friction, vibrations during machine cutting, or
nonlinear feedback control in MEMS devices Chatterjee2007 ; Malas2015 ; Fu2011
.

The example we study here represents a 2-DOF oscillator with nonlinear
elasticity and damping, as shown in FIG. 7. This system is subject to external
excitation along one DOF, which is coupled to the second DOF via a spring with
periodically varying stiffness (see, e.g., Warminski and
SzabelskiSzabelski1997 ; Warminski2DOF1995 ). The three types of excitation
seen in this example (external, parametric and self-excitation) are present in
a number of engineered structures, including coupled vibrations of unbalanced
shafts with time-varying rigidity, whirling in slide bearings, and metal
machine cutting. To study the system dynamics, the cited references construct
approximate analytical solutions by assuming a steady state in the form of a
sum of trigonometric functions. Assuming that the oscillation amplitudes vary
slowly in time and, therefore, dropping higher-order derivatives, the authors
arrive at an expression for the amplitudes, which they solve numerically.

Figure 7: Model of a self-excited oscillator system with external and
parametric excitation. The stiffness of the coupling spring varies
periodically as $k(t)=4k(1+\mu\cos(2\Omega t))$. The negative linear damping
coefficient $c$ introduces self-excitation that is stabilized by the nonlinear
damping with coefficient $\gamma$. The left spring is characterized by linear
and nonlinear stiffness coefficients $k$ and $\kappa$.

The equations of motion can be written in the dimensionless form Szabelski1997

$\displaystyle\ddot{y}_{1}+\tilde{k}(1-\tilde{\mu}\cos{2\Omega
t})(y_{1}-y_{2})+y_{1}+\tilde{\kappa}y_{1}^{3}=q\cos{\Omega t},$ (73)
$\displaystyle\ddot{y}_{2}-(\tilde{c}-\tilde{\gamma}\dot{y}_{2}^{2})\dot{y}_{2}-\tilde{k}M(1-\tilde{\mu}\cos{2\Omega
t})(y_{1}-y_{2})=0.$ (74)

The parameters values are chosen as
$\tilde{k}=4,\tilde{\mu}=0.02,\tilde{\kappa}=0.1,\tilde{c}=0.01,\tilde{\gamma}=0.05,M=0.5$
and $q=0.02$ and $\epsilon=1$. Using SSM-theory to obtain the FRC of this
dynamical system involves the following steps.

(a) (b)

Figure 8: Forced response curve of the self-excited oscillator. The forced
response obtained from a fifth-order approximation of the dynamics on the non-
autonomous SSM accurately describes the full system response. Results are
verified using po-toolbox of coco.

  1. 1.

System setup: Using the standardized form (1), the system (73)-(74) can be
rewritten as

$\displaystyle\bm{M}=\begin{pmatrix}1&0\\\
0&1\end{pmatrix},\quad\bm{K}=\begin{pmatrix}\tilde{k}+1&-\tilde{k}\\\
-\tilde{k}M&\tilde{k}M\end{pmatrix},$ (75)
$\displaystyle\bm{Q}=\tilde{k}\begin{pmatrix}-1&1\\\
M&-M\end{pmatrix},\quad\bm{C}=\begin{pmatrix}0&0\\\
0&-\tilde{c}\end{pmatrix},$ (76)
$\displaystyle\bm{f}(\bm{y},\dot{\bm{y}})=\begin{pmatrix}\tilde{\kappa}y_{1}^{3}\\\
\tilde{\gamma}\dot{y}_{2}^{3}\end{pmatrix},\quad\bm{g}(\Omega
t)=q\begin{pmatrix}\cos{\Omega t}\\\ 0\end{pmatrix}.$ (77)

To solve the invariance equations, we expand the nonlinearity and forcing in
Taylor and Taylor-Fourier series. Details on this procedure can be found in
Appendix A.1.

  2. 2.

Choose master subspace: Linear analysis of the corresponding first-order ODE
(6) provides the eigenvalue pairs $\lambda_{1},\overline{\lambda_{1}}\approx
0.0013\pm 2.5887i$ and $\lambda_{2},\overline{\lambda_{2}}\approx 0.0037\pm
0.5463i$. We choose to analyze the resonant excitation of the first mode with
excitation frequency values $\Omega\approx\operatorname{Im}{\lambda_{1}}$ and
construct a time-periodic SSM over the spectral subspace
$\mathcal{E}=\\{\bm{v}_{1},\overline{\bm{v}_{1}}\\}$ .

  3. 3.

Autonomous SSM and reduced dynamics: The leading order autonomous invariance
equation is solved by choosing

$\displaystyle\bm{W}_{\bm{e}_{1}}=\bm{v}_{1},\
\bm{W}_{\bm{e}_{2}}=\overline{\bm{v}_{1}},$ (78)
$\displaystyle\bm{R}_{\bm{e}_{1}}=\begin{pmatrix}\lambda_{1}\\\
0\end{pmatrix},\ \bm{R}_{\bm{e}_{1}}=\begin{pmatrix}0\\\
\overline{\lambda_{1}}\end{pmatrix}.$ (79)

Subsequently, the autonomous invariance equation (1) can be solved recursively
to the desired order. Near-resonances between the master modes of the type
(27) are detected automatically and the autonomous reduced dynamics are
transformed to their normal form accordingly.

  4. 4.

Assemble non-autonomous invariance equations: The autonomous SSM and reduced
dynamics coefficients do not depend on the external excitation frequency. As
this is not the case for non-autonomous coefficients, the nonautonomous
invariance equation has to be solved for every forcing frequency $\Omega$.
Given a set of forcing frequencies $\\{\Omega_{1},\cdots,\Omega_{z}\\}$, the
nonautonomous invariance equation is solved up to the desired order
independently for each frequency. Near-resonances in the master subspace of
the form (35) are automatically detected at each order of computation, which
determines the nontrivial coefficients of the reduced vector field according
to a normal form style of parametrization. The coefficients of the reduced
dynamics are used to construct the function $\bm{r}(\rho,\psi,\Omega)$ defined
in eq. (61), whose zero level set provides us the forced response according to
Lemma 5. The resulting reduced dynamics up to cubic order at frequency
$\Omega=2.6$ rad/s can be written as

$\displaystyle\dot{p}$
$\displaystyle=(1328.6+25887.3i)p+(-35.9+63.2i)p^{2}\bar{p}$ $\displaystyle\
+(7.4+1152.5i)e^{2i\phi}\bar{p}+(0.2-0.8i)e^{-i\phi}p^{2}$ $\displaystyle\
+(-0.4+1.6i)e^{i\phi}p\bar{p}+(1-1.4i)e^{-2i\phi}p^{3},$ (80)

where we have used a rescaled time variable $s=10^{-4}t$, and defined the
phase variable $\phi=\Omega t$. The dynamics along the second reduced
coordinate is simply given by the complex conjugate of (80).

FIG. 8 shows the FRC obtained from this procedure, where we observe a
interesting self-intersecting loop in the FRC. Once again, we verified the
SSM-based reduced results against the full system using the po-toolbox of
coco, where the relevant continuation parameters were set according to TABLE
LABEL:tab:CocoParametersFRC in Appendix C. The computation time for this self-
excited oscillator using a sequential implementation of the non-autonomous SSM
computation is 3 minutes and 48 seconds. Meanwhile, the verification via
continuation on the full system takes 3 minutes and 32 seconds. We observe
that in this sequential implementation on a low-dimensional example, the ROM
does not bring any computational savings relative to the full system in this
low-dimensional example.

Similar oscillator systems featuring self-excitation or parametric
amplification have been shown to exhibit complicated forced response curves,
including isolated branches of response Warminski2DOF1995 ; Li2020 . To
demonstrate the effectiveness of SSM-based ROMs in reproducing such
complicated forced response curves, we have implemented further examples of
self-excited oscillators with isolas Warminski2DOF1995 and parametric
amplifiers with nonlinear damping Li2020 . These examples are openly available
in SSMTool SSMTool2 . Next, we discuss finite-element-type examples under
parametric and external excitation.

####  VI.2.2 Bernoulli beam with parametric and external excitation

We now return to the nonlinear Bernoulli beam example with the parameters
listed in TABLE LABEL:tab:BB_model. In addition to the parametric excitation,
however, the free end is now forced externally as well. The resulting
equations of motion for the discretised variables are now given by

$\displaystyle\bm{M}\ddot{\bm{y}}+\bm{C}\dot{\bm{y}}+\bm{K}\bm{y}+\epsilon\mu\cos
2\Omega
t\bm{Q}\bm{y}+\bm{f}(\bm{y},\dot{\bm{y}})=\epsilon\bm{f}^{0}_{ext}\cos(\Omega
t).$ (81)

Here $\bm{f}(\bm{y},\dot{\bm{y}})$ again includes the nonlinear damping
coefficient $\gamma$ of the damper, as well as the nonlinear stiffness
$\kappa$ of the spring attached to the tip of the beam. $\bm{Q}$ denotes a
projection onto the tip of the beam along the transverse direction. The
harmonic external excitation vector $\bm{f}^{0}_{ext}$ acts on the tip of the
beam, with a unit magnitude in the transverse direction. The damping
coefficient is set to $\sigma=1$. The external excitation amplitude and the
relative parametric excitation amplitude are taken to be $\epsilon=2\times
10^{-3}$ and $\mu=45[N]$, respectively. The system is forced in resonance with
the first mode with the eigenvalue pair

$\displaystyle\lambda_{1},\overline{\lambda_{1}}\approx-0.1238\pm 6.9995i.$
(82)

(a)

(b)

Figure 9: (a) Forced response obtained for the externally excited Bernoulli
beam with periodically varying stiffness actuation at the tip for $n=40$ DOFs.
The nonlinear spring attached to the beam leads to a hardening-type response.
Due to the interplay of external and parametric excitation, an isola arises
for sufficiently large values of the parametric excitation amplitude $\mu$
($\mu=45[N]$ for this simulation). The isola is independently verified using
an initial condition provided by the results of the SSM computation. (b)
Close-up view of the isolated response branch verified against the full system
in a separate continuation run.

The excitation frequencies for parametric and external excitation are assumed
to have a $2:1$ relationship. Therefore, the direct forcing is resonant with
the first eigenmode of the beam, whereas the parametric term is operated in
the regime of principal parametric resonance. FIG.9a shows the forced response
for a 20 element model with $n=40$ DOFs, where we observe an isolated branch
of forced response under the main branch for the chosen parameter values. The
emergence of such an isola with an increasing parametric excitation amplitude
$\mu$ has been studied in systems where parametric resonance and external
excitation interactRhoads2009 ; Aghamohammadi2022DynamicExcitations . Indeed,
this isola is associated to the coexistence of the two resonances, i.e., the
parametric 2:1 resonance and the direct 1:1 resonance of the external
excitation with the eigenvalues of the first mode.

We verify the FRC obtained from the SSM-based ROM against the full system
using the po toolbox of coco with parameters set according to TABLE
LABEL:tab:CocoParametersFRC. The isolated branch of the forced response, which
appears inside the main branch of the FRC, is not detected by the continuation
run. We use an initial solution on the isola provided by the SSM-based ROM as
an input for coco to independently verify the isolated branch for the full
system. Verification of the main branch took a total of 4 hours and 40
minutes, while the isolated response needed 15 hours and 17 minutes. At the
same time, using Lemma 5, the SSM-based ROM takes only 48 seconds for FRC
extraction with a sequential implementation.

Next, we repeat our analysis for larger system sizes resulting from higher-
order discretizations of the beam. The corresponding computation times are
depicted in FIG. 10. All computations for this example were performed on the
Euler Cluster of ETH Zürich due to the long runtimes of the collocation
routine. We increased the system size up to 25,000 elements, resulting in up
to $n=50000$ DOFs. The sequential computation of the FRC, including the
repeated computation of the non-autonomous SSM and its reduced dynamics for
each forcing frequency $\\{\Omega_{1},\cdots,\Omega_{z}\\}$, now takes 1 hour
and 19 minutes.

Figure 10: Computation times for the FRC of the externally excited Bernoulli
beam example. Results from the SSM-based ROMs were verified against the full
system using the po-toolbox of coco (see FIG. 9). Computation for isolas in
the forced response are reported separately as they were verified in
independent continuation runs.

##  VII Conclusion

We have developed an algorithm for computing non-autonomous SSMs and their
reduced dynamics up to any desired order of accuracy in externally and
parametrically forced systems. In doing so, we have taken into account prior
developments and extended them to the case of parametrically excited
mechanical systems with nonlinearities of arbitrary order and magnitude,
including nonlinear damping. We have put forward automated routines for
computing SSMs for second-order mechanical systems in particular and first-
order dynamical systems in general. In our derivations and computations, we
have made use of the multi-index notation for the expansions and carried out
computations in physical coordinates using minimal eigenvectors of the full
system. All this optimizes memory requirements and avoids modal
transformations that are computationally infeasible for high-dimensional
systems.

Using the SSM-reduced models as the basis for parameter continuation, we have
shown how to compute the stability diagrams for the principal resonance of
periodic, parametrically excited systems. Furthermore, we have shown that in
systems under additional external excitation, the forced response can be
obtained by simply locating the zeros of a two-dimensional function, i.e.,
without numerical continuation. As a result, it is now possible to parallelize
the computation of the periodic responses attached to two-dimensional,
nonautonomous SSMs for systems subject to periodic external and parametric
excitation over any desired range of forcing frequencies. All these algorithms
and their applications to our examples are available in the open-source MATLAB
package SSMTool 2.4 SSMTool2 .

With examples treating coupled Mathieu equations, a parametrically excited
Bernoulli beam, and an axially forced prismatic beam, we have illustrated how
SSM-based model reduction can be used to obtain accurate stability diagrams
and FRCs in parametrically excited systems at a significantly reduced
computational cost. Indeed, we analyzed a number of typical low-dimensional
models studied previously in the literature and showed that SSM-based ROMs
provide accurate predictions even if they do not necessarily yield
computational savings. In contrast, SSM-reduced models also provide major
computational advantages over the full system in more complex examples such as
a Bernoulli beam attached to a nonlinear damper and spring, and subject to
external and parametric excitation. In that example, we have uncovered an
internal isola, which is an effect that has been observed for one-dimensional
models subject to parametric resonance and external excitation Rhoads2009 ;
Aghamohammadi2022DynamicExcitations . Ours appears to be the first
demonstration of this phenomenon in finite element models. We showcased the
power of this mathematically grounded SSM-reduction by increasing system sizes
beyond the capabilities of collocation-based continuation methods when applied
to the full system. We also demonstrated how parallelization can be used to
reduce the computation time by an order of magnitude in SSM-based FRC
calculations.

Future work may target systems with even higher dimensionality, where the
beneficial effects of model reduction should weigh in even more stunningly.
Areas of related applications include various types of parametrically driven
MEMS devices, such as comb-driven micromirrors, among others Frangi2017 ;
Polunin2017 ; Welte2013 . Additionally, the developments presented here open
the door for an in-depth analysis of the interaction of more complex time-
dependent phenomena, such as parametric resonance in MEMS systems with
internal resonances Gobat2021ReducedResonance .

## Data availability

The code used to generate the numerical results included in this paper is
available as part of the open-source MATLAB package SSMTool, available at
https://doi.org/10.5281/zenodo.4614201.

###### Acknowledgements.

We declare no conflict of interest. We thank Mingwu Li for his helpful
comments and advice on the continuation toolbox coco.

##  Appendix A SSM Computation in Multi-Index format

In this section of the appendix, we provide the mathematical details for SSM
computations in physical coordinates with the use of multi-indices, including
the proofs for Lemmas 1 \- 4. We begin by reviewing the multi-index formalism
and its advantages. We then give an in-depth derivation of the expressions
that determine the SSM and its reduced dynamics parametrizations. We present
these results for second- and first-order ODEs

###  A.1 Multi-Index calculus

A multi-index of order $m\in\mathbb{N}$ is a vector $\bm{m}\in\mathbb{N}^{M}$
such that $||\bm{m}||_{1}=m$. The unit multi-indices with a single nonzero
entry in the $j^{\mathrm{th}}$ position are denoted by $\bm{e}_{j}$. The
addition and subtraction of multi-indices are defined element-wise. Using
multi-indices, multivariate monomials can be conveniently expressed as

$\bm{p^{m}}=p_{1}^{m_{1}}\cdots p_{M}^{m_{M}}.$ (83)

Likewise, a time derivative of such monomials are conveniently expressed as

$\frac{d\bm{p^{m}}}{dt}=\sum_{j=1}^{M}\partial_{p_{j}}\bm{p^{m}}\dot{p}_{j}=\sum_{j=1}^{M}m_{j}\bm{p}^{\bm{m-e}_{j}}\dot{p}_{j}.$
(84)

#### Composition Coefficients

In SSM computations, composition of nonlinear forces with the nonlinear SSM
parametrization has to be performed. For this, we need to compute scalar
powers $s$ of power series written in multi-index notation, i.e.,
$(\sum_{\bm{m}}W^{i}_{\bm{m}}\bm{p^{m}})^{s}$. A convenient method, which
offers a recursive definition of the resulting coefficients, uses the concept
of radial derivative discussed by Haro et al. Haro2016 . Ponsioen et al.
Ponsioen2020 use this method to express the coefficients $H_{s,\bm{h}}$ in the
expansion

$\big{(}\sum_{\bm{m}}W^{i}_{\bm{m}}\bm{p^{m}}\big{)}^{s}=\sum_{\bm{h}}H^{i}_{s,\bm{h}}\bm{p^{h}}$
(85)

using the recursive definition

$H^{i}_{s,\bm{h}}=\frac{s}{h_{j}}\sum_{\begin{subarray}{c}\bm{m\leq h}\\\
\end{subarray}}m_{j}W^{i}_{\bm{m}}H^{i}_{s-1,\bm{h-m}},$ (86)

where $j$ is the index of the smallest nonzero element in $\bm{h}$. To
simplify the calculation of this recursion, it is useful to look at some
special cases under the assumption that the coefficients $W^{i}_{\bm{m}}$ are
the autonomous coefficients of the $i^{\mathrm{th}}$ coordinate of the SSM
expansion. The SSM parametrization does not contain constant terms, which
means that all composition coefficients are zero for $|\bm{m}|=0$. The
properties of $H^{i}_{s,\bm{h}}$ are then

  * •

$H^{i}_{1,\bm{0}}=0$,

  * •

$H^{i}_{1,\bm{h}}=W^{i}_{\bm{h}}$,

  * •

$H^{i}_{0,\bm{h}}=0$ for all $\bm{h}\neq\bm{0}$,

  * •

$H^{i}_{0,0}=1$.

We generalize this concept to vector-valued power series by collecting all
terms corresponding to the same multi-index $\bm{h}\in\mathbb{N}^{M}$ and
defining a new composition coefficient $H_{\bm{n},\bm{h}}$ such that

$\displaystyle(\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{W}_{\bm{m}}\bm{p}^{\bm{m}})^{\bm{n}}$
$\displaystyle=\prod_{i=1}^{N}(\sum_{\bm{m}\in\mathbb{N}^{M}}W^{i}_{\bm{m}}\bm{p}^{\bm{m}})^{n_{i}}$
(87)
$\displaystyle=\prod_{i=1}^{N}(\sum_{\bm{u}\in\mathbb{N}^{M}}H^{i}_{n_{i},\bm{u}}\bm{p}^{\bm{u}})$
(88)
$\displaystyle=\sum_{\bm{h}\in\mathbb{N}^{M}}H_{\bm{n},\bm{h}}\bm{p}^{\bm{h}}.$
(89)

In more detail:

$\displaystyle
H_{\bm{n},\bm{h}}:=\sum_{\bm{u}_{1}\in\mathbb{N}^{M}}\cdots\sum_{\bm{u}_{2n}\in\mathbb{N}^{M}}H^{1}_{n_{1},\bm{u}_{1}}\cdots
H^{2n}_{n_{2n},\bm{u}_{2n}}\bigg{|}_{\sum\bm{u}_{i}=\bm{h}}.$ (90)

Although a vectorization of this computation provides a very efficient routine
for the composition of power series, it also poses a major bottleneck in terms
of memory and time requirements in the SSM computation proposed here.

#### Redundancy in tensor notation

We now discuss the computational disadvantages of a tensor-based notation
compared to the use of multi-indices. Let $\bm{p}\in\mathbb{C}^{2}$ and
consider the Taylor-expansion of a nonlinear function $\bm{W}(\bm{p})$ in the
tensor notation as

$\displaystyle\bm{W}(\bm{p})=\sum_{i\in\mathbb{N}}\bm{W}_{i}\bm{p}^{\otimes
i}.$ (91)

Specifically, the second term in the expansion (91) can be represented by a
$2\times 4$ matrix $\bm{W}_{2}$ acting on the Kronecker product
$\bm{p}\otimes\bm{p}$ as

$\displaystyle\bm{W}_{2}\bm{p}\otimes\bm{p}=\begin{bmatrix}W_{11}&W_{12}&W_{13}&W_{14}\\\
W_{21}&W_{22}&W_{23}&W_{24}\end{bmatrix}\begin{bmatrix}p_{1}p_{1}\\\
p_{1}p_{2}\\\ p_{2}p_{1}\\\ p_{2}p_{2}\end{bmatrix}$ (92)
$\displaystyle=\begin{bmatrix}W_{11}p_{1}^{2}+(W_{12}+W_{13})p_{1}p_{2}+W_{14}p_{2}^{2}\\\
W_{21}p_{1}^{2}+(W_{22}+W_{23})p_{1}p_{2}+W_{24}p_{2}^{2}\end{bmatrix}$ (93)

Examining Eq. (93), we note that only 6 of the 8 coefficients contained in the
matrix $\bm{W}_{2}$ are needed to uniquely represent
$\bm{W}_{2}\bm{p}\otimes\bm{p}$. Indeed, the same term in multi-index
notation,

$\displaystyle\bm{W}(\bm{p})=\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{W}_{\bm{m}}\bm{p}^{\bm{m}},$
(94)

is characterised by 6 unique multi-index coefficients given as

$\displaystyle\bm{W}_{\big{(}\begin{subarray}{c}2\\\
0\end{subarray}\big{)}}=\begin{bmatrix}W_{11}\\\
W_{21}\end{bmatrix},\quad\bm{W}_{\big{(}\begin{subarray}{c}1\\\
1\end{subarray}\big{)}}=\begin{bmatrix}W_{12}+W_{13}\\\
W_{22}+W_{23}\end{bmatrix},\quad\bm{W}_{\big{(}\begin{subarray}{c}0\\\
2\end{subarray}\big{)}}=\begin{bmatrix}W_{13}\\\ W_{23}\end{bmatrix}$ (95)

This redundancy in the number of coefficients utilized by the tensor notation
relative to the multi-index notations grows with the order of expansion.
Consider, for instance, the case of a 2-dimensional SSM, where, at any order
$\Gamma$, the expansion is characterized by a tensor with $N2^{\Gamma+1}$
coefficients. However, we have only $\Gamma+1$ unique monomials associated
with these coefficients. Hence, while the number of unique monomials grows
linearly with the order $\Gamma$, the corresponding number of tensor
coefficients grows exponentially in $\Gamma$.

We remark, however, that the tensors can be specially constructed to contain
only a minimal number of nonzero coefficients equal to that in the multi-index
notation. This is possible by exploiting the inherent symmetries in monomials
associated with the expansion, e.g., by using a sparse or symmetric definition
of tensor arrays. In practice, the SSM coefficients which are determined
numerically by solving invariance equations may still lead to densely
populated or assymmetric tensor arrays. At each order of computation, the
computed coefficients would, therefore, require recasting to this minimal
form, resulting in unnecessary computational overhead.

Furthermore, it turns out that the invariance equations (1), (34) that need to
be solved for SSM expansion coefficients are decoupled for individual multi-
indices. For computing a two-dimensional autonomous SSM at order $\Gamma$ in a
vectorized tensor-based notation adopted in Ref. Jain2021, , a linear system
of dimension $N2^{\Gamma+1}$ needs to be solved. Using multi-indices instead,
only $\Gamma+1$ linear systems of dimension $n$ each need to be solved. One
can also solve the linear systems associated with each tensor index
independently, as remarked in Ref. Jain2021,  and as done for second-order
systems in Refs. Vizzaccaro2022, ; Opreni2022, . In that case, a total of
$2^{\Gamma+1}$ systems of linear equations of dimension $n$ must be solved in
the tensor-index notation. Hence, for two-dimensional SSMs, the use of multi-
indices leads to an exponential decrease in the number of linear systems that
needs to be solved relative to the tensor notation as the expansion order
increases.

Thus, we have shown the tensor notation result in redundant computational
costs relative to the multi-index notation.

#### Examples of nonlinear force expansion

We expand the excitation and nonlinearity in the multi-index notation as

$\displaystyle\bm{f}(\bm{y},\dot{\bm{y}})$
$\displaystyle=\sum_{\bm{a,b}\in\mathbb{N}^{n}}\bm{f}_{\bm{a,b}}\bm{y^{a}}\dot{\bm{y}}^{\bm{b}},$
(96) $\displaystyle\bm{g}(\bm{y},\dot{\bm{y}},\bm{\phi})$
$\displaystyle=\sum_{\bm{a,b}\in\mathbb{N}^{n}}\sum_{\bm{\kappa}\in\mathbb{Z}^{k}}\bm{g}_{\bm{a,b},\bm{\kappa}}e^{i\langle\bm{\kappa},\bm{\phi}\rangle}\bm{y^{a}}\dot{\bm{y}}^{\bm{b}},$
(97)

where the coefficients
$\bm{f}_{\bm{a,b}},~{}\bm{g}_{\bm{a,b},\bm{\kappa}}\in\mathbb{C}^{n}$. Note
that these coefficients can be recast in a first-order form simply by
associating a multi-index $\bm{n}\in\mathbb{N}^{N}$ to each pair $\bm{a,b}$
such that in the phase space, $\bm{z^{n}}=\bm{y^{a}}\dot{\bm{y}}^{\bm{b}}$.
The first-order version of the two expansions (96), (97) can then be written
as

$\displaystyle\bm{F}(\bm{z})$
$\displaystyle=\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{F}_{\bm{n}}\bm{z}^{\bm{n}},\quad\bm{F}_{\big{(}\begin{subarray}{c}\bm{a}\\\
\bm{b}\end{subarray}\big{)}}=\begin{bmatrix}\bm{f}_{\bm{a},\bm{b}}\\\
\bm{0}_{n}\end{bmatrix},$ (98) $\displaystyle\bm{G}(\bm{z},\bm{\phi})$
$\displaystyle=\sum_{\bm{n}\in\mathbb{N}^{N}}\sum_{\bm{\kappa}\in\mathbb{Z}^{k}}\bm{G}_{\bm{n},\bm{\kappa}}e^{i\langle\bm{\kappa},\bm{\phi}\rangle}\bm{z}^{\bm{n}},\quad\bm{G}_{\big{(}\begin{subarray}{c}\bm{a}\\\
\bm{b}\end{subarray}\big{)},\bm{\kappa}}=\begin{bmatrix}\bm{g}_{\bm{a,b},\bm{\kappa}}\\\
\bm{0}_{n}\end{bmatrix}.$ (99)

As an example of the above notation, we consider the forcing in the self-
excited oscillator system (73)-(74). The forcing term contains external and
linear parametric excitation given by

$\displaystyle\bm{g}(\Omega t,\bm{y})$ $\displaystyle=\cos{\Omega
t}\begin{pmatrix}q\\\ 0\end{pmatrix}+\lambda\cos{2\Omega
t}\begin{pmatrix}y_{1}-y_{2}\\\ -M(y_{1}-y_{2})\end{pmatrix}$ (100)
$\displaystyle=(e^{i\Omega t}+e^{-i\Omega t})\begin{pmatrix}q/2\\\
0\end{pmatrix}$ $\displaystyle+(e^{i2\Omega t}+e^{-i2\Omega
t})\frac{\lambda}{2}\begin{pmatrix}y_{1}-y_{2}\\\
-M(y_{1}-y_{2})\end{pmatrix}.$ (101)

Here, the set of harmonics is given as $K=\\{1,-1,2,-2\\}$. As the system
contains only 2 DOFs, the spatial multi-indices are two dimensional and the
vectors $\bm{g}_{\bm{n},\kappa}$ that characterise the expansion (97) are
given as

$\displaystyle\bm{g}_{\bm{0},1}$
$\displaystyle=\bm{g}_{\bm{0},-1}=\begin{pmatrix}q/2\\\ 0\end{pmatrix},$ (102)
$\displaystyle\bm{g}_{\bm{e}_{1},2}$
$\displaystyle=\bm{g}_{\bm{e}_{1},-2}=\begin{pmatrix}\lambda/2\\\
-M\lambda/2\end{pmatrix},$ (103) $\displaystyle\bm{g}_{\bm{e}_{2},2}$
$\displaystyle=\bm{g}_{\bm{e}_{2},-2}=\begin{pmatrix}-\lambda/2\\\
M\lambda/2\end{pmatrix}.$ (104)

###  A.2 Second order computation

#### Autonomous SSM Computation

In order to compute the SSMs and their reduced dynamics parametrizations for
second-order mechanical systems, we express them as Taylor series (24),(25)
and (26). Substituting these expansions into the autonomous invariance
equation (23) and collecting unique monomial terms, we obtain the system (1)
of the linear equations for each multi-index $\bm{m}$ as we detail below. At
the leading order, i.e, for monomial multi-index $\bm{p}^{\bm{e}_{i}}$, we
obtain

$\displaystyle\begin{bmatrix}\bm{C}&\bm{M}\\\
\bm{M}&\bm{0}\end{bmatrix}\sum_{j=1}^{M}\begin{bmatrix}\bm{w}_{\bm{e}_{j}}\\\
\bm{\dot{w}}_{\bm{e}_{j}}\end{bmatrix}R^{j}_{\bm{e}_{i}}=\begin{bmatrix}-\bm{K}&\bm{0}\\\
\bm{0}&\bm{M}\end{bmatrix}\begin{bmatrix}\bm{w}_{\bm{e}_{i}}\\\
\bm{\dot{w}}_{\bm{e}_{i}}\end{bmatrix}.$ (105)

Eliminating the velocity variables $\bm{\dot{w}}_{\bm{e}_{i}}$ from system
(A.2), we obtain

$\displaystyle\bm{M}\sum_{j=1}^{M}\sum_{k=1}^{M}\bm{w}_{\bm{e}_{k}}R^{k}_{\bm{e}_{j}}R^{j}_{\bm{e}_{i}}+\bm{C}\sum_{j=1}^{M}\bm{w}_{\bm{e}_{j}}R^{j}_{\bm{e}_{i}}+\bm{K}\bm{w}_{\bm{e}_{i}}=\bm{0}.$
(106)

From the eigenvalue problem (12), we note that the choice
$\bm{w}_{\bm{e}_{j}}=\bm{\phi}_{j}$ and
$R^{j}_{\bm{e}_{j}}=\delta_{ij}\lambda_{i}$ solves equation (106). Together
with eq. (105) this solution choice for the leading-order displacement
parametrization implies that
$\bm{\dot{w}}_{\bm{e}_{i}}=\bm{\phi}_{i}\lambda_{i}$. Higher-order terms of
the invariance equation (23) for a multi-index ${\bm{m}}$ are given as

$\displaystyle\bigg{(}\begin{bmatrix}-\bm{K}&\bm{0}\\\
\bm{0}&\bm{M}\end{bmatrix}$
$\displaystyle-\Lambda_{\bm{m}}\begin{bmatrix}\bm{C}&\bm{M}\\\
\bm{M}&\bm{0}\end{bmatrix}\bigg{)}\begin{bmatrix}\bm{w}_{\bm{m}}\\\
\bm{\dot{w}}_{\bm{m}}\end{bmatrix}$
$\displaystyle=\sum_{j=1}^{M}\begin{bmatrix}\bm{C}&\bm{M}\\\
\bm{M}&\bm{0}\end{bmatrix}\begin{bmatrix}\bm{\phi}_{j}\\\
\bm{\phi}_{j}\lambda_{j}\end{bmatrix}R^{j}_{\bm{m}}+\begin{bmatrix}\bm{1}&\bm{0}\\\
\bm{0}&\bm{M}\end{bmatrix}\begin{bmatrix}{\bm{Y}}_{\bm{m}}\\\
{\bm{V}}_{\bm{m}}\end{bmatrix},$ (107)

where we define

$\displaystyle\Lambda_{\bm{m}}$ $\displaystyle:=\bm{\Lambda}\cdot\bm{m},$
(108) $\displaystyle{\bm{Y}}_{\bm{m}}$
$\displaystyle:=\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u,k}\in\mathbb{N}^{M}\\\
1<u<m\\\
\bm{u+k}-\bm{e}_{j}=\bm{m}\end{subarray}}\big{(}\bm{M}\bm{\dot{w}}_{\bm{u}}\big{)}u_{j}R^{j}_{\bm{k}}+\bm{C}\bm{V}_{\bm{m}}-\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{f}_{\bm{n}}H_{\bm{n},\bm{m}},$
(109) $\displaystyle{\bm{V}}_{\bm{m}}$
$\displaystyle:=\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u,k}\in\mathbb{N}^{M}\\\
1<u<m\\\
\bm{u+k}-\bm{e}_{j}=\bm{m}\end{subarray}}\bm{w}_{\bm{u}}u_{j}R^{j}_{\bm{k}}.$
(110)

We rewrite the second set of equations in system (A.2) to obtain an explicit
expression for the velocity parametrization as

$\displaystyle\dot{\bm{w}}_{\bm{m}}=\Lambda_{\bm{m}}\bm{w}_{\bm{m}}+\sum_{j=1}^{M}\bm{\phi}_{j}R^{j}_{\bm{m}}+{\bm{V}}_{\bm{m}},$
(111)

The velocity parametrization depends purely on multiplicative terms of the
displacement parametrization and the reduced dynamics. Indeed, the expression
(111) is identical to the general result of Lemma 2, which we prove later. Now
we utilize relation (111) to rewrite $\bm{Y}_{\bm{m}}$ in terms of the
displacement parametrization and substitute the resulting expression into the
first set of equations in system (A.2) to obtain

$\displaystyle\overbrace{\bigg{(}\bm{K}+\Lambda_{\bm{m}}\bm{C}+(\Lambda_{\bm{m}})^{2}\bm{M}\bigg{)}}^{=:\bm{L}_{\bm{m}}}\bm{w}_{\bm{m}}=\bm{D}_{\bm{m}}\bm{R}_{\bm{m}}+\bm{C_{m}},$
(112)

where $\bm{C_{m}}:=-\Lambda_{\bm{m}}\bm{M}{\bm{V}}_{\bm{m}}-{\bm{Y}}_{\bm{m}}$
and $\bm{D}_{\bm{m}}\in\mathbb{C}^{n\times M}$ is a matrix whose
$j^{\mathrm{th}}$ column is given by
$(\bm{D}_{\bm{m}})_{j}:=-((\Lambda_{\bm{m}}+\lambda_{j})\bm{M}+\bm{C})\bm{\phi}_{j}$
for every $j=1,\dots,M$. This concludes the proof of Lemma 1.

Now we prove Lemma 2 by applying the chain rule to compute the velocity
parametrization and using the fact that the autonomous reduced dynamics
trajectories satisfy $\dot{\bm{p}}=\bm{R(p)}$. Finally, we collect the terms
that sum up to the same multi-index $\bm{m}$ as

$\displaystyle\dot{\bm{w}}_{\bm{m}}=\frac{d\bm{w}}{dt}\Bigg{|}_{\bm{m}}=[D\bm{w}(\bm{p})\dot{\bm{p}}]_{\bm{m}}$
$\displaystyle=[D\bm{w}(\bm{p})\bm{R}(\bm{p})]_{\bm{m}}$ (113)
$\displaystyle=\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u,k}\in\mathbb{N}^{M}\\\
\bm{u+k}-\bm{e}_{j}=\bm{m}\end{subarray}}\bm{w}_{\bm{u}}u_{j}R^{j}_{\bm{k}}$
(114)

Indeed, equation (114) can be verified against the standard invariance
equation given in eq. (111).

##### Reduced dynamics

The linear system (1) contains $n$ equations for each multi-index $\bm{m}$
with $n$ unknowns associated with the displacement parametrization
$\bm{w}_{\bm{m}}$ and additional $M$ unknowns associated with the reduced
dynamics parametrization $\bm{R}_{\bm{m}}$. This underdeterminacy of system
(1) allows different solution choices for the reduced dynamics
parametrization. Specifically, a normal-form style parametrization is
appealing when (near)-resonances are present among the master subspace
eigenvalues. Indeed, for an eigenvalue $\lambda_{i}$ associated to the master
subspace, a near-resonance of the form

$\displaystyle\lambda_{i}\approx\Lambda_{\bm{m}}$ (115)

results in the coefficient matrix $\bm{L_{m}}$ being nearly singular. To
ensure (robust) solvability of system (1), the coefficients $\bm{R}_{\bm{m}}$
must be chosen such that the right-hand side of system (1) is orthogonal to
the (near) kernel of $\bm{L_{m}}$. In particular, we require (see Jain and
HallerJain2021 )

$\displaystyle\bm{\theta}^{*}_{i}\bigg{(}\bm{K}+\Lambda_{\bm{m}}\bm{C}+(\Lambda_{\bm{m}})^{2}\bm{M}\bigg{)}\bm{w}_{\bm{m}}=0,$
(116)

which directly provides us the essentially nonzero reduced dynamics
coefficients $R^{i}_{\bm{m}}$ in the normal-form parametrization style as

$\displaystyle\bm{\theta}^{*}_{i}((\Lambda_{\bm{m}}+\lambda_{i})\bm{M}+\bm{C})\bm{\phi}_{i}R^{i}_{\bm{m}}=-\bm{\theta}^{*}_{i}(\Lambda_{\bm{m}}\bm{M}{\bm{V}}_{\bm{m}}+{\bm{Y}}_{\bm{m}}).$
(117)

#### Non-autonomous SSM

Collecting terms at $\mathcal{O}(\epsilon)$ in the invariance equation (17)
corresponding to a given monomial multi-index $\bm{m}$ and Fourier harmonic
$\bm{\kappa}$, we obtain

$\displaystyle\bigg{(}\begin{bmatrix}-\bm{K}&\bm{0}\\\
\bm{0}&\bm{M}\end{bmatrix}$
$\displaystyle-\Lambda_{\bm{m},\bm{\kappa}}\begin{bmatrix}\bm{C}&\bm{M}\\\
\bm{M}&\bm{0}\end{bmatrix}\bigg{)}\begin{bmatrix}\bm{x}_{\bm{m},\bm{\kappa}}\\\
\bm{\dot{x}}_{\bm{m},\bm{\kappa}}\end{bmatrix}$
$\displaystyle=\sum_{j=1}^{M}\begin{bmatrix}\bm{C}&\bm{M}\\\
\bm{M}&\bm{0}\end{bmatrix}\begin{bmatrix}\bm{\phi}_{j}\\\
\bm{\phi}_{j}\lambda_{j}\end{bmatrix}S^{j}_{\bm{m},\bm{\kappa}}+\begin{bmatrix}\bm{1}&\bm{0}\\\
\bm{0}&\bm{M}\end{bmatrix}\begin{bmatrix}{\bm{Y}}_{\bm{m},\bm{\kappa}}\\\
{\bm{V}}_{\bm{m},\bm{\kappa}}\end{bmatrix},$ (118)

where were we define

$\displaystyle\Lambda_{\bm{m},\bm{\kappa}}:=\Lambda_{\bm{m}}+i\langle\bm{\Omega},\bm{\kappa}\rangle,$
(119)
$\displaystyle{\bm{Y}}_{\bm{m},\bm{\kappa}}:=\bm{M}\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u},\bm{k}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{k}-\bm{e}_{j}=\bm{m}\\\
|\bm{k}|<m\end{subarray}}\big{(}\bm{\dot{w}}_{\bm{u}}u_{j}S^{j}_{\bm{k},\bm{\kappa}}+\bm{\dot{x}}_{\bm{k},\bm{\kappa}}k_{j}R^{j}_{\bm{u}}\big{)}+\bm{C}\bm{V}_{\bm{m},\bm{\kappa}}$
$\displaystyle-\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{f}_{\bm{n}}\sum_{j=1}^{2n}n_{j}\bigg{(}\sum_{\begin{subarray}{c}\bm{u},\bm{h}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{h}=\bm{m}\end{subarray}}H_{\bm{n-\bm{e}_{j}},\bm{h}}X^{j}_{\bm{u},\bm{\kappa}}\bigg{)}-\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{g}_{\bm{n},\bm{\kappa}}H_{\bm{n},\bm{m}},$
(120)
$\displaystyle{\bm{V}}_{\bm{m},\bm{\kappa}}:=\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u},\bm{k}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{k}-\bm{e}_{j}=\bm{m}\\\
|\bm{k}|<m\end{subarray}}\bm{w}_{\bm{u}}u_{j}S^{j}_{\bm{k},\bm{\kappa}}+\bm{x}_{\bm{k},\bm{\kappa}}k_{j}R^{j}_{\bm{u}}.$
(121)

Analogous to the autonomous case (see eq. (111)), an explicit relation between
displacement and velocity parametrization can be established from the
invariance equation (A.2) as

$\displaystyle\bm{\dot{x}}_{\bm{m},\bm{\kappa}}=\Lambda_{\bm{m},\bm{\kappa}}\bm{x}_{\bm{m},\bm{\kappa}}+\sum_{j=1}^{M}\bm{\phi}_{j}S^{j}_{\bm{m},\bm{\kappa}}+{\bm{V}}_{\bm{m},\bm{\kappa}},$
(122)

which is identical to the general result of Lemma 4, as we show later. We
utilize relation (122) to rewrite ${\bm{Y}}_{\bm{m},\bm{\kappa}}$ in terms of
the displacement parametrization and substitute the resulting expression into
the first set of equations in system (A.2) to obtain

$\displaystyle\overbrace{\big{(}\bm{K}+\Lambda_{\bm{m},\bm{\kappa}}\bm{C}+(\Lambda_{\bm{m},\bm{\kappa}})^{2}\bm{M}\big{)}}^{\bm{L}_{\bm{m},\bm{\kappa}}}\bm{x}_{\bm{m},\bm{\kappa}}=\bm{D}_{\bm{m},\bm{\kappa}}\bm{S}_{\bm{m},\bm{\kappa}}+\bm{C}_{\bm{m},\bm{\kappa}},$
(123)

where
$\bm{C}_{\bm{m},\bm{\kappa}}=-\Lambda_{\bm{m},\bm{\kappa}}\bm{M}{\bm{V}}_{\bm{m},\bm{\kappa}}-{\bm{Y}}_{\bm{m},\bm{\kappa}}$
and $\bm{D}_{\bm{m},\bm{\kappa}}\in\mathbb{C}^{n\times M}$ is a matrix whose
$j^{\mathrm{th}}$ column is given by
$(\bm{D}_{\bm{m},\bm{\kappa}})_{j}:=-(\bm{C}+\bm{M}(\lambda_{j}+\Lambda_{\bm{m},\bm{\kappa}}))\bm{\phi}_{j}$
for every $j=1,\dots,M$. This concludes the proof of Lemma 3.

For proving Lemma 4, we use the chain rule using the fact that the non-
autonomous reduced dynamics trajectories satisfy
$\dot{\bm{p}}=\bm{R}_{\epsilon}(\bm{p},\phi)$. To obtain an expression for
$\dot{\bm{x}}$, we collect $\mathcal{O}(\epsilon)$ terms in the time-
derivative of the displacement parametrization of the SSM as

$\displaystyle\dot{\bm{x}}(\bm{p},\bm{\phi})$
$\displaystyle=\frac{d\bm{w}_{\epsilon}}{dt}\Bigg{|}_{\epsilon}=D\bm{w}(\bm{p})\cdot\bm{S}+D_{\bm{p}}\bm{x}\cdot\bm{R}+D_{\bm{\phi}}\bm{x}\cdot\bm{\Omega}.$
(124)

Finally, we restrict this expression to the multi-index $\bm{m}$ and harmonic
$\bm{\kappa}$ to obtain

$\displaystyle\dot{\bm{x}}_{\bm{m},\bm{\kappa}}=\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u},\bm{k}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{k}-\bm{e}_{j}=\bm{m}\end{subarray}}\bm{w}_{\bm{u}}u_{j}S^{j}_{\bm{k},\bm{\kappa}}+\bm{x}_{\bm{k},\bm{\kappa}}k_{j}R^{j}_{\bm{u}}.$
(125)

Indeed, the expression in equation (125) is identical to the one based on
invariance equation computation given in (122).

##### Nonautonomous reduced dynamics

Similarly to the autonomous case, the underdeterminacy of the linear system
(123) allows for different solution choices for the reduced dynamics
parametrization. Specifically, a normal-form style parametrization is
appealing when (near)-resonances are present among the master subspace
eigenvalues and a given harmonic of the forcing frequency. Indeed, for an
eigenvalue $\lambda_{i}$ associated to the master subspace, a near-resonance
of the form

$\displaystyle\lambda_{i}\approx\Lambda_{\bm{m},\bm{\kappa}}$ (126)

results in the coefficient matrix ${\bm{L}_{\bm{m},\bm{\kappa}}}$ being nearly
singular. To ensure (robust) solvability of system (123), the coefficients
$\bm{S}_{\bm{m},\bm{\kappa}}$ must be chosen such that right-hand side of
system (123) is orthogonal to the (near) kernel of
${\bm{L}_{\bm{m},\bm{\kappa}}}$. In particular, we require (cf. Jain and
HallerJain2021 )

$\displaystyle(\bm{\theta}_{i})^{*}\big{(}\bm{K}+\Lambda_{\bm{m},\bm{\kappa}}\bm{C}+(\Lambda_{\bm{m},\bm{\kappa}})^{2}\bm{M}\big{)}\bm{x}_{\bm{m},\bm{\kappa}}=0,$
(127)

which directly provides the essentially non-zero reduced dynamics coefficients
$S^{i}_{\bm{m},\bm{\kappa}}$ in the normal-form parametrization style as

$\displaystyle(\bm{\theta}_{i})^{*}(\bm{C}$
$\displaystyle+\bm{M}(\lambda_{i}+\Lambda_{\bm{m},\bm{\kappa}}))\bm{\phi}_{i}S^{i}_{\bm{m},\bm{\kappa}}$
$\displaystyle=-(\bm{\theta}_{i})^{*}(\Lambda_{\bm{m},\bm{\kappa}}\bm{M}{\bm{V}}_{\bm{m},\bm{\kappa}}+{\bm{Y}}_{\bm{m},\bm{\kappa}})$
(128)

###  A.3 First-order computation

For general first-order dynamical systems the equations do not exhibit the
internal structure (19) seen for second-order mechanical systems (Lemmas 2 &
4). Hence, we develop general expressions for SSM computation in the first-
order system (3). Once again, we operate in physical coordinates using the
multi-index notation.

#### Autonomous SSM

We substitute the autonomous first-order SSM expansion

$\displaystyle\bm{W}(\bm{p})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{W}_{\bm{m}}\bm{p}^{\bm{m}},$
(129)

where $\bm{W}_{\bm{m}}\in\mathbb{C}^{N}$, and its reduced dynamics expansion
(26) into the autonomous invariance equation (23) to obtain

$\displaystyle\bm{B}(\text{D}_{\bm{p}}\sum_{\bm{u}\in\mathbb{N}^{M}}\bm{W}_{\bm{u}}$
$\displaystyle\bm{p}^{\bm{u}})\sum_{\bm{k}\in\mathbb{N}^{M}}\bm{R}_{\bm{k}}\bm{p}^{\bm{k}}=\bm{A}\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{W}_{\bm{m}}\bm{p}^{\bm{m}}$
$\displaystyle+\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{F}_{\bm{n}}(\sum_{\bm{u}\in\mathbb{N}^{M}}\bm{W}_{\bm{u}}\bm{p}^{\bm{u}})^{\bm{n}},$
(130)

where $\bm{F}_{\bm{n}}\in\mathbb{R}^{N}$ are the coefficients in the
nonlinearity expansion (98) in the first-order form. Collecting coefficients
of every unique monomial in system (130), we obtain a linear system for each
monomial multi-index $\bm{m}$, as we detail below. At the leading-order, i.e,
for monomial multi-index $\bm{p}^{\bm{e}_{i}}$, we obtain

$\displaystyle\bm{B}(\sum_{j=1}^{M}\bm{W}_{\bm{e}_{j}})\bm{R}_{\bm{e}_{i}}=\bm{A}\bm{W}_{\bm{e}_{i}}.$
(131)

From the eigenvalue problem (7), we see that the choice
$\bm{W}_{\bm{e}_{i}}=\bm{v}_{i}$ and
$R^{j}_{\bm{e}_{j}}=\delta_{ij}\lambda_{i}$ solves equation (131). It is
noteworthy that due to this diagonal choice of the leading-order autonomous
reduced dynamics, the equations (130) get decoupled between distinct multi-
indices at each order. Hence, collecting the coefficients of higher-order
terms in system (130) for a multi-index ${\bm{m}}$, we obtain

$\displaystyle\underbrace{\bigg{(}\bm{A}-\Lambda_{\bm{m}}\bm{B}\bigg{)}}_{:=\bm{\mathcal{L}}_{\bm{m}}}\bm{W}_{\bm{m}}=\sum_{j=1}^{M}\bm{B}\bm{v}_{j}R^{j}_{\bm{m}}+\bm{\mathcal{C}}_{\bm{m}},$
(132)

where we define

$\displaystyle\bm{\mathcal{C}}_{\bm{m}}$
$\displaystyle:=\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u},\bm{k}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{k}-\bm{e}_{j}=\bm{m}\\\
1<|\bm{u}|<m\end{subarray}}\bm{B}\bm{W}_{\bm{u}}u_{j}R^{j}_{\bm{k}}-\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{F}_{\bm{n}}H_{\bm{n},\bm{m}}.$
(133)

##### Autonomous reduced dynamics

Analogous to the discussion for second-order systems, we observe that linear
system (132) is underdetermined, which allows for choices of parametrizations
for the reduced dynamics. Indeed, for an eigenvalue $\lambda_{i}$ associated
to the master subspace, a near-resonance of the form (115) results in the
coefficient matrix $\bm{\mathcal{L}_{m}}$ being nearly singular. To ensure the
solvability of the system (132), we require (cf. Jain and HallerJain2021 )

$\displaystyle(\bm{u}_{i})^{*}\bm{\mathcal{L}}_{\bm{m}}\bm{W}_{\bm{m}}=\bm{0},$
(134)

which provides the nonzero reduced dynamics coefficients $\bm{R}^{i}_{\bm{m}}$
in the normal-form parametrization style as

$\displaystyle R^{i}_{\bm{m}}=-(\bm{u}_{i})^{*}\bm{\mathcal{C}}_{\bm{m}},$
(135)

where we have chosen a normalization for the eigenvectors that satisfies
$(\bm{u}_{i})^{*}\bm{B}\bm{v}_{j}=\delta_{ij}$ for every $i,j=1,\dots,M$.

#### Non-autonomous SSM

To derive the non-autonomous invariance equation in first-order form, we
expand the $\mathcal{O}(\epsilon)$-terms in the non-autonomous SSM
parametrization (20) as

$\displaystyle\bm{X}(\bm{p},\bm{\phi})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{M}}\sum_{\bm{\kappa}\in\mathbb{Z}^{k}}\bm{X}_{\bm{m},\bm{\kappa}}e^{i\langle\bm{\kappa},\bm{\phi}\rangle}\bm{p}^{\bm{m}},$
(136)

where $\bm{X}_{\bm{m},\bm{\kappa}}\in\mathbb{C}^{N}$. Substituting the
expansions (136) and (33) into the nonautonomous invariance equation (30), and
collecting terms corresponding to any given Fourier-harmonic $\bm{\kappa}$, we
obtain

$\displaystyle\bm{B}\bigg{(}\text{D}_{\bm{p}}(\sum_{\bm{u}\in\mathbb{N}^{M}}\bm{W}_{\bm{u}}\bm{p}^{\bm{u}})\sum_{\bm{k}\in\mathbb{N}^{M}}\bm{S}_{\bm{k},\bm{\kappa}}\bm{p}^{\bm{k}}+(\partial_{\bm{p}}\sum_{\bm{u}\in\mathbb{N}^{M}}\bm{X}_{\bm{u},\bm{\kappa}}\bm{p}^{\bm{u}})\sum_{\bm{k}\in\mathbb{N}^{M}}\bm{R}_{\bm{k}}\bm{p}^{\bm{k}}+i\langle\bm{\kappa},\bm{\Omega}\rangle\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{X}_{\bm{m},\bm{\kappa}}\bm{p}^{\bm{m}}\bigg{)}=$
$\displaystyle\bm{A}\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{X}_{\bm{m},\bm{\kappa}}\bm{p}^{\bm{m}}+\big{[}\text{D}_{\bm{z}}\bm{F}_{\bm{n}}(\sum_{\bm{u}\in\mathbb{N}^{M}}\bm{W}_{\bm{u}}\bm{p}^{\bm{u}})^{\bm{n}}\big{]}\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{X}_{\bm{m},\bm{\kappa}}\bm{p}^{\bm{m}}+\sum_{\bm{n}\in\mathbb{N}^{N}}\sum_{\bm{m}\in\mathbb{N}^{M}}\bm{G}_{\bm{n},\bm{\kappa}}H_{\bm{n},\bm{m}},$
(137)

where $\bm{G}_{\bm{n},\bm{\kappa}}\in\mathbb{C}^{N}$ are the coefficients in
the forcing expansion (99) in the first-order form. Finally, collecting the
coefficients corresponding to any given monomial multi-index $\bm{m}$ in
system (137), we obtain (cf. Poinsioen et al.Ponsioen2020 )

$\displaystyle\underbrace{\bigg{(}\bm{A}-(\Lambda_{\bm{m}}+i\langle\bm{\kappa},\bm{\Omega}\rangle)\bm{B}\bigg{)}}_{:=\bm{\mathcal{L}}_{\bm{m},\bm{\kappa}}}\bm{X}_{\bm{m},\bm{\kappa}}=\sum_{j=1}^{M}\bm{B}\bm{v}_{j}S^{j}_{\bm{m},\bm{\kappa}}+\bm{\mathcal{C}}_{\bm{m},\bm{\kappa}},$
(138)

where

$\displaystyle\bm{\mathcal{C}}_{\bm{m},\bm{\kappa}}$
$\displaystyle:=\bm{B}\bigg{(}\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u},\bm{k}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{k}-\bm{e}_{j}=\bm{m}\\\
|\bm{k}|<m\end{subarray}}\bm{W}_{\bm{u}}u_{j}S^{j}_{\bm{k},\bm{\kappa}}$
$\displaystyle+\sum_{j=1}^{M}\sum_{\begin{subarray}{c}\bm{u},\bm{k}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{k}-\bm{e}_{j}=\bm{m}\\\
|\bm{u}|<m\end{subarray}}\bm{X}_{\bm{u},\bm{\kappa}}u_{j}R^{j}_{\bm{k}}\bigg{)}-\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{G}_{\bm{n},\bm{\kappa}}H_{\bm{n},\bm{m}}$
$\displaystyle-\sum_{\bm{n}\in\mathbb{N}^{N}}\bm{F}_{\bm{n}}\sum_{j=1}^{2n}n_{j}\bigg{(}\sum_{\begin{subarray}{c}\bm{u},\bm{h}\in\mathbb{N}^{M}\\\
\bm{u}+\bm{h}=\bm{m}\end{subarray}}H_{\bm{n-\bm{e}_{j}},\bm{h}}X^{j}_{\bm{u},\bm{\kappa}}\bigg{)}.$
(139)

Recursively solving the linear system (138) provides us the coefficients
$\bm{X}_{\bm{m},\bm{\kappa}},\bm{S}_{\bm{m},\bm{\kappa}}$ defining the
$\mathcal{O}(\epsilon)$-terms in the parametrizations for the nonautonomous
SSM and its reduced dynamics. At the leading order, system (138) simplifies as
(cf. Jain and Haller Jain2021 )

$\displaystyle\bigg{(}\bm{A}-i\langle\bm{\kappa},\bm{\Omega}\rangle\bm{B}\bigg{)}\bm{X}_{\bm{0},\bm{\kappa}}=\sum_{j=1}^{M}\bm{B}\bm{v}_{j}S^{j}_{\bm{0},\bm{\kappa}}+\bm{G}_{\bm{0},\bm{\kappa}}.$
(140)

At any higher order $m=|\bm{m}|>1$, we have $z_{m}=\binom{m+M-1}{M-1}$
distinct multi-indices. An important consequence of choosing the linear
autonomous reduced dynamics to be diagonal, i.e.,
$R^{i}_{\bm{e}_{j}}=\delta_{ij}\lambda_{j}$, is that the system (138) gets
decoupled between distinct multi-indices at any given order $m$ and for any
given harmonic $\bm{\kappa}$. Thus, instead of solving $z_{m}\times N$ coupled
linear equations in general, we need to solve only $z_{m}$ linear systems of
dimension $N$ each.

##### Reduced Dynamics

Analogous to the discussion for second-order nonautonomous setting, we observe
that the linear system (138) is underdetermined, which allows for choices of
parametrizations for the reduced dynamics. Indeed, for an eigenvalue
$\lambda_{i}$ associated with the master subspace, a near-resonance of the
form (35) results in the coefficient matrix
${\bm{\mathcal{L}}_{\bm{m},\bm{\kappa}}}$ being nearly singular.

To ensure (robust) solvability of system (138), the coefficients coefficients
$\bm{S}_{\bm{m},\bm{\kappa}}$ must be chosen such that right-hand side of
system (123) is orthogonal to the (near) kernel of
${\bm{\mathcal{L}}_{\bm{m},\bm{\kappa}}}$. In particular, we require

$\displaystyle(\bm{u}_{i})^{*}\bm{\mathcal{L}}_{\bm{m},\bm{\kappa}}\bm{X}_{\bm{m},\bm{\kappa}}=\bm{0},$
(141)

which directly provides us the essentially non-zero reduced dynamics
coefficients $S^{i}_{\bm{m},\bm{\kappa}}$ in the normal-form parametrization
style as

$\displaystyle
S^{i}_{\bm{m},\bm{\kappa}}=(\bm{u}_{i})^{*}\bm{\mathcal{C}}_{\bm{m},\bm{\kappa}}.$
(142)

Specifically, at the leading-order, we have
$S^{i}_{\bm{0},\bm{\kappa}}=(\bm{u}_{i})^{*}\bm{G}_{\bm{0},\bm{\kappa}}$ (cf.
Jain and HallerJain2021 ) .

##  Appendix B Proof of Lemma 5

In this section we provide the proofs for the statements in Lemma 5. We use a
polar transformation of the parametrization coordinates
$\bm{p}=[p,\overline{p}]=[\rho e^{i\theta},\rho e^{-i\theta}]$. Thus,
$\dot{p}=(\dot{\rho}+\rho i\dot{\theta})e^{i\theta}$. Specifically, the
reduced dynamics (39) over the two-dimensional nonautonomous SSM ($M=2$)
transforms as

$\displaystyle(\dot{\rho}+\rho i\dot{\theta})$
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{l}}R^{1}_{\bm{m}}(\rho
e^{i\theta})^{m_{1}}(\rho e^{-i\theta})^{m_{2}}e^{-i\theta}$ $\displaystyle\ \
\ \ \ +\epsilon\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{l}\\\
\kappa\in\mathbb{Z}\end{subarray}}S^{1}_{\bm{m},\kappa}e^{i\kappa\phi}(\rho
e^{i\theta})^{m_{1}}(\rho e^{-i\theta})^{m_{2}}e^{-i\theta}$ (143)
$\displaystyle=\sum_{\bm{m}\in\mathbb{N}^{l}}R^{1}_{\bm{m}}\rho^{|\bm{m}|}e^{i\theta(m_{1}-m_{2}-1)}$
$\displaystyle\ \ \ \ \
+\epsilon\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{l}\\\
\kappa\in\mathbb{Z}^{*}\end{subarray}}S^{1}_{\bm{m},\kappa}\rho^{|\bm{m}|}e^{i\kappa\phi}e^{i\theta(m_{1}-m_{2}-1)}.$
(144)

Choosing a normal-form style of parametrization (117) for the autonomous
reduced dynamics, we obtain nonzero coefficients $R^{1}_{\bm{m}}$ only for
multi-indices that satisfy the near resonance condition (27). Hence, we have
$m_{1}-m_{2}-1=0$ for every nonzero $R^{1}_{\bm{m}}$. Thus, the angular part
vanishes for all autonomous terms in eq. B.

Furthermore, choosing a normal-form style of parametrization (128) for the
nonautonomous reduced dynamics, we obtain nonzero coefficients
$S^{1}_{\bm{m},{\kappa}}$ only for the pairs $(\bm{m},{\kappa})$ that satisfy
the near-resonance condition (35). But under the external near-resonance
condition (60), for any pair ($\bm{m}$,$\kappa$) that fulfills the near-
resonance conditions (35), we have $\kappa=(1-m_{1}+m_{2})\kappa_{0}$. This
allows us to rewrite
$e^{i\kappa\phi}e^{i\theta(m_{1}-m_{2}-1)}=e^{i(\kappa/\kappa_{0})(\kappa_{0}\phi-\theta)}$
in eq. (B). Now, introducing the phase shift $\psi=\theta-\kappa_{0}\phi$ in
eq. (B) and using the fact that $\dot{\phi}=\Omega$, we obtain

$\displaystyle(\dot{\rho}+\rho
i(\dot{\psi}+\kappa_{0}\Omega))=\sum_{\bm{m}\in\mathbb{N}^{l}}R^{1}_{\bm{m}}$
$\displaystyle\rho^{|\bm{m}|}$ $\displaystyle+\epsilon
S^{1}_{\bm{0},\kappa_{0}}e^{-i\psi}$
$\displaystyle+\epsilon\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{l}\\\
\kappa\in\mathbb{Z}^{*}\end{subarray}}S^{1}_{\bm{m},\kappa}e^{-i(\kappa/\kappa_{0})\psi}.$
(145)

Finally, separating the real and imaginary parts of eq. (145) results in Eq.
(61), which concludes the proof of statement (i). The zeros of the function

$\displaystyle\bm{r}(\rho,\psi,\Omega)=\bm{0}$ (146)

correspond to fixed points of the polar dynamics (61). For the original non-
autonomous reduced dynamics (39), such a fixed point represents a trajectory
with a constant phase shift $\psi$ relative to the forcing as well as a
constant polar amplitude $\rho$. Hence, this fixed point describes a periodic
orbit of the reduced dynamics. As the SSM parametrization maps the orbits of
the reduced system onto the orbits of the full system, we obtain a periodic
steady-state for the full system which concludes the proof of statement (ii).
For computing the stability of these fixed points, we rewrite the polar
dynamics (145) in the form

$\displaystyle\begin{bmatrix}\dot{\rho}\\\
\dot{\psi}\end{bmatrix}=\begin{bmatrix}a(\rho)\\\ b(\rho)/\rho\end{bmatrix}$
$\displaystyle+\epsilon\bigg{(}\frac{\bm{Q}(\psi)}{\rho}\begin{bmatrix}\operatorname{Re}(S^{1}_{\bm{0},\kappa_{0}})\\\
\operatorname{Im}(S^{1}_{\bm{0},\kappa_{0}})\end{bmatrix}$
$\displaystyle+\sum_{\begin{subarray}{c}\bm{m}\in\mathbb{N}^{2}\\\
\kappa\in\mathbb{Z}\end{subarray}}$
$\displaystyle\rho^{|\bm{m}|-1}\bm{Q}((\kappa/\kappa_{0})\psi)\begin{bmatrix}\operatorname{Re}(S^{1}_{\bm{m},\kappa})\\\
\operatorname{Im}(S^{1}_{\bm{m},\kappa})\end{bmatrix}\bigg{)}.$ (147)

Evaluating the Jacobian of the right-hand side of eq. (147) at a hyperbolic
fixed point $(\rho,\psi)$ concludes its stability via linearization. This
proves statement (iii).

##  Appendix C Details on the continuation routine

Within SSMTool 2.4 as well as for the full system analysis, we make use of the
po toolbox provided by the continuation toolbox coco Dankowicz2013 . In this
section, we give details on performing parameter continuation over the SSM-
based ROMs. We also list the continuation parameters used for each example in
this paper.

###  C.1 Sensitivity coefficients

For an efficient numerical continuation of periodic orbits in coco, the
Jacobian of the reduced vector field with respect to the continuation
parameters $\epsilon$ and $\Omega$ and with respect to the parametrization
coordinates $\bm{p}$ should be provided. The Jacobian with respect to the
parametrization coordinates is straightforward, as the vector field is
polynomial in $\bm{p}$. For the derivative with respect to the forcing
amplitude $\epsilon$, only the $\mathcal{O}(\epsilon)$-terms need to be
considered, which yields

$\displaystyle\partial_{\epsilon}\bm{R}_{\epsilon}(\bm{p},\phi)=\bm{S}(\bm{p},\phi).$
(148)

As the nonautonomous coefficients are obtained by solving system (34), they
implicitly depend on the forcing frequency, i.e.,
$\bm{S}_{\mathbf{m},\kappa}=\bm{S}_{\mathbf{m},\kappa}(\Omega)$. Accounting
for this frequency-dependence of the reduced dynamics coefficients, we write
the derivative with respect to $\Omega$ as

$\displaystyle\partial_{\Omega}\bm{R}_{\epsilon}(\bm{p},\bm{\phi})$
$\displaystyle=\epsilon\sum_{\bm{m}\in\mathbb{N}^{M}}\sum_{\kappa\in\mathbb{Z}}(\partial_{\Omega}\bm{S}_{\bm{m},\kappa}e^{i\kappa\Omega
t}+i\kappa t\bm{S}_{\bm{m},\kappa}e^{i\kappa\Omega t}){\bm{p}}^{\bm{m}},$
(149)

where the computation of $\partial_{\Omega}\bm{S}_{\bm{m},\kappa}$ is not
straightforward and we choose to ignore these coefficients in Eq. (149). This
choice may lead to an inaccurate computation of the Jacobian, which may reduce
the rate of convergence to the solution during numerical continuation.
However, this choice does not compromise the solution accuracy, which is
controlled by separate numerical tolerances (see Table
LABEL:tab:CocoParametersDescription). These tolerances are related to the
solution of the invariance equations and to the residual in the Newton
iteration, which we evaluate up to machine precision. Opreni et al. Opreni2022
extend this choice to the SSM and its reduced dynamics parametrizations by
assuming the coefficients $\bm{W}_{\bm{m},\kappa},\bm{S}_{\bm{m},\kappa}$ to
be independent of the forcing frequency $\Omega$, and computing them only once
across the entire frequency range. Such an assumption cannot guarantee
solution accuracy, but offers additional computational advantages. Hence, in
contrast to Opreni et al. Opreni2022 , we choose to solve the nonautonomous
invariance equations for each value of the forcing frequency here.

###  C.2 Period-doubling bifurcations and stability boundary

The po-toolbox of coco, which is used for periodic-orbit continuation in our
SSM-based ROMs as well as in the full system, uses Floquet multipliers to
determine the stability of a periodic orbit and detect bifurcations (if any).
Specifically, a change in the stability of a periodic orbit occurs when a
Floquet multiplier crosses the unit circle in the complex plane, resulting in
a period-doubling bifurcation Dankowicz2013 ; MingwuPart2 . In a system
subject to periodic, parametric excitation, switching to a branch of period-
doubling bifurcations and performing continuation along this branch using coco
then provides us the stability boundary for the trivial periodic response.

###  C.3 Parameters for collocation

The parameters for the continuation of the SSM-based ROMs are presented in
TABLE LABEL:tab:CocoParametersFRCSSM. TABLE LABEL:tab:CocoParametersFRC lists
the collocation and continuation parameters used for various simulations of
the full system. Note that since the level set method of SSMTool 2.4 does not
require continuation and relies solely on the evaluation of the zero level set
of a function, there are no continuation parameters listed for the examples in
Section VI.2. TABLE LABEL:tab:CocoParametersDescription gives a short
description of each parameter along with its default value in coco.

Table 5: Parameters for continuation with po Dankowicz2013 Parameter | Default | Description
---|---|---
NAdapt | 0 | Period for refinement of coll. mesh
h_0 | 0.1 | Initial continuation stepsize
h_max | 0.5 | Max. continuation stepsize
h_min | 0.01 | Min. continuation stepsize
h_fac_max | 2 | Max. factor for apdapting stepsize
h_fac_min | 0.5 | Min. factor for adapting stepsize
MaxRes | 0.1 | Maximal residual norm
bi_direct | 1 | Bidirectional continuation (boolean)
PtMX | 100 | Max. number of continuation steps
al_max | 7 | Max. angle btw. consecutive tangents
ItMX | 10 | Iterations before adapting stepsize
TOL | 1e-6 | Tolerance for Newton Iteration
NTST | 10 | Sub-intervals in coll. discretization
NCOL | 4 | Number of collocation points
Table 6: Collocation and continuation parameters set in the po-toolbox of coco for the SSM-based ROMs in different examples. Figure | 3a | 3b | 5a | 5b | 6a | 6b
---|---|---|---|---|---|---
Parameter | Value |  |  |  |  |
NAdapt | 0 | 0 | 0 | 0 | 0 | 1
h_0 | 0.1 | 1e-4 | 0.1 | 1e-4 | 0.1 | 1e-4
h_max | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 50
h_min | 0.01 | 1e-4 | 0.01 | 1e-4 | 0.01 | 1e-4
bi_direct | 1 | 0 | 1 | 0 | 1 | 0
PtMX | 50 | 35 | 50 | 40 | 50 | 100
Table 7: Collocation and continuation parameters set in the po-toolbox of coco for the full system in different examples. Figure | 3a | 3b | 5a | 5b | 6a | 6b | 8 | 9 |  9 Isola
---|---|---|---|---|---|---|---|---|---
Parameter | Value |  |  |  |  |  |  |  |
NAdapt | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 1
h_0 | 0.1 | 1e-4 | 0.1 | 0.1 | 0.1 | 1e-4 | 0.1 | 0.1 | 0.1
h_max | 0.5 | 0.5 | 0.5 | 0.5 | 0.5 | 50 | 0.5 | 0.5 | 0.5
MaxRes | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1 | 0.1
bi_direct | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 0
PtMX | 10 | 60 | 100 | 60 | 100 | 500 | 350 | 200 | 90
ItMX | 10 | 15 | 10 | 15 | 10 | 15 | 10 | 10 | 10
NTST | 10 | 30 | 10 | 30 | 10 | 70 | 70 | 20 | 80

##  Appendix D Modal equations for the prismatic beam

In this section, we derive the ODEs governing the evolution of the modal
coordinates of a clamped-hinged prismatic beam as a function of the physical
parameters, i.e., elastic modulus $E$, density $\rho$, radius of gyration $r$,
area $A$, and moment of inertia of the cross section $I$, and a characteristic
length $L$. The characteristic length may be the length of the beam or a
typical wavelength for oscillatory modes of the beam. We generally assume that
the beam length is obtained by multiplying the characteristic length $L$ by a
dimensionless constant $l$ as $lL$. At the hinged end, the beam is subject to
an axial load $p_{a}(t)$. The PDE governing the behaviour of the transverse
displacement $w$ of the beam in a nondimensional form is given as
Nayfeh1974NonlinearElements ; Barari2011NonlinearBeams

$\displaystyle\frac{\partial^{4}w}{\partial
y^{4}}+\frac{\partial^{2}w}{\partial t^{2}}$
$\displaystyle=+\epsilon\big{(}H\frac{\partial^{2}w}{\partial
y^{2}}-2c\frac{\partial w}{\partial
t}p_{t}(y,t)-p_{a}(t)\frac{\partial^{2}w}{\partial y^{2}}\big{)},$ (150)
$\displaystyle w(0)$ $\displaystyle=w^{\prime\prime}(0)=w(l)=w^{\prime}(l)=0,$
(151)

where the following dimensionless quantities are used
Nayfeh1974NonlinearElements

$\displaystyle y$ $\displaystyle=\frac{y^{\prime}}{L},$ $\displaystyle t$
$\displaystyle=\sqrt{\frac{Er^{2}}{\rho L^{4}}}t^{\prime},$ $\displaystyle w$
$\displaystyle=\frac{L}{r^{2}}w^{\prime},$ (152) $\displaystyle p$
$\displaystyle=\frac{L^{7}}{r6EA}p^{\prime},$ $\displaystyle c$
$\displaystyle=\frac{L^{4}}{2r^{3}\sqrt{\rho
E}}\frac{r^{2}}{L^{2}}c^{\prime},$ (153)

with $(\bullet)^{\prime}$ denoting the corresponding quantities in physical
units. Specifically, $p^{\prime}$ denotes the distributed load applied along
the beam; $c^{\prime}$ is the distributed damping coefficient;
$y^{\prime},t^{\prime},w^{\prime}$ represent the axial coordinate, the time,
and the transverse displacement. The axial stretching forces due to the
nonlinear bending-stretching coupling are accounted by the term $H$ defined as

$\displaystyle H=\frac{r^{2}}{2L^{2}}\int_{0}^{l}\big{(}\frac{\partial
w}{\partial y}\big{)}^{2}dy.$ (154)

In this example, we choose the slenderness ratio to be
$\frac{r^{2}}{L^{2}}=1\times 10^{-4}$ in line with the original treatment of
Nayfeh et al.Nayfeh1974NonlinearElements , who analyze the prismatic beam
under weak nonlinearities. A modal expansion of the displacement field $w$
yields

$\displaystyle w(y,t)=\sum_{i}\psi_{i}(y)z_{i}(t),$ (155)

where $\psi_{i}(y)$ denotes the $i^{\mathrm{th}}$ eigenfunction associated to
the eigenfrequency $\omega_{i}$ of the linearized version of the PDE (150)
around the trivial solution. These eigenfunctions are solutions of the
eigenvalue problem

$\displaystyle\frac{\partial^{4}\psi_{i}}{\partial y^{4}}$
$\displaystyle=\omega_{i}^{2}\psi_{i},$ (156) $\displaystyle\psi(0)$
$\displaystyle=\psi^{\prime\prime}(0)=\psi(l)=\psi^{\prime}(l)=0,$ (157)

where we choose an orthonormal representation for these eigenfunctions, i.e.,

$\displaystyle\int_{0}^{l}\psi_{i}(y)\psi_{j}(y)=\delta_{ij}dy.$ (158)

Nayfeh et al. Nayfeh1974NonlinearElements have computed explicit expressions
for these eigenfunctions in terms of trigonometric and hyperbolic functions.
Substituting the modal expansion (155) into system (150), and performing a
Galerkin projection, we obtain the set of ODEs governing the evolution of the
modal coordinates $z_{j}(t)$ as

$\displaystyle\omega_{j}^{2}z_{j}+\ddot{z}_{j}+2c\epsilon\dot{z}_{j}=\sum_{i,k,s}\epsilon\alpha_{jiks}z_{i}z_{k}z_{s}+\epsilon\sum_{i}p_{a}(t)z_{i}a_{ji},$
(159)

where

$\displaystyle a_{ji}$ $\displaystyle=\int_{0}^{l}\psi_{j}\
\psi^{\prime\prime}_{i}dy,$ (160) $\displaystyle\alpha_{jiks}$
$\displaystyle=\bigg{(}\frac{r^{2}}{2L^{2}}\int_{0}^{l}\psi^{\prime}_{s}\psi^{\prime}_{k}dy\bigg{)}\bigg{(}\int_{0}^{l}\psi_{j}\psi^{\prime\prime}_{i}dy\bigg{)}.$
(161)

We have assumed the external load to be applied at the hinged tip of the beam
and set the transversal excitation to zero. Hence, the axial forcing
effectively leads to a parametric excitation of the transverse modal
coordinates in system (159).

## References

## References

  * (1) Champneys, A.  Dynamics of parametric excitation.  In _Encyclopedia of Complexity and Systems Science_ (Springer, New York, NY, 2013).  URL https://doi.org/10.1007/978-3-642-27737-5_144-3.
  * (2) Sanches, L., Michon, G., Berlioz, A. & Alazard, D.  Parametrically excited helicopter ground resonance dynamics with high blade asymmetries.  _Journal of Sound and Vibration_ 331, 3897–3913 (2012).
  * (3) Coleman, R. P. & Feingold, A. M.  Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades (1956).  URL https://digital.library.unt.edu/ark:/67531/metadc60767/m1/3/.
  * (4) Ecker, H.  Beneficial Effects of Parametric Excitation in Rotor Systems.  In _IUTAM Bookseries_ , vol. 25, 361–371 (Springer, Dordrecht, 2011).  URL https://link.springer.com/chapter/10.1007/978-94-007-0020-8_31.
  * (5) Tehrani, M. G., Bonisoli, E. & Scapolan, M.  Energy Harvesting Perspectives from Parametric Resonant Systems.  In _Conference Proceedings of the Society for Experimental Mechanics Series_ , vol. 9, 223–232 (Springer, Cham, 2015).  URL https://link.springer.com/chapter/10.1007/978-3-319-15233-2_23.
  * (6) Rhoads, J. F. _et al._ Generalized parametric resonance in electrostatically actuated microelectromechanical oscillators.  _Journal of Sound and Vibration_ 296, 797–829 (2006).
  * (7) Li, L., Hiller, T., Bamieh, B. & Turner, K.  Amplitude control of parametric resonances for mass sensing.  In _Proceedings of IEEE Sensors_ , vol. 2014, 198–201 (Institute of Electrical and Electronics Engineers Inc., 2014).