Robust computation of higher-dimensional invariant tori from individual trajectories

TL;DR

The paper addresses the challenge of computing invariant tori of dimension $d>1$ from a single, short trajectory without continuation or initial guesses and without fixing a coordinate system. It develops a four-step pipeline: (i) Birkhoff reduced rank extrapolation to identify multiple frequencies, (ii) Bayesian MAP estimation to infer a valid rotation vector, (iii) Korkine-Zolatarev lattice reduction to select an optimal homology basis, and (iv) a least-squares Fourier parameterization of the torus. The authors demonstrate robustness by recovering many $d=2$ tori in a standard map and extending to $d=3$ tori in the Earth–Moon ER3BP, including island structures, while diagnosing near-resonant and filamentary cases as primary failure modes. The approach is coordinate-agnostic and non-perturbative, enabling practical torus computation in physical systems with limited a priori knowledge, and the accompanying SymplecticMapTools.jl package provides public access to the implementation.

Abstract

We present a method for computing invariant tori of dimension greater than one. The method uses a single short trajectory of a dynamical system without any continuation or initial guesses. No preferred coordinate system is required, meaning the method is practical for physical systems where the user does not have much \textit{a priori} knowledge. Three main tools are used to obtain the rotation vector of the invariant torus: the reduced rank extrapolation method, Bayesian maximum a posteriori estimation, and a Korkine-Zolatarev lattice basis reduction. The parameterization of the torus is found via a least-squares approach. The robustness of the algorithm is demonstrated by accurately computing many two-dimensional invariant tori of a standard map example. Examples of islands and three-dimensional invariant tori are shown as well.

Figures (10)

• Figure 1: A trojan invariant torus for the ER3BP plotted in the $(\xi,\eta)$, $(\dot \xi, \dot \eta)$, $(\xi,\dot \xi)$, and $(\eta, \dot \eta)$ projections. The length $601$ trajectory is given in red, while the computed invariant torus with lines that are constant in coordinates are given in gray.
• Figure 2: The validation error \ref{['eq:validation-error']} for the (left) projection and (right) least-squares methods, for both the unweighted and weighted versions of each method. The red path on the unweighted least-squares method indicates the path of the adaptive method (see Sec. \ref{['subsec:fourier-adaptive']}).
• Figure 3: A schematic of how the Fourier coefficients for the ER3BP example shear under the transformation $L_1$.
• Figure 4: Four projections of a trajectory of the coupled standard map \ref{['eq:coupled-sm']}.
• Figure 5: A $C^\infty$-windowed discrete Fourier transform of a length $10000$ trajectory of the 2D standard map \ref{['eq:coupled-sm']} shown in Fig. \ref{['fig:wonky-standard-map']} with the observable function in Eq. \ref{['eq:sm-observable']}. The 60 highest peaks $(\Omega_j, H_j)$ for $1\leq j \leq 60$ are plotted in black dots, cf. Table \ref{['tab:wonky-coefficients']}. The only two peaks in the 15 largest that form a valid rotation vector are marked with red stars ($j=12,14$). The optimal rotation vector for the anisotropy of the invariant torus is marked in magenta 'x's ($j=31,60$).