Low-dimensional reduced models of high-dimensional nonlinear dynamical systems are critically needed in various branches of applied science and engineering. Such simplified models would significantly reduce computational costs and enable physical interpretability, design optimization and efficient controllability. As of yet, however, no generally applicable procedure has emerged for the reliable and robust identification of nonlinear reduced models.
Instead, the most broadly used approach to reducing nonlinear dynamical systems has been a fundamentally linear technique, the proper orthogonal decomposition (POD), followed by a Galerkin projection [1, 2, 3]
. Projecting the full dynamics to the most energetic linear modes, POD requires the knowledge of the governing equations of the system and hence is inapplicable when only data is available. As purely data-based alternatives, machine learning methods are broadly considered and tested in various fields[4, 5, 6, 7]
. While the black-box approach of machine learning might often seem preferable to a detailed nonlinear analysis, the resulting neural network models require extensive tuning, lack physical interpretability, generally perform poorly outside their training range and tend to be unnecessarily complex. This has inspired a number of approaches that seek a blend of machine learning with a priori information about the underlying physics [9, 10]. Still within the realm of machine learning, sparse regression has also shown promise in approximating the right-hand sides of low-dimensional, simple dynamical systems with functions taken from a preselected library . Another recent approach is cluster-based network modeling, which uses the toolkit of network science and statistical physics for modeling nonlinear dynamics .
A popular alternative to POD and machine learning is the dynamic mode decomposition (DMD) , which approximates directly the observed system dynamics. The original DMD and its later variants fit a linear dynamical system to temporally evolving data, possibly including further functions of the original data, over a given finite time interval 
. DMD provides an appealingly simple yet powerful algorithm to infer a local model near steady states where the nonlinear dynamics is always approximately linear. This linear model is also more globally valid if constructed over observables lying in a span of some eigenfunctions of the Koopman operator, which maps observables evaluated over initial states into their evaluations over current states[14, 15, 16]. This relationship between DMD and the Koopman operator has motivated an effort to machine-learn Koopman eigenfunctions from data in order to linearize nonlinear dynamical systems globally on the space of their observables [17, 18, 19].
Finding physically realizable observables that fall in a Koopman eigenspace is, however, often described as challenging or difficult
. A more precise assessment would be that such a find is highly unlikely, given that the probability of any countable set of a priori selected observables falling in any Koopman eigenspace is zero. In addition, those eigenspaces can only be determined explicitly in simple, low-dimensional systems. In practice, therefore, DMD can only provide a justifiable model near an attracting fixed point of a dynamical system. While Koopman modes still have the potential to linearize the observer dynamics on larger domains, those domains cannot include more than one attracting or repelling fixed point[20, 21, 19]. Indeed, DMD and Koopman mode expansions fail to converge outside neighborhoods of fixed points even in the simplest, one-dimensional nonlinear systems with two fixed points [22, 20]. In summary, while these data-driven model reduction methods are powerful and continue to inspire ongoing research, their applicability is limited to locally linearized systems and globally linearizable nonlinear systems, such as the Burgers equation .
The focus of this paper is the development of data-driven, simple and predictive reduced-order models for essentially nonlinear dynamical systems, i.e., non-linearizable systems. Determining exact linearizability conclusively from data is beyond reach. In contrast, establishing that a dynamical system is non-linearizable in a domain of interest is substantially simpler: one just needs to find an indication of coexisting isolated stationary states in the data. By an isolated stationary state, we mean here a compact and connected invariant set with an open neighborhood that contains no other compact and connected invariant set. Examples of such stationary states include hyperbolic fixed points, periodic orbits, invariant spheres and quasiperiodic tori; closures of homoclinic orbits and heteroclinic cycles; and chaotic attractors and repellers. If a data set indicates the coexistence of any two sets from the above list, then the system is conclusively non-linearizable in the range of the available data. Specifically, there will be no homeomorphism (continuous transformation with a continuous inverse) that transforms the orbits of the underlying dynamical system into those of a linear dynamical system. While this is a priori clear from dynamical systems theory, several studies have specifically confirmed a lack of convergence of Koopman-mode expansions already for the simplest case of two coexisting fixed points, even over subsets of their domain of attraction or repulsion [22, 20].
Non-linearizable systems are ubiquitous in science, technology and nature. Beyond the well-known examples of chaotic dynamical systems and turbulent fluid flows , any bifurcation phenomenon, by definition, involves coexisting steady states and hence is automatically non-linearizable. Indeed, aerodynamic flutter , buckling of beams and shells , bistable microelectromechanical systems , traffic jams  or even tipping points in climate change  are all fundamentally non-linearizable, just to name a few. Figure 1 shows some examples of non-linearizable systems emerging in technology, nature and scientific modeling.
We will show here that a collection of classic and recent mathematical results from nonlinear dynamical systems theory enables surprisingly accurate and predictive low-dimensional modeling from data for a number of non-linearizable phenomena. Our construct relies on the recent theory of spectral submanifolds (SSMs), the smoothest invariant manifolds that act as nonlinear continuations of non-resonant eigenspaces from the linearization of a system at a stationary state (fixed point, periodic orbit or quasiperiodic orbit; ). Using appropriate SSM embeddings [30, 31, 32] and an extended form of the classic normal form theory , we obtain sparse dynamical systems describing the reduced dynamics on the slowest SSMs of the system, which are normally hyperbolic and hence robust under perturbations .
We construct the extended normal form within the slowest SSM as if the eigenvalues of the linearized dynamics within the SSM had zero real parts, although that is not the case. As a result, our normalization procedure will not render the simplest possible (linear) normal form for the SSM dynamics, valid only near the underlying isolated stationary state. Instead, our procedure yields a sparsified nonlinear, polynomial normal form on a larger domain of the SSM that can also capture nearby coexisting stationary states. This fully data-driven normalization algorithm learns the normal form transformation and the coefficients of the normal form simultaneously by minimizing an appropriately defined conjugacy error between the unnormalized and normalized SSM dynamics.
For a generic observable of an oscillatory dynamical system without an internal resonance, a two-dimensional data-driven model calculated on the slowest SSM of the system turns out to capture the correct asymptotic dynamics. Such an SSM-reduced model is valid on domains in which the nonlinearity and any possible external forcing are strong enough to create non-linearizable dynamics, yet are still moderate enough to render the eigenspace of the linear system relevant.
More generally, oscillatory systems with independent internal resonances in their spectrum can be described by reduced models on -dimensional SSMs. In both the resonant and the non-resonant cases, the models can be refined by increasing the degree of their nonlinearity rather than by increasing their dimension. As we show in examples, the resulting SSM-based models are explicit, deterministic and even have the potential to predict system behavior outside the range of the training data away from bifurcations. Most importantly, we find that the models also accurately predict forced response, even though they are only trained on data collected from unforced systems.
We illustrate the power of data-driven SSM-reduced models on high-dimensional numerically generated data sets and on experimental data. These and further examples are also available as MATLAB®
live scripts, which are part of a general open-source package,SSMLearn, that performs this type of model reduction and prediction for arbitrary data sets.
2.1 Spectral submanifolds and their reduced dynamics
A recent result in dynamical systems is that all eigenspaces (or spectral subspaces) of linearized systems admit unique nonlinear continuations under well-defined mathematical conditions. Specifically, spectral submanifolds (SSMs), as defined by , are the unique smoothest invariant manifolds that serve as nonlinear extensions of spectral subspaces under the addition of nonlinearities to a linear system. The SSM formulation and terminology we use here is due to ; the Methods section A.1 discusses the history of these results and further technical details.
We consider -dimensional dynamical systems of the form
with a constant matrix and with class functions and , where is the
-dimensional torus. The elements of the frequency vectorare rationally independent, and hence the function is quasiperiodic in time. The assumed degree of smoothness for the right-hand side of (1) is , with referring to analytic. The small parameter signals that the forcing in system (1) is moderate so that the structure of the autonomous part is still relevant for the full system dynamics. Rigorous mathematical results on SSMs are proven for small enough , but continue to hold in practice for larger values of as well, as we will see in examples. Note that eq. (1) describes equations of motions of physical oscillatory systems. It does not cover phenomenological models of phase oscillators, such as the Kuramoto model .
The eigenvalues of , with multiplicities counted, are ordered based on their real parts, , as
Their corresponding real modal subspaces (or eigenspaces),
, are spanned by the imaginary and real parts of the corresponding eigenvectors and generalized eigenvectors of. To analyze typical systems, we assume that holds for all eigenvalues, i.e., is a hyperbolic fixed point for .
A spectral subspace is a direct sum
of an arbitrary collection of modal subspaces, which is always an invariant subspace for the linear part of the dynamics in (1). Classic examples of spectral subspaces are the stable and unstable subspaces, comprising all modal subspaces with and , respectively. Projections of the linearized system onto the nested hierarchy of slow spectral subspaces,
provide exact reduced-order models for the linearized dynamics over an increasing number of time scales under increasing , as sketched in Fig. 2a. This is why a Galerkin projection onto is an exact model reduction procedure for linear systems, whose accuracy can be increased by increasing . A fundamental question is whether nonlinear analogues of spectral subspaces continue to organize the dynamics under the addition of nonlinear and time-dependent terms in the full system (1).
Let us fix a specific spectral subspace within either the stable or the unstable subspace. If is non-resonant (i.e., no nonnegative, low-order, integer linear combination of the spectrum of is contained in the spectrum of outside ), then has infinitely many nonlinear continuations in system (1) for small enough . These invariant manifolds are of smoothness class , with the spectral quotient measuring the ration of the fastest decay exponent outside to the slowest decay exponent inside (see eq. ((13)) of the Methods section A.1). All such manifolds are tangent to for , have the same quasiperiodic time dependence as does and have a dimension equal to that of .
Of these infinitely may invariant manifolds, however, there will be a unique smoothest one, the spectral submanifold (SSM) of , denoted . This manifold is smooth if and can therefore be approximated more accurately than the other infinitely many nonlinear continuations of . In particular, SSMs have convergent Taylor expansions if the dynamical system (1) is analytic (). Then the reduced dynamics on a slow SSM, , can be approximated with arbitrarily high accuracy using arbitrarily high-order Taylor expansions, without ever increasing the dimension of (see Fig. 2b). Such an approximation for dynamical systems with known governing equations is now available for any required order of accuracy via the open source MATLAB® package SSMTool . In contrast, reduced models obtained from projection-based procedures can only be improved by increasing their dimensions.
The nearby coexisting stationary states in Fig. 2 happen to be contained in the SSM. In specific examples, however, these states may also be off the SSM, contained instead in one of the infinitely many additional nonlinear continuations, , of the spectral subspace . The Taylor expansion of the dynamics on and are, however, identical up to order . Therefore, the reduced models we will compute on the SSM also correctly capture the nearby stationary states on , as long as the polynomial order of the model stays below . In large physical systems, this represents no limitation, given that .
2.2 Embedding SSMs via generic observables
If at least some of the real parts of the eigenvalues in (2) are negative, then longer-term trajectory data for system (1) will be close to an attracting SSM, as illustrated in Fig. 2b. This is certainly the case for data from experiments that are run until a nontrivial, attracting steady state emerges (see, e.g., Fig. 1e). Measurements of trajectories in the full phase space, however, are seldom available from such experiments. Hence, if data about system (1) is only available from observables, the construction of SSMs and their reduced dynamics has to be carried out in the space of those observables.
An extended version of Whitney’s embedding theorem guarantees that almost all (in the sense of prevalence) smooth observable vectors provide an embedding of a compact subset of a -dimensional SSM, , into the observable space for high enough . Specifically, if we have simultaneous and independent continuous measurements, , of the observables, then almost all maps are embeddings of , and hence the top right plot of Fig. 3 is applicable with probability one.
In practice, we may not have access to independent observables and hence cannot invoke Whitney’s theorem. In that case, we invoke the Takens delay embedding theorem , which covers observable vectors built from uniformly sampled, consecutive measured instances of a single observable. More precisely, if is a generic scalar quantity measured at times apart, then the observable vector for delay-embedding is formed as . We discuss the embedding, , of an autonomous SSM, , in the observable space in more detail in the Methods section A.2.
2.3 Data-driven extended normal forms on SSMs
Once the embedded SSM, , is identified in the observable space, we seek to learn the reduced dynamics on . An emerging requirement for learning nonlinear models from data has been model sparsity , without which the learning process would be highly sensitive. The dynamics on , however, is inherently non-sparse, which suggests that we learn its Poincaré normal form  instead. This classic normal form is the simplest polynomial form to which the dynamics can be brought via successive, near-identity polynomial transformations of increasing order.
Near the origin on a slow SSM, however, this simplest polynomial form is just the restriction of the linear part of system (1) to , as long as infinitely many non-resonance conditions are satisfied for the operator . The Poincaré normal form on would, therefore, only capture the low-amplitude, linearized part of the slow SSM dynamics.
To construct an SSM-reduced model for non-linearizable dynamics, we use extended normal forms. This idea is motivated by normal forms used in the study of bifurcations of equilibria on center manifolds depending on parameters [33, 41]. In that setting, the normal form transformation is constructed at the bifurcation point where the system is non-linearizable by definition. The same transformation is then used away from bifurcations, even though the normal form of the system would be linear there. One, therefore, gives up the maximal possible simplicity of the normal form but gains a larger domain on which the normal form transformation is invertible and hence captures truly nonlinear dynamics. In our setting, there is no bifurcation at , but we nevertheless construct our normal form transformation as if the eigenvalues corresponding to the slow subspace were purely imaginary. This procedure leaves additional, near-resonant terms in the SSM-reduced normal form, enhancing the domain on which the transformation is invertible and hence the normal form is valid.
We determine the normal form coefficients directly from data via the minimization of a conjugacy error (see the Methods section). This least-square minimization procedure renders simultaneously the best fitting normal form coefficients and the best fitting normal form transformation. As we will find in a specific example, this data-driven procedure can yield accurate reduced models even beyond the formal domain of convergence of equation-driven normal forms.
The simplest extended normal form on a slow SSM of an oscillatory system arises when the underlying spectral subspace corresponds to a pair of complex conjugate eigenvalues. Writing in polar coordinates and truncating at cubic order,  finds this normal form on the corresponding two-dimensional, autonomous SSM, , to be
The dynamics of (5) is characteristically non-linearizable when , given that a limit cycle coexists with the fixed point in that case. Further coexisting steady states will arise when forcing is added to the system, as we discuss in the next section. We note that the cubic normal form on two-dimensional SSMs has also been approximated from data in . That non-sparse procedure fits the full observer dynamics to a low-dimensional, discrete polynomial dynamical system, then performs an analytic SSM reduction and a classic normal form transformation on the SSM.
For higher accuracy, the extended normal form on an oscillatory SSM of dimension is of the form
If the linearized frequencies are non-resonant, then the functions and only depend on . Our numerical procedure determines these functions up to the necessary order that ensures a required accuracy for the reduced-order model on the SSM. This is illustrated schematically for a four-dimensional slow SSM ( in the bottom right plot of Fig. 3.
2.4 Predicting forced dynamics from unforced data
With the normalized reduced dynamics (6) on the embedded SSM, , at hand, we can also make predictions for the dynamics of the embedded quasiperiodic SSM, , of the full system (1). This forced SSM is guaranteed to be an -close perturbation of for moderate external forcing amplitudes. A strict proof of this fact is available for small enough , but as our examples will illustrate, the smooth persistence of the SSM, , generally holds for all moderate values in practice. Such moderate forcing is highly relevant in a number of technological settings, including system identification in structural dynamics and fluid-structure interactions, where the forcing must be moderate to preserve the integrity of the structure.
We discuss the general extended normal form on in the Methods section A.3. In the simplest and most frequent special case, the external forcing is periodic () and is the embedding of the slowest, two-dimensional SSM corresponding to a pair of complex conjugate eigenvalues. Using the modal forcing amplitude and modal phase shift in the general normal form (25),  introduces the new phase coordinate and lets , , to obtain the planar, autonomous, extended normal form on as
at leading order in . All stable and unstable periodic responses on the SSM are fixed points of system (7), with their amplitudes and phases satisfying the equations
The first analytic formula in (8) predicts the forced response curve (FRC) of system (1), i.e., the relationship between response amplitude, forcing amplitude and forcing frequency, from the terms and of the extended normal form of the autonomous SSM, . These terms are constructed from trajectories of the unforced system, thus eq. (8) predicts the behavior of a non-linearizable dynamical system under forcing based solely on unforced training data. The stability of the predicted periodic response follows from a simple linear analysis at the corresponding fixed point of the ODE (7). The first formula in (8) also contains another frequently used notion of nonlinear vibration analysis, the dissipative backbone curve , which describes the instantaneous amplitude-frequency relation along freely decaying vibrations within the SSM.
As we will also show in examples, our unforced model-based predictions for forced periodic response (see the Methods section A.4) are confirmed by numerical simulation or dedicated laboratory experiments on forced systems.
We now illustrate data-driven, SSM-based modeling and prediction on several numerical and experimental data sets describing non-linearizable physical systems. Further applications are described in . Both the numerical and the experimental data sets were initialized without knowledge of the exact SSM. All our computations have been carried out by the publicly available MATLAB® package, SSMLearn, whose repository also contains further examples not discussed here. The main algorithm behind SSMLearn is illustrated in Fig. 3, with more detail given in the Methods section A.5.
To quantify the errors of an SSM-based reduced model, we use the normalized mean-trajectory-error (). For observations of the observable vector and their model-based reconstructions, , this modeling error is defined as
Here is a relevant normalization vector, such as the data point with the largest norm. When validating the reduced dynamics for a given testing trajectory, we run the reduced model from the same initial condition for the comparison. Increasing the order of the truncated normal form polynomials in eq. (6) generally reduces the error to any required level but excessively small errors can lead to overfitting. In our examples, we will be allowing model errors in the order of to avoid overfitting.
2.5.1 Finite-element model of a damped-forced beam
(a). In contrast to the classic Euler-Bernoulli beam, the von Kármán model captures moderate deformations by including a nonlinear, quadratic term in the kinematics. We first construct a 33 degree-of-freedom, damped, unforced finite element model (i.e.,and in eq. (1)) for an aluminum beam of length 1 [m], width 5 [cm], thickness 2 [cm] and material damping modulus 10 [ ] (see the Supplementary Information for more detail).
Our objective is to learn from numerically generated trajectory data the reduced dynamics on the slowest, two-dimensional SSM, , of the system, defined over the slowest two-dimensional () eigenspace of the linear part. To do so, we generate two trajectories starting from initial beam deflections caused by static loading of 12 [kN] and 14 [kN] at the midpoint, as shown in Fig. 4(a). The latter trajectory, shown in Fig. 4(b), is used for training, the other for testing. Along the trajectories, we select our single observable to be the midpoint displacement of the beam.
The beam equations are analytic , and hence the SSM, , admits a convergent Taylor expansion near the origin. The minimal embedding dimension for the two-dimensional, , as required by Whitney’s theorem, is , which is not satisfied by our single scalar observable We therefore employ delay-embedding using with [ms]. By Takens’s theorem, this delayed observable embeds the SSM in with probability one.
A projection of the embedded SSM, onto three coordinates is shown in Fig. 4(c). On , SSMLearn returns the -order extended normal form
to achieve our preset reconstruction error bar of on the test trajectory (), shown in Fig. 4(d).
We now use the model (10), trained on a single decaying trajectory, to predict the forced response of the beam for various forcing amplitudes and frequencies in closed form. We will then compare these predictions with analytic forced response computations for the forced SSM, , obtained from SSMTool  and with numerical simulations of the damped-forced beam. The periodic forcing is applied at the midpoint node; the Taylor expansion order in SSMTool for the analytically computed dynamics on is set to , as in (10). Figure 4(e) shows the FRCs (green) and the backbone curve (blue) predicted by SSMLearn based on formula (8) from the single unforced trajectory in Fig. 4(b). To obtain the relevant forcing amplitudes in the delay-observable space, we have followed the calibration procedure described in the Methods section A.4 for the forcing values [N] at the single forcing frequency [Hz]. Recall that coexisting stable (solid lines) and unstable (dashed lines) periodic orbits along the same FRC are hallmarks of non-linearizable dynamics and hence cannot be captured by the model reduction techniques we reviewed in the Introduction for linearizable systems.
The data-based prediction for the FRCs agrees with the analytic FRCs for low forcing amplitudes but departs from it for higher amplitudes. Remarkably, as the numerical simulations (red) confirm, the data-based FRC is the correct one. The discrepancy between the two FRCs for large amplitudes only starts decreasing under substantially higher-order Taylor series approximations used in SSMTool (see the Supplementary Information). This suggests the use of the data-based approach for this class of problems even if the exact equations of motion are available.
2.5.2 Vortex-shedding behind a cylinder
Next we consider the classic problem of vortex shedding behind a cylinder . Our input data for SSM-based reduced modeling are the velocity and pressure fields over a planar, open fluid domain with a hole representing the cylinder section, as shown in Fig. 5(a). The boundary conditions are no-slip on the circular inner boundary, standard outflow on the outer boundary at the right side, and fixed horizontal velocity on the three remaining sides . The Reynolds number for this problem is the ratio between the cylinder diameter times the inflow velocity and the kinematic viscosity of the fluid.
Available studies [50, 51, 8] report that, at low Reynolds number, the two-dimensional unstable manifold, , of the wake-type steady solution, , in Fig. 5(b) connects to the limit cycle shown in Fig. 5(c). Here we evaluate the performance of SSMLearn on learning this unstable manifold as an SSM, along with its reduced dynamics, from trajectory data at Reynolds number equal to 70. For this SSM, we again have and , as in our previous example. There is no external forcing in this problem, and hence we have in eq. (1). In contrast to prior studies that often consider a limited number of observables [51, 52, 8], here we select the full phase space of the discretized Navier-Stokes simulation to be the observable space for illustration, which yields in eq. (1). We generate nine trajectories numerically, eight of which will be used for training and one for testing the SSM-based model.
The nine initial conditions of our input trajectory data are small perturbations from the wake-type steady solution along its unstable directions, equally spaced on a small amplitude circle on this unstable plane. All nine trajectories quickly converge to the unstable manifold and then to the limit cycle representing periodic vortex shedding.
We choose to parametrize the SSM, , with two leading POD modes of the limit cycle, which have been used in earlier studies for this problem. The training trajectories projected onto these two POD modes are shown in Fig. 5(d). To limit the modeling error (9) to less than , SSMLearn requires a polynomial order of 18 in the SSM computations. For this order, our approach can accommodate the strong mode deformation observed for this problem , manifested by a fold of the SSM over the unstable eigenspace in Fig. 5(f). Figure 5(g) shows the strongly nonlinear geometry of projected to the observable subspace formed by the velocities and the pressure of a probe point in the wake.
To capture the SSM-reduced dynamics with acceptable accuracy, we need to compute the extended normal form up to order 11, obtaining
To describe a transition qualitatively from a fixed point to a limit cycle, the reduced-order dynamical model should be at least of cubic order . Capturing the qualitative behavior (i.e., the unstable fixed point and the stable limit cycle), however, does not imply a low error for the model. Indeed, the data-driven cubic normal form for this example gives a reconstruction error of normalized over the limit cycle amplitude, mainly arising from an out-of-phase convergence to the limit cycle along the testing trajectory. In contrast, the normal form in eq. (11) reduced this error drastically to on the testing trajectory, as shown in Fig. 5(e).
We show in Section 1.2.3 of the Supplementary Information that for comparable accuracy, the Sparse Identification of Nonlinear DYnamics (SINDy) approach of  returns non-sparse nonlinear models for this example. Similarly, while the DMD  can achieve highly accurate curve-fitting on the available training trajectories with a high-dimensional linear model, that model only captures linearizable dynamics near the origin. As a consequence, its trajectories grow without bound over longer integration times and hence fail to capture the limit cycle.
2.5.3 Fluid sloshing experiments in a tank
Fluid oscillations in a tank exhibit highly nonlinear characteristics . To describe such non-linearizable softening effects observed in the sloshing motion of surface waves, Duffing-type models have been proposed . While amplitude variations observed in forced experiments can be fitted to forced softening Duffing equations, nonlinear damping remains a challenge to identify .
The experiments we use to construct an SSM-reduced nonlinear model for sloshing were performed in a rectangular tank of width 500 [mm] and depth 50 [mm], partially filled with water up to a height of 400 [mm], as shown in Fig. 6(a). The tank was mounted on a platform excited harmonically by a motor. The surface level was detected via image processing from a monochrome camera. As an observable we used the horizontal position of the computed center of mass of the water at each time instant, normalized by the tank width. This physically meaningful scalar is robust with respect to image evaluation errors .
We identify the unforced nonlinear behavior of the system from data obtained in resonance decay experiments . In those experiments (as in Fig. 6a, but with a shaker instead of a motor), once a periodic steady state is reached under periodic horizontal shaking of the tank, the shaker is turned off and the decaying sloshing is recorded. We show such a decaying observable trajectory (orange line) in Fig. 6(b), with the shaker switched off slightly before zero time. This damped oscillation is close, by construction, to the two-dimensional, slowest SSM of the system. We use three such decaying observer trajectories (two for training and one for model testing) for the construction of a two-dimensional (), autonomous, SSM-based reduced order model for . For delay embedding dimension, we again pick , the minimal value guaranteed to be generically correct for embedding the SSM by Takens’s theorem. The delay used in sampling is [s]. For this input and for a maximal reconstruction error of , SSMLearn identifies a nearly flat SSM in the delayed observable space–see Fig. 6(c)–with a cubic extended normal form
This lowest-order, Stuart–Landau-type normal form (see (5)) already constitutes an accurate reduced order model with on the testing data set (see Fig. 6(b)). The amplitude-dependent nonlinear damping, , provided by this model is plotted in Fig. 6(d) with respect to the physical amplitude.
In another set of experiments with the setup of Fig. 6a, steady states of periodically forced sloshing were measured in sweeps over a range of forcing frequencies under three different shaker amplitudes. As in the previous beam example, we identify the corresponding forcing amplitude, , in (7) at the maximal amplitude response of each frequency sweep. Shown in Fig. 6(e,f), the closed-form predictions for FRCs from eq. (8) (solid lines) match closely the experimental FRCs (dots). Given the strong nonlinearity of the FRC, any prediction of this curve from a DMD-based model is bound to be vastly inaccurate, as we indeed show in Section 1.3 of the Supplementary information.
The phase of the forced response relative to the forcing has been found difficult to fit to forced Duffing-type models , but the present modeling methodology also predicts this phase accurately using the second expression in (8). The blue curve in Fig. 6(d) shows the backbone curve of decaying vibrations, which terminates at the highest amplitude occurring in the training data set. This plot therefore shows that the closed-form FRC predictions obtained from the SSM-based reduced model are also effective for response amplitudes outside the training range of the reduced model.
We have described a data-driven model reduction procedure for non-linearizable dynamical systems with coexisting isolated stationary states. Our approach is based on the recent theory of spectral submanifolds (SSMs), which are the smoothest nonlinear continuations of spectral subspaces of the linearized dynamics. Slow SSMs form a nested hierarchy of attractors and hence the dynamics on them provide a hierarchy of reduced order models with which generic trajectories synchronize exponentially fast. These SSMs and their reduced models smoothly persist under moderate external forcing, yielding low-dimensional, mathematically exact reduced order models for forced versions of the same dynamical system. The normal hyperbolicity of SSMs also ensures their robustness under small noise.
All these results have been implemented in the open source MATLAB® package, SSMLearn, which we have illustrated on data sets arising from forced nonlinear beam oscillations, vortex shedding behind a cylinder and water sloshing in a vibrating tank. For all three examples, we have found that two-dimensional data-driven extended normal forms on the slowest SSMs provide sparse yet accurate models of non-linearizable dynamics in the space of the chosen observables. Beyond matching training and testing data, SSM-reduced models prove their intrinsic, qualitative meaning by predicting non-linearizable, forced steady states purely from decaying, unforced data.
In this brief report, examples of higher-dimensional SSMs and multi-harmonic forcing have not been considered, even though SSMLearn is equipped to handle them. Higher-dimensional SSMs are required in the presence of internal resonances or in non-resonant problems in which initial transients also need to be captured more accurately. A limitation of our approach for non-autonomous systems is the assumption of quasiperiodic external forcing. Note, however, that even specific realizations of stochastic forcing signals can be approximated arbitrarily closely with quasiperiodic functions over any finite time interval of interest. A further limitation in our work is the assumption of smooth system dynamics. For data from non-smooth systems, SSMLearn will nevertheless return an equivalent smooth reduced-order model whose accuracy is a priori known from the available mean-squared error of the SSM fitting and conjugacy error of the normal form construction. We are addressing these challenges in ongoing work to be reported elsewhere. Further applications of SSMLearn to physical problems including higher-dimensional coexisting steady states (see, e.g., ) are also underway.
Appendix A Methods
All data discussed in the results presented here is publicly available in the SSMLearn repository at github.com/haller-group/SSMLearn.
The code supporting the results presented here is publicly available in the SSMLearn repository at github.com/haller-group/SSMLearn.
a.1 Existence of SSMs
In the context of rigid body dynamics, invariant manifolds providing generalizations of invariant spectral subspaces to nonlinear systems were first envisioned and formally constructed as nonlinear normal modes by  (see  for a recent review of related work). Later studies, however, pointed out the non-uniqueness of nonlinear normal modes in specific examples ([60, 61]).
In the mathematics literature, 
obtained general results on the existence, smoothness and degree of uniqueness of such invariant manifolds for mappings on Banach spaces. These results use a special parameterization method to construct the manifolds even in evolutionary partial differential equations that admit a well-posed flow map in both time directions (see for a mechanics application). The results have been extended to a form applicable to dynamical systems with quasiperiodic time dependence . An extensive account of the numerical implementation of the parametrization method with a focus on computing invariant tori and their whiskers in Hamiltonian systems is also available .
 discussed the existence of the SSM, , depending on its absolute spectral quotient,
where is the stable (unstable) spectrum of if the SSM is stable (unstable). For a stable SSM, is the integer part of the quotient of the minimal real part in the spectrum of and the maximal real part of the spectrum of restricted to .
Based on , we call a -dimensional spectral subspace non-resonant if for any set of nonnegative integers satisfying , the eigenvalues, , of satisfy
This condition only needs to be verified for resonance orders between and . In particular, a resonance between and is allowed if , in which case each strongly resonant spectral subspace gives rise to a unique nearby spectral submanifold.
If violates the non-resonance condition (14), then can be enlarged to a higher-dimensional spectral subspace until the non-resonance relationship (14) is satisfied. In the absence of external forcing (), the non-resonance condition (14) can also be relaxed with the help of the relative spectral quotient,
to the form
This is indeed a relaxation because condition (16) is only violated if both the real and the imaginary parts of eigenvalues involved are in the exact same resonance with each other. In contrast, (14) is already violated when the real parts are in resonance with each other.
If in eq. (2) and all subspaces are non-resonant, then the nested set of slow spectral submanifolds,
gives a hierarchy of local attractors. All solutions in a vicinity of approach the reduced dynamics on one of these attractors exponentially fast, as sketched in Fig. 2b for the limit. As we will see, non-linearizable dynamics tend to emerge on due to near-resonance between the linearized frequencies within and the forcing frequencies . The specific location of nontrivial steady states in is then determined by a balance between the nonlinearities, damping and forcing.
A resonant subspace can be enlarged by adding the next modal subspaces to it until in the hierarchy (4) becomes non-resonant and hence admits an SSM, . This technical enlargement is also in agreement with the physical expectation that all interacting modes have to be included in an accurate reduced-order model. Finally, we note that SSMs are robust features of dynamical systems: they inherit smooth dependence of the vector field in (1) on parameters .
For discrete-time dynamical systems of the form
We close by noting that in a neighborhood of an SSM, an invariant family of surfaces resembling the role of coordinate planes in a linear system exists . This invariant spectral foliation (ISF) can, in principle, be used to generate a nonlinear analogue of linear modal superposition in a vicinity of a fixed point. Constructing the ISF from data has shown both initial promise and challenges to be addressed.
a.2 Embedding the SSM in the observable space
Originally conceived for autonomous systems, the Takens delay embedding theorem  has been strengthened and generalized to externally forced dynamics . By these results, the embedding for a -dimensional compact SSM subset, , in the delay observable space as is guaranteed for almost all choices of the observable if and some generic assumptions regarding periodic motions on are satisfied .
Of highest importance in technological applications is the case of time-periodic forcing (), with frequency and period . In this case, the Whitney and Takens embedding theorems can be applied to the associated period- sampling map (or Poincaré map) of the system based at time . This map is autonomous and has a time-independent SSM that coincides with the -dimensional SSM, , of the full system (1). In this case, by direct application of the embedding theorems to the discrete dynamical system generated by
, the typically sufficient embedding dimension estimate is improved tofor Whitney’s and Takens’s theorem.
Technically speaking, the available data will never be exactly on an SSM, as these embedding theorems assume. By the smoothness of the embeddings, however, points close enough to the SSM in the phase space will be close to in the observable space under the embeddings. Moreover, as slow SSMs attract nearby trajectories exponentially, the distance of observable data from the embedded slow SSM will shrink exponentially fast. Therefore, even under uncorrelated noise in the measurements, mean-squared estimators are suitable for learning slow SSMs from data in the observable space, as we illustrate in the Supplementary Information.
After a possible coordinate shift, the trivial fixed point of the autonomous limit of system (1) will be mapped into the origin of the observable space. To find an embedded, -dimensional SSM, , attached to this origin for , we focus on observable domains in which is a graph over its tangent space at the origin . Such domains always exist and are generally large enough to capture non-linearizable dynamics in most applications (but see below). Note that coincides with the image of the spectral subspace in the observable space.
To learn such a graph-style parametrization for from data, we define a matrix
with columns that are orthonormal vectors spanning the yet unknown. The reduced coordinates for a point are then defined as the orthogonal projection . We week a Taylor-expansion for near the origin, denoting by the family of all monomials of variables from degree 2 to . For example, if and , then . As a graph over , the manifold is approximated as , where the matrices and contain coefficients for the -variate linear and nonlinear monomials, respectively. Learning from a data set of observations then amounts to finding the matrices that minimize the mean-square reconstruction error along the training data:
The simplest solution to this problem is with the additional constraint
, which represents a basic nonlinear extension of the principal component analysis.
The above graph-style parametrization of the SSM breaks down for larger values if develops a fold over . That creates an issue for model reduction if a nontrivial steady state on falls outside the fold, as the limit cycle does in our vortex shedding example. In that case, alternative parametrization methods for can be used to enhance the domain of the SSM-reduced model. These methods include selecting the columns of to be the leading POD modes of the nontrivial steady state, or enlarging the embedding space with (further) delayed observations. In these cases, the columns of are still orthonormal vectors spanning
In both Fig. 4 and Fig. 4, the SSM, , is nearly flat in the delay-embedding space. This turns out to be a universal property of delay embedding for small delays and low embedding dimensions (see the Supplementary Information).
For small (i.e., for moderate forcing), the autonomous SSM, , already captures the bulk nonlinear behavior of system (1). Indeed, for this forcing range, the reduced dynamics on the corresponding SSM can simply be computed as an additive perturbation of the autonomous dynamics on [47, 68, 69] (see section 2.4).
a.3 SSM dynamics via extended normal forms
For an autonomous SSM , the reduced dynamics is governed by a vector field
with a flow map . We can generically assume that the Jacobian is semisimple, i.e., , where is a diagonal matrix containing the eigenvalues of . Classic normal form theory would seek to simplify the reduced dynamics (19) in a vicinity of via a nonlinear change of coordinates, , so that the transformed vector field with flow map has a diagonal linear part and has as few nonlinear terms in its Taylor expansion as possible. In our present setting, the origin is assumed hyperbolic, in which case the classic normal form is simply under appropriate non-resonance conditions that generically hold . The corresponding normal form transformation , however, is only valid on a small enough domain in which the dynamics is linearizable.
To capture non-linearizable behavior, we employ extended normal forms motivated by those used to unfold bifurcations . In this approach, we construct normal forms that do not remove those polynomial terms from (19) whose removal would result in small denominators in the Taylor coefficients and hence decrease its domain of convergence. Instead, we seek a normal form for (19) of the form
where the matrices , and contain the coefficients for the appropriate -variate monomials. To identify near-resonances, we let be the matrix of integers whose columns are the powers of the -variate monomials from order 2 to . We then define a matrix containing all relevant integer linear combinations of eigenvalues as follows:
Following the approach used in universal unfolding principles , we collect in a set the row and column indices of the entries of for which near-resonances occur, i.e., for which the corresponding entry of is smaller in norm than a small, preselected threshold. (The default threshold is in SSMLearn.) The entries of and with indices contained in are then set to zero but the corresponding monomial terms are retained in . Conversely, coefficients of non-near-resonant entries of and are selected in a way so that the corresponding non–near-resonant monomials vanish from the normal form . As a result, the matrix is sparse, containing only the coefficients of essential, near-resonant monomials.
For example, if , and the eigenvalues of form a complex pair with , then we have
Only two elements of are (near-) zero, and hence the reduced dynamics in extended normal form will require learning the following coefficients:
The corresponding cubic polar form (5) is then obtained from the relations and .
For a data-driven construction of the extended normal form (20), we first obtain an estimate for the Jacobian
from linear regression. This determines the matrixand the types of monomials arising in and . Next, we note that the flow map of the SSM-reduced dynamics and the flow map of the extended normal form are connected through the conjugacy relationship . We find the nonzero complex coefficients of and by minimizing the error in this exact conjugacy over the available data points, represented in the coordinates. Specifically, we determine the nonzero elements of and as
Once is known, we obtain the coefficients of via regression.
As initial condition for the minimization problem (24), we set all unknown coefficients to zero. This initial guess assumes linear dynamics, which the minimization corrects as needed. We can compute the time derivative in (24) reliably using finite differences, provided that the sampling time of the trajectory data is small compared to the fastest timescale of the SSM dynamics. For larger sampling times, one should use the discrete formulation of SSM theory, as discussed in section A.1 and . In that formulation, the conjugacy error must be formulated for the 1-step prediction error of the normal form flow map . The matrix defined in eq. (21) also carries over to the discrete time setting, with defined as the diagonal matrix of the logarithms of the eigenvalues of .
a.4 Prediction of forced response from unforced training data
Forced SSMs continue to be embedded in our observable space, provided that we also include the phase of the forcing among our observables . (In the simplest case of periodic forcing, this inclusion is not necessary, as we pointed out Section 2.2). The quasiperiodic SSM-reduced normal form of system (1) in the observable space takes the general form
where the terms and are the forcing amplitudes and phases for each mode of the SSM and for each forcing harmonic , while are the set containing the indexes of the resonant forcing frequencies for mode (see the Supplementary Information). The normal form (25) will capture non-linearizable dynamics arising from resonant interactions between the eigenfrequencies of the spectral subspace (which may also contain internal resonances) and the external forcing frequencies in . One can use numerical continuation  to find nontrivial co-existing steady states (such as periodic orbits and invariant tori) in eq. (25) under varying forcing amplitudes and forcing frequencies.
To predict forced response from the SSM-based model trained on unforced data, the forcing amplitude relevant for eq. (7) in the observable space needs to be related to the forcing amplitude relevant for system (1) in the physical phase space. This involves (1) employing a single forcing amplitude-frequency pair in the experiment (2) measuring the periodic observable response (3) computing the corresponding normalized reduced and normalized response amplitude (4) substituting into the first formula in (8) and (5) solving for in closed form. This can then be used to make a prediction for the full FRC and response phase via (8) in the experiment for arbitrary forcing frequencies. The predicted FRC may have several connected components, including isolated responses (isolas) that are notoriously difficult to detect by numerical or experimental continuation .
a.5 Summary of the algorithm
The data-driven model reduction method used in this paper is available in the open-source MATLAB® package SSMLearn. User input is the measured trajectory data of the autonomous dynamical system (), the SSM dimension , the polynomial orders or approximation for the SSM and for the extended normal form, as well as the type of the dynamical system (discrete or continuous). If the number of observables is not sufficient for manifold embedding, the data is automatically augmented with delays to reach the minimum embedding dimension . If the manifold learning returns poor results (due to, e.g., insufficient closeness of the data to the SSM), then the starting value of can be increased until a good embedding is found. Then, the algorithm learns the SSM geometry in observable space and, after unsupervised detection of the required normal form, identifies the extended normal form of the reduced dynamics. The level of accuracy can be increased with larger polynomial orders, keeping in mind that excessive orders may lead to overfitting.
SSMLearn also offers all the tools we have used in this paper to analyze the reduced dynamics and make predictions for forced response from unforced training data. In particular, it contains the MATLAB®-based numerical continuation core coco . which can compute steady state and help with the design of nonlinear control strategies. In principle, there are no restrictions on the dimensions of the reduced order model, yet the larger the SSM is, the more computationally expensive the problem becomes.
Qualitative or partial a priori knowledge of the linearized dynamics (e.g., some linearized modes and frequencies) helps in finding good initial conditions for trajectories to be used in SSMLearn. For example, the resonance decay method  (which we exploited in our sloshing example), targets a specific 2-dimensional, stable SSM in laboratory experiments. This method consists of empirically isolating a resonant periodic motion on the SSM based on its locally maximal amplitude response under a forcing frequency sweep. Discontinuing the forcing will then generate transient decay towards the equilibrium in a close proximity of the SSM. For noisy data, filtering or dimensionality reduction can efficiently de-noise the data , provided that the polynomial orders used for the description of the SSM and its reduced dynamics are not excessively large (see the Supplementary Information). For higher-dimensional SSMs, it is desirable to collect diverse trajectories to avoid bias towards specific motions. Good practice requires splitting the data sets into training, testing and validation parts.
-  P.J. Holmes, J.L. Lumley, G. Berkooz, and C.W. Rowley. Turbulence, Coherent Structures, Dynamical Systems and Symmetry. Cambridge Monographs on Mechanics. Cambridge University Press, 2 edition, 2012.
-  J. Awrejcewicz, V.A. Krys’ko, and Vakakis A.F. Order Reduction by Proper Orthogonal Decomposition (POD) Analysis, pages 279–320. Springer, Berlin, Heidelberg, 2004.
-  K. Lu, Y. Jin, Y. Cehn, Y. Yang, L. Hou, Z. Zhang, Z. Li, and C. Fu. Review for order reduction based on proper orthogonal decomposition and outlooks of applications in mechanical systems. Mech. Sys. and Signal Proc., 123:264–297, 2019.
-  S.L. Brunton, J.L. Proctor, and J.N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
-  K.S. Mohamed. Machine Learning for Model Order Reduction. Springer, Cham, 2018.
-  T. Daniel, F. Casenave, N. Akkari, and D. Ryckelynck. Model order reduction assisted by deep neural networks (ROM-net). Adv. Model. and Simul. in Eng. Sci., 7:105786, 2020.
-  M. Calka, P. Perrier, J. Ohayon, C. Grivot-Boichon, M. Rochette, and Y. Payan. Machine-learning based model order reduction of a biomechanical model of the human tongue. Computer Methodes and Programs in Biomedicine, 198:105786, 2021.
-  J.-C. Loiseau, S.L. Brunton, and B.R. Noack. From the POD-Galerkin method to sparse manifold models, pages 279–320. De Gruyter, Berlin, 2020.
-  G.E. Karniadakis, I.G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Physics-informed machine learning. Nature Rev. Phys., 123:422–440, 2021.
S. Li and Y. Yang.
Data-driven identification of nonlinear normal modes via physics-integrated deep learning.Nonlinear Dyn., 106:3231–3246, 2021.
-  D. Fernex, B.R. Noack, and R. Semaan. Cluster-based network modeling–From snapshots to complex dynamical systems. Science Advances, 7(25):eabf5006, 2021.
-  P.J. Schmid. Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech., 656:5–28, 2010.
-  J.N. Kutz, S.L. Brunton, B.W. Brunton, and J.L. Proctor. Dynamic Mode Decomposition. SIAM, Philadelphia, PA, 2016.
-  C.W. Rowley, I. Mezić, S. Bagheri, P. Schlachter, and D.S. Henningson. Spectral analysis of nonlinear flows. J. Fluid Mech., 641:115–127, 2009.
-  I. Mezić. Analysis of fluid flows via spectral properties of the Koopman operator. Ann. Rev. Fluid Mech., 45(1):357–378, 2013.
-  A. Mauroy, I. Mezić, and Y. Susuki. The Koopman Operator in Systems and Control Concepts, Methodologies, and Applications: Concepts, Methodologies, and Applications. Springer, New York, 2020.
-  B. Lusch, J.N. Kutz, and S.L. Brunton. Deep learning for universal linear embeddings of nonlinear dynamics. Nature Commun., 9:4950:1–10, 2018.
S.E. Otto and C.W. Rowley.
Linearly recurrent autoencoder networks for learning dynamics.SIAM J. Applied Dynamical Systems, 18(1):558–593, 2019.
-  E. Kaiser, J.N. Kutz, and S.L. Brunton. Data-driven discovery of koopman eigenfunctions for control. Mach. Learn.: Sci. Technol., 2:035023, 2021.
-  J. Page and R.R. Kerswell. Koopman mode expansions between simple invariant solutions. J. Fluid Mech., 879:1–27, 2019.
-  S.L. Brunton, B.W. Brunton, J.L. Proctor, and J.N. Kutz. Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control. PLOS ONE, 11:1–19, 2016.
-  S. Bagheri. Koopman-mode decomposition of the cylinder wake. J. Fluid Mech., 726:596–623, 2013.
-  J. Page and R.R. Kerswell. Koopman analysis of Burgers equation. Phys. Rev. Fluids, 3:071901, 2018.
-  E.H. Dowell. Panel flutter - a review of the aeroelastic stability of plates and shells. AIAA Journal, 8(3):385–399, 1970.
-  A. Abramian, E. Virot, E. Lozano, S.M. Rubinstein, and T.M. Schneider. Nondestructive prediction of the buckling load of imperfect shells. Phys. Rev. Lett., 125:225504, 2020.
-  P. Podder, D. Mallick, A. Amann, and S. Roy. Influence of combined fundamental potentials in a nonlinear vibration energy harvester. Sci Rep, 6:37292, 2016.
-  G. Orosz and G. Stépán. Subcritical hopf bifurcations in a car-following model with reaction-time delay. Proc. Royal Soc. A, 462(2073):2643–2670, 2006.
-  P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox. Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system. Phil. Trans. Royal Soc. A, 370(1962):1166–1184, 2012.
-  G. Haller and S. Ponsioen. Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction. Nonlinear Dyn., 86(3):1493–1534, 2016.
-  H. Whitney. The self-intersections of a smooth n-manifold in 2n-space. Ann. of Math., 45(2):220–246, 1944.
-  J. Stark, D.S. Broomhead, M.E. Davies, and J. Huke. Takens embedding theorems for forced and stochastic systems. Nonlinear Analysis: Theory, Methods and Applications, 30(8):5303–5314, 1997.
-  J. Stark. Delay embeddings for forced systems. I. Deterministic forcing. J. Nonlinear Sci., 9:255–332, 1999.
-  J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields. Springer, New York, 1983.
-  N. Fenichel. Persistence and smoothness of invariant manifolds for flows. Indiana University Mathematics Journal, 21(3):193–226, 1971.
-  Y. Kuramoto. Chemical Oscillations, Waves and Turbulence. Springer, Berlin, 1984.
-  S. Jain and G. Haller. How to compute invariant manifolds and their reduced dynamics in high-dimensional finite-element models? Nonlinear Dyn., 2021.
-  T. Sauer, J.A. Yorke, and M. Casdagli. Embedology. J. Stat. Phys., 65:579–616, 1997.
-  F. Takens. Detecting strange attractors in turbulence. In D. Rand and L. Young, editors, Dynamical Systems and Turbulence, Warwick 1980, pages 366–381. Springer Berlin Heidelberg, 1981.
-  H. Poincaré. Les Méthodes Nouvelles de la Mécanique Céleste. Gauthier-Villars et Fils, Paris, 1892.
-  S. Sternberg. On the structure of local homeomorphisms of euclidean n-space, II. Am. J. Math., 80(3):623–631, 1958.
-  J. Murdock. Normal Forms and Unfoldings for Local Dynamical Systems. Springer Monographs in Mathematics. Springer-Verlag New York, 2003.
-  S. Ponsioen, T. Pedergnana, and G. Haller. Automated computation of autonomous spectral submanifolds for nonlinear modal analysis. J. Sound and Vibration, 420:269–295, 2018.
-  L.D. Landau. On the problem of turbulence. Dokl. Akad. Nauk SSSR, 44(8):339–349, 1944.
-  J.T. Stuart. On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. Part 1. The basic behaviour in plane Poiseuille flow. Journal of Fluid Mechanics, 9(3):353–370, 1960.
-  K. Fujimura. Centre manifold reduction and the Stuart-Landau equation for fluid motions. Proceedings: Mathematical, Physical and Engineering Sciences, 453(1956):181–203, 1997.
-  R. Szalai, D. Ehrhardt, and G. Haller. Nonlinear model identification and spectral submanifolds for multi-degree-of-freedom mechanical vibrations. Proc. Royal Society A, 473(2202):20160759, 2017.
-  T. Breunung and G. Haller. Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems. Proc. Royal Soc. A, 474:20180083, 2018.
-  M. Cenedese, J. Axås, H. Yang, M. Eriten, and G. Haller. Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems. arXiv:2110.01929, 2021.
-  S. Jain, P. Tiso, and G. Haller. Exact nonlinear model reduction for a von Kármán beam: slow-fast decomposition and spectral submanifolds. Journal of Sound and Vibration, 423:195–211, 2018.
-  D. Barkley and R.D. Henderson. Three-dimensional Floquet stability analysis of the wake of a circular cylinder. Journal of Fluid Mechanics, 322:215–241, 1996.
-  B.R. Noack, K. Afanasiev, M. Morzyński, G. Tadmor, and F. Thiele. A hierarchy of low-dimensional models for the transient and post-transient cylinder wake. Journal of Fluid Mechanics, 497:335–363, 2003.
-  C.W. Rowley and S.T.M. Dawson. Model reduction for flow analysis and control. Annual Rev. Fluid Mech., 49(1):387–417, 2017.
-  G.I. Taylor. An experimental study of standing waves. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 218(1132):44–59, 1953.
-  J.R. Ockendon and H. Ockendon. Resonant surface waves. Journal of Fluid Mechanics, 59(2):397–413, 1973.
-  B. Bäuerlein and K. Avila. Phase lag predicts nonlinear response maxima in liquid-sloshing experiments. J. Fluid Mech., 925, 2021.
-  M. Peeters, G. Kerschen, and J.C. Golinval. Dynamic testing of nonlinear vibrating structures using nonlinear normal modes. J. Sound and Vibration, 330(3):486–509, 2011.
-  N. Deng, B.R. Noack, M. Morzyński, and L.R. Pastur. Low-order model for successive bifurcations of the fluidic pinball. Journal of Fluid Mechanics, 884:A37, 2020.
-  S.W. Shaw and C. Pierre. Normal modes for non-linear vibratory systems. J. Sound and Vibration, 164(1):85–124, 1993.
-  L. Renson, G. Kerschen, and B. Cochelin. Numerical computation of nonlinear normal modes in mechanical engineering. J. Sound and Vibration, 364:177–206, 2016.
-  S.A. Neild, A.R. Champneys, D.J. Wagg, T.L. Hill, and A. Cammarano. The use of normal forms for analysing nonlinear mechanical vibrations. Phil. Trans. Royal Society A, 373(2051):20140404, 2015.
-  G.I. Cirillo, A. Mauroy, L. Renson, G. Kerschen, and R. Sepulchre. A spectral characterization of nonlinear normal modes. J. Sound and Vibration, 377:284–301, 2016.
-  X. Cabré, E. Fontich, and R. de la Llave. The parameterization method for invariant manifolds I: manifolds associated to non-resonant subspaces. Indiana Univ. Math. J., 52(2):283–328, 2003.
-  F. Kogelbauer and G. Haller. Rigorous model reduction for a damped-forced nonlinear beam model: an infinite-dimensional analysis. J. Nonlinear Sci., 28:1109–1150, 2018.
-  A. Haro and R. de la Llave. A parameterization method for the computation of invariant tori and their whiskers in quasi-periodic maps: rigorous results. J. Differential Eqs., 228(2):530–579, 2006.
-  A. Haro, M. Canadell, J.-L. Figueras, A. Luque, and J.M. Mondelo. The Parameterization Method for Invariant Manifolds: from Rigorous Results to Effective Computations. Springer, New York, 2016.
-  R. Szalai. Invariant spectral foliations with applications to model order reduction and synthesis. Nonlinear Dyn., 101:2645–2669, 2020.
-  C.M. Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics. Springer-Verlag New York, 2006.
-  S. Ponsioen, T. Pedergnana, and G. Haller. Analytic prediction of isolated forced response curves from spectral submanifolds. Nonlinear Dyn., 98:2755–2773, 2019.
-  S. Ponsioen, S. Jain, and G. Haller. Model reduction to spectral submanifolds and forced-response calculation in high-dimensional mechanical systems. J. Sound and Vibration, 488:115640, 2020.
-  H. Dankowicz and F. Schilder. Recipes for Continuation. Society for Industrial and Applied Mathematics, 2013.