Direct Parametrisation of Invariant Manifolds.

We write the governing equations as an $\mathfrak{n}$-th order polynomial ODE for the physical state $\mathbf{u}(t) \in \mathbb{C}^{N}$:

with linear coefficient matrices $\mathbf{B}_k \in \mathbb{C}^{N\times N}$ and an equilibrium at the origin (with $\mathbf{r}=\mathbf{0}$). The right‑hand side $\mathbf{F}$ encodes all nonlinearities and external couplings in polynomial form. The external state $\mathbf{r}\in\mathbb{C}^{\widehat{n}}$ obeys an autonomous polynomial ODE $\dot{\mathbf{r}} = \mathbf{E}(\mathbf{r})$ with equilibrium at the origin. We solve our cohomological equations directly on the $N$-dimensional $\mathfrak{n}$-th order system; we never convert to a $\mathfrak{n}N$-dimensional first-order system.

Most reduced‑order models constrain the solution to a linear subspace (e.g. a modal subspace) and project the dynamics onto a linear subspace (the same one in Galerkin methods), neglecting the orthogonal complement. We instead search for an invariant manifold tangent to a chosen master (spectral) subspace and seek a smooth parametrisation $\mathbf{W}: \mathbb{C}^{n} \to \mathbb{C}^{N}$:

where $\mathbf{z} = (\widetilde{\mathbf{z}}, \mathbf{r})\in\mathbb{C}^{n}$ combines the $\widetilde{n}$ master‑mode coordinates and the $\widehat{n}$ external states, and $n = \widetilde{n}+\widehat{n}$. We expand both $\mathbf{W}$ and $\mathbf{R}$ as power series in $\mathbf{z}$:

where $\mathbf{z}^{\bm{\alpha}} = z_1^{\alpha_1} z_2^{\alpha_2} \cdots z_n^{\alpha_n}$ for multiindices $\bm{\alpha}\in\mathbb{N}^{n}$. We solve the invariance condition monomial by monomial in graded‑lexicographic order; at each monomial $\mathbf{z}^{\bm{\gamma}}$, the equations involve only previously computed coefficients, yielding a causal hierarchy.

We linearise at the origin and solve the polynomial eigenvalue problem

Selecting a set of $\widetilde{n}$ eigensolutions defines the master spectral subspace $E$. We also support defective cases: if a Jordan chain arises, we include the corresponding generalised eigenvectors. The invariant manifold that we parametrise is tangent to $E$ at the origin.

At first order, the parametrisation $\mathbf{W}(\mathbf{z}) = \mathbf{\Phi} \mathbf{z} + \mathcal{O}(\|\mathbf{z}\|^2)$ expands along the generalised eigenvectors, while the reduced dynamics $\mathbf{R}(\mathbf{z}) = \mathbf{\Lambda} \mathbf{z} + \mathcal{O}(\|\mathbf{z}\|^2)$ inherits the associated eigenvalues (or Jordan blocks).

At higher orders, nonlinearities curve the manifold away from the flat spectral subspace; the parametrisation and the reduced dynamics encode this geometric deformation monomial by monomial.

At each monomial $\mathbf{z}^{\bm{\gamma}} = z_1^{\gamma_1} z_2^{\gamma_2} \cdots z_n^{\gamma_n}$ we have a linear system for the unknown coefficients $\mathbf{W}_{\bm{\gamma}}$ and $\mathbf{R}_{\bm{\gamma}}$:

with $\mathbf{L} = \mathbf{B}_{0} + \mathbf{B}_{1} s + \dotsm + \mathbf{B}_{\mathfrak{n}} s^{\mathfrak{n}} $ given the superharmonic $s = \langle \bm{\gamma}, \bm{\lambda} \rangle = \gamma_1 \lambda_1 + \gamma_2 \lambda_2 + \dotsm + \gamma_n \lambda_n$. The operator $\mathbf{C}$ built from $s$, the $\mathbf{B}_k$ matrices and the eigenvectors. We then impose orthogonality conditions to fix the non‑unique resonant master components: for each resonant master mode with left eigenvector $\bm{\ell}$, we set $\bm{\ell}^\dag \mathbf{W}_{\bm{\gamma}}^{(0)} = 0$ (or a user‑supplied value) by adding the corresponding rows. The non‑resonant rows of $\mathbf{R}_{\bm{\gamma}}$ we set to zero, leaving a square system of $N + \widetilde{n}_\mathrm{res}$ equations that we solve directly.

The operator $\mathbf{L}(s)$ becomes ill‑conditioned when $s$ approaches an eigenvalue $\lambda$ of the linearised problem. We detect such inner resonances with the master modes and outer resonances with the neglected modes. For inner resonances, we absorb the corresponding monomial into the reduced dynamics $\mathbf{R}$ — this is precisely what yields the normal form. For outer resonances, we require the right‑hand side to be orthogonal to the left nullspace of $\mathbf{L}(\lambda)$; if not, we must augment the master subspace with those outer modes. MORFE offers four resonance‑classification strategies:

• Graph – every monomial of degree $\ge 2$ is resonant with all master modes (conservative but may inflate the ROM).
• Complex normal form – flags a mode when $|s-\lambda| < \varepsilon$ for a user‑supplied tolerance.
• Real normal form – extends to conjugate pairs: resonant if $|s-\lambda| < \varepsilon$ or $|s-\lambda^*| < \varepsilon$.
• Condition‑number estimate – flags when $\rho\,\kappa(\lambda)/|s-\lambda| > \kappa_{\max}$, with $\rho$ the spectral radius and $\kappa(\lambda)$ the eigenvalue condition number, giving an adaptive, physics‑aware criterion.

Users may also manually mark any monomial as resonant or non‑resonant, a crucial override for known internal resonances (e.g. 1:2 or 1:3) that automatic tests might miss due to numerical tolerances.

We treat external forcing as an autonomous ODE in the augmented state $\mathbf{z} = (\widetilde{\mathbf{z}}, \widehat{\mathbf{z}})$:

where $\mathbf{Q}$ is the Schur factor that makes $\mathbf{Q}^{-1}\widehat{\mathbf{\Lambda}}\mathbf{Q}$ upper triangular. The eigenvalues of $\widehat{\mathbf{\Lambda}}$ become part of the superharmonic $s = \bm{\lambda}\cdot\bm{\gamma}$, so monomials involving external coordinates enter the cohomological equations on the same footing as internal ones. Resonances involving $\widehat{\mathbf{z}}$ are automatically detected and absorbed into the reduced dynamics — no averaging, no Floquet transformation, no separate treatment of periodic, quasi‑periodic, or even chaotic forcing. The external system can be as simple as a harmonic oscillator ($\dot{\widehat{\mathbf{z}}} = \mathrm{i}\Omega\,\widehat{\mathbf{z}}$) or as complex as a Lorenz attractor; we handle polynomial nonlinearities in $\mathbf{g}$ as well, enabling fully coupled reduced‑order models of driven nonlinear structures.

We never assemble or store the high‑dimensional multilinear operators $\mathbf{G}_{ijk} \in \mathbb{R}^{N\times N\times N}$ or $\mathbf{H}_{ijkl} \in \mathbb{R}^{N\times N\times N\times N}$. Instead, each nonlinear term is a FEMMultilinearMap that exposes an element‑level interface: scatter_qp! interpolates a $\mathbf{W}$-column to quadrature‑point gradients, accumulate_qp! evaluates the multilinear integrand at each quadrature point, and assemble_element! scatters the element residual to global DOFs. During the DPIM solve, a single element loop per monomial deduplicates scatter operations across terms: if multiple terms share a $\mathbf{W}$-column, we interpolate it to quadrature‑point gradients only once per element. The element‑local working buffer scales as $\mathcal{O}(n_\text{unique} \times n_\text{qp})$ — independent of $N^3$ or $N^4$. This batched assembly strategy, implemented with Horner‑style recurrence for $\mathbf{L}$, $\mathbf{\Xi}$, and $\mathbf{C}$, makes MORFE equally efficient for models with tens of thousands of DOFs and for small analytical benchmarks.

The DPIM core accepts multilinear terms $\mathbf{T}_j$ in two complementary forms. For analytically defined systems, we supply closed‑form expressions — a Duffing‑type cubic stiffness $\beta x_1 x_2 x_3$, velocity‑squared damping $\gamma \dot{x}_1\dot{x}_2$, or mixed parametric coupling $x_1 x_2 \dot{x}\,\widehat{z}$ — and the algebra evaluates them directly. For large‑scale finite‑element models, we delegate the physics to the FEM backend: the module evaluates the weak‑form integrand at quadrature points for given discrete field values and their derivatives, while the DPIM solver drives the monomial loop and assembles the global right‑hand side. This clean separation keeps the DPIM algebra entirely independent of spatial discretisation; any FEM library and any polynomial physical model slot in without modifying the solver core. We have demonstrated this modularity with Gridap.jl, Ferrite.jl, and the legacy MORFE2.0 backend, covering geometrically nonlinear structural mechanics, piezoelectric coupling, and thermomechanical interaction with identical solver logic.