🔢 mathematics

Numerical computation of high-order expansions of invariant manifolds of high-dimensional tori

This paper presents a highly parallel, efficient numerical algorithm with O(N log N) complexity that computes high-order Taylor-Fourier expansions of invariant manifolds for high-dimensional, reducible tori in stroboscopic Poincaré maps, utilizing a two-step quadratically convergent scheme and multiple shooting strategies to handle highly unstable cases.

Original authors: Joan Gimeno, Àngel Jorba, Begoña Nicolás, Estrella Olmedo

Published 2026-09-07

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Joan Gimeno, Àngel Jorba, Begoña Nicolás, Estrella Olmedo

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the vast landscape of physics and engineering, predicting the future motion of an object often relies on finding a few special paths that act as the skeleton for the entire system. Just as a mountain range defines the flow of rivers, these special paths—known as equilibrium points, periodic orbits, and quasi-periodic solutions—organize the chaotic behavior of everything from a swinging pendulum to a spacecraft navigating the solar system. When a system is simple, these paths are easy to spot. But when the system is subjected to multiple, competing rhythms, such as the gravitational tug of the Sun, Earth, and Moon all acting at once, the motion becomes quasi-periodic. This means the object never exactly repeats its path, yet it stays confined to a specific, doughnut-shaped surface in space called a torus. Understanding the shape of these surfaces and the invisible tunnels that lead into or away from them is crucial for designing stable satellite orbits or understanding the long-term stability of our solar system. However, when these surfaces exist in high-dimensional spaces with many competing frequencies, calculating their precise shape and the paths that connect to them becomes a task so computationally heavy that it has historically been nearly impossible to do with high accuracy.

A team of researchers has now developed a new computational procedure that successfully maps these complex, high-dimensional surfaces and the invisible tunnels, known as invariant manifolds, that branch off from them. The team focused on systems where the motion is driven by several different frequencies that do not line up in a simple ratio, creating a quasi-periodic rhythm. Their method works in two distinct phases. First, they calculate the precise shape of the torus itself, along with a mathematical tool that simplifies the system's behavior around it, allowing them to understand how the system reacts to small disturbances. If the torus has directions where things naturally drift away or fall in, the second phase of their method calculates the detailed, high-order expansions of these drifting paths. These paths are the "highways" of the system, guiding objects toward stability or flinging them into chaos. The researchers demonstrated that their approach can handle systems with up to five different frequencies, a level of complexity that previously required such massive amounts of computer memory and time that it was often unfeasible.

To achieve this, the researchers built a method that breaks the problem into smaller, manageable pieces that can be solved simultaneously. Instead of trying to compute the entire shape of the torus and its connecting paths in one giant, slow step, their algorithm divides the work among many processors in a computer. They use a technique called "jet transport," which allows the computer to track not just the position of an object, but also how that position changes when slightly tweaked, all the way up to very high levels of detail. This is essential because the paths they are looking for are extremely sensitive; a tiny error in calculation can lead to a completely wrong prediction of where a spacecraft might end up. By using this high-precision tracking combined with a strategy called multiple shooting, which breaks a long, difficult journey into shorter, safer segments, they can compute these paths even when the system is violently unstable. In their tests, they successfully computed the stable and unstable paths for a model of a pendulum being pushed by four different rhythms, and for a sophisticated model of the Earth-Moon system influenced by the Sun and five other natural frequencies.

The results of their work show that this new approach is not only accurate but also incredibly efficient when run on modern computers with many cores. In one experiment involving the Earth-Moon model, the researchers computed the unstable paths for a point near the Moon's orbit that is notoriously difficult to study because it is so unstable. When they used a single processor, the calculation took nearly six hours. By spreading the work across sixteen processors, they reduced the time to just twenty-six minutes, a speed-up that proves the method scales well as more computing power is added. The team also verified the accuracy of their results using several rigorous tests, confirming that the calculated paths satisfy the fundamental laws of motion to a precision of one part in ten billion. This level of precision is vital for real-world applications, such as planning the trajectory of a satellite that needs to stay in a specific orbit for years without crashing or drifting away.

The significance of this work lies in its ability to handle the "curse of dimensionality," where adding more frequencies to a system usually makes the calculation time explode exponentially. The researchers found that their method grows much more slowly, making it possible to study systems that were previously out of reach. They explicitly noted that while their current code works for systems with one unstable direction, the underlying logic can be extended to handle systems with multiple unstable directions, opening the door to even more complex models. They also highlighted that their approach is ready to be adapted for graphics processing units, the specialized chips found in modern computers that are designed for massive parallel processing, which could make these calculations even faster in the future. By providing a clear, efficient, and highly parallel way to map these invisible structures, the researchers have given scientists and engineers a powerful new tool to navigate the complex, multi-rhythmic dance of the universe.

Technical Summary: Numerical Computation of High-Order Expansions of Invariant Manifolds of High-Dimensional Tori

Problem Statement
The paper addresses the challenge of accurately computing reducible invariant tori and their associated stable and unstable invariant manifolds within stroboscopic Poincaré maps. These structures are critical for understanding the global dynamics of quasi-periodically forced systems, such as spacecraft motion in the Earth-Moon-Sun system. The specific difficulties addressed include:

High Dimensionality: As the number of incommensurable frequencies increases, the dimension of the torus grows, requiring a large number of Fourier modes ( $N$ ) for accurate representation.
High-Order Accuracy: Standard linear approximations are insufficient for globalizing manifolds in highly unstable regimes; high-order Taylor-Fourier expansions are required.
Computational Cost: Traditional methods for computing these objects scale poorly with the number of Fourier modes and the dimension of the phase space, often becoming infeasible for high-dimensional tori.
Strong Instability: In cases where the torus is highly unstable (large hyperbolic eigenvalues), direct integration of the flow over the Poincaré map period leads to significant numerical errors, necessitating specialized strategies like multiple shooting.

Methodology
The authors propose a two-step numerical procedure that combines the parametrization method with jet transport and parallel computing.

Computation of the Invariant Torus and Floquet Data:
- The method utilizes a quadratically convergent Newton scheme to simultaneously compute the Fourier series of the torus $\phi(\theta)$ , the Floquet transformation $C(\theta)$ , and the Floquet matrix $B$ .
- This step relies on the assumption that the torus is reducible, allowing the linearized dynamics around the torus to be transformed into a constant coefficient system.
- The algorithm solves cohomological equations for Fourier coefficients independently for each mode, facilitating parallelization.
Computation of Invariant Manifolds:
- Parametrization: The stable and unstable manifolds are represented as Taylor-Fourier expansions $W(\theta, \sigma) = \sum a_k(\theta)\sigma^k$ . The invariance equation $P(W(\theta, \sigma), \theta) = W(\theta + \rho, \lambda\sigma)$ is solved order by order.
- Jet Transport: To compute the high-order derivatives of the Poincaré map required for the Taylor expansion coefficients, the authors employ "jet transport." This technique uses automatic differentiation to propagate truncated power series (jets) through the numerical integration of the ODE, avoiding the tedious manual derivation of high-order variational equations.
- Stable Manifold Strategy: For stable manifolds, the method utilizes the inverse Poincaré map to avoid the amplification of numerical errors caused by unstable directions during forward integration.
- Multiple Shooting: For highly unstable tori, the integration interval is split into $r$ sections. The invariance equations are modified to link these sections, effectively reducing the integration time per step and controlling error propagation. The linear reducibility is extended to this multi-section framework using block matrices.
Parallel Implementation:
- The algorithms are implemented in C using OpenMP for shared-memory parallelism.
- Parallelism is exploited at two levels:
  - Map Evaluation: The evaluation of the Poincaré map and its derivatives (via ODE integration) for different angle values in the mesh is performed concurrently across CPU cores.
  - Linear Solvers: The solution of the decoupled linear systems for Fourier coefficients in the cohomological equations is parallelized.
- The use of FFTW3 for Fast Fourier Transforms is integrated, though the primary speedup comes from the ODE integration steps.

Key Contributions

Extension of Parametrization Methods: The paper extends the reducible torus computation methods (previously limited to the torus and Floquet data in [JO09]) to include high-order parametrizations of stable and unstable invariant manifolds.
Jet Transport Integration: It demonstrates the application of jet transport to compute high-order derivatives of Poincaré maps for invariant manifolds, a technique previously applied mainly to fixed points.
Multiple Shooting for Manifolds: The authors adapt the multiple shooting strategy to the computation of high-order manifold expansions, addressing the issue of strong instability that prevents accurate single-shooting integration.
Scalability and Efficiency: The proposed algorithm exhibits a computational complexity of $O(N \log N)$ and memory usage of $O(N)$ , where $N$ is the number of Fourier modes. The high degree of parallelism allows the method to scale efficiently with the number of processors.

Results
The authors validate the method through two applications:

Quasi-Periodically Forced Pendulum: A system with a 4-dimensional frequency vector (resulting in a 5-dimensional flow torus, reduced to 4D in the map). The method successfully computed the torus and invariant manifolds up to order 10.
- Performance: Using 16 cores, the computation of the torus achieved a speed-up of 12.38, and the manifold computation achieved a speed-up of 13.74 (unstable) and 12.25 (stable).
- Accuracy: Tests confirmed that the invariance equations were satisfied with residuals around $10^{-11}$ to $10^{-13}$ , and the decay of Fourier coefficients confirmed the analyticity of the solutions.
Earth-Moon System (SSSM Model): A celestial mechanics model involving five natural frequencies. The authors computed invariant tori replacing the $L_1$ $L_{1}$ and $L_2$ $L_{2}$ equilibrium points in the presence of solar perturbations.
- Instability Handling: Due to extreme instability (eigenvalues $\approx 10^2$ ), multiple shooting with 3 and 4 sections was employed.
- Performance: Significant speed-ups were observed (e.g., 13.09 for the $L_2$ unstable manifold on 16 cores), reducing computation time from nearly 6 hours to 26 minutes.

Significance and Claims
The paper claims that the computation of high-order Taylor-Fourier expansions for invariant manifolds of high-dimensional tori is now feasible, even for systems with strong instability. The significance lies in:

Feasibility: Demonstrating that complex, high-dimensional dynamical structures can be computed with high precision using modern parallel hardware.
Robustness: The combination of jet transport and multiple shooting allows for the accurate computation of manifolds in regimes where standard integration fails.
Efficiency: The algorithm's $O(N \log N)$ scaling and high parallelism make it suitable for modern multi-core architectures, overcoming the bottlenecks of previous distributed-memory implementations.

The authors conclude that while the current work focuses on fixed rotation dynamics and known frequency vectors, the framework is adaptable to autonomous systems and unknown frequencies, and future work will explore GPU parallelization to further enhance performance.

Technical Summary: Numerical Computation of High-Order Expansions of Invariant Manifolds of High-Dimensional Tori

More like this