Source-linked AI summary
How to Compute Invariant Manifolds and their Reduced Dynamics in High-Dimensional Finite-Element Models
Shobhit Jain, George Haller
TL;DR
High-dimensional finite-element models make invariant-manifold computation challenging because established continuation approaches are limited as system dimension grows. The paper develops a scalable parametrization-based methodology and demonstrates fast, reliable reduced-dynamics computations in large mechanical structures.
Problem
Existing invariant-manifold computations and continuation packages primarily target low-dimensional systems, limiting their applicability to mechanics problems with thousands to millions of degrees of freedom.
Method
The paper reformulates the parametrization method in physical coordinates, using only master modes to compute local invariant manifolds and their reduced dynamics.
Results
The methodology computes Lyapunov subcenter manifolds, conservative backbone curves, spectral submanifolds, and forced response curves, including fast and reliable computations for very large structures.
Takeaways & Limitations
The approach facilitates mathematically rigorous nonlinear model reduction for very high-dimensional finite-element mechanical systems.
Takeaways & Limitations
Certain outer resonances prevent the invariant manifold from existing, while aircraft-wing coefficient assembly can require substantial memory.
Abstract
from arXiv · showhide
Invariant manifolds are important constructs for the quantitative and qualitative understanding of nonlinear phenomena in dynamical systems. In nonlinear damped mechanical systems, for instance, spectral submanifolds have emerged as useful tools for the computation of forced response curves, backbone curves, detached resonance curves (isolas) via exact reduced-order models. For conservative nonlinear mechanical systems, Lyapunov subcenter manifolds and their reduced dynamics provide a way to identify nonlinear amplitude-frequency relationships in the form of conservative backbone curves. Despite these powerful predictions offered by invariant manifolds, their use has largely been limited to low-dimensional academic examples. This is because several challenges render their computation unfeasible for realistic engineering structures described by finite-element models. In this work, we address these computational challenges and develop methods for computing invariant manifolds and their reduced dynamics in very high-dimensional nonlinear systems arising from spatial discretization of the governing partial differential equations. We illustrate our computational algorithms on finite-element models of mechanical structures that range from a simple beam containing tens of degrees of freedom to an aircraft wing containing more than a hundred-thousand degrees of freedom.
1 Introduction
Computing invariant manifolds in high-dimensional mechanics is hindered by the computational cost of global continuation and by diagonalization and coordinate changes in local parametrization methods. The paper develops a scalable physical-coordinate methodology using only the minimum number of eigenvectors, with broad applicability and an open-source implementation.
- Computational challenges: Global collocation and spectral approaches become computationally intractable as the system dimension N grows in PDE-discretized mechanical systems.These methods require solving O(...) nonlinear algebraic equations at each continuation step, while existing packages mainly demonstrate automated continuation on low-dimensional examples.
- Computational challenges: Diagonalization-based parametrization methods are unfeasible for high-dimensional mechanics because they require all N eigenvalues and destroy sparsity in physical coordinates.The associated linear coordinate change also removes the inherent sparsity of nonlinear coefficients in mechanics applications.
- Proposed methodology: The paper develops a scalable local parametrization methodology that uses physical coordinates and only the minimum number of eigenvectors.The method derives Taylor-series coefficient expressions for invariant manifolds attached to fixed points in the autonomous setting.
- Scope and implementation: The procedure applies to smooth first-order systems with periodic or quasiperiodic forcing and computes stable, unstable, center, Lyapunov subcenter, and spectral submanifolds.It supports arbitrarily high orders of accuracy while retaining physical coordinates.
- Scope and implementation: The computational techniques are implemented in the open-source MATLAB package SSMTool 2.0.The implementation is intended to automate computation across the invariant-manifold types covered by the methodology.
2 General setup
The paper formulates finite-element mechanical models as first-order systems and defines invariant structures through generalized eigenspaces and spectral subspaces. Although modal diagonalization supports theoretical invariant-manifold analysis, it is computationally intractable for high-dimensional finite-element models.
- Mechanical model: Finite-element discretization yields second-order mechanical equations with mass, damping, stiffness, nonlinear internal forces, and autonomous, periodic, or quasiperiodic external forcing.The forcing is characterized by a frequency vector with K=0, K=1, or K>1 rationally incommensurate frequencies.
- First-order formulation: The methodology applies to any first-order system, including mechanical systems of dimension N=2n, with equivalent formulations determined by a nonsingular matrix N.When the mechanical matrices are symmetric, suitable choices of N make the first-order matrices symmetric.
- Invariant structures: Generalized eigenspaces of the linearized system generate invariant spectral subspaces, including stable, unstable, and center subspaces, while SSMs provide smoothest nonlinear continuations of stable spectral subspaces.SSM existence and uniqueness under appropriate conditions are obtained using the parametrization method [52].
- Computational limitation: Modal diagonalization is useful for proving invariant-manifold existence and uniqueness but is computationally intractable for actual computation in high-dimensional finite-element problems.The parametrization framework uses right and left eigenvectors spanning a master spectral subspace and its adjoint.
3 Pitfalls of the diagonalized form (14)
The diagonalized form (14) becomes impractical for high-dimensional finite-element systems because complete spectral data are costly to compute and coordinate transformation destroys sparsity, causing enormous coefficient-storage requirements.
- 3 Pitfalls of the diagonalized form (14): The physical-coordinate nonlinearity is represented with Kronecker products, where z⊗k contains N^k degree-k monomials and F_k stores their coefficients.For vectors, z⊗z produces zz^T with entries z_i z_j; higher-order tensor products generalize this construction.
- 3 Pitfalls of the diagonalized form (14): The same coefficient-based representation is used for the modal-coordinate nonlinearity T in q∈C^N.The diagonalized formulation expresses nonlinear terms in modal coordinates using the corresponding polynomial coefficient arrays.
- 3.1 Eigenvalue and eigenvector computation: Complete eigenvector bases remain infeasible for very high-dimensional systems, although small subsets of eigenvalues and eigenvectors can be computed and stored efficiently.Invariant-manifold computations typically require all eigenvectors to remain in active memory.
- 3.1 Eigenvalue and eigenvector computation: Nearly a year is required to compute all eigenvalues as plate degrees of freedom grow to approximately 100,000, whereas sparse iterative solvers compute only small subsets efficiently.Figure 2 compares full eig computations with MATLAB eig against five smallest-magnitude eigenpairs computed using eigs.
- 3.2 Unfeasible memory requirements due to coordinate-change: Local invariant-manifold approximations transform the system into modal coordinates, represent the manifold over a master subspace, and solve Taylor-series invariance equations recursively.The procedure substitutes a polynomial ansatz into the diagonalized system at each order.
- 3.2 Unfeasible memory requirements due to coordinate-change: Transforming z=Vq destroys the governing equations’ physical-coordinate sparsity, increasing polynomial coefficient counts by orders of magnitude and making memory use unfeasible.Physical sparsity arises because numerical meshes couple primarily neighboring degrees of freedom.
- 3.2 Unfeasible memory requirements due to coordinate-change: Several terabytes to petabytes of memory are required to store transformed nonlinear coefficients even for moderately sized square-plate meshes.Figure 3 compares coefficient storage in physical and modal coordinates as phase-space dimension increases.
4 Computing invariant manifolds of fixed points in physical coordinates
The method computes invariant manifolds directly in physical coordinates using only eigenvectors and eigenvalues of the master subspace, together with their reduced dynamics. It solves polynomial invariance equations recursively while exploiting parametrization freedom, including choices suited to resonances and real-valued computations.
- 4 Computing invariant manifolds of fixed points in physical coordinates: Invariant manifolds are computed in physical coordinates using only the eigenvectors and eigenvalues associated with the master subspace E, unlike commonly used approaches [37] [41].
- 4 Computing invariant manifolds of fixed points in physical coordinates: The parametrization W maps M-dimensional coordinates p to the full phase space, while R defines reduced dynamics whose trajectories map onto full-system trajectories on the manifold.
- 4 Computing invariant manifolds of fixed points in physical coordinates: W and R are represented by multivariate polynomial expansions, with sparse physical-coordinate nonlinearity coefficients that remain straightforward to store despite their large size.
- 4 Computing invariant manifolds of fixed points in physical coordinates: At each degree i≥2, the polynomial invariance equation is solved recursively by vectorizing the cohomological equation for the unknown manifold coefficients W_i.The order-i cohomological operator is determined entirely by the linear full and reduced systems, through A, B, and R_1.
- 4 Computing invariant manifolds of fixed points in physical coordinates: Parametrization underdeterminacy permits context-dependent choices: real master-subspace bases halve memory requirements, while suitable reduced dynamics can preserve solvability when resonances make the operator singular.
4.1 Choice of parametrization
The parametrization choice is governed by the resonance structure of the cohomological operator: inner resonances can be handled by adjusting the reduced dynamics, whereas outer resonances obstruct solvability. The authors derive physical-coordinate procedures for normal-form and graph-style parametrizations using only master-subspace eigenvectors, with the graph style avoiding a priori resonance identification.
- Choice of parametrization: Inner resonances can be accommodated by adjusting the reduced dynamics, but outer resonances make the cohomological equation generally unsolvable because its right-hand side cannot be placed in the operator range.Both resonance types render Li singular; inner resonances permit parametrization adjustment, while outer resonances do not.
- Normal form parametrization: The normal-form procedure computes reduced dynamics directly in physical coordinates by projecting onto the left kernel of Li, using only eigenvalues and eigenvectors associated with the master spectral subspace.Nontrivial reduced coefficients are retained precisely for inner-resonant monomials, providing the normal-form coefficients without transforming to diagonal coordinates.
- Normal form parametrization: In lightly damped systems, near-resonances between eigenvalue imaginary parts should be included in the normal-form parametrization because exact resonance conditions are seldom satisfied.The relevant resonance relations require simultaneous agreement of real and imaginary parts, so near-resonant modes can matter when real parts are small.
- Graph style parametrization: The normal-form parametrization retains only resonant monomials, whereas the graph style generally includes more nontrivial monomials and therefore produces more complex reduced dynamics.The graph style directly provides reduced-dynamics coefficients without evaluating Gi and does not require specific inner resonances to be identified beforehand.
4.2 Applications
The methodology computes parameter-dependent center manifolds and normal-form reduced dynamics without diagonalization, recovering the Lorenz pitchfork bifurcation. Applied to conservative mechanical systems, it computes Lyapunov subcenter manifolds and directly yields conservative backbone curves under non-resonance conditions.
- Lorenz-system application: The Lorenz computation detects inner resonances within the two-dimensional center subspace but no outer resonances with the remaining eigenvalues, enabling a nontrivial normal-form reduced dynamics.The reduced dynamics are obtained recursively to higher order after selecting the normal-form parametrization.
- Lorenz-system application: The computed parameter-dependent center-manifold dynamics recover the Lorenz system’s pitchfork bifurcation with respect to µ = ρ − 1.The construction uses the two-dimensional center subspace associated with the zero eigenvalues and a normal-form parametrization.
- Lyapunov subcenter manifolds: An LSM is a unique analytic two-dimensional nonlinear extension of a modal subspace spanned by a distinct eigenvalue pair ±iω_m when the master frequency is non-resonant with the remaining frequencies.The LSM is computed from the invariance equation in physical coordinates using only the master mode ϕ_m.
- Lyapunov subcenter manifolds: For a conservative mechanical system, the method directly provides the backbone curve relating periodic-orbit amplitude ρ and frequency ω on the Lyapunov subcenter manifold.This result holds for the LSM associated with mode ϕ_m under the non-resonance condition ω_m ≠ kω_i for all integer k and remaining modes i.
5 Invariant manifolds and their reduced dynamics under non-autonomous forcing
Under non-autonomous forcing, invariant manifolds become whiskers attached to invariant tori or periodic orbits and require angular-variable parameterization. Their leading-order Fourier formulation handles near-resonances, while the reduced dynamics yields forced-response curves and isolas for high-dimensional systems.
- Forced invariant manifolds: The parametrization method computes an (M + K)-dimensional whisker attached to a forced invariant torus, with manifold and reduced dynamics parameterized by forcing phases.The torus is periodic for one forcing frequency and quasiperiodic for multiple frequencies; the construction reduces to the autonomous manifold when ϵ = 0.
- Fourier computation: Leading-order non-autonomous manifold and reduced-dynamics terms are computed through Fourier expansions, producing independent linear systems for each forcing multi-index κ.The independent coefficient solves enable parallel computation with minimal cross-process communication.
- Resonance treatment: Near-resonant forcing creates small divisors, so resonant terms are incorporated into the reduced dynamics through normal-form or graph parametrizations.If resonance involves modes outside the master subspace, the Taylor approximation’s convergence domain is reduced; including those modes can mitigate the issue.
- Forced-response curves: The generalized reduced dynamics extracts forced-response curves near resonance for periodically forced systems, extending prior SSM formulas beyond proportional damping and synchronous mono-harmonic forcing.Forcing-response curves arise from fixed points of the polar reduced system, corresponding to periodic orbits with frequency ηΩ.
- Applications and scope: The resulting zero-level-set formulation identifies detached resonance curves, or isolas, and remains applicable to high-dimensional systems with inner resonances at arbitrarily high accuracy.The formulation also supports choosing equivalent scalar functions to avoid numerical difficulties.
6 Numerical examples
The numerical examples demonstrate that local SSM computations recover nonlinear forced-response curves and stable or unstable periodic orbits from reduced dynamics across finite-element mechanical systems. They agree with global continuation results while remaining computationally feasible for substantially larger models.
- Computational method: The method computes FRCs by solving local invariance equations, simulating two-dimensional SSM reduced dynamics, and mapping polar-coordinate trajectories onto the manifold.The resulting reduced dynamics on SSMs attached to periodic orbits provide the FRCs through Lemma 2.
- Mechanical examples: The examples span finite-element-like oscillator chains, beams, shallow arches, and an aircraft wing, with responses showing softening in shallow arches and expected hardening in the wing.For the wing, stable and unstable FRC periodic orbits are obtained from hyperbolic fixed points of the reduced dynamics.
- FRC validation: Local SSM computations closely reproduce FRCs from collocation and harmonic-balance continuation and identify stable and unstable periodic orbits through hyperbolic reduced-dynamics fixed points.This agreement is reported for both the oscillator-chain and beam examples, with NLvib used for harmonic balance and coco for collocation.
- Computational scale: For the aircraft-wing-scale computation, model generation and coefficient assembly required 1 hour 21 minutes 38 seconds without parallelization, with memory consumption the main bottleneck.Computing the first 10 eigenvalues took approximately 5 seconds; coefficient assembly caused the peak memory demand.
7 Conclusions
The work reformulates invariant-manifold and reduced-dynamics computation for high-dimensional nonlinear mechanics in physical coordinates, avoiding infeasible full linear system diagonalization. It demonstrates applications to center manifolds, LSMs and backbone curves, and SSMs and forced response curves, while identifying extensions to higher-dimensional manifolds and parallel computation.
- Computational methodology: The parametrization method enables invariant manifolds and reduced dynamics to be computed in physical coordinates when classical linear system diagonalization is infeasible.The reformulation targets high-dimensional nonlinear mechanics problems.
- Engineering applications: The methodology covers parameter-dependent center manifolds, Lyapunov subcenter manifolds with conservative backbone curves, and spectral submanifolds with forced response curves.It also supports fast and reliable forced-response-curve computations through a normal-form-style parametrization of spectral submanifolds in very large systems.
- Extensions: Although the examples focused on two-dimensional spectral submanifolds, the automated procedure and implementation can treat higher-dimensional invariant manifolds.Their reduced dynamics can directly compute forced response curves in internally resonant systems with energy transfer among multiple modes, as planned in Li et al. [73] and Li & Haller [74].
- Future improvements: Parallelization is a planned improvement because the invariance equations can be solved independently for each monomial/Fourier multi-index.This development was underway and reserved for future reporting.
Funding
The authors report that no specific funding was received for this work.
- No specific funding was received for this work.
A Basic propositions
The section establishes Kronecker-product spectral properties for semisimple matrices and generalized eigenvalue matrix pairs. These results also characterize the associated tensor-product eigenvectors and, when the eigenvalue product equals one, a resulting singular matrix and its kernels.
- Proposition 1: Proposition 1 shows that the i-fold Kronecker power of a semisimple matrix is semisimple, with eigenvalues formed by products of the original eigenvalues and tensor-product eigenvectors.Linear independence of the original eigenvectors and Kronecker-product arguments establish semisimplicity.
- Proposition 2: Proposition 2 shows that λµ is a generalized eigenvalue of the Kronecker matrix pair (C ⊗ A),(D ⊗ B), with left and right eigenvectors e ⊗ u and f ⊗ v.The result follows by applying the generalized eigenvalue relations for the two original matrix pairs and Kronecker-product identities.
- Proposition 2: When λµ = 1, Proposition 2 further states that the resulting matrix is singular, with f ⊗ v and e ⊗ u belonging to its right and left kernels.The kernel characterization is obtained from the corresponding Kronecker-product eigenvector relations.
B Proof of Lemma 1
The proof shows that reduced dynamics on the Lyapunov subcenter manifold has constant-amplitude steady states and amplitude-dependent angular frequency. Because the manifold is foliated by periodic orbits, this frequency-amplitude relation defines the conservative backbone around mode m.
- B Proof of Lemma 1: The scalar amplitude equation decouples from phase, so its only steady states have constant amplitude ρ = ρ(0).
- B Proof of Lemma 1: The phase equation ˙θ = ω(ρ) therefore gives a constant angular frequency determined by the steady-state amplitude.The reduced-dynamics functions a and ω are analytic because the LSM is analytic.
- B Proof of Lemma 1: The periodic-orbit foliation of the LSM makes ω(ρ) the oscillation frequency at each constant amplitude, yielding the conservative backbone around mode m through relation (69).
C Proof of Lemma 2
The proof derives the polar reduced dynamics, characterizes fixed points as periodic reduced-system orbits whose projection yields the FRC, obtains the phase shift, and establishes stability through the Jacobian eigenvalues.
- Polar reduced dynamics: Comparing real and imaginary parts after a phase shift yields the polar reduced dynamics, proving statement (i).The phase shift is introduced as ψ = θ − ηφ after dividing by e^iθ.
- Fixed points and FRC: Fixed points of the polar reduced dynamics represent periodic reduced-system orbits, and their projection onto (ρ, Ω) provides the FRC.The fixed points form a one-dimensional zero submanifold when solving the fixed-point equations; eliminating ψ gives the corresponding (ρ, Ω) relation.
- Phase shift: Solving the fixed-point equations for cos ψ and sin ψ determines the phase shift, proving statement (iii).The phase shift follows from the resulting trigonometric expressions.
- Stability: The Jacobian evaluated at each fixed point provides eigenvalues that determine the stability of hyperbolic fixed points by linearized stability analysis.This standard-form formulation establishes statement (iv).