Globalizing manifold-based reduced models for equations and data

Abstract

One of the very few mathematically rigorous nonlinear model reduction methods is the restriction of a dynamical system to a low-dimensional, sufficiently smooth, attracting invariant manifold. Such manifolds are usually found using local polynomial approximations and, hence, are limited by the unknown domains of convergence of their Taylor expansions. To address this limitation, we extend local expansions for invariant manifolds via Padé approximants, which re-express the Taylor expansions as rational functions for broader utility. This approach significantly expands the range of applicability of manifold-reduced models, enabling reduced modeling of global phenomena, such as large-scale oscillations and chaotic attractors of finite element models. We illustrate the power of globalized manifold-based model reduction on several equation-driven and data-driven examples from solid mechanics and fluid mechanics.

Manifold-based model reduction often uses power series expansions that fail to converge on larger domains. Here, authors extend the reach of such approximations using global rational approximants. This enables the reduced models to capture previously inaccessible nonlinear dynamics.

Introduction

Reduced-order modeling is a widespread technique that seeks to simplify high-dimensional nonlinear systems while retaining their essential dynamical features. Among reduced-order modeling procedures, manifold-based methods have been steadily gaining momentum. This can be largely attributed to the prevalence of data-driven approaches that successfully build on the fundamental manifold hypothesis^{1, 2–3} of machine learning.

Center manifold reduction⁴, geometric singular perturbation theory^5,6 and inertial manifold theory^7,8 all rely on the existence of low-dimensional attracting invariant manifolds in the phase space of a dynamical system. These methods constitute mathematically rigorous examples of nonlinear model reduction and yield truly predictive models. However, for systems encountered in practice that are not close to bifurcations, invariant manifolds can only be realistically constructed when they emanate from a known robust stationary state, such as a hyperbolic fixed point.

In those cases, seeking the invariant manifolds perturbatively and expressing them as local Taylor expansions at the known stationary state is justified. Traditionally, only stable and unstable manifolds, as continuations of the stable and unstable subspaces of the linearization, were approximated in this fashion⁴. The recent theory of spectral submanifolds (SSMs) extends this approach to arbitrary non-resonant spectral subspaces of the linearized system⁹. In particular, SSMs are now known to exist as smooth continuations of stable subspaces (like-mode SSMs) and of subspaces spanned by stable and unstable modes (mixed-mode SSMs)¹⁰. This model reduction approach has been used in a broad range of physical settings to deduce very low-dimensional, mathematically justified polynomial models^9,10.

SSM reduction has been successfully applied to obtain accurate reduced models of nonlinear vibrations observed in high-dimensional finite element models^11,12 and experiments^13,14, multistable fluid flows^15,16, chaotic systems¹⁷, fluid-structure interaction problems^18,19 and control of soft robots²⁰. In an equation-driven setting, SSM reduction starts with the solution of an invariance equation through local Taylor expansions^11,21,22. Data-driven SSM reduction¹³ also uses polynomial expansions to obtain an approximation for SSMs.

SSMs are ideal tools for model reduction because their existence, uniqueness, and smoothness properties are precisely known. Specifically, if the governing equations are analytic, SSMs of attracting fixed points are guaranteed to also possess convergent Taylor series near the fixed point. However, the domain of convergence is generally unknown.

This has been a general limitation of invariant manifold-based reduction methods, restricting the range of applicability of the resulting reduced-order model to an a priori unknown domain. Importantly, the convergence-limiting singularity of the Taylor expansion is not a physical singularity of the dynamical system, and hence, invariant manifolds continue to exist even outside the domain of convergence of the Taylor expansions for those manifolds. Therefore, the breakdown of convergence arises without any prior indication. This major limitation of local approximation methods of invariant manifolds restricts the user to potentially small physical domains.

Here, we overcome this limitation by extending the local information contained in the Taylor series to considerably larger domains via the process generally known as analytic continuation²³. Among commonly used methods of analytic continuation, we focus on Padé approximants²⁴. Padé approximants are rational functions whose Taylor expansion around a point coincides with that of the original function up to a given order, but can represent the original function more efficiently and globally.

Padé approximants have been used in theoretical physics and applied mathematics for summing divergent series²⁵ with applications in cosmology²⁶, quantum electrodynamics²⁷, fluid dynamics^28,29 and solving the Helmholtz equation³⁰. Closest in spirit to our present work is the use of Padé approximants for center manifold reduction^31,32. The latter use is, however, restricted to equation-driven model reduction near Hopf bifurcations.

In the context of SSM-reduced models, we must consider generalizations of Padé approximants to account for multivariate functions describing the parametrization and the reduced dynamics. Since two-dimensional non-resonant SSMs are typical in applications, we focus on the bivariate and univariate cases.

In “Results”, we present our results on using Padé approximants to construct reduced models on global SSMs (gSSMs). We discuss four examples, including the Kolmogorov flow³³, a nonlinear von Kármán beam in periodic and chaotic regimes^11,17, and a data-driven model of an inverted flag experiment¹⁸. The mathematical details of SSMs and Padé approximants are discussed in “Methods”. The Supplementary Information contains further applications and examples.

Results

Spectral submanifolds and Padé approximants

We consider an n − dimensional nonlinear dynamical system

1

$\stackrel{°}{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{f}\left(\mathbf{x}\right)+\epsilon {\mathbf{f}}_{\mathrm{ext}}\left(\mathbf{x},t\right),\mathbf{x}\in {\mathbb{R}}^n,\mathbf{A}\in {\mathbb{R}}^{n\times n},\mathbf{f}\in \mathcal{O}\left(\mid \mathbf{x}{\mid }^2\right),n\ge 1,0\le \epsilon \ll 1$where f_ext(x, t) = f_ext(x, t + T) represents time-periodic external forcing. We assume that the nonlinearity f and the forcing f_ext(x, t) are also analytic functions of x. We assume that, for ε = 0, the fixed point x = 0 is hyperbolic and the spectrum of A is non-resonant. The slow spectral subspace is denoted E, and is defined as the d-dimensional real subspace spanned by the eigenvectors of A associated to its d eigenvalues closest to zero and hence is an attracting, slow invariant manifold for the linearized dynamics.

We focus here on SSMs, which are the nonlinear continuations of spectral subspaces under the addition of the nonlinear terms in (1). Although there are multiple invariant manifolds tangent to the spectral subspace E, there is a unique smoothest one, called the primary spectral submanifold⁹ denoted $\mathcal{W}\left(E\right)$. If the fixed point x = 0 is stable, then $\mathcal{W}\left(E\right)$ is even known to be analytic. Due to the slowness of the subspace E, the SSM, $\mathcal{W}\left(E\right)$, is an attracting slow manifold for the autonomous system (1). The members of the invariant manifold family with reduced smoothness are called fractional (or secondary) SSMs¹⁰. For completeness, the necessary assumptions for the existence and uniqueness of primary SSMs^9,10,34 are recalled in “Spectral submanifolds”.

In the autonomous limit with ε = 0, the d − dimensional (d ≤ n) primary SSM, $\mathcal{W}\left(E\right)$, can locally be represented as the image of a parametrization $\mathbf{W}:U\subset {\mathbb{R}}^d\to {\mathbb{R}}^n$, over some open set $U\subset {\mathbb{R}}^d$ as

2

$\mathcal{W}\left(E\right)=\left\{\mathbf{x}=\mathbf{W}\left(\mathbf{p}\right)\mid \mathbf{p}\in U\right\}\subset {\mathbb{R}}^n.$

The reduced dynamics $\stackrel{°}{\mathbf{p}}=\mathbf{R}\left(\mathbf{p}\right)$, with $\mathbf{R}:U\to {\mathbb{R}}^d$ are conjugate to (1), i.e., $\mathcal{W}\left(E\right)$ is invariant under the time evolution of (1) and therefore satisfies the invariance equation

3

$\mathbf{A}\mathbf{W}\left(\mathbf{p}\right)+\mathbf{f}\left(\mathbf{W}\left(\mathbf{p}\right)\right)=D\mathbf{W}\left(\mathbf{p}\right)\stackrel{°}{\mathbf{p}}.$

We refer to “Spectral submanifolds” for a discussion on SSMs of the non-autonomous system with ε > 0.

We solve Eq. (3) by representing the parametrization of $\mathcal{W}\left(E\right)$ and its reduced dynamics as a power series truncated to some order N, i.e.,

4

$\begin{array}{ccc}\hfill {\mathbf{W}}^N\left(\mathbf{p}\right)& \hfill =\hfill & {\sum }_{\mid \mathbf{k}\mid =0}^N{\mathbf{W}}_{\mathbf{k}}{\mathbf{p}}^{\mathbf{k}},\hfill \\ \hfill {\mathbf{R}}^N\left(\mathbf{p}\right)& \hfill =\hfill & {\sum }_{\mid \mathbf{k}\mid =0}^N{\mathbf{R}}_{\mathbf{k}}{\mathbf{p}}^{\mathbf{k}}.\hfill \end{array}$

We define the multi-index k = (k₁, …, k_d) and ∣k∣ = k₁ + k₂ + … + k_d, so that ${\mathbf{p}}^{\mathbf{k}}=p_1^{k_1}{p}_2^{k_2}\dots p_d^{k_d}$ refers to a scalar monomial of the components of p with total order ∣k∣. The coefficients W_k and R_k are vectors in ${\mathbb{R}}^n$ and ${\mathbb{R}}^d$, respectively, for all k.

The coefficients R_k depend on the style of parametrization used. In the graph-style parametrization, the reduced coordinates are obtained as projections onto the spectral subspace E, while in the normal form style parametrization, non-resonant terms are set to zero. The difference between these two choices is explained in more detail by, e.g., refs. ^11,22,35.

Since the primary autonomous SSM is analytic, there is a domain of convergence denoted as $U_0\subset U\subset {\mathbb{C}}^d$, for which the N → ∞ limit exists, i.e.,

5

${lim}_{N\to \infty }{\mathbf{W}}^N\left(\mathbf{p}\right)=\mathbf{W}\left(\mathbf{p}\right),\forall \mathbf{p}\in U_0.$

For a system whose slowest mode is oscillatory and not in resonance with higher modes, the optimal model reduction is achieved by a two-dimensional SSM tangent to a single oscillatory eigenspace of the autonomous problem. In that case, we can select $\mathbf{p}={\left(p,\overline{p}\right)}^T$ with $p\in \mathbb{C}$. With the normal form style parametrization, the reduced dynamics only contain near-resonant terms of the form $p^{k+1}{\overline{p}}^k$ and $p^k{\overline{p}}^{k+1}$, and it is conveniently expressed in polar coordinates¹¹. We introduce

6

$p=\rho e^{i\theta },\overline{p}=\rho e^{-i\theta },$which allows us to write the SSM-reduced dynamics as

7

$\begin{array}{ccc}\hfill \stackrel{°}{\rho }& \hfill =\hfill & \kappa \left(\rho \right)\rho ,\hfill \\ \hfill \stackrel{°}{\theta }& \hfill =\hfill & \omega \left(\rho \right).\hfill \end{array}$

The functions κ(ρ) and ω(ρ) represent the instantaneous dependence of the damping rate on the amplitude and the frequency on the amplitude, respectively. These functions are obtained from the coefficients of the autonomous reduced dynamics, R_k. Their Taylor expansions are of the form

8

$\kappa \left(\rho \right)={\sum }_{k=0}^{\infty }{\kappa }_k{\rho }^{2k},\omega \left(\rho \right)={\sum }_{k=0}^{\infty }{\omega }_k{\rho }^{2k}.$

The software package SSMTool^11,36 returns the Taylor coefficients W_k, κ_k and ω_k up to arbitrary orders. However, the expansions (8) only converge as long as the amplitude ρ corresponding to the physical response of the system is inside the domain of convergence.

For the expansions in (8), the domain of convergence is the interior of a disk of radius R, the radius of convergence. As a corollary of the analyticity of holomorphic functions²³, that disk contains the nearest singularity of the complex extension of the function. This result also generalizes to the multivariate setting, although defining the domain of convergence is more complicated³⁷. In addition, in contrast to the scalar case, singularities of multivariate functions are never isolated.

The a priori unknown domain of convergence of their local Taylor expansions represents a limitation of SSM-reduction approaches. Although solutions of the invariance equation (21) for SSMs exist up to any order N, their formal sum W^N has little to do with the primary SSM outside the domain of convergence U₀, even if W is well-defined on U\U₀.

This fundamental limitation has forced most invariant manifold studies to focus on deriving reduced-order models under small perturbations near a fixed point. This, however, impedes predicting the system’s response to large perturbations or the discovery of steady states far away from the known fixed points.

However, SSMs, as invariant manifolds, are known to extend globally in the phase space, as long as the flow map of the dynamical system (1) remains well-defined on them for all times. A straightforward extension of the locally known parametrization is to evolve a set of initial conditions from inside the domain of convergence globally under the flow map. This, however, assumes that the flow map of the full high-dimensional system is known in detail, making the reduced-order model redundant. This envisioned globalization is only partly achieved by the data-driven construction of SSMs¹³, which starts from a limited number of trajectories of (1) and finds the observed invariant manifold using regression.

In this work, we propose a different approach to extend the range of applicability of SSM-reduced models. As long as an analytic function is known on some open domain, fundamental results in complex analysis guarantee that this knowledge can be extended to the entire domain of analyticity of the function, possibly using a different representation of the function. Switching to such a representation, other than the Taylor series (4), is known as analytic continuation²³, a powerful technique to globalize the maps W(p) and R(p).

Globalization of invariant manifolds via Padé approximation

A well-known method for analytic continuation is the Padé approximation²⁴, which has had numerous applications in theoretical physics and engineering. To describe this procedure, we introduce a multivariate rational function of the form, using the same notation as in Eq. (4),

9

$\left[N/M\right]\left(\mathbf{z}\right)=\frac{{\sum }_{\mid \mathbf{k}\mid =0}^N{a}_{\mathbf{k}}{\mathbf{z}}^{\mathbf{k}}}{{\sum }_{\mid \mathbf{k}\mid =0}^M{b}_{\mathbf{k}}{\mathbf{z}}^{\mathbf{k}}},b_{\mathbf{0}}=1,\mathbf{z}\in {\mathbb{R}}^{\ell },$where the orders of the numerator and denominator are fixed constants N, M. This formulation covers the cases of univariate (ℓ = 1) and multivariate (ℓ ≥ 2) functions as well.

We represent solutions of the invariance equation (21) as rational functions of the form (9), which is achieved by requiring that the Taylor series of (9) matches that of the function to be approximated. Rational functions of the form (9) have the advantage that they can incorporate singularities that would otherwise limit the convergence of Taylor series. Based on these observations, we seek to extend SSM-reduced models using Padé approximants. We call the extended representations of an SSM the gSSM.

Increasing the orders of the numerator and the denominator in (9), one expects that the approximants provide better approximations for the functions W(p), R(p), κ(ρ), and ω(ρ). In most practical applications, diagonal approximants (i.e., those with numerators and denominators of the same order) have proven to be the most effective. Moreover, diagonal Padé approximants of a univariate function are related to the continued-fraction representation of the function.

For meromorphic functions with an a priori unknown number of poles, strict convergence is only guaranteed in measure by the theorems of refs. ³⁸ and ³⁹. These theorems cover the case of diagonal approximants and state that the sequence [M/M](z) converges to the given function as M → ∞ on bounded compact subsets of $\mathbb{C}$, except for sets of measure zero. Similar theorems also exist for the multivariate case⁴⁰. These exceptional sets correspond to zero sets of the denominator in (9).

Data-driven global reduced-order models

If the equations of motion of a dynamical system are known, we will rely on the established theory of Padé approximants and invariant manifolds and use highly optimized computational methods to solve the invariance equation (21) and construct the gSSM-reduced models.

In many practical applications, however, the governing equations are only partially known or completely unknown, and yet a predictive reduced-order model is required. In such cases, we must rely directly on data-driven methods to approximate the gSSMs. In particular, the SSMLearn algorithm of ref. ¹³ works on observations of trajectories to approximate the slow SSM $\mathcal{W}\left(E\right)$ locally.

A common use case is when a single scalar observation $y\left(t\right)\in \mathbb{R}$ is recorded. In that case, by the Takens embedding theorem⁴¹, a d − dimensional attracting SSM can be embedded generically in the space of delayed measurements,

10

$\mathbf{y}\left(t\right)={\left(y\left(t\right),y\left(t-\tau \right),\dots ,y\left(t-\left(p-1\right)\tau \right)\right)}^T,$for some time-lag τ as long as the number of delays p is more than twice the dimension of the SSM d. The slow SSM can then be parametrized as a graph over its tangent space at the fixed point from which it emanates. The reduced coordinates η on the SSM are obtained by projecting the delay-embedded measurements onto that d − dimensional tangent space, i.e., letting

11

$\mathit{\eta }={\mathbf{V}}^T\mathbf{y},$where the columns of V span the tangent space. The tangent space can be efficiently approximated using the leading principal components of the delay-embedded measurements after discarding initial transients^18,19.

The parametrization of the SSM is obtained by regression using the observed trajectories and their reduced coordinates. The most straightforward choice is a polynomial basis, which is justified by the existence of a locally convergent Taylor expansion of the SSM and by the universal approximation property of polynomials²³.

Motivated by the success of Padé approximants in enhancing the convergence properties of equation-driven models, we generalize the regression step of the data-driven SSM-reduction to rational approximants. Related approaches are rational interpolation and multi-point Padé approximation. While the former requires the approximant to fit the data exactly, the latter matches the asymptotic behavior of the function at multiple expansion points.

Rational function regression⁴² seeks a rational approximant of the form (9). In addition, we enforce that all components of the vector function share the same denominator and that the denominator is never zero on the training data. Having a common denominator for all components of the vector function makes it simpler to avoid spurious singularities. We elaborate on the steps of this regression task in “Rational function regression”.

Example 1: connecting orbit in Kolmogorov flow

Our first example is the 2D Kolmogorov flow, governed by the Navier-Stokes equations in a periodic domain subject to spatially periodic forcing. The flow domain is given by x, y ∈ [0, 2π] and the time evolution of the vorticity ω(x, y) = ∇ × u is governed by the non-dimensionalized equation

12

$\frac{\partial \omega }{\partial t}=-\left(\mathbf{u}\cdot \nabla \right)\omega +\frac{1}{\text{Re}}\mathrm{\Delta }\omega -L\cos Ly,$where L denotes the forcing wave number. This influences the size of the large-scale flow structures, the bifurcations observed in (12) and the properties of the turbulent dynamics at high Re⁴³. The laminar solution is a fixed point of (12) for all Reynolds numbers that can be written as

13

${\omega }_0\left(x,y\right)=-\frac{\text{Re}}{L}\cos Ly.$

Following³³, we set L = 4 and discretize the system using 576 = 24 × 24 Fourier modes. The numerical implementation is based on ref. ⁴⁴ and results in a system of ODEs in the form (1) for the Fourier amplitudes $\widehat{\omega }\left(k_x,k_y\right)$ with n = 576. Further details on the implementation can be found in the Supplementary Information.

The laminar flow ω₀ is already unstable for low Reynolds numbers, losing stability in a bifurcation around Re = 9.1. The state ω₀ is connected to two new stable fixed points ω_1,2 by its 1D unstable manifold. This unstable manifold coincides with the 1D slow SSMs of ω_1,2, forming two heteroclinic orbits.

We fix Re = 11 and consider the stable fixed point ω₁. We compute the parametrized slow SSM of ω₁ as i.e., $\widehat{\omega }=\mathbf{W}\left(\xi \right)$. This computation is carried out automatically by SSMTool, which returns the coefficients of the Taylor expansion of W(ξ) and the reduced dynamics on the SSM as

14

$\widehat{\omega }\left(k_x,k_y\right)=\mathbf{W}\left(\xi \right)={\sum }_{k=0}^N{\mathbf{W}}_k{\xi }^k,$

15

$\stackrel{°}{\xi }=R\left(\xi \right)={\sum }_{k=0}^N{R}_k{\xi }^k.$

We visualize the connecting orbit, i.e., the 1D slow SSM in a 3D slice of the 576-D phase space in Fig. 1a. The stable fixed points ω_1,2 are marked with black and blue dots, and the unstable fixed point ω₀ is red. The heteroclinic orbit, which is obtained by direct numerical integration of (12) (black curve), is approximated by the Taylor expansion of the SSM up to order-16. The domain of convergence is clearly limited, and it does not contain the unstable fixed point.

[See PDF for image]

Fig. 1

Model reduction for the Kolmogorov flow.

Heteroclinic orbits (black) connect ω_1,2 and ω₀. a Projection of the phase space onto three dominant Fourier modes (1, 4), (0, 4) and (2, 4). We show the slow SSM $\mathcal{W}\left(E\right)$ (black), which is tangent to the spectral subspace E (gray), its order-16 Taylor expansion (red) near the fixed point ω₁, and the order [5/5] gSSM approximation (orange). The curves are also projected onto the horizontal axes. The vorticity fields corresponding to the three fixed points are shown in the insets. b A trajectory on the heteroclinic orbit obtained by backwards integration and its SSM-reduced and gSSM-reduced counterparts. c SSM-reduced and gSSM-reduced dynamics.

To construct a globalized slow SSM, we compute the Padé approximant of each component (Fourier mode) of the vector W(ξ). Although, in principle, spurious poles could arise for each component separately, we find that the diagonal approximants are well-behaved. An alternative approach would be to construct the vector Padé approximants⁴⁵, that is, to find a common denominator for all components of W(ξ)⁴⁶.

The componentwise computed [5/5](ξ) Padé approximant approximates the heteroclinic orbit connecting ω₁ and ω₀ remarkably well, as shown by the orange curves in Fig. 1a. We also observe that the manifold $\mathcal{W}\left(E\right)$ can only be represented as a graph over E for this segment, since the derivative $\frac{\partial }{\partial \xi }\mathbf{W}\left(\xi \right)$ diverges at a fold point near ω₀. The parametrization, therefore, cannot be continued to capture ω₂. Note, however, that the Taylor approximation diverges well before encountering this unremovable singularity of the graph-style parametrization.

To verify the validity of the reduced-order models, predictions should be compared to trajectories of the full system. However, since the fixed point is stable, a nearby initial condition will leave its neighbourhood along the heteroclinic orbit only in backward time. Therefore, we integrate the initial condition ξ(0) = 10⁻⁵, close to the stable fixed point, backward under the SSM-reduced and the gSSM-reduced dynamics.

Figure 1b shows the SSM-reduced trajectory, which exhibits finite-time blowup and is only a reliable model near the stable fixed point. In contrast, the gSSM-reduced trajectory converges to ω₀ in backward time and hence captures the global dynamics accurately.

Figure 1c shows the reduced dynamics $\stackrel{°}{\xi }$ on the SSM and the gSSM. In addition to the trivial fixed point at ξ = 0, the gSSM-reduced model contains the unstable fixed point ω₀, given by the intersection with $\stackrel{°}{\xi }=0$ and hence provides a robust reduced model of the system. In contrast, based on the Taylor approximated SSM, the existence of an unstable fixed point cannot be concluded.

Example 2: Von Kármán beam

We now consider the model of a nonlinear von Kármán beam with clamped-free boundary conditions⁴⁷, shown in Fig. 2a. The beam has length L = 1 m, height h = 1 mm, width b = 0.1 m, Young’s modulus E = 70 GPA, viscous damping rate α = 10⁷Pa s and density 2700 kg/m³. We use a finite element discretization with 10 elements, resulting in 30 degrees of freedom, with a phase space of dimension n = 60. The discretization is implemented in the SSMTool³⁶ repository.

[See PDF for image]

Fig. 2

Model reduction of a von Kármán beam.

a von Kármán beam with clamped-free boundary conditions. b Trajectory of the unforced system (black) with its order-16 SSM approximation (red) and gSSM-approximation (orange) with a [5/5] and [5/4] Padé approximant. c, d Representation of the end point displacement q_end and the full-order trajectory shown in (b) with the SSM and gSSM-reduced trajectories. e, f Forced response defined as the maximal end point displacement due to a forcing amplitude ε = 0.5, 1.7, 2.8. e Shows the SSM-prediction and (f) shows the gSSM prediction. A supplementary animation showing the gSSM-prediction is also available.

We construct the slowest SSM of the undeformed configuration, which is tangent to the spectral subspace corresponding to the eigenvalues λ_1,2 = −0.0019 ± 5.1681i. Due to the Taylor approximation, the autonomous SSM can only capture the decaying oscillations of the beam up to amplitudes of around 5 mm, which is verified by backward integration of the SSM-reduced and the full models. This can be seen in Fig. 2b–c, which shows the autonomous trajectory exiting the domain of convergence of the Taylor series. This is evident for the high-order Taylor approximant, which exhibits finite-time blowup.

In contrast, the gSSM model, globalized using a [5/5] Padé approximant for the reduced dynamics and a [5/4] approximant for the parametrization, remains well-behaved for even larger amplitudes. Due to spurious singularities in the parametrization, we chose the [5/4] approximant instead of [5/5]. Since the convergence of Padé-approximants is only guaranteed in measure, singularities coinciding with zero sets of the denominator of (9) must be actively avoided before deploying the reduced-order model, as we demonstrate in the Supplementary Information.

The backbone curve, i.e., the instantaneous relationship between the normal form amplitude ρ and the frequency, is given by ω(ρ) in (7). The Taylor expansion of ω(ρ) and its [5/5] Padé approximant are given as

16

$\omega \left(\rho \right)=5.16+9.3\cdot {10}^4{\rho }^2-4.5\cdot 10^9{\rho }^4+4.2\cdot 10^{14}{\rho }^6+O\left({\rho }^8\right),\left[5/5\right]\left(\rho \right)=\frac{5.16+10^6{\rho }^2+3.7\cdot 10^{10}{\rho }^4}{1+1.8\cdot 10^5{\rho }^2+4.7\cdot {\rho }^4}.$

The Taylor coefficients in (16) are growing rapidly, and the alternating sign pattern of the coefficients suggests the convergence-limiting singularity is along the imaginary axis. This can be inferred using the method of refs. ^48,49, which we specialize to our examples in the Supplementary Information.

When an external forcing $\epsilon {\mathbf{f}}_{ext}\cos \left(\mathrm{\Omega }t\right)$ is applied to the beam with a frequency Ω almost in resonance with the slowest eigenfrequency Im(λ_1,2), we can use the unforced SSM to make predictions about the forced response. To leading order, the reduced dynamics of the forced system are simply a perturbed version of those of the unforced one (see “Spectral submanifolds”). Furthermore, the forced response corresponds to the periodic orbits of the reduced model. For 2D SSMs, these are directly given by the equation (26) in “Spectral submanifolds”.

We compare the SSM-reduced forced response to the forced response of the full system, obtained by direct numerical continuation using COCO⁵⁰. Predictions with a high-order Taylor approximation of the SSM are shown in Fig. 2e. They accurately capture the forced response for small forcing amplitudes up to ε = 0.5, as initially reported by ref. ¹¹. For higher amplitudes, the Taylor approximated reduced trajectory reaches its boundary of convergence, and the model breaks down. In contrast, the gSSM prediction using the [5/5] Padé approximant remains accurate for larger amplitudes as well. Note that the small errors in the peaks of the forced response curve are due to our initial assumption of a small forcing amplitude. Indeed, (26) technically holds only for small ε.

Example 3: chaotic von Kármán beam

In our next example, we construct a 2D mixed-mode SSM to characterize the chaotic behavior of a periodically forced buckled von Kármán beam, shown in Fig. 3. We adopt the same finite element code used in “Example 2: Von Kármán beam”, developed by ref. ¹¹, with pinned-pinned boundary conditions. To induce buckling, a compressive force is applied at the rightmost element, equal to 145% of the critical value f_crit = 1.5 kN. This gives rise to a pair of stable fixed points. The same system, much closer to the critical point, was also analyzed using a data-driven model by ref. ¹⁷.

[See PDF for image]

Fig. 3

Reduced model of the buckled beam.

a Sketch of the beam in the buckled configuration, with no external forcing. b Trajectories in the unstable manifold of the unstable fixed point (black). Their order-18 SSM-reduced (red) and order [6/6] gSSM-reduced (orange) approximations are also shown. c, d The SSM and gSSM in the physical space with the direction field of the reduced dynamics indicated on the surface of the manifold. The predicted trajectories connecting the unstable fixed point to the stable fixed points are shown in black.

The spectral subspace of the linear part associated with the buckling instability is two-dimensional and has the corresponding eigenvalues λ₁ = 23.48 and λ₂ = −23.48. The other modes are all stable and correspond to oscillatory dynamics. This eigenvalue configuration indicates a considerable departure from criticality and center manifold theory.

We use SSMTool to find an order-18 approximation of the 2D SSM of the unforced beam and its reduced dynamics. Since the SSM is tangent to a spectral subspace with real eigenvalues, we denote the reduced coordinates as $\mathit{\eta }=\left({\eta }_1,{\eta }_2\right)\in {\mathbb{R}}^2$. The reduced dynamics are given by

17

${\stackrel{°}{\eta }}_1=23.48{\eta }_1-2800{\eta }_1^3-1760{\eta }_2^3-7297{\eta }_1^2{\eta }_2-6268{\eta }_1{\eta }_2^2+O\left({\mathit{\mid \eta \mid }}^4\right),{\stackrel{°}{\eta }}_2=-23.52{\eta }_2+2790{\eta }_1^3+1760{\eta }_2^3+7297{\eta }_1^2{\eta }_2+6268{\eta }_1{\eta }_2^2+O\left({\mathit{\mid \eta \mid }}^4\right).$

Analysis of the reduced dynamics shows, however, that the fixed points born due to the buckling instability lie outside the domain of convergence. We show in Fig. 3b the time series of a trajectory initialized on the unstable manifold of the fixed point. The SSM-reduced trajectory based on Taylor expansion blows up once the reduced trajectory exits the domain of convergence. In contrast, a comparable, [6/6] Padé approximant and the gSSM-reduced trajectory correctly capture the convergence to both of the buckled states.

This is even more apparent in Fig. 3c, d, showing the image of the parametrization and the direction field of the reduced dynamics. The SSM-reduced model in Fig. 3b predicts diverging, unphysical displacements for the mid-point of the beam. The gSSM-model in Fig. 3c correctly identifies all fixed points and the orbits connecting them. We have, therefore, extended the SSM-reduced model obtained from local information around the unstable fixed point to a globally valid one.

We now extend the autonomous gSSM model to account for periodic forcing. A periodic force on the middle node of the beam, as shown in Fig. 4a, can make the dynamics of the beam chaotic. This has been observed and reported in the data-driven model of ref. ¹⁷, who used simulation data of the full system. We now characterize the chaotic behavior without relying on any full-order simulations.

[See PDF for image]

Fig. 4

Chaotic response of the buckled beam.

a Buckled von Kármán beam with periodic external forcing. b Time series of the reduced coordinate η₁ on a chaotic trajectory of the full system (black). Also shown are the SSM-reduced forced model, which diverges immediately (red), and the gSSM-reduced trajectory (orange). c Sampling of the Poincaré map of the gSSM-reduced model (orange) and the true system (black). d, e The autonomous SSM and gSSM with the chaotic trajectory of the full model. A supplementary animation comparing the SSM- and gSSM-predictions is also available.

To compare the full-order and the reduced-order dynamics, we take one of the buckled fixed points as an initial condition and simulate the beam under the influence of periodic forcing acting on the mid-point with Ω = 25.3 rad/s and ∣εf_ext∣ = 21.1 N. Since this initial condition is known to be on the SSM, we run the same trajectory with the reduced dynamics after adding the leading-order contribution of the forcing as in (23).

The full-order trajectory is shown in Fig. 4b. As expected from the previous analysis and Fig. 3, the forced SSM-reduced model blows up within a fraction of a second. In contrast, the globally valid gSSM model exhibits sustained chaotic behavior, closely matching the full-order trajectory for short times. Figure 4d–e shows that the chaotic attractor appearing as a result of the periodic forcing extends way outside the domain of convergence of the Taylor series of the SSM.

We also construct the Poincaré-map of the full model and the gSSM model by sampling the trajectories at multiples of the driving period $T=\frac{2\pi }{\mathrm{\Omega }}$. Because the gSSM-model is a simple 2D ODE, we can sample the Poincaré-map with a fine resolution to obtain the structure in the reduced phase space shown in Fig. 4c. Overlaying the Poincaré-map obtained from the full system, we see a close correspondence with the predicted attractor.

In addition, we estimate the leading Lyapunov exponents based on the exponential rate of divergence of initially close trajectories⁵¹ as λ_gSSM = (3.0 ± 0.02)1/s and λ_full = (3.1 ± 0.05)1/s. The reduced model is in close agreement with the full model, even though the forced dynamics were approximated by simply projecting the forcing term onto the tangent space of the SSM according to (23). Incorporating higher-order corrections of the forced dynamics, as in (22), further improves the model, leading to more accurate short-time predictions. We present these comparisons, along with additional properties of the chaotic gSSM model, in the Supplementary Information.

Example 4: data-driven model of an inverted flag experiment

We now consider, as a data-driven example, the dynamics of an inverted flag, which is a flexible elastic sheet in the counterflow of a water tunnel. This configuration has generated recent interest due to its applications in energy harvesting⁵² and vegetation⁵³. From a modeling perspective, the inverted flag is a complicated fluid-structure interaction problem that benefits from reduced-order modeling. Here, we rely on experimental data obtained by ref. ¹⁸, whose experimental configuration is shown in Fig. 5a–b. The elastic sheet is mounted in a water tunnel and is recorded from below with a video camera. Using classical image processing tools, the deflection of the tip of the flag y(t) is recorded during the experiment.

[See PDF for image]

Fig. 5

Reduced model of the inverted flag experiment.

a Experimental setup of the inverted flag and snapshot of the experiment (b), courtesy of Giovanni Berti. The geometric parameters are H = 150 mm, L = 150 mm, U = 1m/s. c Phase portrait of the gSSM-reduced dynamics of the inverted flag obtained from the [5/5](η) approximation. The unstable fixed points are marked with colored dots. Blue curves denote the stable manifold of the undeflected state, which connects to the two coexisting deflected fixed points. The red curve denotes the unstable manifold of the saddle, which wraps around the stable limit cycle. d, e Predictions of the reduced model on test trajectories. The black curve is the true trajectory, and the dotted orange curve is the gSSM prediction.

The main parameters of the system are the bending stiffness K_B and the Reynolds number of the incoming flow, governed by the mean velocity U. By tracking the displacement of the tip of the flag, a data-driven SSM-reduced model was obtained¹⁸ for both the large-amplitude periodic flapping regime and the chaotic flapping regime. In these regimes, the undeflected state of the flag is an unstable fixed point, which has a low-dimensional attracting mixed-mode SSM.

We focus here on the large-amplitude flapping regime with K_B = 0.21 and Re = 6 × 10⁴. The slow SSM is tangent to a spectral subspace associated with an unstable real eigenvalue and a stable real eigenvalue, and is two-dimensional. In addition to the saddle-type fixed point corresponding to the undeflected flag, the slow SSM contains two additional fixed points, which correspond to the deflected, but still stationary flag. The periodic flapping is a stable limit cycle within this slow SSM.

We now construct a data-driven gSSM-reduced model based on the tip-deflection data of ref. ¹⁸. A total of 16 trajectories are used for training. To reconstruct the slow SSM, we embed the tip displacement data y(t) trajectories using p = 25 time delays to form the observable $\mathbf{y}\in {\mathbb{R}}^{25}$. As in ref. ¹⁸, we approximate the tangent space of the SSM at the fixed point using the two leading principal components of the delay-embedded trajectories. Since we use a moderate number of time delays and a short delay time, a leading-order (linear) approximation for the SSM suffices, i.e., we let

18

$\left(\begin{array}{c}\hfill {\eta }_1\hfill \\ \hfill {\eta }_2\hfill \end{array}\right)=\mathit{\eta }={\mathbf{V}}^T\mathbf{y},\mathbf{y}=\mathbf{V}\mathit{\eta }.$

The reduced dynamics are now approximated using rational functions, as detailed in “Rational function regression”, by solving the minimization problem (44). We find that a [5/5] approximant gives the optimal reconstruction error, as computed on a validation trajectory, which was not used in the training. The resulting reduced vector field and the predictions on test trajectories are shown in Fig. 5c–e. A diagonal approximant also ensures that the reduced dynamics remain well-behaved outside the domain of the training data. As opposed to classical polynomial regression, the values of the approximants remain bounded or only grow mildly for large ∣η∣.

To obtain a comparable test accuracy to that shown in Fig. 5d–e, previously, an order-11 polynomial approximation was used for the reduced dynamics. This required determining a total of 154 coefficients R_k. The data-driven rational approximation, on the other hand, requires only 60 coefficients, which is a significant reduction, for the cost of a slightly increased computational burden. In addition, the rational approximants extrapolate to larger domains in a more controlled way.

As we show in the Supplementary Information, preventing the rational functions from becoming singular is essential for an accurate approximation. Other available methods, such as the rational function extension of the Sparse Identification of Nonlinear Dynamics (SINDy) algorithm^54,55 enforce no such constraints.

Discussion

We have presented a method to extend the range of validity of invariant manifold-based reduced-order models by applying Padé approximation, a classic analytic continuation technique. The Taylor coefficients of the parametrization of the SSM and the reduced dynamics were obtained with the robust numerical routines of SSMTool. We then extend their range of validity globally to obtain the gSSM-reduced model, by applying Padé approximation to these mappings that agree with the local Taylor expansions up to a prescribed order.

We have demonstrated the method on high-dimensional examples of dynamical systems exhibiting global nonlinear behaviour, such as transitions between steady states, large-amplitude oscillations, or even chaotic behaviour. In all cases, gSSM models obtained with Padé approximation significantly extend the domain of validity of the reduced model, reaching well beyond the domain of convergence of the classical Taylor series. We have also shown that the sign pattern of the univariate Taylor coefficients can be used to infer the location of the singularity.

To apply our approach, unexpected singularities in the approximants must be checked, and the optimal approximant of the Padé table must be determined according to the steps we have outlined in “Padé approximants”. Luckily, the numerical effort required to compute Padé approximants is negligible compared to the computation of the original Taylor coefficients of the SSM and the reduced dynamics. Therefore, even if no singularity-free Padé approximants can be found up to a given order, the benefits of a gSSM model far outweigh the costs of deriving it.

We have shown that a data-driven analog of Padé approximation, rational function regression, can be directly applied to experimental data. Although rational function regression requires more computational power than the standard polynomial approximation implemented in SSMLearn, it can produce more accurate reduced models with fewer unknown coefficients, resulting in ultimately simpler models.

We have demonstrated the method on the inverted flag example of ref. ¹⁸, where the rational approximants of the SSM and the reduced dynamics were able to predict the coexisting deflected and undeflected fixed points and the stable limit cycle. As a main benefit, we note that the data-driven approach does not suffer from the same singularity issues as the equation-driven one.

As we have demonstrated, Padé approximants and their data-driven extensions are particularly well-suited for approximating the global reduced dynamics on invariant manifolds. In addition to being able to represent more complex functions, they retain physical interpretability. In contrast to modern machine learning methods, the coefficients of the gSSM-reduced dynamics can be directly related to the underlying physics of the system. Specifically, the coefficients of the rational function can be interpreted as nonlinear frequencies or damping rates, as we have shown in “Example 2: Von Kármán beam”.

Finally, we note that although we have focused here on SSM-reduced models, the globalization method applies to other perturbative methods for dynamical systems, such as Poincaré-Lindstedt series⁵⁶, geometric singular perturbation theory^5,6, or model reduction based on local linearization results⁵⁷.

Methods

Spectral submanifolds

Let us assume that x = 0 is a hyperbolic fixed point of the system (1) with ε = 0 and ${\lambda }_1,\dots ,{\lambda }_n\in \mathbb{C}$ are the eigenvalues of A with the corresponding eigenvectors denoted ${\mathbf{e}}_1,\dots ,{\mathbf{e}}_n\in {\mathbb{C}}^n$. We assume that a spectral gap condition holds for some d < n, i.e.,

19

$\text{Re}{\lambda }_n\le \dots \le \text{Re}{\lambda }_{d+1}<\text{Re}{\lambda }_d\le \dots \le \text{Re}{\lambda }_1.$

Let us denote the d – dimensional slow spectral subspace of A as E, which is defined as the span of the real and imaginary parts of the eigenvectors corresponding to the d eigenvalues closest to zero. If the nonresonance condition

20

$\begin{array}{c}{\lambda }_k\ne {\sum }_{j=1}^n{m}_j{\lambda }_j,k=1,\dots ,n,\hfill \\ m_j\in \mathbb{N},{\sum }_{j=1}^n{m}_j>1,\hfill \end{array}$holds, a family of SSMs exists tangent to E, as discussed in refs. ^9,10. The primary SSM is the unique, smoothest member of the family of d − dimensional invariant manifolds tangent to E at the origin and is denoted as $\mathcal{W}\left(E\right)$^9,10.

SSMs also exist for the non-autonomous system, i.e., for ε > 0. In that case, SSMs are slow invariant manifolds attached to an anchor trajectory³⁴. For simplicity, we focus here on the case of periodic forcing with a single harmonic, where the anchor trajectory is a periodic orbit, and the original formulation of ref. ¹¹ applies. Denoting the parametrization of the slow SSM as W_ε(p, Φ), where Φ is the phase of the periodic forcing, the invariance equation reads as

21

$\begin{array}{c}\mathbf{A}{\mathbf{W}}_{\epsilon }\left(\mathbf{p},\mathrm{\Phi }\right)+\mathbf{f}\left({\mathbf{W}}_{\epsilon }\left(\mathbf{p},\mathrm{\Phi }\right)\right)+\epsilon {\mathbf{f}}_{\mathrm{ext}}\cos \mathrm{\Phi }\hfill \\ =D{\mathbf{W}}_{\epsilon }\left(\mathbf{p},\mathrm{\Phi }\right)\stackrel{°}{\mathbf{p}}+D_{\mathrm{\Phi }}{\mathbf{W}}_{\epsilon }\stackrel{°}{\mathrm{\Phi }},\hfill \end{array}$where $\stackrel{°}{\mathrm{\Phi }}=\mathrm{\Omega }$ is the forcing frequency. The reduced dynamics is denoted as $\stackrel{°}{\mathbf{p}}={\mathbf{R}}_{\epsilon }\left(\mathbf{p},\mathrm{\Phi }\right)$.

Due to the guaranteed smoothness of primary SSMs, the invariance equation can be solved using a Taylor expansion in p and a Fourier expansion in Φ. The coefficients are then obtained by imposing the invariance equation order-by-order in the reduced coordinates. As shown by refs. ^{58, 59–60}, to first order in ε and to order-$\widehat{N}$ in the reduced coordinates, the resulting expression is

22

${\mathbf{R}}_{\epsilon }\left(\mathbf{p},\mathrm{\Phi }\right)=\mathbf{R}\left(\mathbf{p}\right)+\epsilon {\sum }_{\mid \mathbf{k}\mid =0}^{\widehat{N}}{\mathbf{S}}_{\mathbf{k}}\left(\mathrm{\Phi }\right){\mathbf{p}}^{\mathbf{k}}+O\left({\epsilon }^2\right),$where the coefficients S_k(Φ) can also be Fourier-expanded. Accounting for the phase-dependence of the nonautonomous coefficients in the reduced dynamics increases the accuracy of the reduced models for higher forcing amplitudes, at the cost of computing the coefficients S_k(Φ) for each forcing frequency of interest. These computations are already implemented in SSMTool⁶⁰.

However, the phase-dependence of the parametrization and the reduced dynamics are often only small effects. It is, therefore, common to keep only the leading-order phase-dependent term in the reduced dynamics, resulting in

23

${\mathbf{R}}_{\epsilon }\left(\mathbf{p}\right)=\mathbf{R}\left(\mathbf{p}\right)+\epsilon {\mathbf{V}}^{*}{\mathbf{f}}_{ext}\cos \left(\mathrm{\Omega }t\right),$where the operator V^* projects onto the tangent space and can be computed using the left eigenvectors of A. Specifically, for a non-resonant two-dimensional SSM, we have

24

$\stackrel{°}{\rho }=\kappa \left(\rho \right)\rho +\epsilon f\sin \psi$

25

$\stackrel{°}{\psi }=\omega \left(\rho \right)-\mathrm{\Omega }+\epsilon f\frac{1}{\rho }\cos \psi ,$where we have introduced the phase lag ψ = θ − Ωt, and the forcing amplitude f is the projection of the external forcing amplitude¹¹.

The forced response of the system is then obtained by seeking fixed points of the reduced dynamics, which are given by the solutions of the equation $\stackrel{°}{\rho }=\stackrel{°}{\psi }=0$. The amplitude ρ* of the response satisfies the implicit equation¹¹

26

${\left(\mathrm{\Omega }-\omega \left({\rho }_{*}\right)\right)}^2-\frac{{\epsilon }^2{f}^2}{{\rho }_{*}^2}+\kappa {\left({\rho }_{*}\right)}^2=0.$

Padé approximants

We now seek to improve the convergence properties of Taylor series approximations by summing the series outside the domain of convergence. One of the most popular methods of analytic continuation is Padé approximation, which is often used to sum divergent perturbative series. Padé approximation is primarily carried out on functions of a single complex variable, and hence it is directly applicable to 1D SSMs and to the functions κ(ρ) and ω(ρ).

Consider the function

27

$f\left(z\right):\mathbb{C}\to \mathbb{C},$which could represent the parametrization of an invariant manifold, the reduced dynamics or κ(ρ) and ω(ρ) in (26). Let us denote its Taylor series representation around z = 0 as

28

$f\left(z\right)={\sum }_{k=0}^{\infty }{c}_k{z}^k.$

The Padé approximant of type (N, M) is defined as the rational function

29

$\left[N/M\right]\left(z\right)=\frac{{\sum }_{k=0}^N{a}_k{z}^k}{{\sum }_{k=0}^M{b}_k{z}^k},$where a₀ = f(0) and b₀ can be chosen as b₀ = 1 without loss of generality. The Padé approximant is the best rational approximation of f around 0²⁴ in the sense that its Taylor expansion around 0 is the same as that of f up to order N + M, i.e.,

30

$f\left(z\right){\sum }_{k=0}^M{b}_k{z}^k-{\sum }_{k=0}^N{a}_k{z}^k=O\left(z^{N+M+1}\right).$

Therefore, the coefficients b_k and a_k can be determined by requiring

31

${\sum }_{k=0}^{\infty }{c}_k{z}^k{\sum }_{m=0}^M{b}_m{z}^m={\sum }_{k=0}^N{a}_k{z}^k$for orders up to N + M in z. This results in the linear equations for the coefficients a_k, b_k,

32

$a_k={\sum }_{m=0}^k{c}_{k-m}{b}_m,k=0,1,\dots ,N+M.$

We follow the robust approach of ref. ⁶¹, who first solve the homogeneous equations

33

${\sum }_{m=0}^k{c}_{k-m}{b}_m=0,k=N+1,\dots ,N+M,$with the convention of b_j = 0 for j < 0 or j > M. Using the singular value decomposition (SVD) of the Toeplitz matrix with elements [c_k−m], the coefficients b_k can be computed. The remaining unknown coefficients a_k are then given as

34

$a_k={\sum }_{m=0}^k{c}_{k-m}{b}_m,k=0,1,\dots ,N.$This method is robust against numerical errors in the Taylor coefficients. In addition, many of the spurious poles of the approximant are removed, although not all of them⁶².

To apply Padé approximation to higher-dimensional SSMs, we need to generalize the method to multivariate power series. The multivariate generalization is not as straightforward as the univariate case, but multiple definitions exist. The most commonly used are Chisholm approximants⁶³ and the homogeneous approximants^24,40,64. Let us consider the multivariate function

35

$f:{\mathbb{C}}^d\to \mathbb{C}.$given as a convergent Taylor series

36

$f\left(\mathbf{z}\right)={\sum }_{\mid \mathbf{k}\mid =0}^{\infty }{c}_{\mathbf{k}}{\mathbf{z}}^{\mathbf{k}}.$

We adopt the homogeneous approximants (9), as defined originally by refs. ^40,64 and require that the Taylor expansion of [N/M](z) around z = 0 coincides with that of f(z) up to order N + M, as

37

$f\left(\mathbf{z}\right){\sum }_{\mid \mathbf{k}\mid =0}^M{b}_{\mathbf{k}}{\mathbf{z}}^{\mathbf{k}}-{\sum }_{\mid \mathbf{k}\mid =0}^N{a}_{\mathbf{k}}{\mathbf{z}}^{\mathbf{k}}=O\left({\mathbf{z}}^{N+M+1}\right)$

38

${\sum }_{\mid \mathit{\ell }\mid =0}^{\mid \mathbf{k}\mid }{c}_{\mathbf{k}-\mathit{\ell }}{b}_{\mathit{\ell }}=a_{\mathbf{k}}\text{for}\mid \mathbf{k}\mid =0,1,\dots ,N+M.$

Since two-dimensional, oscillatory SSMs are the relevant objects for most systems, we consider the bivariate case as an illustration, i.e., with $\mathbf{z}={\left(z_1,z_2\right)}^T$, where the coefficients are indexed on a lattice ${\mathcal{L}}_{N_1,N_2}$ defined as ${\mathcal{L}}_{N_1,N_2}=\left\{\left(i,j\right)\in {\mathbb{N}}^2:N_1\le i+j\le N_2\right\}\subset {\mathbb{N}}^2$.

The conditions (38) become

39

${\sum }_{k=0}^{\alpha }{\sum }_{\ell =0}^k{c}_{\beta -\ell ,\alpha -k-\beta +\ell }{b}_{\ell ,k-\ell }=a_{\beta ,\alpha -\beta },$for $\left(\alpha ,\beta \right)\in {\mathcal{L}}_{0,N}$ and

40

${\sum }_{k=0}^{\alpha }{\sum }_{\ell =0}^k{c}_{\beta -\ell ,\alpha -k-\beta +\ell }{b}_{\ell ,k-\ell }=0,$for $\left(\alpha ,\beta \right)\in {\mathcal{L}}_{N,N+M}$. Note, however, that the total number of unknowns and equations in (39) and (40) are not equal. Therefore, we seek a least-squares solution to (40)⁶⁵.

To summarize, gSSM-reduction consists of the following steps to compute the [N/M] Padé approximants of the parametrization and the reduced dynamics of a slow SSM.

1. Compute the Taylor series expansion of the SSM up to order N + M using either the normal form style parametrization or the graph-style parametrization. This returns W(p), ω(ρ) and κ(ρ).

2. Solve the homogeneous linear equations (40) and (33) for the coefficients of the denominators of the parametrization and of ω and κ using the robust SVD-based method⁶¹.

3. Check the zero sets of the denominator functions. If they contain points in the region of interest near the origin, adjust the orders N and M. We found that diagonal and near diagonal [N ± 2/N ± 2] approximants tend to work the best.

4. Compute the coefficients of the numerators of the parametrization and of ω and κ by evaluating (39) and (34).

Although Padé approximation is a well-researched topic, convergence results are generally limited. For meromorphic functions f(z) with a finite number of poles at z₁, …, z_k, the theorem of de Montessus de Ballore⁶⁶ guarantees the convergence of the approximants

41

$\left[M/k\right]\left(z\right)\text{as}M\to \infty$globally. A stronger result is available for Stieltjes functions, i.e., for functions of the form

42

$f\left(z\right)={\int }_0^{\infty }\frac{d\mu \left(t\right)}{1+zt},$for some positive measure μ. In this case, the Padé approximants converge to f(z) for all z outside the negative real axis²⁴.

In the Supplementary Information, we discuss the Stiltjes-type center manifold of Euler’s system^{9,67, 68–69}, which, although non-analytic, can be described globally using Padé approximants.

Rational function regression

Let us denote the function to be approximated as $\mathbf{f}:{\mathbb{R}}^d\to {\mathbb{R}}^{\ell }$. This could represent f = W with ℓ = n or f = R with ℓ = d. We assume that the value of f, denoted ζ_i = f(η_i), is known at points η_i, …, η_K in the domain of interest. We then approximate f as

43

$\mathbf{f}\left(\mathit{\eta }\right)\approx \left[N/M\right]\left(\mathit{\eta }\right)=\frac{{\sum }_{\mid \mathbf{k}\mid =0}^N{\mathbf{a}}_{\mathbf{k}}{\mathit{\eta }}^{\mathbf{k}}}{{\sum }_{\mid \mathbf{k}\mid =0}^M{b}_{\mathbf{k}}{\mathit{\eta }}^{\mathbf{k}}},$where we have chosen a common denominator with coefficients b_k for all components of f, similarly to vector Padé approximants^45,46. To avoid introducing singularities for the approximants, we require that the denominator is non-zero at all points η_i. This is equivalent to requiring that the denominator is strictly positive.

The coefficients are determined by minimizing the error

44

${\mathcal{E}}_r={\sum }_{i=1}^K{∣{\mathit{\zeta }}_i-\frac{{\sum }_{\mid \mathbf{k}\mid =0}^N{\mathbf{a}}_{\mathbf{k}}{\mathit{\eta }}_i^{\mathbf{k}}}{{\sum }_{\mid \mathbf{k}\mid =0}^M{b}_{\mathbf{k}}{\mathit{\eta }}_i^{\mathbf{k}}}∣}^2,\mathrm{such}\mathrm{that}$

45

${\sum }_{\mid \mathbf{k}\mid =0}^M{b}_{\mathbf{k}}{\mathit{\eta }}_i^{\mathbf{k}}\ge \delta \text{for}i=1,\dots ,K,$for some small δ > 0. The constrained minimization problem is solved by gradient-based optimization methods⁷⁰. The initial guess for the coefficients is obtained by solving the linearized problem, minimizing

46

${\sum }_{i=1}^K{∣\left({\sum }_{\mid \mathbf{k}\mid =0}^M{b}_{\mathbf{k}}{\mathit{\eta }}_i^{\mathbf{k}}\right){\mathit{\zeta }}_i-{\sum }_{\mid \mathbf{k}\mid =0}^N{\mathbf{a}}_{\mathbf{k}}{\mathit{\eta }}_i^{\mathbf{k}}∣}^2,$subject to the same positivity constraint. This is a linear least-squares problem, which can be solved efficiently using the method of ref. ⁴², based on robust Padé approximation⁶¹.

Acknowledgements

We are grateful to Mohammad Farazmand for sharing his code to simulate the Kolmogorov flow.

Author contributions

B.K. and G.H. designed the research, analyzed the examples, and edited the manuscript. B.K. carried out the formal analysis, developed the numerical algorithm, and wrote the manuscript.

Peer review

Peer review information

Nature Communications thanks the anonymous reviewers for their contribution to the peer review of this work. A peer review file is available.

Funding

Open access funding provided by Swiss Federal Institute of Technology Zurich.

Data availability

The data discussed in the manuscript is available publicly in the globalized-SSM repository at https://github.com/haller-group/globalized-SSM.

Code availability

The code for the numerical implementation of the methods is available publicly in the globalized-SSM repository at https://github.com/haller-group/globalized-SSM.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-61252-9.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Baldi, P; Hornik, K. Neural networks and principal component analysis: learning from examples without local minima. Neural Netw.; 1989; 2, pp. 53-58. [DOI: https://dx.doi.org/10.1016/0893-6080(89)90014-2]

2. Gorban, AN; Tyukin, IY. Blessing of dimensionality: mathematical foundations of the statistical physics of data. Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.; 2018; 376, 20170237.3797498 [DOI: https://dx.doi.org/10.1098/rsta.2017.0237]

3. Floryan, D; Graham, MD. Data-driven discovery of intrinsic dynamics. Nat. Mach. Intell.; 2022; 4, pp. 1113-1120. [DOI: https://dx.doi.org/10.1038/s42256-022-00575-4]

4. Guckenheimer, J. & Holmes, P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, vol. 42 (Springer New York, New York, NY, 1983).

5. Fenichel, N. Persistence and smoothness of invariant manifolds for flows. Indiana Univ. Math. J.; 1971; 21, pp. 193-226.287106 [DOI: https://dx.doi.org/10.1512/iumj.1972.21.21017]

6. Kuehn, C. Multiple Time Scale Dynamics, vol. 191 (Springer International Publishing, Cham, 2015).

7. Foias, C; Sell, GR; Temam, R. Inertial manifolds for nonlinear evolutionary equations. J. Differ. Equ.; 1988; 73, pp. 309-353.943945 [DOI: https://dx.doi.org/10.1016/0022-0396(88)90110-6]

8. Koronaki, ED; Evangelou, N; Martin-Linares, CP; Titi, ES; Kevrekidis, IG. Nonlinear dimensionality reduction then and now: aims for dissipative PDEs in the ML era. J. Comput. Phys.; 2024; 506, 112910.4720994 [DOI: https://dx.doi.org/10.1016/j.jcp.2024.112910]

9. Haller, G; Ponsioen, S. Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction. Nonlinear Dyn.; 2016; 86, pp. 1493-1534.3562433 [DOI: https://dx.doi.org/10.1007/s11071-016-2974-z]

10. Haller, G; Kaszás, B; Liu, A; Axås, J. Nonlinear model reduction to fractional and mixed-mode spectral submanifolds. Chaos: Interdiscip. J. Nonlinear Sci.; 2023; 33, 063138.4599208 [DOI: https://dx.doi.org/10.1063/5.0143936]

11. Jain, S; Haller, G. How to compute invariant manifolds and their reduced dynamics in high-dimensional finite element models. Nonlinear Dyn.; 2022; 107, pp. 1417-1450. [DOI: https://dx.doi.org/10.1007/s11071-021-06957-4]

12. Bettini, L; Cenedese, M; Haller, G. Model reduction to spectral submanifolds in piecewise smooth dynamical systems. Int. J. Non-Linear Mech.; 2024; 163, 104753. [DOI: https://dx.doi.org/10.1016/j.ijnonlinmec.2024.104753]

13. Cenedese, M; Axås, J; Bäuerlein, B; Avila, K; Haller, G. Data-driven modeling and prediction of non-linearizable dynamics via spectral submanifolds. Nat. Commun.; 2022; 13, 1:CAS:528:DC%2BB38Xht1CmtLbJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35169152][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8847615][DOI: https://dx.doi.org/10.1038/s41467-022-28518-y] 872.

14. Cenedese, M; Axås, J; Yang, H; Eriten, M; Haller, G. Data-driven nonlinear model reduction to spectral submanifolds in mechanical systems. Philos. Trans. R. Soc. A; 2022; 380, 20210194.44604891:STN:280:DC%2BB2MjmvFGgug%3D%3D [DOI: https://dx.doi.org/10.1098/rsta.2021.0194]

15. Kaszás, B; Cenedese, M; Haller, G. Dynamics-based machine learning of transitions in Couette flow. Phys. Rev. Fluids; 2022; 7, L082402. [DOI: https://dx.doi.org/10.1103/PhysRevFluids.7.L082402]

16. Kaszás, B; Haller, G. Capturing the edge of chaos as a spectral submanifold in pipe flows. J. Fluid Mech.; 2024; 979, A48.4695480 [DOI: https://dx.doi.org/10.1017/jfm.2023.956]

17. Liu, A; Axås, J; Haller, G. Data-driven modeling and forecasting of chaotic dynamics on inertial manifolds constructed as spectral submanifolds. Chaos: Interdiscip. J. Nonlinear Sci.; 2024; 34, 033140.4722159 [DOI: https://dx.doi.org/10.1063/5.0179741]

18. Xu, Z et al. Data-driven modelling of the regular and chaotic dynamics of an inverted flag from experiments. J. Fluid Mech.; 2024; 987, R7.47479631:CAS:528:DC%2BB2cXis1Cms7rI [DOI: https://dx.doi.org/10.1017/jfm.2024.411]

19. Axås, J; Cenedese, M; Haller, G. Fast data-driven model reduction for nonlinear dynamical systems. Nonlinear Dyn.; 2023; 111, pp. 7941-7957. [DOI: https://dx.doi.org/10.1007/s11071-022-08014-0]

20. Alora, JI; Cenedese, M; Schmerling, E; Haller, G; Pavone, M. Practical deployment of spectral submanifold reduction for optimal control of high-dimensional systems. IFAC-PapersOnLine; 2023; 56, pp. 4074-4081. [DOI: https://dx.doi.org/10.1016/j.ifacol.2023.10.1734]

21. Cabre, X; Fontich, E; de la Llave, R. The parameterization method for invariant manifolds I: manifolds associated to non-resonant subspaces. Indiana Univ. Math. J.; 2003; 52, pp. 283-328.1976079 [DOI: https://dx.doi.org/10.1512/iumj.2003.52.2245]

22. Haro, À., Canadell, M., Figueras, J. L., Luque, A. & Mondelo, J. M. The Parameterization Method for Invariant Manifolds: From Rigorous Results to Effective Computations, vol. 195 (Springer International Publishing, Cham, 2016).

23. Rudin, W. Real and complex analysis, 3rd edn. (McGraw-Hill, New York, 1987).

24. Baker, G. A. & Graves-Morris, P. Padé Approximants, 2nd edn. volume 59 (Cambridge University Press, Cambridge, 1996). https://doi.org/10.1017/CBO9780511530074.

25. Bender, C. M. & Orszag, S. A. Advanced Mathematical Methods for Scientists and Engineers I (Springer New York, New York, NY, 1999). https://doi.org/10.1007/978-1-4757-3069-2.

26. Alho, A; Uggla, C. Global dynamics and inflationary center manifold and slow-roll approximants. J. Math. Phys.; 2015; 56, 012502.3390832 [DOI: https://dx.doi.org/10.1063/1.4906081]

27. Zinn-Justin, J. Convergence of Padé approximants in the general case. Rocky Mt. J. Math. 4, https://doi.org/10.1216/RMJ-1974-4-2-325 (1974).

28. Razafindralandy, D; Hamdouni, A; Allery, C. Numerical divergent series resummation in fluid flow simulation. Eur. J. Comput. Mech.; 2008; 17, pp. 431-451. [DOI: https://dx.doi.org/10.3166/remn.17.431-451]

29. Andrianov, I; Shatrov, A. Padé approximants, their properties, and applications to hydrodynamic problems. Symmetry; 2021; 13, 1869. [DOI: https://dx.doi.org/10.3390/sym13101869]

30. Bonizzoni, F; Nobile, F; Perugia, I; Pradovera, D. Least-Squares Padé approximation of parametric and stochastic Helmholtz maps. Adv. Comput. Math.; 2020; 46, 46. [DOI: https://dx.doi.org/10.1007/s10444-020-09749-3]

31. Sinou, JJ; Thouverez, F; Jézéquel, L. Methods to reduce non-linear mechanical systems for instability computation. Arch. Comput. Methods Eng.; 2004; 11, pp. 257-344.2109050 [DOI: https://dx.doi.org/10.1007/BF02736228]

32. Sinou, JJ; Thouverez, F; Jezequel, L. Non-linear stability analysis of a complex rotor/stator contact system. J. Sound Vib.; 2004; 278, pp. 1095-1129. [DOI: https://dx.doi.org/10.1016/j.jsv.2003.10.048]

33. Chandler, GJ; Kerswell, RR. Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow. J. Fluid Mech.; 2013; 722, pp. 554-595.3137827 [DOI: https://dx.doi.org/10.1017/jfm.2013.122]

34. Haller, G; Kaundinya, RS. Nonlinear model reduction to temporally aperiodic spectral submanifolds. Chaos: Interdiscip. J. Nonlinear Sci.; 2024; 34, 043152.4736865 [DOI: https://dx.doi.org/10.1063/5.0187080]

35. Stoychev, AK; Römer, UJ. Failing parametrizations: what can go wrong when approximating spectral submanifolds. Nonlinear Dyn.; 2023; 111, pp. 5963-6000. [DOI: https://dx.doi.org/10.1007/s11071-022-08154-3]

36. Jain, S., Li, M., Thurnher, T. & Haller, G. SSMTool 2.6: computation of invariant manifolds in high-dimensional mechanics problems https://doi.org/10.5281/zenodo.13909795 (2024).

37. Scheidemann, V. Introduction to Complex Analysis In Several Variables (Springer Nature Switzerland, Cham, 2023).

38. Nuttall, J. The convergence of Padé approximants of meromorphic functions. J. Math. Anal. Appl.; 1970; 31, pp. 147-153.288279 [DOI: https://dx.doi.org/10.1016/0022-247X(70)90126-5]

39. Pommerenke, C. Padé approximants and convergence in capacity. J. Math. Anal. Appl.; 1973; 41, pp. 775-780.328090 [DOI: https://dx.doi.org/10.1016/0022-247X(73)90248-5]

40. Gončar, AA. A local condition for the single-valuedness of analytic functions of several variables. Math. USSR-Sb.; 1974; 22, pp. 305-322. [DOI: https://dx.doi.org/10.1070/SM1974v022n02ABEH001694]

41. Takens, F. In Handbook of Dynamical Systems, vol. 3, pp. 345–377 (Elsevier, 2010).

42. Austin, AP et al. Practical algorithms for multivariate rational approximation. Comput. Phys. Commun.; 2021; 261, 107663.42109431:CAS:528:DC%2BB3cXit12gtr%2FI [DOI: https://dx.doi.org/10.1016/j.cpc.2020.107663]

43. Platt, N; Sirovich, L; Fitzmaurice, N. An investigation of chaotic Kolmogorov flows. Phys. Fluids A: Fluid Dyn.; 1991; 3, pp. 681-696.1097332 [DOI: https://dx.doi.org/10.1063/1.858074]

44. Farazmand, M; Saibaba, AK. Tensor-based flow reconstruction from optimally located sensor measurements. J. Fluid Mech.; 2023; 962, A27.45838311:CAS:528:DC%2BB3sXht1GitLnL [DOI: https://dx.doi.org/10.1017/jfm.2023.269]

45. Sidi, A. Vector Extrapolation Methods with Applications (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2017). https://doi.org/10.1137/1.9781611974966.

46. Lane, J. S. & Akers, B. F. Steady Boussinesq convection: parametric analyticity and computation. Stud. Appl. Math. p. e12740 https://doi.org/10.1111/sapm.12740 (2024).

47. Jain, S; Tiso, P; Haller, G. Exact nonlinear model reduction for a von Kármán beam: Slow-fast decomposition and spectral submanifolds. J. Sound Vib.; 2018; 423, pp. 195-211. [DOI: https://dx.doi.org/10.1016/j.jsv.2018.01.049]

48. Mercer, GN; Roberts, AJ. A centre manifold description of contaminant dispersion in channels with varying flow properties. SIAM J. Appl. Math.; 1990; 50, pp. 1547-1565.1080507 [DOI: https://dx.doi.org/10.1137/0150091]

49. Domb, C; Skyes, MF. On the susceptibility of a ferromagnetic above the Curie point. Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci.; 1957; 240, pp. 214-228.1:CAS:528:DyaG2sXosVWmsQ%3D%3D

50. Dankowicz, H. & Schilder, F. Recipes for Continuation (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2013). https://doi.org/10.1137/1.9781611972573.

51. Ott, E. Chaos in Dynamical Systems, 2nd edn. (Cambridge University Press, Cambridge, 2002).

52. Hsu, MC; Bazilevs, Y. Fluid–structure interaction modeling of wind turbines: simulating the full machine. Comput. Mech.; 2012; 50, pp. 821-833. [DOI: https://dx.doi.org/10.1007/s00466-012-0772-0]

53. de Langre, E. Effects of wind on plants. Annu. Rev. Fluid Mech.; 2008; 40, pp. 141-168.2402390 [DOI: https://dx.doi.org/10.1146/annurev.fluid.40.111406.102135]

54. Mangan, NM; Brunton, SL; Proctor, JL; Kutz, JN. Inferring biological networks by sparse identification of nonlinear dynamics. IEEE Trans. Mol. Biol. Multi-Scale Commun.; 2016; 2, pp. 52-63. [DOI: https://dx.doi.org/10.1109/TMBMC.2016.2633265]

55. Brunton, SL; Proctor, JL; Kutz, JN. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proc. Natl. Acad. Sci.; 2016; 113, pp. 3932-3937.34940811:CAS:528:DC%2BC28XkvFGis7w%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27035946][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4839439][DOI: https://dx.doi.org/10.1073/pnas.1517384113]

56. Falcolini, C; De La Llave, R. Numerical calculation of domains of analyticity for perturbation theories in the presence of small divisors. J. Stat. Phys.; 1992; 67, pp. 645-666.1171147 [DOI: https://dx.doi.org/10.1007/BF01049723]

57. Haller, G; Kaszás, B. Data-driven linearization of dynamical systems. Nonlinear Dyn.; 2024; 112, pp. 18639-18663. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39219721][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11362512][DOI: https://dx.doi.org/10.1007/s11071-024-10026-x]

58. Breunung, T; Haller, G. Explicit backbone curves from spectral submanifolds of forced-damped nonlinear mechanical systems. Proc. R. Soc. A: Math. Phys. Eng. Sci.; 2018; 474, 20180083.3832468 [DOI: https://dx.doi.org/10.1098/rspa.2018.0083]

59. Ponsioen, S. Exact nonlinear model reduction in structural dynamics. Ph.D. thesis, ETH Zurich https://doi.org/10.3929/ETHZ-B-000401859 (2019).

60. Thurnher, T; Haller, G; Jain, S. Nonautonomous spectral submanifolds for model reduction of nonlinear mechanical systems under parametric resonance. Chaos; 2024; 34, 073127.4778569 [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38995987][DOI: https://dx.doi.org/10.1063/5.0168431]

61. Gonnet, P; Güttel, S; Trefethen, LN. Robust Padé approximation via SVD. SIAM Rev.; 2013; 55, pp. 101-117.3089442 [DOI: https://dx.doi.org/10.1137/110853236]

62. Mascarenhas, WF. Robust Padé approximants may have spurious poles. J. Approx. Theory; 2015; 189, pp. 76-80.3280671 [DOI: https://dx.doi.org/10.1016/j.jat.2014.10.002]

63. Chisholm, JSR. Rational approximants defined from power series in N variables. Proc. R. Soc. Lond. A. Math. Phys. Sci.; 1974; 336, pp. 421-452.361572 [DOI: https://dx.doi.org/10.1098/rspa.1974.0028]

64. Cuyt, A. How well can the concept of Padé approximant be generalized to the multivariate case? J. Comput. Appl. Math.105, 25–50 (1999).

65. Guillaume, P. & Huard, A. Multivariate Padé approximation. J. Comput. Appl. Math.121, 197–219 (2000).

66. de Montessus de Ballore, R. Sur les fractions continues algébriques. Rend. Palermo; 1905; 19, pp. 1-73. [DOI: https://dx.doi.org/10.1007/BF03014011]

67. Euler, L. De seriebus divergentibus. Novi Comment. Acad. Sci. Imp.Petropol.; 1760; 5, pp. 205-237.

68. Bromwich, T. J. I. An Introduction to the Theory of Infinite Series, 2nd edn. (Macmillan and Company, London, 1947).

69. Arnold, V. I. Geometrical Methods in the Theory of Ordinary Differential Equations, vol. 250 (Springer New York, New York, NY, 1988). https://doi.org/10.1007/978-1-4612-1037-5

70. Bates, D. M. & Watts, D. G. Nonlinear Regression Analysis and Its Applications (Wiley, New York, 1988).