We present a data-driven and interpretable approach for reducing the dimensionality of chaotic systems using spectral submanifolds (SSMs). Emanating from fixed points or periodic orbits, these SSMs are low-dimensional inertial manifolds containing the chaotic attractor of the underlying high-dimensional system. The reduced dynamics on the SSMs turn out to predict chaotic dynamics accurately over a few Lyapunov times and also reproduce long-term statistical features, such as the largest Lyapunov exponents and probability distributions, of the chaotic attractor. We illustrate this methodology on numerical data sets including delay-embedded Lorenz and Rössler attractors, a nine-dimensional Lorenz model, a periodically forced Duffing oscillator chain, and the Kuramoto–Sivashinsky equation. We also demonstrate the predictive power of our approach by constructing an SSM-reduced model from unforced trajectories of a buckling beam and then predicting its periodically forced chaotic response without using data from the forced beam.

In high-dimensional systems, inertial manifolds containing the low-dimensional chaotic attractor play a significant role in understanding and predicting chaotic behaviors. Such low-dimensional, attracting inertial manifolds contain the essential dynamics of the system and, therefore, offer a tool for model reduction. We construct such inertial manifolds as spectral submanifolds (SSMs) emanating from an unstable steady state and develop SSM-reduced models on them. SSM-reduced models of chaotic systems give accurate short-term predictions and simultaneously reproduce long-term statistical properties. We further demonstrate the predictability of such SSM-reduced models on several numerical data sets including the delay-embedded Lorenz and Rössler attractors, forced Duffing oscillator chain, the Kuramoto–Sivashinsky equation, and a periodically forced buckling beam.

Chaos is ubiquitous in a variety of systems, varying from the double pendulum to turbulent flows.1–3 Most of these chaotic systems are high-dimensional but their chaotic dynamics take place in low-dimensional attractors. The dynamics within these strange attractors provides an ideal reduced model for the full system, as long as appropriate localized coordinates can be constructed along the attractor. This typically requires the construction of a smooth, low-dimensional manifold into which the attractor is properly embedded.

Based on Whitney’s4 and Takens’s5 embedding theorems, several algorithms have been developed for attractor reconstruction and for learning appropriate coordinates near those attractors from delay-embedded data.6–8 In general, a higher embedding dimension preserves more information at the cost of increasing system dimensionality. To infer the optimal embedding dimension and delay embedding time lag, using average mutual information (AMI)9 and the false nearest neighbor (FNN) method10 has become a standard practice. For high-dimensional time series, the dimensionality of the measurements is typically higher than that of the intrinsic system dynamics, requiring dimensional reduction for system identification and modeling.

For this purpose, it is customary to first perform model reduction and then fit a low-dimensional chaotic model in the reduced space. Frequently used linear model reduction techniques, such as the principal component analysis (PCA)11 or the dynamic mode decomposition (DMD) and its variants,12 cannot capture intrinsically non-linearizable phenomena such as chaos. To counter this limitation, several nonlinear model reduction techniques have been developed, including kernel PCA and local linear embedding.13,14 Other popular methods for data-driven model reduction include autoencoders, deep learning, and neural networks. These methods aim to extract hidden information with multiple levels of representation from nonlinear and complex dynamical systems. Various algorithms for black-box models have been developed, compared, and applied in various fields, including fluid dynamics and neural science.15,16

To improve the physical interpretability of models produced by these approaches, previous studies have utilized manifold learning17 or the concept of inertial manifolds.18 Inertial manifolds are smooth, finite-dimensional invariant manifolds that contain the global attractor and attract all solutions exponentially. In principle, therefore, inertial manifolds are the perfect tool for reducing complex systems. To utilize this tool, a recent method has been developed involving training neural networks to learn an atlas of low-dimensional inertial manifolds in high-dimensional spatiotemporal chaotic systems,19,20 which show significant promise.

Subsequent to model reduction, modeling chaos in the reduced coordinates is a further challenge. Such reduced models can be classified into local models and global models.21,22 Local models focus on capturing the dynamics within a small region of the state space for short-term prediction. An example of local chaos modeling is the k-th nearest neighbor (kNN) method,23 which only relies on neighboring trajectory information to make forecasts. Although this method is inherently discrete and contains a large number of parameters, it has been proven to be practically useful when a significant amount of data is available.21 In contrast, global models aim to capture the overall behavior across the whole phase space. Popular global modeling techniques include neural ODEs24 and recurrent neural networks, such as long short-term memory networks (LSTM).25,26 While these neural network models can make forecasts over short Lyapunov times and reproduce statistical properties of the system, they often require a large time domain for training data, extensive parameter tuning, and trial-and-error model selection.

A data-driven modeling alternative that offers simplicity and interpretability is the Sparse Identification of Nonlinear Dynamical Systems (SINDy) approach, which applies sparse regression to fit the right-hand side of a system of ODEs to the time-history of select observables that are assumed to form a closed, reduced dynamical system.27 While this methodology has been successful in several settings,28,29 its results are dependent on the choice of the selected function library, require an a priori knowledge of a reduced observable set, and suffer from sensitivity to noise. For increased interpretability, recent developments combine DMD or SINDy with physical principles.30,31

In this work, we introduce a data-driven model reduction method for chaotic dynamics using spectral submanifolds. A spectral submanifold (SSM) is the smoothest nonlinear continuation of a non-resonant spectral subspace of the linear part of a system at a stationary state.32 This stationary state can be a fixed point, periodic orbit, or quasiperiodic orbit. SSM-based model reduction has so far been used in the data-driven modeling of non-linearizable and deterministic dynamics in mechanical systems.33 SSM-reduction has also been successfully combined with model-predictive control of soft robots34,35 and augmented with a multivariate delay-embedding technique.36 In recent work, the mathematical theory behind SSM-reduction has been extended to fractional and mixed-mode SSMs,37 with the latter type of SSM also containing unstable modes in addition to the slowest stable modes.

This paper aims to further extend data-driven SSM-reduction to chaotic systems. Given an observable space that has a dimension higher than the Lyapunov dimension of the strange attractor, we seek to reconstruct a mixed-mode SSM that contains the chaotic attractor. Such an SSM then constitutes an explicitly constructed, very low-dimensional inertial manifold that is anchored to a stationary state. We illustrate examples that mixed-mode, SSM-reduced models of chaotic systems give accurate short-term predictions while they simultaneously reproduce long-term statistical properties, such as the Lyapunov exponent and probability distribution density.

We consider a dynamical system of the form
x ˙ = F ( x ) = A x + f ( x ) , x R n , A R n × n , f = O ( | x | 2 ) C .
(1)
We assume that A R n × n is diagonalizable and x = 0 is a hyperbolic fixed point, i.e.,
0 Re [ spect ( A ) ] ,
(2)
where spect ( A ) = { λ 1 , λ 2 , , λ n } denotes the spectrum of the linear part of the dynamics at the origin. We further assume the following non-resonance conditions within spect ( A ):
λ j k = 1 n m k λ k , m k N , k = 1 n m k 2 , j = 1 , , n .
(3)
We then select a d-dimensional spectral subspace E R n that is the direct sum of a family of real eigenspaces of A. Without loss of generality, we assume that the spectrum of A within E, spect ( A | E ), contains p purely real eigenvalues and q pairs of complex conjugate eigenvalues, where p + 2 q = d. After a linear change of coordinates, the linear part of the system (1) within the spectral subspace E can be assumed diagonal. We denote the coordinates along the real eigenvectors by u R p and their corresponding Jordan block by
A = diag [ λ 1 , , λ p ] R p × p .
(4)
Similarly, we denote the coordinates along the complex conjugate q pairs of eigenvectors by z C q, and their real Jordan block by
B i = ( α i ω i ω i α i ) R 2 × 2 , i = 1 , , q .
(5)
We finally assume that the rest of the spectrum of A outside spect ( A | E ) comprises r purely real eigenvalues and s pairs of complex conjugate eigenvalues that correspond to the coordinates
v R r , w C s ,
(6)
respectively, within our final set of coordinates ( u , z , v , w ).
Spectral submanifolds (SSMs) are invariant manifolds W ( E ) that serve as nonlinear continuations of a ( p + 2 q )-dimensional spectral subspace
E = { ( u , z , v , w ) : v = 0 , w = 0 }
(7)
of the linearized dynamics at a fixed point. When the linear part of system (1) is asymptotically stable and the non-resonance conditions (3) hold, then there is a unique, primary W ( E ) that is as smooth as the full system.32,38 Recent work by Haller et al.37 extended this result to general hyperbolic fixed points with general (like-mode or mixed-mode) spectral subspaces. (Like-mode spectral subspaces contain only eigenspaces of the same stability type, in contrast to their mixed-mode counterparts.) In this paper, we focus only on the unique, C -smooth primary SSMs, W ( E ), that are tangent to mixed-mode spectral subspaces E of hyperbolic fixed points.
By the existence theory of such primary SSMs,37 for any positive integer K 2, there exists a unique set of coefficients stored in a matrix G such that in the ( u , z , v , w ) coordinates introduced above, the SSM, W ( E ), tangent to a spectral subspace E can locally be written near x = 0 as a C graph over the ( u , z ) variables
( v , w ) T = G ( u , z ) = G ( u , z , z ¯ ) 1 : K + O ( | ( u , z ) | K ) , G = [ G 1 , G 2 , , G K ] C ( n d ) × k = 1 K d k , G k C ( n d ) × d k .
(8)
Here, the superscript ( ) l : r denotes a vector of all monomials of the variables listed in ( ) from order l to order r and d k is the total number of d-variate monomials composed of those variables at order k. For instance, if p = 2 , q = 0 in Eqs. (4) and (5), and the polynomial order we select to approximate the SSM is K = 3, then we have ( u , z , z ¯ ) = [ u 1 , u 2 ] T, d 2 = 3 and d 3 = 4, generating the set of monomials
( u , z , z ¯ ) 2 : 3 = [ u 1 2 , u 1 u 2 , u 2 2 , u 1 3 , u 1 2 u 2 , u 1 u 2 2 , u 2 3 ] T ,
(9)
and the coefficient matrix G 2 : 3 R ( n 2 ) × ( 3 + 4 ) is real.
The reduced dynamics on the SSM, W ( E ), are given by
( u ˙ z ˙ ) = ( A 0 0 0 0 B 1 0 0 0 0 0 0 0 0 B q ) ( u z ) + f u , z ( u , z , G ( u , z ) ) ,
(10)
where f u , z denotes the ( u , z ) coordinate component of f = ( f u , z , f v , f w ).

The open source Matlab software SSMTool39 can compute a Taylor series approximation of the SSM, W ( E ), up to any required order when system (1) is explicitly known. When system (1) is only known from data, then the open source SSMLearn38 algorithm can approximate W ( E ) purely from observed trajectories.

In practice, physical systems in which all state variables can be simultaneously monitored are rare. This necessitates the use of invariant manifold reconstruction methods relying only on a limited number of observables. Specifically, using the Takens delay-embedding theorem,5,40 it is possible to embed an invariant manifold of dimension d in the space of 2 d + 1 (or higher) time-shifted samples of a single, generic scalar observable. Recently, Ref. 36 discussed the use of delay embedding to recover specifically SSMs and their tangent spaces in such observable spaces.

To recall this methodology, we let W be a d-dimensional spectral submanifold of (1) and denote the reduced dynamics on W by η ˙ = R ( η ), η R d, which generates the reduced flow map R t : η 0 η ( t , η 0 ). Then, there exists a smooth SSM parametrization M : R d R n, such that the reduced flow map R t is conjugate to the restriction of the full flow map F t : R n R n of system (1) to W R n,
F t ° M = M ° R t .
(11)
We then consider time-resolved, scalar-valued observations s ( t ) = ζ ( x ( t ) ) of system (1) from a generic observable ζ : R n R. We define a delay coordinate map with m delays: Ψ : R n R m by stacking m subsequent entries of the time series s ( t ), sampled uniformly in time with sampling time T s, in a vector
y ( t ) = Ψ ( x ( t ) ) = [ s ( t ) , s ( t + T s ) , , s ( t + ( m 1 ) T s ) ] T R m .
(12)

The map Ψ is thus a mapping from the phase space of system (1) to the delay-embedded space R m. Let F Ψ : R m R m be the flow map induced by the flow of system (1) on this delay embedding space. By construction, if the observation map ζ maps 0 to 0, the x = 0 equilibrium point of system (1) is mapped into a fixed point y = 0, i.e., F Ψ ( 0 ) = 0. According to the Takens embedding theorem,5,40 for a generic observable function ζ and under certain nondegeneracy conditions on the sampling time T s, if m > 2 d, then the delay coordinate map Ψ restricted to W provides an embedding of W into R m. As a consequence, W ~ = Ψ ( W ) is diffeomorphic to the manifold W.

The delay-embedded dynamics F Ψ on W ~ and the reduced dynamics restricted to the invariant manifold W are conjugate, i.e.,
F Ψ ° Ψ = Ψ ° F t | W ,
(13)
which provides the theoretical basis for identifying an SSM attached to the y = 0 origin of the m-dimensional delay-embedding space.

For dynamical systems with (generally non-smooth) global attractors, inertial manifolds have been defined as finite-dimensional, forward-invariant manifolds that are at least Lipschitz and contain the global attractor.18,41 Therefore, attracting SSMs with an attractor in their reduced dynamics function as inertial manifolds, at least locally near a stationary state. If the attractor within the SSM is globally unique, then an SSM is, in fact, an inertial manifold. Unlike general inertial manifolds, however, SSMs can be precisely located and computed up to arbitrary order of accuracy, which makes them the perfect tools for a mathematically precise reduction of high-dimensional systems to very low-dimensional models.

In order to contain the attractor, the dimension of the SSM should be larger than the Hausdorff dimension of the attractor.42 In a data-driven setting with no prior information on the attractor dimension, the necessary embedding dimension is often estimated via the false nearest neighbors (FNN) method.1,43 This method is based on the idea that if the embedding dimension is too low, neighboring points in the embedding space may appear closer than they actually are, leading to false neighbors. Specifically, the FNN algorithm examines data points in the delay-embedding space and checks if they are indeed close to each other in all directions. The percentage of the false nearest neighbors goes down when the embedding dimension is increased. When the percentage becomes lower than a predefined value, the data disentangle in the embedding space, thus giving a correct approximation for the attractor dimension. The SSMs we construct will have a dimension of d that is at least as large as the estimate obtained from the FNN algorithm.

According to the SSM theory discussed in Sec. II A, the primary SSM, W ( E ), tangent to the spectral subspace E at the origin, can locally be written as a polynomial expansion over E. The open source Matlab and Python implementations of SSMLearn44 compute the SSM parametrization and its reduced dynamics purely from observed data.33,38 In this section, we discuss the computation of the SSM-reduced dynamics and predictions from such reduced models of the full system dynamics.

To provide a clearer explanation, in this section, we only consider q = 0, i.e., the chosen spectral subspace E only contains p real eigenvalues. Our procedure, however, extends directly to spectral subspaces with complex eigenvalues, such as those arising in Ref. 37. The applications of numerical datasets in Sec. IV contain examples with general choices of spectral subspaces, i.e., with both q = 0 and q 0.

Let the coordinate y R ρ denote data obtained from system (1) either via direct observations or delay-embedding (see Sec. II B). Also, let η R d denote coordinates along the corresponding spectral subspace E of D Ψ ( E ), respectively. As a graph over this spectral subspace, the embedded SSM can be approximated by a polynomial expansion. According to the C version of the SSM theory outlined in Sec. II A, for any positive integer order K 2, the coefficients of this polynomial expansion are unique. We, therefore, seek a parametrization M : R d R ρ of the embedded SSM, W ~ ( E ) R ρ, as a K-th order polynomial in the form
M ( η ) = V η 1 : K = V 1 η + V 2 : K η 2 : K , V = [ V 1 , V 2 , , V K ] R ρ × k = 1 K d k , V k R ρ × d k ,
(14)
with the reduced coordinates η = V 1 T y, where the matrix V 1 R ρ × d has orthonormal columns that span the tangent space of the manifold W ~ ( E ) R ρ at y = 0.
To find the constant in the above parametrization from data, we use SSMLearn to solve the constrained minimization problem
V = [ V 1 , V 2 : K ] = arg min V 1 , V 2 : K j y j V 1 V 1 T y j V 2 : K η 2 : K , V 1 T V 1 = I , V 1 T V 2 : K = 0 ,
(15)
where the last constraint represents a basic nonlinear extension of the principal component analysis.
To quantify the accuracy of this SSM approximation, we define the invariance error as follows. From a specific observation y i 0, we compute the corresponding reduced coordinate η i = V 1 T y i as the projection of y i onto the tangent space V 1. By lifting η i to the full embedding space via the manifold parametrization (14) and comparing it with the original data point y i, we obtain a measurement of the manifold fitting accuracy for a given choice of the approximation order K. For a total of P observations, we denote the lifted data points by M ( η ) = M ( V 1 T y i ) and define the normalized invariance error as
Invariance Error = 1 | | y _ | | 1 P i = 1 P | | y i M ( V 1 T y i ) | | , | | y _ | | = max i | | y i | | .
(16)
A larger manifold parametrization order K will generally decrease the invariance error, but excessively large K values will lead to overfitting. In the examples from this paper, K is often chosen at the local minimum of the invariance error. The value of the invariance error depends on specific problems but keeping it below 0.5 % suffices in our experience.

We obtain a model from the reduced dynamics on the SSM by fitting the right-hand side of an ODE to the numerically determined, SSM-reduced vector field in the reduced coordinates. In this paper, we employ two alternative approximations for the reduced dynamics models: a polynomial fit and a nearest neighbor interpolation.

In our first prediction method, we approximate the reduced dynamics in the form
η ˙ = R 1 η + R 2 : K η 2 : K ,
(17)
where the coefficient matrix R = [ R 1 , R 2 : K ] R d × k = 1 K d k is obtained via regression as
R = arg min R j | | η ˙ j R η j 1 : K | | .
(18)
We then diagonalize the linear part of the reduced dynamics (17) by a linear transformation ξ = W 1 η, where ξ R d denotes modal coordinates and W R d × d is the matrix of the eigenvectors of R 1. In those new coordinates, the ODE (17) becomes
ξ ˙ = Λ ξ + N 2 : K ξ 2 : K ,
(19)
where Λ R d × d is a diagonal matrix containing the eigenvalues of R , and N 2 : K R d × k = 2 K d k is the coefficient matrix of nonlinear monomials of ξ from order 2 to K. The coefficient matrix N 2 : K can be specified via a recursive sequence of normal form transformations that preserve the dynamics, as explained in Ref. 38. Here, we do not perform this step because subsequent normal form transformations are generally defined on increasingly smaller neighborhoods of y = 0. This would limit our ability to capture larger-sized attractors on the SSM.

A global polynomial representation of the reduced dynamics may require very high-order polynomials, and yet fail to reproduce the chaotic dynamics.21 Also, the lack of training data outside of the attractor may lead to limited data structures in monomial spaces, which makes it difficult to fit a globally stable polynomial model to vector fields with chaotic trajectories. As an alternative, for complicated strange attractors, we turn to the simplest and earliest method of local forecasting:45,46 we fit maps locally to data points, and then use the resulting collection of local polynomials to form a global model.

This idea brings us to our second prediction method, the local k-th nearest neighbor (kNN) approach. The kNN prediction method23,47 identifies the k nearest neighbors in the training data set and then, by using various interpolation or regression methods, makes local predictions based on those neighbors.48 This method is suitable for predicting highly nonlinear or even chaotic behavior.49 In this paper, for a prediction based at a point η, we first search for its k-th nearest neighbor η j, j = 1 , , k in the training data set that minimizes the Euclidean norm η η j 2. We then use a zeroth-order linear fit to give the prediction η at the next time step, i.e., define the weight w j as
w j = η η j 2 j = 1 k η η j 2
(20)
and take the weighted sum
η = j = 1 k w j η η j 2 .
(21)

Although the kNN method, like other local chaos prediction methods, has the limitation of discontinuities, large numbers of adjustable parameters, and the ability to predict only within the training data range, as long as the training data is sufficiently dense on the attractor, the method still has good predictive power.

In this paper, we only use polynomial regression and the kNN method to fit the reduced dynamics. After obtaining the data sets in reduced coordinates, however, it is also possible to combine SSM-based model reduction with other forecasting techniques. We discuss this further in the Conclusions.

For non-chaotic dynamics, the prediction errors from SSM-reduced models are quantified by the normalized mean-trajectory-error (NMTE),38 which is the normalized difference between the test trajectory and the predicted trajectory. Chaotic dynamics, however, have sensitive dependence on initial conditions, which implies that the predictions for individual trajectories will invariably fail for longer times. Indeed, even a well-constructed reduced model can only be expected to capture statistical measures of the chaotic dynamics in long-term predictions. Accordingly, in this paper, we will evaluate the ability of SSM-reduced models to reproduce both the Lyapunov exponent and the probability density distribution of the attractors.

First, we consider the Lyapunov characteristic exponents,50,51 the most often used indicators of chaotic flows. Let x ( t ; x 0 ) denote a trajectory of system (1) starting from the initial condition x 0 at t = 0. An initial perturbation ε ( 0 ) = ε 0 to x 0 will evolve along x ( t ; x 0 ) as ε t = Φ ( t ; x 0 ) ε 0, where Φ ( t ; x 0 ) is the fundamental matrix solution of the linearized system (1) along the trajectory x ( t ; x 0 ). The maximal exponential rate of expansion along x ( t ; x 0 ) is then the maximal Lyapunov exponent (MLE)
λ max ( x 0 ) = lim t 1 2 t log λ max [ Φ T ( t ; x 0 ) Φ ( t ; x 0 ) ] ,
(22)
if the limit exists.

As the maximum rate of separation of nearby trajectories, λ max ( x 0 ) defines a local measure of predictability for a dynamical system. The Lyapunov exponents are independent of the initial condition x 0 and depend only on the trajectory x ( t ; x 0 ). If x ( t ; x 0 ) is contained in a chaotic (or strange) attractor, then λ max ( x 0 ) is constant for almost all trajectories along the attractor.52 

One way to approximate λ max ( x 0 ) is to compute the exponent of the rate at which trajectories from nearby initial conditions diverge from the trajectory x ( t ; x 0 ), i.e., calculate
λ ( t ; x 0 , ϵ 0 ) = log | x ( t ; x 0 + ϵ 0 ) x ( t ; x 0 ) | | ϵ 0 |
(23)
for a randomly chosen initial perturbation ϵ 0 to x 0 and for large enough times t. For a chaotic system, λ ( t ; x 0 , ϵ 0 ) is approximately a linear function of time.53 When a good linear fit to λ ( t ; x 0 , ϵ 0 ) is not possible, the system is generally not chaotic.54 While λ ( t ; x 0 , ϵ 0 ) / t depends on the initial condition x 0 and the perturbation ϵ 0, this ratio will converge to λ max ( x 0 ) for generic choices of x 0 and ϵ 0.52 Based on this observation, we randomly pick two nearby initial conditions on chaotic attractors and use a linear fit to their temporal separation rate to estimate the Lyapunov exponent λ max of the attractor.
A related measure of predictability is the Lyapunov time, defined as the reciprocal of λ max of the attractor,55 i.e.,
Lyapunov Time = 1 λ max .
(24)
This quantity gives the characteristic time scale of chaotic dynamics. We call a prediction a short-time prediction if it is carried out over times that are low-order multiples of the Lyapunov time.
Along with the Lyapunov exponent, the probability density56,57 is also an invariant property of attractors and hence is often used for comparisons between the original and modeled chaotic attractors.58,59 If a chaotic attractor contains the orbit x ( t ), its probability measure σ is defined as the asymptotic average of Dirac deltas along x ( t ),60 i.e.,
σ = lim T 1 T 0 T δ x ( t ) d t ,
(25)
where the integral of δ x ( t ) over the entire trajectory is equal to one.

On a chaotic attractor, the invariant measure σ is independent of the initial conditions. To approximate the probability density of the chaotic attractors, we will take a sampling time that is long enough for σ to converge. Along each axis of the system, we can either count the number of data points between intervals and plot the histogram or estimate the probability distribution σ based on a normal kernel function.22 In our examples, we will compare the probability density of the original data points with those of the predictions made by SSM-reduced model using the same kernel function. The closeness of the two plots indicates how statistically accurate our long-term SSM-based predictions are.

In summary, we will test SSM-based prediction quality for short times via trajectory error calculations and for long times via Lyapunov exponent and probability density distribution calculations.

In this section, we perform SSM-based model reduction on four numerical data sets containing chaotic attractors.

In our first two examples, we observe only a scalar time series and delay-embed an SSM containing the attractor. In our subsequent examples, we assume that all phase space variables can be fully observed. For the first delay-embedded Lorenz attractor, we use a global polynomial representation for the reduced dynamics. For the other examples, we adopt the kNN nearest neighbor approach.

For autonomous ODE examples, the SSM-reduction outlined in Sec. II A will be directly applicable. We also explore one non-autonomous (time-periodic) ODE example in Sec. IV C, where we will carry out the discrete version of the same procedure for the Poincaré map of the system. Specifically, we perform model reduction on the Poincaré map, identifying an SSM emanating from its unstable fixed point.

We further include the example of the delay-embedded Rössler attractor in Sec. IV D, in which the SSM is attached to the unstable subspace corresponding to a complex conjugate pair of eigenvalues. The applications of SSM-based model reduction range from classic chaotic attractors to low-dimensional attractors in high-dimensional systems, such as the Kuramoto–Sivashinski equation in Sec. IV E and the finite element beam model in Sec. IV F.

As our first example, we consider the Lorenz model61 
x ˙ = σ ( y x ) , y ˙ = x ( ρ z ) y , z ˙ = x y β z ,
(26)

with the classic parameter values σ = 10, ρ = 28, and β = 8 / 3. This system has an unstable fixed point at ( x , y , z ) = 0 and hence falls in the setting of this paper. We will only use a single scalar signal in our analysis and hence our methodology is fully unaware of the dimension and form of the ODE (26). As an SSM containing a chaotic attractor must be at least three-dimensional, Takens’s theorem guarantees success (with probability one) in embedding the SSM in an at least seven-dimensional space. We then seek to reconstruct the Lorenz attractor within the 3D SSM.

To compare the SSM-reduced model results with the broadly used SINDy algorithm mentioned in the Introduction, we will adopt the same setting as in the Matlab code given by Brunton et al.27 Accordingly, we will observe only the x coordinate of system (26) to generate a scalar observable time series. We generate two trajectories (one for training and the other for testing) over the interval ( 0 , 100 ) with a time step of 0.001. The test data starts with the classic initial condition ( 8 , 8 , 27 ) T, while the initial condition of the training data is slightly separated ( 8 , 8 , 27 + ϵ ) T, where ϵ = 10 9. After truncating the data over the first short period of the time interval ( 0 , 1 ), each of these two individual trajectories lies on the low-dimensional attractor of system (26).

We use FNN to estimate the dimension d of the SSM that prevails in the data. The different choices of model reduction parameters and corresponding manifold fitting results are shown in Table I. The high FNN percentage, 31.4 %, in the first row of Table I indicates that the SSM containing the strange attractor has a dimension higher than two. This prompts us to select the SSM dimension as d = 3, Based on the Takens embedding theorem, a 3D SSM can generically be delay-embedded in a ρ = 2 × d + 1 = 7 dimensional space. Table I shows the accuracy of SSM reconstruction in this setting under different polynomial approximation orders for the 3D SSM.

TABLE I.

Different choices of SSM dimension, false nearest neighbor percentage (FNN), order of manifold parametrization, and the corresponding invariance error. Due to the high FNN percentage for d = 2, we did not compute an SSM approximation for that case.

SSM dimension (d
FNN percentage  31.4%  0.0%  0.0%  0.0% 
Delay embedding dimension (ρ
Manifold expansion order ( K … 
Invariance Error  …  3.85%  1.25%  1.45% 
SSM dimension (d
FNN percentage  31.4%  0.0%  0.0%  0.0% 
Delay embedding dimension (ρ
Manifold expansion order ( K … 
Invariance Error  …  3.85%  1.25%  1.45% 

In delayed coordinates, the origin is a fixed point corresponding to the origin of system (26). We therefore seek a 3D SSM attached to the origin in the seven-dimensional embedding space and determine the reduced dynamics on that SSM. By computing the normalized invariance error defined in Eq. (16) and identifying its minimum, we choose the order K = 3 for the polynomial expansion for the SSM. We find that the 3D reduced dynamics on the SSM does admit a chaotic attractor, as shown in Fig. 1.

FIG. 1.

Reconstruction and prediction results of Lorenz attractor in 3D. The inertial manifold is reconstructed by observing only x coordinate as a scalar time serie, delay embedded into 7D, and performing model reduction to 3D. Dynamics on SSM is represented by a fifth-order polynomial.

FIG. 1.

Reconstruction and prediction results of Lorenz attractor in 3D. The inertial manifold is reconstructed by observing only x coordinate as a scalar time serie, delay embedded into 7D, and performing model reduction to 3D. Dynamics on SSM is represented by a fifth-order polynomial.

Close modal
In the modal coordinates (19) that diagonalize the linear part of the SSM-reduced dynamics, we obtain the polynomial reduced-order model on the 3D SSM in the form
u ˙ = 11.0194 u 0.0001 u 2 0.0005 u v 0.0025 u w 0.0005 v 2 + 0.0026 w 2 + 0.0317 u 3 + 0.0622 u 2 v 0.0011 u 2 w + 0.0421 u v 2 + 0.0102 u v w 0.0143 u w 2 + 0.0058 v 3 + 0.0235 v 2 w 0.0349 v w 2 + 0.0108 w 3 , v ˙ = 2.3859 v + 0.0036 u w + 0.0002 v 2 + 0.0006 v w 0.0025 w 2 0.0273 u 3 0.0550 u 2 v 0.0001 u 2 w 0.0391 u v 2 0.0096 u v w + 0.0117 u w 2 0.0059 v 3 0.0231 v 2 w + 0.0357 v w 2 0.0131 w 3 , w ˙ = + 9.3871 w 0.0001 u 2 0.0004 u v + 0.0032 u w + 0.0009 v w 0.0016 w 2 0.0144 u 3 0.0303 u 2 v 0.0001 u 2 w 0.0233 u v 2 0.0054 u v w + 0.0058 u w 2 0.0041 v 3 0.0147 v 2 w + 0.0246 v w 2 0.0112 w 3 .
(27)

One might wonder what we have achieved here, given that system (27) is substantially more complicated than the original system (26) and has the same dimension as that system. One should not forget, however, that the input to our procedure was a single scalar observable and we used no information about the dimension or the form of the Lorenz system (26). The procedure we follow here is the same for any system that is only known from a scalar observable. This procedure will generally yield a smooth manifold, the SSM, that contains the attractor of the system.

With the SSM-reduced model, we can now make predictions for test trajectories not used in constructing this model, and compare it with predictions obtained from the SINDy algorithm. From the same training data, SINDy27 constructs a Hankel matrix by stacking 10 delayed time series as rows, and uses singular value decomposition (SVD) to select the first three dominant modes in these time-delay coordinates. Each column of the data is then trained with a different sparsification parameter ( λ SINDy = 0.01 , 0.2 , 2) to fit a third-order polynomial right-hand side of an ODE to the data.

Figure 2 shows the comparison of both short- and long-term predictions from SINDy and SSM reduction. Both methods are able to forecast for short times but our approach shows a more accurate result [Figs. 2(a) and 2(b)]. We then perform a Lyapunov analysis by plotting the separation of two trajectories against time in a log plot in Fig. 2(c). The predictions given by SSMLearn admit a clear linear fit which indicates our model is chaotic, while the trajectories reconstructed by SINDy are not chaotic. Indeed, as noted in Sec. 4.5 in the Appendix of Ref. 27, although the skeleton of the Lorenz attractor is captured by SINDy, the attractor is reconstructed as a non-chaotic, quasi-periodic orbit. In contrast, the Lyapunov exponent from the SSM-based model (27) is 0.89, which is very close to the real value of 0.90 computed directly from the Lorenz system (26).

FIG. 2.

Time series forecasting results for delay-embedded Lorenz attractor, and its reconstructed Lyapunov exponent. Subplots (a) and (b) show more accurate predictions from SSM-reduction than from SINDy over approximately 2.97 Lyapunov times (NMTE less than 0.1). Subplot (c) plots the separation of nearby trajectories against time. The red prediction made by the SINDy model has the maximal Lyapunov exponent 0.23 compared with the real MLE 0.90 of system (26). A blue linear fit to the SSM prediction confirms that trajectories separate exponentially, and the system is chaotic with the maximal Lyapunov exponent 0.89. At the same time, the red prediction made by SINDy suggests that the system is not chaotic.

FIG. 2.

Time series forecasting results for delay-embedded Lorenz attractor, and its reconstructed Lyapunov exponent. Subplots (a) and (b) show more accurate predictions from SSM-reduction than from SINDy over approximately 2.97 Lyapunov times (NMTE less than 0.1). Subplot (c) plots the separation of nearby trajectories against time. The red prediction made by the SINDy model has the maximal Lyapunov exponent 0.23 compared with the real MLE 0.90 of system (26). A blue linear fit to the SSM prediction confirms that trajectories separate exponentially, and the system is chaotic with the maximal Lyapunov exponent 0.89. At the same time, the red prediction made by SINDy suggests that the system is not chaotic.

Close modal
We consider a 3D viscous fluid layer that is uniformly heated from below, and let t and r be the time and space variables. The density field ρ ( r , t ) describes the mass distribution of the fluid, and the physical state is determined by the pressure p ( r , t ), the temperature T ( r , t ), and the velocity field v ( r , t ). The dynamic behavior is governed by three nonlinear PDEs, namely, the equation of continuity, the Navier–Stokes equation, and the equation of thermal conductivity
v t + ( v ) v = ρ ρ 0 g 1 ρ 0 p + ν 2 v , div v = 0 , T t + ( v ) T = χ 2 T ,
(28)
where g = ( 0 , 0 , g ) denotes the vector of gravity, ρ 0 is the standard density, ν is the kinematic viscosity, and χ is the thermal conductivity.
Lorenz applied a double Fourier expansion in deriving his 3D governing equation (26) in Ref. 61. Reference 62 applied a similar approach but extended the double Fourier expansion to triple expansion. After truncation and selection of the leading-order-terms in the expansion, Ref. 62 obtains the nine-dimensional ODE
C ˙ 1 = σ b 1 C 1 C 2 C 4 + b 4 C 4 2 + b 3 C 3 C 5 σ b 2 C 7 , C ˙ 2 = σ C 2 + C 1 C 4 C 2 C 5 + C 4 C 5 σ 2 C 9 , C ˙ 3 = σ b 1 C 3 + C 2 C 4 b 4 C 2 2 b 3 C 1 C 5 + σ b 2 C 8 , C ˙ 4 = σ C 4 + C 2 C 3 C 2 C 5 + C 4 C 5 + σ 2 C 9 , C ˙ 5 = σ b 5 C 5 + 1 2 C 2 2 1 2 C 4 2 , C ˙ 6 = b 6 C 6 + C 2 C 9 C 4 C 9 , C ˙ 7 = b 1 C 7 r C 1 + 2 C 5 C 8 C 4 C 9 , C ˙ 8 = b 1 C 8 + r C 3 2 C 5 C 7 + C 2 C 9 , C ˙ 9 = C 9 r C 2 + r C 4 2 C 2 C 6 + 2 C 4 C 6 + C 4 C 7 C 2 C 8 ,
(29)
where σ is the Prandtl number, r is the Rayleigh number, and the constant parameters b i represent a measure of the geometry of the square cell.
Depending on r and σ, the dynamics can be stationary, periodic, chaotic, or hyperchaotic, with the latter referring to the case wherein multiple positive Lyapunov exponents exist. For certain parameters, four symmetric low-dimensional chaotic attractors are observed in system (29). From Ref. 62, we choose the system parameters as
σ = 1 2 , r = 14.2 , a = 1 2 , b 1 = 4 ( 1 + a 2 ) 1 + 2 a 2 , b 2 = 1 + 2 a 2 2 ( 1 + a 2 ) , b 3 = 2 ( 1 a 2 ) 1 + a 2 , b 4 = a 2 1 + a 2 , b 5 = 8 a 2 1 + 2 a 2 , b 6 = 4 1 + 2 a 2 .
(30)

The origin is a fixed point of system (29) with the eigenvalues listed in Table II. As in our previous example, no information about these eigenvalues will be used in our data-driven SSM reduction; we only list these eigenvalues for reference. Note that the 1:1 resonances between λ 2 and λ 3 and between λ 8 and λ 9 do not violate the non-resonance condition (3). There is also a 5th-order resonance between λ 4, λ 5, and λ 6, which violate Eq. (3). However, as earlier work on SSMs shows,32,38 a polynomial expansion for an SSM, W ( E ), can still be computed as long as the linearized spectrum within E does not contain λ 4, λ 5, and λ 6. We will choose E in this fashion.

TABLE II.

Eigenvalues of the 9D Lorenz attarctor.

λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9
1.9263  −0.2741  −0.2741  −0.5000  −0.6667  −2.6667  −3.4263  −4.7259  −4.7259 
λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9
1.9263  −0.2741  −0.2741  −0.5000  −0.6667  −2.6667  −3.4263  −4.7259  −4.7259 

As in the 3D Lorenz example, we observe only the C 1 coordinate of system (29) as a scalar time series and delay-embed the signal. We generate two random initial conditions (one for training and the other for testing) within the unit sphere of the origin in the 9D phase space and observe the first coordinate C 1. We then select and truncate the trajectories to make them lie close to one of the four symmetric strange attractors. Assuming no knowledge of the original system or the dimension of the chaotic attractor, we again use FNN to estimate the dimension d of the SSM and pick d = 3 according to the false neighbor percentage in Table III. By the Takens embedding theorem, we delay-embed this 3D SSM into a ρ = 2 × d + 1 = 7 dimensional space. The training and testing trajectories in this seven-dimensional delay-embedded space are shown in Fig. 3.

FIG. 3.

The chaotic attractor of the 9D Lorenz model (29) in the 7D delay-embedded space, projected onto the coordinates y 1, y 3, and y 7.

FIG. 3.

The chaotic attractor of the 9D Lorenz model (29) in the 7D delay-embedded space, projected onto the coordinates y 1, y 3, and y 7.

Close modal
TABLE III.

Different choices of SSM dimension, false nearest neighbor percentage (FNN), order of manifold parametrization, and the corresponding invariance error of the delay-embedded 9D Lorenz model. Due to the high FNN percentage for d = 2, we did not compute an SSM approximation for that case.

SSM dimension (d
FNN percentage  41.6%  0.0%  0.0%  0.0%  0.0%  0.0%  0.0% 
Delay embedding dimension (ρ
Manifold expansion order ( K … 
Invariance Error  …  2.24%  0.59%  0.22%  0.09%  0.05%  0.04% 
SSM dimension (d
FNN percentage  41.6%  0.0%  0.0%  0.0%  0.0%  0.0%  0.0% 
Delay embedding dimension (ρ
Manifold expansion order ( K … 
Invariance Error  …  2.24%  0.59%  0.22%  0.09%  0.05%  0.04% 

Using SSMLearn, we approximate a 3D SSM attached to the origin of the 7D delay-embedded space. We minimize the invariance error over different choices of SSM dimension and polynomial approximation orders. As shown in Table III, these considerations lead us to construct a 3D SSM via a sixth-order polynomial expansion. This construct produces an invariance error of order 0.05 %, which is low enough to guarantee high accuracy for our SSM-reduced model.

In the reduced coordinates, we find that a polynomial representation of the SSM-reduced dynamics fails to capture the complex dynamics with sufficient accuracy. Instead, we turn to the k-th nearest neighbor method discussed in Sec. III B. The number of neighbors in kNN algorithm is often chosen to be larger than the data dimension.23 Here, we select k = 4 and perform the kNN method on a training set containing 1.5 × 10 6 data points. Note here that the training data for SSM fitting ( 10 5 points) is in general much less than the number of points required for reduced dynamics fitting. This is because the kNN prediction algorithm requires data to fully cover the attractor in its phase space. Figure 4(a) shows the prediction results from this nearest neighbor method. Since chaotic dynamics have sensitive dependence on initial conditions, the forecast is only accurate for 5.94 Lyapunov times.

FIG. 4.

Trajectory prediction in the time domain [plot (a)] and in the delay-embedding space [plot (b)] from the reduced models for an initial condition not contained in the training data. The length of the SSM prediction (blue) interval is approximately 5.94 Lyapunov times, while the SINDy model (red) fails for trajectory integration at around 1.82 Lyapunov times. We obtain the prediction of the SSM-reduced model by making predictions on 3D reduced coordinates using the kNN method, then projecting the trajectory back in 7D delay-embedded space. The training set contains 1.5 × 10 6 points.

FIG. 4.

Trajectory prediction in the time domain [plot (a)] and in the delay-embedding space [plot (b)] from the reduced models for an initial condition not contained in the training data. The length of the SSM prediction (blue) interval is approximately 5.94 Lyapunov times, while the SINDy model (red) fails for trajectory integration at around 1.82 Lyapunov times. We obtain the prediction of the SSM-reduced model by making predictions on 3D reduced coordinates using the kNN method, then projecting the trajectory back in 7D delay-embedded space. The training set contains 1.5 × 10 6 points.

Close modal

We attempted to perform the SINDy algorithm on the same training data to identify a seven-dimensional SINDy model and compare it to SSM-reduced model prediction results. We tested several combinations of different model orders (from two to six) and sparsity parameters ( λ SINDy = 0.0001 to 0.1). We also tried to follow the SINDy tutorial for63 to tune the hyperparameter λ SINDy by minimizing RMSE of predicted trajectories for different choices of λ SINDy. However, the resulting SINDy models do not contain a global attractor in general, and their solutions tend to blow up after some time of integration. As an example, Fig. 4 shows results from a third-order SINDy model with λ SINDy = 0.001. One could also try to set different sparsity parameters for each coordinate, as for the delay-embedded Lorenz example in Ref. 27, but we did not pursue that exploration here.

To test the statistical properties of the SSM-reduced model, we further compute the NMTE and Lyapunov exponent of 200 pieces of the test trajectory whose initial conditions are not contained in the training data. After adding a small perturbation, we plot the deviation of the original and the predicted model trajectory against time on a log scale. For 200 trajectories of the original system (29), this analysis yields the Lyapunov exponent 0.032 as the mean of the observed separation exponents. Remarkably, this exponent is closely approximated by the value 0.033 yielded by our SSM-reduced model in Fig. 5(b).

FIG. 5.

Trajectory prediction error (a) and the Lyapunov exponent analysis (b) on test data of the 3D SSM-reduced model. We compute the time average of the NMTE and the nearby trajectory separation rate over 200 test trajectories but only plot 20 of them in this figure. The average trajectory separation rate of system (29) is 0.032, which is closely approximated by the mean Lyapunov exponent 0.033 computed from subplot (b) of the SSM-reduced model.

FIG. 5.

Trajectory prediction error (a) and the Lyapunov exponent analysis (b) on test data of the 3D SSM-reduced model. We compute the time average of the NMTE and the nearby trajectory separation rate over 200 test trajectories but only plot 20 of them in this figure. The average trajectory separation rate of system (29) is 0.032, which is closely approximated by the mean Lyapunov exponent 0.033 computed from subplot (b) of the SSM-reduced model.

Close modal
We now consider the classic forced-damped Duffing oscillator64 given by
x ¨ + δ x ˙ + α x + β x 3 = γ cos ω t ,
(31)
where δ, α, β, and γ are the coefficients for damping, linear stiffness, nonlinear stiffness, and forcing amplitude, respectively. We use this nonlinear oscillator as the third element in an otherwise linear chain of oscillators, as shown in Fig. 6. We denote the displacement of the i-th degree of freedom by the coordinate x i. We then apply periodic forcing along x 3 as illustrated in Fig. 6.
FIG. 6.

Forced harmonic oscillator chain with a single Duffing oscillator added. We apply periodic forcing to this Duffing oscillator.

FIG. 6.

Forced harmonic oscillator chain with a single Duffing oscillator added. We apply periodic forcing to this Duffing oscillator.

Close modal
We choose the third mass as m 3 = 1 kg and the remaining mass as m 1 = m 2 = m 4 = 0.1 kg each; the third linear damping coefficient attached to the ground is c 3 = 0.1 N s/m; the other damping coefficients are c 1 = c 2 = c 3 = c 4 = 0.75 N s/m. The linear spring stiffness is k = 1 N/m and the third nonlinear spring has the force relation F 3 ( x 3 ) = k 3 x 3 + β x 3 3, where k 3 = 3 N/m and β = 0.25 N/ m 3. The periodic forcing that acts on the Duffing oscillator has an amplitude of 2.25 N and a frequency of 2 rad/s. The resulting equations of motion is given by a system of second-order ODEs of the form
0.1 x ¨ 1 + 0.75 x ˙ 1 + 2 x 1 x 2 = 0 , 0.1 x ¨ 2 + 0.75 x ˙ 2 x 1 + 2 x 2 x 3 = 0 , x ¨ 3 + 0.1 x ˙ 3 x 2 x 3 x 4 + 0.25 x 3 3 = 2.25 cos 2 t , 0.1 x ¨ 4 + 0.75 x ˙ 4 x 3 + 3 x 4 = 0.
(32)

When written in first-order form, these equations define an eight-dimensional, non-autonomous dynamical system. For the autonomous part of the dynamics, the origin is an unstable fixed point. This unstable equilibrium has eight eigenvalues, as listed in Table IV. Again, these are listed only for reference and will not be used in our SSM construction.

TABLE IV.

Eigenvalues of the autonomous part of the Duffing oscillator chain.

λ1 λ2 λ3,4 λ5,6 λ7,8
1.2204  −5.8047  −1.5713 ± 0.6269i  −3.6865 ± 4.0004i  −3.7500 ± 3.9922i 
λ1 λ2 λ3,4 λ5,6 λ7,8
1.2204  −5.8047  −1.5713 ± 0.6269i  −3.6865 ± 4.0004i  −3.7500 ± 3.9922i 

When we add the periodic forcing term to the autonomous part of the equation, the unstable fixed point at the origin bifurcates into an unstable periodic orbit. As long as the amplitude of the time-dependent term in system (32) is moderate, the corresponding T = π-periodic Poincaré map will also have an unstable fixed point. This hyperbolic fixed point is the anchor point for our SSM computation. Our strategy will be to locate an SSM for this map that contains an observed chaotic attractor of the system. As in our previous examples, this SSM will act as an inertial manifold whose restricted dynamics provide a reduced model for the dynamics on and near the attractor. The corresponding geometry is illustrated in Fig. 7.

FIG. 7.

SSM construction for the Poincaré map of system (32) in the phase space.

FIG. 7.

SSM construction for the Poincaré map of system (32) in the phase space.

Close modal

For the eight-dimensional chaotic system (32), we consider all eight degrees of freedom as observables. We select seven initial conditions (six for training and one for testing) randomly within the unit sphere around the origin and generate the corresponding trajectories for 10 5 times the forcing period. We track the evolution of the origin under increasing forcing amplitude into a hyperbolic periodic orbit using the numerical continuation package COCO.65 In Fig. 8(a), we plot the coordinates x 3 and x ˙ 3 for the original continuous dynamics in blue, and the unstable periodic orbit computed by COCO in red. When we pass to the Poincaré map, the continuous orbits become discrete dots in the Poincaré section as in Fig. 8(b). Based on the theory of mixed-mode SSMs,37  SSMLearn computes a polynomial expansion for an SSM attached to the hyperbolic fixed point [red dot in Fig. 8(b)]. Note that our construct relies on the existence of the hyperbolic fixed point, which we compute using COCO from system (32). However, even when the exact location of this anchor point is not known, SSMLearn can still provide an approximation for the SSM automatically by regressing a manifold (with the addition of constant terms in its parameterization) to the available data. The location of the fixed point itself will then follow as a prediction from the SSM-reduced dynamics.

FIG. 8.

(a) Chaotic trajectories (blue) and a saddle-type periodic orbit of system (32) projected onto the ( x 3 , x ˙ 3 ) plane. (b) The same trajectories and periodic orbits as they appear as a fixed point and a trajectory on the attractor of the Poincaré map, projected again to the ( x 3 , x ˙ 3 ) plane.

FIG. 8.

(a) Chaotic trajectories (blue) and a saddle-type periodic orbit of system (32) projected onto the ( x 3 , x ˙ 3 ) plane. (b) The same trajectories and periodic orbits as they appear as a fixed point and a trajectory on the attractor of the Poincaré map, projected again to the ( x 3 , x ˙ 3 ) plane.

Close modal

As in the previous examples, we use the invariance error as a criterion for determining the correct SSM dimension. Table V shows different invariance errors obtained for different choices of SSM dimension and expansion order. Based on these results, we choose a 3D SSM construction via a seventh-order polynomial expansion.

TABLE V.

Different choices of SSM dimension, manifold expansion order and the corresponding invariance error for the oscillator chain (32).

SSM dimension (d) 2D SSM 3D SSM
SSM expansion orders ( K ) 
Invariance error  8.26%  4.24%  1.89%  0.061%  0.005%  0.0001% 
SSM dimension (d) 2D SSM 3D SSM
SSM expansion orders ( K ) 
Invariance error  8.26%  4.24%  1.89%  0.061%  0.005%  0.0001% 

After obtaining reduced coordinates on the 3D SSM, we again use the kNN algorithm discussed in Sec. III B to reconstruct the chaotic attractor of the system on this SSM. We set k = 6 and use a training set with the total size of 6 × 10 5 points. The prediction is generated by the 3D reduced model, then lifted to the eight-dimensional state space via the SSM parametrization. In Fig. 9, we visualize the test data (red) and its prediction (blue) by showing them in various combinations of three of the eight coordinates. From Fig. 9, we see that the red dots and the blue dots coincide in the 8D state space, and the shape of the chaotic attractor is accurately reconstructed. To verify the accurate reproduction of long-term statistical features of the SSM-reduced model, we further plot the probability distribution of the original and the reconstructed dynamics in all of the coordinates separately. From Fig. 10, the statistical properties of the reduced-order model are very close to those of the original attractor.

FIG. 9.

Predictions for the strange attractor of the forced Duffing oscillator chain in the eighth-dimensional full state space, using the kNN prediction method on a three-dimensional SSM. Here k = 6 and the training set contains 6 × 10 5 points.

FIG. 9.

Predictions for the strange attractor of the forced Duffing oscillator chain in the eighth-dimensional full state space, using the kNN prediction method on a three-dimensional SSM. Here k = 6 and the training set contains 6 × 10 5 points.

Close modal
FIG. 10.

Probability distribution for the chaotic attractor of the full system (32) (red) and for the attractor of the 3D, SSM-reduced model in various sets of coordinates (blue).

FIG. 10.

Probability distribution for the chaotic attractor of the full system (32) (red) and for the attractor of the 3D, SSM-reduced model in various sets of coordinates (blue).

Close modal

We also attempted to apply the map version of SINDy66 to the same eight-dimensional training data (six trajectories) for maps. However, after testing different combinations of various values of polynomial order (from three to six) and sparsity parameter ( λ SINDy = 0.001 to 0.05), we failed to identify a feasible model using SINDy and the prediction results invariably blew up after some time. This is likely due to the high complexity of the Poincaré map dynamics.

In our previous examples, we considered mixed-mode SSMs tangent to spectral subspaces with only real eigenvalues [i.e., q = 0 in Eq. (5)]. Here, we turn to an example of an SSM tangent to an eigenspace with complex conjugate eigenvalues, the Rössler attractor system67 
x ˙ = y z , y ˙ = x + a y , z ˙ = b + z ( x c ) .
(33)
This system has two unstable fixed points at
( x 0 y 0 z 0 ) = ( 1 2 ( c ± c 2 4 a b ) 1 2 a ( c c 2 4 a b ) 1 2 a ( c ± c 2 4 a b ) ) .
One of these fixed points lies in the center of the Rössler attractor loop and the other one lies relatively far from the attractor. We choose the classical parameter values a = b = 0.2 , c = 5.7 and focus on the unstable equilibrium ( x 0 , y 0 , z 0 ) = ( 0.007 , 0.035 , 0.035 ) located in the center. We build the mixed-mode SSM attached to this fixed point that has one real eigenvalue 5.687 and a complex conjugate pair 0.097 ± 0.995 i.

As in the delay-embedded Lorenz attractor example in Sec. IV A, we observe only the x coordinate of system (33) to generate scalar observable time series. Then, using the information of the fixed point x 0, we shift the origin of the system to this unstable fixed point. The training and test data are generated separately by the initial conditions ( 1 , 1 , 1 ) and ( 1 , 1 , 2 ) over the time interval ( 0 , 7 000 ) with a time step of 0.01. After the first ( 0 , 20 ) time interval, both the training and test trajectories lie very close to the Rössler attractor. Following the same procedure as in Sec. IV A, we use FNN to estimate the SSM dimension d and pick the manifold expansion order according to the invariance error. We select the SSM attached to the unstable fixed point to be d = 3 dimensional, the delay-embedding dimension to be ρ = 7, and the polynomial expansion order to be K = 4. Namely, we delay-embed the scalar observable in 7D and use SSMLearn to perform model reduction onto a 3D reduced space. The shape of the chaotic attractor in three-dimensional reduced coordinates is shown in Fig. 11(a).

FIG. 11.

The delay-embedded Rössler attractor in reduced coordinates and prediction from the SSM-reduced model. Subplot (a) shows the training and test data in three-dimensional reduced coordinates. Subplots (b) and (c) show the time series forecasting results and the reconstructed Lyapunov exponent. The prediction remains accurate for about 3.3 Lyapunov times. The blue linear fit of the reduced-order model in the subplot (c) gives the maximal Lyapunov exponent 0.0946, which is near 0.0938, the actual MLE of system (33).

FIG. 11.

The delay-embedded Rössler attractor in reduced coordinates and prediction from the SSM-reduced model. Subplot (a) shows the training and test data in three-dimensional reduced coordinates. Subplots (b) and (c) show the time series forecasting results and the reconstructed Lyapunov exponent. The prediction remains accurate for about 3.3 Lyapunov times. The blue linear fit of the reduced-order model in the subplot (c) gives the maximal Lyapunov exponent 0.0946, which is near 0.0938, the actual MLE of system (33).

Close modal

In the reduced coordinates, we use the kNN method to make forecasts for the test trajectory. We select k = 4 and perform the kNN prediction method on a training set containing 7 × 10 5 data points. From Fig. 11, the prediction stays accurate for around 3.3 Lyapunov times of system (33), and the reconstructed maximum Lyapunov exponent 0.0946 is close to the actual MLE of the system 0.0938.

As in the other examples, we also tried to apply SINDy to the same seven-dimensional training data. We tested different combinations of SINDy parameters (polynomial order ranging from three to five, and the sparsity parameter λ SINDy ranging from 0.001 to 0.1) but failed to identify a feasible model. The prediction results generated by these SINDy models generally blew up, i.e., the models do not contain a stable attractor. This may be because the dynamics in these delay-embedded coordinates do not have a sparse polynomial form in general.

To test SSM-based model reduction in higher-dimensional chaotic systems, we consider the Kuramoto–Sivashinsky (KS) equation.68,69 This equation is often used as a simplified model to study spatiotemporal chaos and pattern formation in fluid dynamics, particularly in the context of understanding turbulence. Its dynamics in one spatial dimension is given by the PDE
u t + 1 2 ( u 2 ) x + u x x + u x x x x = 0 , x [ L 2 , L 2 ] ,
(34)
where x is the spatial coordinate and t 0 is the time. Equation (34) describes the time evolution of the velocity u = u ( x , t ) on a periodic domain u ( x , t ) = u ( x + L , t ) with length L. When the velocity is 0 everywhere on x, i.e., u ( x ) 0, the system has an unstable fixed point from which SSMs originate.

To generate a data set for the discovery of SSMs, we adapt the Matlab code from Ref. 68, which uses a spectral solver for numerical solutions of the Kuramoto–Sivashinsky equation. We consider the system size L = 22, for which a truncation of the number of Fourier modes to the range 16 N 128 yields accurate solutions for this system size,69 therefore, we pick N = 64. We discretize L into 1024 segments, generate trajectories with the time step 0.25, and discard the first 100 time units to get trajectories close to the attractor. The training set which is generated for the time length 1.25 × 10 5 contains approximately 5 × 10 5 data points. In this example, we only test the short-term predictability of the SSM-reduced model and therefore generate a short test set for a 100 time interval.

For SSMLearn, computing the parametrization of an SSM requires solving the implicit minimization in Eq. (15), which can be computationally demanding for high-dimensional systems. Therefore, we turn to the fastSSM algorithm,70 which significantly reduces the computational cost when the observable space dimension is high. The fastSSM algorithm calculates the tangent space of the SSM by singular value decomposition and the manifold parametrization via explicit polynomial regression. As a result, the computational cost scales linearly with the input data dimension. We refer to Ref. 70 for a detailed discussion.

To estimate the SSM dimension from data, we plot the manifold fitting error as a function of the SSM dimension and SSM polynomial expansion order as shown in Fig. 12. When we increase the SSM dimension to eight, there is a significant drop in manifold fitting error. This indicates the chaotic attractor of the system can be captured by an inertial manifold with a dimensionality of at least eight. Therefore, we fit an eight-dimensional SSM with a third order polynomial expansion to the 1024-dimensional trajectory data. The invariance error of the manifold fitting is 0.35 % for the test data. The estimated dimension d = 8 agrees with the estimations in Refs. 71 and 69 and other model reduction results on inertial manifolds.24,72

FIG. 12.

Manifold fitting error for different SSM dimensions, as functions of the SSM fitting order. From this figure, we fit a d = 8 dimensional SSM with a K = 3rd order expansion to the 1024-dimensional input data, which gives an invariance error 0.35 % for the test set.

FIG. 12.

Manifold fitting error for different SSM dimensions, as functions of the SSM fitting order. From this figure, we fit a d = 8 dimensional SSM with a K = 3rd order expansion to the 1024-dimensional input data, which gives an invariance error 0.35 % for the test set.

Close modal

After obtaining the data in 8D reduced coordinates, we perform the kNN algorithm with k = 9 to make forecasts for the test trajectory. The forecasting quality and the validity range of each prediction depend on the amount of training data used and the individual initial condition of the test trajectory. We only illustrate one of the SSM prediction results in Fig. 13. In general, the normalized trajectory errors reach 0.5 between 10 and 60 time intervals. A detailed analysis of Lyapunov exponents and average predictability for this infinite-dimensional example is beyond the scope of the present paper.

FIG. 13.

Test trajectory and its prediction by the SSM-reduced model with system size L = 22. The top two plots are the velocity u = u ( x , t ) and its prediction u ~ = u ~ ( x , t ) for the time interval ( 0 , 90 ). The SSM prediction resembles the real chaotic flow pattern but gradually diverges from it. The bottom two plots show the pointwise prediction error and the time evolution of the normalized trajectory error.

FIG. 13.

Test trajectory and its prediction by the SSM-reduced model with system size L = 22. The top two plots are the velocity u = u ( x , t ) and its prediction u ~ = u ~ ( x , t ) for the time interval ( 0 , 90 ). The SSM prediction resembles the real chaotic flow pattern but gradually diverges from it. The bottom two plots show the pointwise prediction error and the time evolution of the normalized trajectory error.

Close modal

In our final and most complex example, we use data-driven model reduction to a mixed-mode SSM to identify a low-order model for the chaotic attractor of a periodically forced buckling beam. The nonlinear finite-element model we use for a von Kármán beam is a second-order approximation of the kinematics of a beam, taking into account the deformation of the cross-section for both bending and stretching. We specifically adopt here the finite-element Matlab code developed in Ref. 73 and the same buckling beam setting as in Ref. 37. The boundary condition of the beam is set as pinned–pinned, as shown in Fig. 14, and an axial compressive force is applied to the right node.

FIG. 14.

Dynamic buckling of a pinned–pinned von Kármán beam, with the subsequent addition of a time-periodic vertical load at its midpoint.

FIG. 14.

Dynamic buckling of a pinned–pinned von Kármán beam, with the subsequent addition of a time-periodic vertical load at its midpoint.

Close modal

The finite-element model of the von Kármán beam consists of 13 nodes, each with three degrees of freedom (DOFs): axial direction, transverse direction, and rotation. Under the above boundary conditions, the beam model has 12 elements and 36 degrees of freedom. With the displacement and velocity of each DOF counted, the full model is a 72-dimensional dynamical system in its phase space. We adopt the physical parameters listed in Table VI from Ref. 73.

TABLE VI.

Physical parameters for the von Kármán beam.

Symbol Meaning Value [Unit]
t  Time  …  [s] 
L  Length of the beam  [m] 
h  Height of the beam  1 × 10−2  [m] 
b  Width of the beam  5 × 10−2  [m] 
A  Width of the beam  5 × 10−4  [m2
E  Young’s modulus  190 × 109  [Pa] 
κ  Viscous damping rate of material  7 × 106  [Pa⋅s] 
ρ  Density  7850  [kg/m3
Symbol Meaning Value [Unit]
t  Time  …  [s] 
L  Length of the beam  [m] 
h  Height of the beam  1 × 10−2  [m] 
b  Width of the beam  5 × 10−2  [m] 
A  Width of the beam  5 × 10−4  [m2
E  Young’s modulus  190 × 109  [Pa] 
κ  Viscous damping rate of material  7 × 106  [Pa⋅s] 
ρ  Density  7850  [kg/m3
Euler’s critical load for buckling is given by P n = n 2 π 2 E I l 2. For the first buckling mode n = 1, the critical buckling force is
F buckling = π 2 E I l 2 = π 2 E b h 3 12 l 2 .
(35)

For the buckling force value F = 1.1 F buckling, the system has three fixed points: one unstable fixed point corresponding to purely axial displacement and two stable fixed points corresponding to the upper and lower buckled states. The first 10 eigenvalues of the unstable fixed points are listed in Table VII. Therefore, each unstable steady state has one pair of real eigenvalues 11.06 and 11.10, the other pairs correspond to stable spirals. The positive real eigenvalue 11.06 has the physical meaning of the beam leaving the unstable fixed point and approaching the upper or lower buckled states, and the negative real eigenvalue 11.10 corresponds to a trajectory settling down at the unstable state. The other complex pairs of eigenvalues represent the dynamics of higher modes. According to the mixed-mode SSM theory,37 there exists a two-dimensional SSM tangent to the spectral subspaces corresponding to the two real eigenvalues. As in our earlier examples, these eigenvalues will not be used in our SSM construction.

TABLE VII.

First 10 eigenvalues of the buckled von Kármán beam for reference.

λ1 λ2 λ3,4 λ5,6 λ7,8 λ9,10
11.06  −11.10  −0.36 ± 119.36i  −1.83 ± 295.56i  −5.80 ± 541.50i  −14.19 ± 858.19i 
λ1 λ2 λ3,4 λ5,6 λ7,8 λ9,10
11.06  −11.10  −0.36 ± 119.36i  −1.83 ± 295.56i  −5.80 ± 541.50i  −14.19 ± 858.19i 

We generate four training trajectories from the autonomous beam model (without the vertical forcing shown in Fig. 14) for the non-dimensionalized time interval ( 0 , 7.89 ), which covers 150 times the oscillating period of the slowest mode. Each trajectory contains 1.5 × 10 4 points. The initial conditions are generated either by adding a small perturbation to the unstable stationary state, or by applying a large transverse force ( 10 4 N) at the middle node. After applying this force, the beam is further down than the lower buckled state and will first oscillate around the two stable fixed points with a large amplitude, then converge to one of them. The higher modes corresponding to the complex eigenvalue pairs die out, and the trajectories converge to a 2D mixed-mode SSM, which we locate from the training trajectories using SSMLearn. Converged portions of two training trajectories and the SSM are shown in Fig. 15. Based on invariance error and NMTE for non-chaotic test data (one trajectory) in Table VIII, we construct a seventh-order SSM expansion and a 9th-order polynomial representation of the reduced dynamics in SSMLearn.

FIG. 15.

Non-chaotic training trajectories of the buckling beam model on a 2D, autonomous SSM.

FIG. 15.

Non-chaotic training trajectories of the buckling beam model on a 2D, autonomous SSM.

Close modal
TABLE VIII.

Different choices of SSM expansion order and their corresponding invariance error for non-chaotic test data of the buckling beam. After selecting the manifold expansion order to be K = 7, we choose the SSM-model polynomial order to be ninth according to their NMTE.

Manifold expansion order ( K ) 
Invariance error  0.10%  0.10%  0.14%  0.14%  0.097% 
Reduced model order (SSM expansion order K = 7 10  11 
NMTE for non-chaotic test data  0.24%  0.24%  0.12%  0.12%  0.16% 
Manifold expansion order ( K ) 
Invariance error  0.10%  0.10%  0.14%  0.14%  0.097% 
Reduced model order (SSM expansion order K = 7 10  11 
NMTE for non-chaotic test data  0.24%  0.24%  0.12%  0.12%  0.16% 

Next, we apply transverse periodic forcing on the middle node of the beam, as shown in Fig. 14, which makes the beam oscillate chaotically. For moderate forcing amplitudes, the autonomous SSM, W ( E ), is known to persist in the form of a nearby, time-periodic SSM, W ( E , t ).37 The leading-order forced dynamics on W ( E , t ) can be obtained by simply adding to the reduced dynamics of W ( E ), the projection of the time-periodic forcing on E.37 This fact enables one to make predictions for the forced response based solely on a model trained on unforced data.

This procedure is implemented in SSMLearn38,44 and will be discussed here. The periodic forcing in Fig. 14 has the amplitude F 0 = 3.4 N and the frequency Ω = 15.4 rad/s. We then generate the chaotic response of the buckling beam to this periodic force, and this one chaotic trajectory as the test data for our reduced-order modeling. Figure 16 shows this test trajectory lies entirely on the 2D SSM and can be projected onto the two-dimensional reduced coordinates.

FIG. 16.

(a) Chaotic response of the beam under periodic vertical forcing in the full phase space. This chaotic trajectory lies on the 2D SSM that we computed earlier from the non-chaotic data. (b) Chaotic trajectories of one testing data in two-dimensional reduced coordinates.

FIG. 16.

(a) Chaotic response of the beam under periodic vertical forcing in the full phase space. This chaotic trajectory lies on the 2D SSM that we computed earlier from the non-chaotic data. (b) Chaotic trajectories of one testing data in two-dimensional reduced coordinates.

Close modal
Specifically, up to leading order, the forced, SSM-reduced model on 2D SSM, W ( E , t ), is in the form37 
η ˙ = R ( η ) + ( A cos Ω t sin ϕ A cos Ω t cos ϕ ) .
(36)

After calibrating the forcing amplitude A and forcing phase ϕ to a single forced trajectory (see Ref. 38 for details), we find that the SSM-reduced model computed from the autonomous system is indeed able to forecast the chaotic forced response, as shown in Fig. 17. The Lyapunov exponent of the reconstructed system 1.732 is very close to its counterpart 1.742 calculated for the full, 72-dimensional system. Thus, the statistical properties of the chaotic attractor are successfully predicted by an SSM-reduced model in this example too.

FIG. 17.

(a) Prediction for a chaotic beam trajectory from an SSM-reduced model trained on unforced, non-chaotic data, then appended with the leading-order forcing term. The prediction lasts for about 3.31 Lyapunov time. (b) Separation of nearby trajectories against time. The original 72-dimensional system (red) and the reduced two-dimensional model (blue) exhibit nearly equal Lyapunov exponents (1.742 and 1.732).

FIG. 17.

(a) Prediction for a chaotic beam trajectory from an SSM-reduced model trained on unforced, non-chaotic data, then appended with the leading-order forcing term. The prediction lasts for about 3.31 Lyapunov time. (b) Separation of nearby trajectories against time. The original 72-dimensional system (red) and the reduced two-dimensional model (blue) exhibit nearly equal Lyapunov exponents (1.742 and 1.732).

Close modal

Even though a short-term prediction shown in Fig. 17(a) quickly loses its accuracy, one should not forget: This prediction by the 2D, SSM-reduced model is made entirely based on unforced, non-chaotic trajectories and yet it does approximate both the geometry and the long-term statistics (as described by the maximal Lyapunov exponent) with notable accuracy for the chaotic attractor of the full, 72-dimensional system. This shows the high physical interpretability of SSM-reduced models, which other data-driven model reduction techniques generally do not provide. For instance, the SINDy algorithm we used in some of our previous examples was not realistic to run on this high-dimensional example and the algorithm would not be equipped to make forced predictions from unforced training data.

Using the recent theory of mixed-mode spectral submanifolds,37 we have developed a methodology to derive SSM-reduced models for dynamical systems with chaotic attractors. Specifically, we combined the SSMLearn algorithm with the method of nearest neighbors to obtain an accurate enough representation of the SSM-reduced dynamics on large enough domains that capture the attractors of the system. In order to contain the chaotic attractor and the dominant behavior of the underlying system, the SSM must have a dimension higher than the attractor dimension. In a strictly data-driven setting, the SSM dimension can be approximated using the false nearest neighbor method, then further determined by computing the invariance error.

This methodology is generally applicable to systems in which an inertial manifold containing the chaotic attractor is, in fact, an SSM emanating from an unstable steady state. In the examples treated here, that steady state was either a fixed point or a periodic orbit. Building the inertial manifold specifically as a spectral submanifold enables us to construct mathematically justified, low-dimensional models for the dynamics on and near the attractors. The models are either explicit polynomial ODEs or ODEs with right-hand sides interpolated using nearest neighbor points on the SSM. The core part of our approach is the open-source SSMLearn algorithm,44 a general Matlab (Python) package for SSM-based model reduction. We additionally deploy the false nearest neighbor (FNN) algorithm74 to identify the optimal SSM dimension from data. For SSM-based dimension reduction, no information about the time derivative data is required, and the amount of data needed is relatively low.

We have applied two methods to model the reduced dynamics: a global polynomial fit and the kNN method. Both methods have limitations. The global polynomial fit requires computing the time derivatives of the trajectories and can become unstable for complex systems. The kNN method only uses local information to make predictions and therefore requires a large amount of data that densely covers the whole attractor. To counter these challenges, it is possible to combine SSM-model reduction with other prediction methods, provide the weighting functions,75,76 use other function bases such as radial basis functions,77 or employ machine learning and deep learning methods78,79 such as neural ODEs.24 

Our SSM-based attractor models have yielded accurate local predictions over short Lyapunov times and also reproduced closely the statistical and ergodic properties of chaotic attractors over longer times. Examples supporting these conclusions included the classic Lorenz and Rössler system, the extended 9D Lorenz model, a forced oscillator chain, the Kuramoto–Sivashinsky equation, and a periodically forced buckling beam.

We also compare our approach with the data-driven nonlinear modeling method SINDy.27 While the advantages of SINDy include simplicity, ease of implementation, and versatility, the method scales unfavorably with the input data dimension, and depends on the choice of the function library and the coordinates, as discussed in Sec. IV D. In efforts to address these challenges, SINDy has also been combined with various coordinate-identifying methods and model reduction techniques, such as learning appropriate observables in an assumed latent space. For SSM-based model reductions, the existence of such low-dimensional invariant manifolds is theoretically justified in the presence of a hyperbolic steady state. The requirement of such a steady state is a limitation of our method, but in our view, it is necessary for the localization of invariant manifolds from data to be on firm theoretical ground.

Importantly, since SSM-reduced models are fully dynamics-based, they can predict forced chaotic responses based solely on the knowledge of unforced (non-chaotic) trajectory data as demonstrated by our forced buckling beam example. We are unaware of any other data-driven model reduction approach that has this capacity.

Future applications of our methodology will include other PDE models, such as the 2D Kolmogorov flow,80,81 which also has finite-dimensional attractors and coexisting unstable steady states. Further extensions of mixed-mode SSM theory to more general anchoring steady states (such as invariant tori and other bounded invariant sets) are needed for the present methodology to become applicable to inertial manifolds that emanate from invariant sets other than fixed points or periodic orbits.

We acknowledge helpful discussions with Mattia Cenedese, Mingwu Li, and Shobhit Jain. A.L. was supported by the Swiss National Foundation (SNF) Grant No. 200021_214908.

The authors have no conflicts to disclose.

Aihui Liu: Formal analysis (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Joar Axås: Conceptualization (equal); Formal analysis (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). George Haller: Conceptualization (lead); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing – review & editing (lead).

The code used to generate the numerical results included in this paper is available as part of the open-source Matlab package SSMLearn at https://github.com/haller-group/SSMLearn, Ref. 44.

1.
S. H.
Strogatz
,
Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering
(
Westview Press
,
2000
).
2.
G. L.
Baker
and
J. P.
Gollub
,
Chaotic Dynamics: An Introduction
, 2nd ed. (
Cambridge University Press
,
1996
).
3.
R. L.
Devaney
,
An Introduction to Chaotic Dynamical Systems
(
Chapman and Hall/CRC
,
2021
).
4.
H.
Whitney
, “
Differentiable manifolds
,”
Ann. Math.
37
,
645
680
(
1936
).
5.
F.
Takens
, “Detecting strange attractors in turbulence,” in Lecture Notes in Mathematics (Springer, Berlin, 1981), pp. 366–381.
6.
W.
Gilpin
, “Deep reconstruction of strange attractors from time series,” arXiv:2002.05909 [cs.LG] (2020).
7.
K. H.
Kraemer
,
G.
Datseris
,
J.
Kurths
,
I.
Kiss
,
J.
Ocampo-Espindola
, and
N.
Marwin
, “
A unified and automated approach to attractor reconstruction
,”
New J. Phys.
23
,
033017
(
2021
).
8.
Z.
Lu
,
B. R.
Hunt
, and
E.
Ott
, “
Attractor reconstruction by machine learning
,”
Chaos
28
,
061104
(
2018
).
9.
A. M.
Fraser
and
H. L.
Swinney
, “
Independent coordinates for strange attractors from mutual information
,”
Phys. Rev. A
33
,
1134
1140
(
1986
).
10.
C.
Rhodes
and
M.
Morari
, “
The false nearest neighbors algorithm: An overview
,”
Comput. Chem. Eng.
21
,
S1149
S1154
(
1997
).
11.
I.
Scherl
,
B.
Strom
,
J. K.
Shang
,
O.
Williams
,
B. L.
Polagye
, and
S. L.
Brunton
, “
Robust principal component analysis for modal decomposition of corrupt fluid flows
,”
Phys. Rev. Fluids
5
,
054401
(
2020
).
12.
P. J.
Schmid
, “
Dynamic mode decomposition and its variants
,”
Annu. Rev. Fluid Mech.
54
,
225
254
(
2022
).
13.
L.
van der Maaten
,
E.
Postma
, and
H.
Herik
, “
Dimensionality reduction: A comparative review
,”
J. Mach. Learn. Res.
10
,
66
71
(
2007
).
14.
C.
Touzé
,
A.
Vizzaccaro
, and
O.
Thomas
, “
Model order reduction methods for geometrically nonlinear structures: A review of nonlinear techniques
,”
Nonlinear Dyn.
105
,
1141
1190
(
2021
).
15.
M.
Sangiorgio
,
F.
Dercole
, and
G.
Guariso
,
Deep Learning in Multi-step Prediction of Chaotic Dynamics
(
Springer International Publishing
,
2021
).
16.
S.
Shahi
,
F. H.
Fenton
, and
E. M.
Cherry
, “
Prediction of chaotic time series using recurrent neural networks and reservoir computing techniques: A comparative study
,”
Mach. Learn. Appl.
8
,
100300
(
2022
).
17.
K.
Lee
and
K. T.
Carlberg
, “
Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders
,”
J. Comput. Phys.
404
,
108973
(
2020
).
18.
C.
Foias
,
G. R.
Sell
, and
R.
Temam
, “
Inertial manifolds for nonlinear evolutionary equations
,”
J. Differ. Equ.
73
,
309
353
(
1988
).
19.
A. J.
Linot
and
M. D.
Graham
, “
Deep learning to discover and predict dynamics on an inertial manifold
,”
Phys. Rev. E
101
,
062209
(
2020
).
20.
D.
Floryan
and
M.
Graham
, “
Data-driven discovery of intrinsic dynamics
,”
Nat. Mach. Intell.
4
,
1113
1120
(
2022
).
21.
H. D. I.
Abarbanel
,
Analysis of Observed Chaotic Data
(
Springer New York
,
1996
).
22.
H. D. I.
Abarbanel
,
R.
Brown
,
J. J.
Sidorowich
, and
L. S.
Tsimring
, “
The analysis of observed chaotic data in physical systems
,”
Rev. Mod. Phys.
65
,
1331
1392
(
1993
).
23.
J. D.
Farmer
and
J. J.
Sidorowich
, “
Predicting chaotic time series
,”
Phys. Rev. Lett.
59
,
845
848
(
1987
).
24.
A. J.
Linot
and
M. D.
Graham
, “
Data-driven reduced-order modeling of spatiotemporal chaos with neural ordinary differential equations
,”
Chaos
32
,
073110
(
2022
).
25.
P. R.
Vlachas
,
W.
Byeon
,
Z. Y.
Wan
,
T. P.
Sapsis
, and
P.
Koumoutsakos
, “
Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
474
,
20170844
(
2018
).
26.
H.
Eivazi
,
H.
Veisi
,
M. H.
Naderi
, and
V.
Esfahanian
, “
Deep neural networks for nonlinear model order reduction of unsteady flows
,”
Phys. Fluids
32
,
105104
(
2020
).
27.
S. L.
Brunton
,
J. L.
Proctor
, and
J. N.
Kutz
, “
Discovering governing equations from data by sparse identification of nonlinear dynamical systems
,”
Proc. Natl. Acad. Sci. U.S.A.
113
,
3932
3937
(
2016
).
28.
U.
Fasel
,
J. N.
Kutz
,
B. W.
Brunton
, and
S. L.
Brunton
, “
Ensemble-sindy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
478
,
20210904
(
2022
).
29.
A. A.
Kaptanoglu
,
L.
Zhang
,
Z. G.
Nicolaou
,
U.
Fasel
, and
S. L.
Brunton
, “
Benchmarking sparse system identification with low-dimensional chaos
,”
Nonlinear Dyn.
111
,
13143
13164
(
2023
).
30.
P. J.
Baddoo
,
B.
Herrmann
,
B. J.
McKeon
,
J.
Nathan Kutz
, and
S. L.
Brunton
, “
Physics-informed dynamic mode decomposition
,”
Proc. R. Soc. A: Math. Phys. Eng. Sci.
479
,
20220576
(
2023
).
31.
B.
de Silva
,
D. M.
Higdon
,
S. L.
Brunton
, and
J. N.
Kutz
, “
Discovery of physics from data: Universal laws and discrepancies
,”
Front. Artif. Intell.
3
,
25
(
2019
).
32.
G.
Haller
and
S.
Ponsioen
, “
Nonlinear normal modes and spectral submanifolds: Existence, uniqueness and use in model reduction
,”
Nonlinear Dyn.
86
,
1493
1534
(
2016
).
33.
M.
Cenedese
,
J.
Axås
,
H.
Yang
,
M.
Eriten
, and
G.
Haller
, “
Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems
,”
Philos. Trans. R. Soc. A
380
,
20210194
(
2022
).
34.
J. I.
Alora
,
M.
Cenedese
,
E.
Schmerling
,
G.
Haller
, and
M.
Pavone
, “Data-driven spectral submanifold reduction for nonlinear optimal control of high-dimensional robots,” in 2023 IEEE Int. Conf. Robot. Autom. (IEEE, 2022), pp. 2627–2633.
35.
F.
Mahlknecht
,
J. I.
Alora
,
S.
Jain
,
E.
Schmerling
,
R.
Bonalli
,
G.
Haller
, and
M.
Pavone
, “Using spectral submanifolds for nonlinear periodic control,” in 2022 IEEE 61st Conference on Decision and Control (CDC) (IEEE, 2022), pp. 6548–6555.
36.
J.
Axås
and
G.
Haller
, “
Model reduction for nonlinearizable dynamics via delay-embedded spectral submanifolds
,”
Nonlinear Dyn.
111
,
22079
22099
(
2023
).
37.
G.
Haller
,
B.
Kaszás
,
A.
Liu
, and
J.
Axås
, “
Nonlinear model reduction to fractional and mixed-mode spectral submanifolds
,”
Chaos
33
,
063138
(
2023
).
38.
M.
Cenedese
,
J.
Axås
,
B.
Bäuerlein
,
K.
Avila
, and
G.
Haller
, “
Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds
,”
Nat. Commun.
13
,
872
(
2022
).
39.
S.
Jain
,
M.
Li
,
T.
Thurnher
, and
G.
Haller
(2023). “SSMTool: Computation of invariant manifolds in high-dimensional mechanics problems,”
Zenodo
,
v2.5
.
40.
T.
Sauer
,
J. A.
Yorke
, and
M.
Casdagli
, “
Embedology
,”
J. Stat. Phys.
65
,
579
616
(
1991
).
41.
R.
Temam
, “
Inertial manifolds
,”
Math. Intell.
12
,
68
74
(
1990
).
42.
J. C.
Robinson
, “Dimensions, embeddings, and attractors,” in Cambridge Tracts in Mathematics (Cambridge University Press, 2010).
43.
H.
Kantz
and
T.
Schreiber
,
Nonlinear Time Series Analysis
, 2nd ed. (
Cambridge University Press
,
2003
).
44.
M. Cenedese and J. Axås (2022). “SSMLearn,” GitHub. https://github.com/haller-group/SSMLearn.
45.
D.
Ruppert
, “Local polynomial regression and its applications in environmental statistics,” in
Statistics for the Environment
, Vol. 3 (Wiley, 1997), pp. 155–173.
46.
L.-Y.
Su
, “
Prediction of multivariate chaotic time series with local polynomial fitting
,”
Comput. Math. Appl.
59
,
737
744
(
2010
).
47.
D.
Yankov
,
D.
Decoste
, and
E.
Keogh
,
Ensembles of Nearest Neighbor Forecasts
(
Springer, Berlin
,
2006
), Vol. 4212, pp. 545–556.
48.
F.
Martínez
,
M.
Frías
,
M.
Pérez
, and
A.
Rivera Rivas
, “
A methodology for applying k-nearest neighbor to time series forecasting
,”
Artif. Intell. Rev.
52
,
2019
2037
(
2019
).
49.
B.
Finkenstädt
, “A nearest neighbor approach to forecast nonlinear time series,” in Nonlinear Dynamics in Economics: A Theoretical and Statistical Approach to Agricultural Markets (Springer, Berlin, 1995), pp. 122–147.
50.
A.
Wolf
,
J.
Swift
,
H.
Swinney
, and
J.
Vastano
, “
Determining Lyapunov exponents from a time series
,”
Phys. D
16
,
285
317
(
1985
).
51.
D.
Lathrop
, “
Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering
,”
Phys. Today
68
,
54
55
(
2015
).
52.
T.
Tél
,
M.
Gruiz
, and
K.
Kulacsy
,
Chaotic Dynamics: An Introduction Based on Classical Mechanics
(
Cambridge University Press
,
2006
), pp. 1–408.
53.
M. T.
Rosenstein
,
J. J.
Collins
, and
C. J. D.
Luca
, “
A practical method for calculating largest lyapunov exponents from small data sets
,”
Phys. D: Nonlinear Phenom.
65
,
117
134
(
1993
).
54.
J. B.
Dingwell
, “Lyapunov exponents,” in Wiley Encyclopedia of Biomedical Engineering (John Wiley & Sons, Ltd, 2006).
55.
P.
Gaspard
, “Chaos, scattering and statistical mechanics,” in Cambridge Nonlinear Science Series (Cambridge University Press, 1998).
56.
D.
Ruelle
, “Chaotic evolution and strange attractors,” in Lezioni Lincee (Cambridge University Press, 1989).
57.
L. M.
Berliner
, “
Statistics, probability and chaos
,”
Stat. Sci.
7
,
69
90
(
1992
).
58.
M.
Dörfle
and
R.
Graham
, “
Probability density of the Lorenz model
,”
Phys. Rev. A
27
,
1096
1105
(
1983
).
59.
R.
Graham
and
H. J.
Scholz
, “
Analytic approximation of the Lorenz attractor by invariant manifolds
,”
Phys. Rev. A
22
,
1198
1204
(
1980
).
60.
J.-P.
Eckmann
and
D.
Ruelle
, “
Ergodic theory of chaos and strange attractors
,”
Rev. Mod. Phys.
57
,
617
656
(
1985
).
61.
E. N.
Lorenz
, “
Deterministic nonperiodic flow
,”
J. Atmos. Sci.
20
,
130
148
(
1963
).
62.
P.
Reiterer
,
C.
Lainscsek
,
F.
Schuerrer
,
C.
Letellier
, and
J.
Maquet
, “
A nine-dimensional Lorenz system to study high-dimensional chaos
,”
J. Phys. A
31
,
7121
7139
(
1998
).
63.
A.
Kaptanoglu
,
B.
de Silva
,
U.
Fasel
,
K.
Kaheman
,
A.
Goldschmidt
,
J.
Callaham
,
C.
Delahunt
,
Z.
Nicolaou
,
K.
Champion
,
J.-C.
Loiseau
,
J.
Kutz
, and
S.
Brunton
, “
PySINDy: A comprehensive python package for robust sparse system identification
,”
J. Open Source Softw.
7
,
3994
(
2022
).
64.
I.
Kovacic
and
M. J.
Brennan
, “Background: On georg duffing and the duffing equation,” in The Duffing Equation (John Wiley & Sons, Ltd, 2011), Chap. 1, pp. 1–23.
65.
Z.
Ahsan
,
H.
Dankowicz
,
M.
Li
, and
J.
Sieber
, “
Methods of continuation and their implementation in the COCO software platform with application to delay differential equations
,”
Nonlinear Dyn.
107
,
3181
3243
(
2022
).
66.
J. J.
Bramburger
and
J. N.
Kutz
, “
Poincaré maps for multiscale physics discovery and nonlinear floquet theory
,”
Phys. D: Nonlinear Phenom.
408
,
132479
(
2020
).
67.
O.
Rössler
, “
An equation for continuous chaos
,”
Phys. Lett. A
57
,
397
398
(
1976
).
68.
P.
Cvitanović
, “Turbulence,” in Chaos: Classical and Quantum (Niels Bohr Inst., Copenhagen, 2016), Chap. 30.
69.
P.
Cvitanović
,
R. L.
Davidchack
, and
E.
Siminos
, “
On the state space geometry of the Kuramoto–Sivashinsky flow in a periodic domain
,”
SIAM J. Appl. Dyn. Syst.
9
,
1
33
(
2010
).
70.
J.
Axås
,
M.
Cenedese
, and
G.
Haller
, “
Fast data-driven model reduction for nonlinear dynamical systems
,”
Nonlinear Dyn.
111
,
7941
7957
(
2022
).
71.
X.
Ding
,
H.
Chaté
,
P.
Cvitanović
,
E.
Siminos
, and
K. A.
Takeuchi
, “
Estimating the dimension of an inertial manifold from unstable periodic orbits
,”
Phys. Rev. Lett.
117
,
024101
(
2016
).
72.
P. R.
Vlachas
,
G.
Arampatzis
,
C.
Uhler
, and
P.
Koumoutsakos
, “
Multiscale simulations of complex systems by learning their effective dynamics
,”
Nat. Mach. Intell.
4
,
359
366
(
2022
).
73.
S.
Jain
,
P.
Tiso
, and
G.
Haller
, “
Exact nonlinear model reduction for a von Kármán beam: Slow-fast decomposition and spectral submanifolds
,”
J. Sound Vib.
423
,
195
211
(
2018
).
74.
H. D. I.
Abarbanel
and
M. B.
Kennel
, “
Local false nearest neighbors and dynamical dimensions from observed chaotic data
,”
Phys. Rev. E
47
,
3057
3068
(
1993
).
75.
J. P.
Crutchfield
and
B. S.
McNamara
, “
Equations of motion from a data series
,”
Complex Syst.
1
,
417
452
(
1987
).
76.
M.
Giona
,
F.
Lentini
, and
V.
Cimagalli
, “
Functional reconstruction and local prediction of chaotic time series
,”
Phys. Rev. A
44
,
3496
3502
(
1991
).
77.
M.
Casdagli
, “
Nonlinear prediction of chaotic time series
,”
Phys. D: Nonlinear Phenom.
35
,
335
356
(
1989
).
78.
D. S.
Karunasinghe
and
S.-Y.
Liong
, “
Chaotic time series prediction with a global model: Artificial neural network
,”
J. Hydrol.
323
,
92
105
(
2006
).
79.
M.
Sangiorgio
, “Deep learning in multi-step forecasting of chaotic dynamics,” in Special Topics in Information Technology, edited by L. Piroddi (Springer International Publishing, Cham, 2022), pp. 3–14.
80.
N.
Platt
,
L.
Sirovich
, and
N.
Fitzmaurice
, “
An investigation of chaotic Kolmogorov flows
,”
Phys. Fluids A
3
,
681
696
(
1991
).
81.
G. J.
Chandler
and
R. R.
Kerswell
, “
Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow
,”
J. Fluid Mech.
722
,
554
595
(
2013
).
Published open access through an agreement with ETH Zurich