Variational Continuation for Double Pendulum Periodic Orbits Leo Yao,∗ Ziming Liu,† and Max Tegmark‡
arXiv:2609.05337v1 [cs.LG] 4 Sep 2026
Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts and The NSF AI Institute for Artificial Intelligence and Fundamental Interactions (Dated: September 7, 2026) We present a Hessian-based approach to numerically continue periodic orbits in dynamical systems. A loop (periodic orbit candidate) is parametrized as a Fourier series; a loss function is defined based on the deviation of the loop from the physical differential equations. Unlike previous work relying on hand-derived Jacobians, our method automates the process by leveraging automatic differentiation, a common machine learning technique. The continuation direction can be determined by the flat directions of the loss landscapes (directions with zero eigenvalues), making the search of periodic orbits efficient and guided. Our method is integrator-free, precisely initializes oscillations around unstable fixed points, and efficiently detects orbit family intersections and subharmonic bifurcations. As a demonstration, we present full continuations of periodic double pendulum oscillations from fixed points, showing bifurcations along orbit families and categorizing branches of periodic orbits. In particular, we find periodic orbits where both pendulum masses are never simultaneously at rest, which to our knowledge has been missing in the literature.
I.
INTRODUCTION
Finding and analyzing periodic orbits plays a key role in understanding dynamical systems [1–3]. Periodic orbits can not only facilitate calculations of statistical properties of trajectories (e.g., via cycle expansions) [4–6] but also are important due to the inherent simplicity of recurring trajectories. Through both statistical methods and intuition, periodic orbits offer a window through which we can start to understand the dynamics of otherwise difficult-to-interpret physical systems. Although various numerical methods have been developed to search for periodic orbits [7–10], there does not exist an efficient method to map out the full orbit space. Our contribution is a continuation method that leverages automatic differentiation and Hessian eigendecompositions to efficiently navigate through the space of periodic orbit solutions. We use the double pendulum to demonstrate our method, as shown in Figure 1. Although the double pendulum is one of the most elementary systems in physics, to the best of our knowledge, we are the first to survey bifurcations in the periodic orbit manifold and visualize them on a Poincaré section. This diagram not only shows bifurcation structures, but also reveals more complicated orbits, particularly those without a moment when both masses are at rest. This simple fact, albeit easily understood, has been missing from the literature. We briefly review the two basic numerical methods for finding periodic orbits. The shooting method first extracts close recurrences from integrated initial conditions and then converges a periodic orbit via Newton-Raphson descent or other methods [11]. However, a fundamental property of chaotic systems is the exponential diver-
∗ [email protected] † [email protected] ‡ [email protected]
gence of close trajectories, requiring an increasingly precise initial guess for longer orbits or higher-dimensional systems, even close to fixed points. Converging longer orbits becomes exponentially more difficult, or even impossible due to accumulating numerical precision errors. Variations, such as multiple shooting [12], alleviate the convergence issue partially but not fully, and introduce additional complexities of their own. Another orbit-finding method that avoids convergence issues is the variational method. Instead of varying initial conditions and evolving until obtaining a closed orbit, the initialization is a closed loop that does not necessarily satisfy the equations of motion. The loop is then varied until satisfying the evolution equations, resulting in a periodic orbit. As the loop is adapted to the equations of motion locally, no integration of conditions and subsequent exponential divergences occurs. Variations on this method include choice of loop representation, method of optimization, and method of loop initialization. When introduced by Lan and Cvitanović, the orbit was represented as evenly spaced points in time, with the first derivative given by a four-point approximation [13]. Representations with different spacings, such as a weighting between phase space and time spacing [14], and representations involving Fourier coefficients in time have also been implemented since. Methods of optimization range from Newton descent and other Jacobian-based methods to adjoint-free and other linear algebra methods [15]. In most periodic orbit studies, initializations are made either by perturbing from known fixed points or by extracting close recurrences from a grid of initial conditions. These methods do not incorporate the underlying orbit structure into additional initializations as it is probed, finding solutions first and classifying families of solutions later with human input. Numerical continuation methods have been previously applied in periodic orbit searches, but the majority of attempts have used older integrator methods. Computing the tangent vector in pseudo-arclength continuation has involved integrat-
2
OD1
IB ID2 OD1B
In-Phase
ID1 OD2
OB Up-Down
OD1 Out-Phase
Down-Up
FIG. 1. Converged double pendulum periodic orbits in the Poincaré section with θ̇1 = 0; bifurcating branches intersect and extend off continuous orbit families between fixed points. Fixed points circled in black. In-Phase, Out-Phase, Down-Up, and Up-Down start from fixed points and run to period divergences. Bifurcation branches abbreviated with I and O for base branch, B for branching at same period, D for period-doubling, and numbered for uniqueness.
ing the Jacobian [16], a computationally expensive and numerically unstable task for longer orbits, and there has not been recent development on other propagation techniques. It is our motivation to revisit this continuation problem, given great recent advancement of computational tools, especially the technique of automatic differentiation. The ability to map out the whole periodic orbit spectrum can shed great light on the properties of dynamical systems. The two simplest physical systems demonstrating chaotic behavior and a rich periodic orbit spectrum are the double pendulum and plane-circular restricted three-body problem (PCR3BP). Two Hamiltonian systems each with a four-dimensional state space, both have multiple known fixed points, stable and unstable, and small perturbations from them forming periodic orbits. From periodic orbit family theory [17–20], we expect periodic orbits to be organized into continuous, single-
parameter families. An individual family of orbits cannot bifurcate or arbitrarily conclude; orbit families can only terminate at fixed points, or if the period or a phase space coordinate diverges. Bifurcations in the orbit spectrum are caused by crossings of families at the same and higher multiples of periods, and these crossings also correspond to changes in orbit stability. These connected, one-dimensional families of orbits make dimension-four Hamiltonian systems the perfect application for numerical continuation techniques. While orbits in the PCR3BP have been extensively studied, most work on the double pendulum has focused on systems with external driving and damping forces. These are studied to model various real-world applications, but lead to non-Hamiltonian systems with perturbed steady states and limit cycles [21–24]. For the pure double pendulum, fixed points (with the pendulum bobs vertically up or down) and small angle peri-
3
FIG. 2. An example of a periodic orbit spectrum consisting of multiple orbit families. At period T1 , a crossing orbit family creates a bifurcation in the orbit spectrum, shown in purple. At period T2 , a period-doubling bifurcation occurs due to a crossing of an orbit family with twice the period, shown in dashed blue.
odic orbits are known, as well as basic numerical investigations, such as the Lyapunov exponent and autocorrelation function [25, 26]. Periodic orbits are occasionally found as a byproduct of other explorations, such as auto-discovering conserved quantities [27], but no systematic search has been conducted. More recently, there has been some research on homoclinic and heteroclinic orbits, which has also yielded orbits around unstable fixed points [28]. However, a full numerical survey of double pendulum periodic orbits has not been conducted in the literature. Orbits are theoretically predicted to form continuous families, but extensions of small angle oscillations to higher amplitudes has not yet been attempted. While individual examples of more complicated periodic orbits exist, no methodical investigations of the orbit spectrum and bifurcations leading to non-stationary periodic orbits have been conducted. In this paper, we introduce a novel Hessian-based variational continuation method for periodic orbits. In Section II, we represent a periodic orbit loop under a Fourier parametrization and study Hessian eigenvectors of the variational loss. In Section III, we demonstrate the automatic discovery of small oscillations around unstable fixed points, continuation of orbit families with precise loop initializations, and detection of bifurcations in the orbit spectrum. We present full extensions of symmetric periodic orbits from small oscillations to period divergence, examples of orbit family crossings and spectrum bifurcations, and discuss integrating the Hessian method with other orbit-finding techniques. We summarize our conclusions in Section IV. II. A.
METHOD
Variational Orbit Setup
We consider a general dynamical system parametrized by some state vector z in a phase space M, where time evolution z(t) is determined by the ordinary differential equation: dz = f (z) dt
(1)
A periodic orbit is defined by an initial condition z0 in phase space and a time T > 0 such that z(T ) = z0 . To set up our variational method, we take a closed loop L, defined as a period T and a trajectory z(t) in phase space, with 0 ≤ t < T , z ∈ M, and z(T ) = z(0). At each point z(t) on the trajectory, we can quantify the deviation of the loop trajectory from physical evolution by taking the squared norm of the difference: ℓ(t) =
dz(t) − f (z(t)) dt
2
(2)
We then average over the entire trajectory to obtain a loss function for the loop: ℓ(L) =
1 T
Z T ℓ(t) dt = 0
1 T
Z T 0
dz(t) − f (z(t)) dt
2
dt
(3) This loss function quantifies the overall deviation of the loop from a physical trajectory. Although the loop space is in principle infinite-dimensional, we need to choose a finite-dimensional parameterization. We will parameterize the loop as a Fourier series (as detailed in Section II D), but for now, it suffices to denote the parameters as a vector θ, so ℓ(L) becomes ℓ(θ). Zeroes of ℓ correspond to loops that always locally match physical evolution, and therefore are periodic orbits. As the loss function is always nonnegative, at periodic orbits the loss ℓ and gradients ∇θ ℓ both vanish.
FIG. 3. For an arbitrary trajectory z(t) in phase space, the may not necessarily match the physical loop tangent dz(t) dt evolution f (z(t)). For a periodic orbit, these vectors must match across the entire loop.
B.
Hessian Analysis
Once we have a single periodic orbit, how can we use it to obtain a new nearby orbit? Traditional continuation methods start from an initial condition, then integrate the entire Jacobian matrix [29]. This integration is both numerically unstable and ignores the originally discovered orbit structure by just using an initial condition. By instead considering the entire orbit loop and leveraging automatic differentiation, we can more easily obtain such a perturbation direction by analyzing the Hessian of the variational loss landscape. The Hessian H of the loss function ℓ describes the local curvature of the loss around a parametrized loop θ within a parameter space. For a periodic orbit, l is at a local minimum, so the Hessian H(θ) must be positive
4 semi-definite. The eigenvalues provide information about curvature magnitudes and associated directions. For a zero eigenvalue, the associated eigenvector gives a direction in which perturbing θ maintains zero loss, and therefore maintains a periodic orbit. For a nontrivial orbit (not a fixed point), the full Hessian must have at least one such flat direction, corresponding to shifting the phase of the loop in time. Additional vanishing eigenvalues correspond to additional perturbation directions for the loop, creating a subspace of possible perturbations. The dimension of this nullspace, and therefore the dimensionality of the connected periodic orbit space at θ, is given by the number of zero eigenvalues of the Hessian H(θ). We can also consider a subspace of the loop parameter space, restricting the possible perturbations to the loop. Taking the subspace Hessian then gives only local curvatures corresponding to possible loop perturbations within this subspace. For example, enforcing a phase condition removes time evolution from the possible perturbations, removing the corresponding zero eigenvalue. By restricting trivial or known perturbations from the subspace until there is a unique zero eigenvalue, we can extract the direction of chosen adjacent connected orbits by reading off the remaining unique flat direction. An adjacent orbit initialization can be made by taking a finite step in loop space along the flat direction. The orbit family may not be perfectly linear, but for a sufficiently small step size, the initialization will be close to the orbit family, and can be converged to a periodic orbit on the family. This directional step is the same tangent propagation step used in pseudo-arclength continuation, but made in loop space. Further, as the new initialization is an entire loop, the subsequent periodic orbit optimization benefits from the much larger radius of convergence of the variational method.
C.
is given by the corresponding eigenvector. By sweeping values of T , all possible small oscillations around fixed points, which are starting points for periodic orbit branches, can be systematically discovered. Continuation of periodic orbit branches can be improved by the Hessian method. The larger radius of convergence in loop space allows for larger steps to be taken and longer orbits to be continued, allowing exploration further along orbit families even when integrator methods fail due to precision issues. The eigenvalue direction also gives an orthogonal orbit constraint in loop space, allowing for regular steps in pseudo-arclength to be taken while better reflecting the overall structure of the orbit compared to just an initial condition. The Hessian method also offers a means to detect bifurcations of periodic orbit families. If an orbit family crosses another, there will be two unique flat directions at the crossing point, one along each family. This gives a two-dimensional nullspace at the bifurcation point, and can be detected as a second eigenvalue drops to zero as it is approached. Constraining the subspace to eliminate the known orbit family at the bifurcation point, the direction of the crossing orbit family can be determined and used to obtain an initialization for it. Subharmonic bifurcations can be detected with a modification to the loop. For a period-doubling bifurcation, a new loop can be constructed with two windings of the original orbit. This doubly-wound loop can still be propagated along the original family, though with all components at doubled (even) frequencies. However, subharmonic (odd frequency) components can now also be perturbed, allowing a second flat direction corresponding to the period-multiplied bifurcation. By transforming a loop representation to a multiply-winded version and testing for zero eigenvalues, other subharmonic bifurcations in the orbit spectrum can also be detected and initialized.
Hessian Applications
The method of Hessian analysis allows brute force methods for finding periodic orbits, such as grid search, to be replaced by a more methodical approach leveraging the underlying geometric structure of periodic solutions. Fixed points are (trivial) periodic orbits for any period T , and the starting point for periodic orbit searches. While small amplitude periodic perturbations around stable fixed points can be found by linearization, orbits around unstable fixed points require manual converging, which can be increasingly difficult for chaotic system dynamics and may collapse back down to the fixed point. The Hessian method allows for precise initializations of small amplitude periodic perturbations, even in highly chaotic systems. Loops can be initialized to the fixed point for various periods T , and the existence of periodic perturbations can be tested by computing the Hessian. A drop of the minimal eigenvalue to zero indicates a period T with a periodic perturbation, and this perturbation
D.
Implementation
We consider the ideal double pendulum with equal point masses and equal-length massless arms, setting all parameters to 1 (m1 = m2 = l1 = l2 = g = 1). A system state is specified by arm angles θ1 , θ2 and angular velocities θ̇1 , θ̇2 . A loop L is represented by a period T and a closed trajectory z(t). We parametrize z(t) with a Fourier decomposition in time with a finite maximum frequency cutoff K: z(t) = a0 +
K X k=1
ak cos
2kπt T
+ bk sin
2kπt T
,
where 0 ≤ t ≤ T. (4) We implement computations in PyTorch [30] with float64 precision, leveraging autograd for Hessian computation. As parameters and gradients vary by orders of magnitude between low and high frequency components,
5
FIG. 4. On left, a branch of orbits L emerging from a fixed point at period T0 (above) is found by sweeping period T at the fixed point, and detecting a drop in the minimal eigenvalue λ (below). At center, the increased radius of convergence (shaded) of the variational method allows for exploration further along orbit families, with an orthogonal loop space constraint. On right, a crossing of orbit families at L′ (above) is found by detecting a drop in the off-branch eigenvalue λ (below).
we use the Rprop algorithm [31], which adapts step sizes per parameter based solely on the signs of the gradients. We optimize to a minimum of √ the integrated error, and use a convergence condition of ℓ < 10−10 . Higher energy, longer period orbits require higher frequency components to be suitably represented by the Fourier parametrization. As we propagate along a family of orbits, we monitor the final converged value of ℓ, increasing the frequency cutoff K if the convergence condition is not reached. Symmetries in the orbit are also exploited to reduce the number of parameters optimized. For spatially symmetric orbits of the double pendulum, we evolve (require gradient on) only cosine positional components and sine momentum components, enforcing the symmetry constraint and halving the parameter count. This constraint both enforces a phase condition and avoids accidentally propagating to perturbed orbits.
III. A.
RESULTS
Fixed Point Analysis
We start by using a sweep of period T to analyze the fixed points of the double pendulum. We choose a step size of ∆T = 0.01, which controls the precision of the detected period T and the accuracy of the corresponding orbit initialization. As the loop representations of fixed points can be immediately initialized (all non-constant Fourier components are zero), the only computation necessary is computing the Hessian, and the entire sweep takes a few minutes to run for K = 16 Fourier terms. Plotting the minimum Hessian eigenvalue against pe-
riod T gives sharp drops to zero at T ≈ 3.39 and T ≈ 8.21 for the stable fixed point (both masses down), and T ≈ 5.28 around both unstable saddle points (one mass down, one mass up). These T values match the theoretically known periods of oscillation for normal modes. No eigenvalue drop is noticed for the unstable fixed point (both masses up), as small-amplitude periodic orbits do not exist around it.
The corresponding eigenvectors give precise initializations for periodic orbits. Taking a step of size ∆θ = 0.5◦ in phase space, we can compute an initial condition from the loop L at time t = 0, then integrate for a time T given by our detected period. Comparing the final result of the integration with the starting point, we obtain relative deviations from the perfect periodicity of within 2% compared to the phase space step size (Table I). Our method is robust to local stability, providing precise initial conditions around both stable and unstable fixed points without the need for a computationally expensive search.
If we take the Hessian of all free parameters for the stable fixed point, we also obtain an additional eigenvalue dip at T ≈ 6.78. This corresponds to the same normal mode at the T ≈ 3.39 eigenvalue zero, but with two oscillations in a single loop instead of one. Constraining the free parameters of the Hessian to require antisymmetry in the two halves of the orbit (only odd frequency components) removes this period-doubled version of the oscillation from the eigenvalue zeroes (Figure 5).
6 TABLE I. For each of the detected normal modes, a step in loop space corresponding to a phase space deviation of 0.5◦ is taken, an initial condition is obtained by evaluating the Fourier parametrization at t = 0, and this phase space condition is integrated for the detected period T . The deviation of the integration result from the original condition is small compared to the step size of 0.5◦ , implying a precise initialization corresponding to a periodic orbit. Normal Mode
Integrated Deviation (degrees)
Down-Down, In-Phase Down-Down, Out-Phase Down-Up Up-Down
1.76 × 10−4 1.02 × 10−3 0.0982 0.0976
FIG. 5. Minimum eigenvalue for the period T with both masses down (θ1 = 0, θ2 = 0); eigenvalue minimums give out of phase and in-phase normal modes and correct corresponding periods T .
B.
Symmetric Orbit Propagation
From each normal mode, we propagate a family of symmetric orbits by repeatedly converging an initialization to a periodic orbit, taking the Hessian, and taking a loop space step in the zero-eigenvalue direction. During convergence, we constrain the training dynamics perpendicular to the step direction, ensuring that propagation proceeds along the orbit family. We use a step size of 0.5◦ for all orbit families. We constrain parameters both to enforce a phase condition and to limit to odd frequency components, each optimization halving the number of trainable parameters. We start with a frequency cutoff of K = 16 from the fixed point analysis,√doubling the cutoff K if the convergence condition of ℓ < 10−10 is not reached. The frequency cutoff rose to K = 1024 for the longest converged orbits; converging an orbit took from seconds for K = 16 to a few hours for K = 1024, and memory requirements were primarily due to Hessian size. From small oscillations, we were able to successfully propagate all four orbit families to orbits approaching a period divergence, with masses within 5◦ of vertical. Through the entire energy spectrum, loop loss remained below the convergence condition, while once-
around phase space integrator error scaled exponentially with orbit period T (Figure 6). For orbits approaching vertical, the error after integrating once around diverges, despite the precision of the initial condition obtained from the loop. If using scipy.integrate at double precision, attempting to use an integrator method to converge longer period orbits will fail, as accumulating error from integration destroys the orbit. The variational method is robust to orbit period and chaotic phase space dynamics, demonstrating an advantage in finding orbits that are otherwise impossible to compute without extended-precision numerics.
FIG. 6. Loop and integrator losses for in-phase oscillations; variational convergence succeeds even for longer periods T as masses approach vertical, while integrator error diverges.
C.
Bifurcation Detection
To detect crossing orbit families, we look for additional zero eigenvalues along an orbit family, after removing the zeroes corresponding to time evolution and branch direction; we call these “off-branch eigenvalues”. Plotting the minimum off-branch eigenvalue along an orbit family shows sharp drops to zero, indicating discrete crossing points as expected. To test for period-doubling bifurcations, we first generate a representation of a twice-wound loop by doubling the frequencies of all components and the period. We fill zeroes into the odd frequency components and take the Hessian, again extracting the minimum off-branch eigenvalue. In addition to the bifurcation points of the original loop, additional eigenvalue minimums corresponding to subharmonic perturbations are also seen (Figure 7). Similar to fixed points, eigenvectors give loop perturbations that will maintain zero loss and a periodic orbit. As the dimensionality of the null subspace is greater, an orthogonalization is theoretically required to separate the perturbation subspace from the original orbit branch. In practice, a difference in numerical precision separated the orbit family and perturbation eigenvalues.
7 turbation. By building up a collection of branching orbit families, we can trace a path to each orbit through the bifurcation structure, starting from fixed points and normal modes, then following orbit families and crossings. Motions can be systematically categorized in terms of base and branching families, and complicated motions can be intuitively broken down into sums of simpler oscillations.
FIG. 7. Minimum off-branch eigenvalue for out-of-phase oscillations; additional eigenvalue minimums give same period and period-doubling bifurcations leading to new orbit families.
D.
Interpreting Orbits
Eigenvector directions at bifurcation points lead to branching orbit families and additional connected orbits. By discovering continuous orbit families instead of individual disconnected orbits, we can classify orbits in terms of the families that they lie on, and the orbit families themselves can be classified based on their structure. Along orbit families, higher energy periodic solutions can be understood as extensions of lower energy orbits. In the double pendulum, branches from fixed points start off with small oscillations, continuing to orbits where masses are very close to vertical. These higher energy periodic solutions can be interpreted as extensions of lowenergy normal modes, and the shape of the oscillations also reflects this continuation (Figure 8).
FIG. 8. Out-of-phase oscillations at various energies along the continuous family; higher energy oscillations smoothly extend and deform the shapes of the low-energy normal mode.
Through intersections of orbit families, we can also interpret motions that are not continuations of normal modes. If we consider a bifurcation point and the corresponding orbit, the eigenvector direction corresponds to a small perturbation leading to the branching orbit family. Other orbits on this branching family are continuations of this perturbation. We can understand this orbit family as a combination of the base orbit, and an oscillation perturbing from it (Figure 9). In the case of period-doubling bifurcations, multiple oscillations of the base orbit may correspond to a single period of the per-
FIG. 9. Bifurcating orbits are created by combinations of perturbations at same and higher multiples of period, which can be interpreted by following propagated branches.
E.
Extending Bifurcation Branches
From each bifurcation point, we propagate orbits in the directions of the off-branch perturbations. We maintain the constraint on training dynamics, but relax the constraints on symmetry and phase. When close to bifurcation points, unconstrained loss landscape geometry favors convergence towards the simpler orbit branch, instead of the bifurcation branch. The training dynamics constraint is especially important to guide initial convergence when starting to propagate a bifurcation, in order to maintain the direction of the initial perturbation and not jump back to the already known orbit family. Due to the gradient signal being weaker in our constrained direction, final convergence is also notably slower in the vicinity of bifurcation points, despite an exponentially decreasing loss function. To accelerate a slow rate of exponential convergence, linear extrapolation was employed during the training process. The direction of convergence was estimated based on previous epochs and a linear search was conducted along the estimated direction, with the new initialization made at the training loss minimum along the line. Doing so sped up convergence significantly, with significant drops in losses seen at each extrapolation step (Figure 10).
F.
Phase Space Exploration
Many examples of successfully extended branches were obtained across the phase space, with a majority of
8 extended to search for additional orbit families up to numerical constraints.
IV.
FIG. 10. Comparison of loop and integrator losses during training with and without extrapolation; linear search greatly speeds up slow exponential convergence.
branches approaching period divergences or another termination condition. A variety of rich behaviors in the periodic orbit spectrum are revealed, such as an orbit family forming a closed loop and an additional bifurcation point. A plot of all converged periodic orbits is shown in Figure 1, with branches labeled by base branch and type of bifurcation. The structure of our converged branches and bifurcations matches orbit family theory. Both in-phase doublings show an orbit family passing continuously through the bifurcation point, with the first (ID1) looping back in on itself and the second (ID2) proceeding towards a period divergence. The first out-phase doubling (OD1) crosses the out-of-phase family (Out-Phase); the points on each side represent the same orbits (up to time reversal), but show initial conditions at different phases. The out-phase bifurcation branch (OB) continues past the bifurcation point to negative θ̇2 , following the same trace of initial angles but with opposite initial angular velocity. The additional bifurcation point and branch (OD1B) is found along an out-of-phase period-doubling (OD1). Approaching a bifurcation point, the Hessian nullspace increases in dimension; this causes the computed zeroeigenvalue directions to vary depending on numerical precision. This bifurcation was found from two separate runs initially which followed the same orbit family before diverging. The in-phase (IB, very short) and out-phase bifurcation branches (OB) cut off mid-propagation due to convergence failures. This is also likely due to a bifurcation point, but the convergence dynamics did not lead to a successful continued propagation along either orbit family. This failure mode is not present in the original branches, due to the phase and symmetry conditions. As the orbit branches are otherwise smooth, a constraint based on the recent orbit family direction could be used to continue propagation in future work. We otherwise did not explicitly search for bifurcation points along non-symmetric orbit families. Theoretically, the periodic orbit structure is infinite; by continuing to compute off-branch eigenvalues and propagating new branches along perturbations, our method can be
CONCLUSION
We have introduced a loss landscape perspective to searching for periodic orbits in dynamical systems. By probing loss curvatures using the Hessian, we can systematically explore the structure of orbit families, determining connected classes of orbits and bifurcations. We obtain a numerical continuation method in loop space, with variational initialization and convergence properties. By combining the Rprop optimizer and linear extrapolation, we are able to converge orbits of varying complexity across the phase space of the ideal double pendulum. To our knowledge, this is the first survey of orbits that both discovers orbits off normal-mode branches and systematically classifies them into continuous families. The success of this optimization routine also demonstrates the applicability of machine learning optimizers to problems in physical systems. While Rprop quickly gives satisfactory results without much tuning, taking advantage of different hyperparameters and algorithms may further improve performance. On an implementation level, constraints such as the propagation direction can be easily expressed using developed tools and methods from machine learning such as parametrizations. A limitation of using the Rprop optimizer is that it is a gradient-based algorithm, which scales unfavorably in complexity for higher-dimensional systems compared to recent Jacobian-free methods. However, the Hessian probing of loss landscapes, which takes an orbit and finds eigenvector directions and new loop initializations, can be decoupled from methods to converge the initialization to a new orbit. Only a single Hessian evaluation is needed to determine the local loss landscape structure, and our propagation method is compatible with other variational methods in the literature. If computing the full Hessian is too expensive, approximations can be made by using reduced parameter loop representations or Hessian approximation algorithms. Instead of finding individual, disconnected orbits, methodically searching orbits allows for solutions to be understood as part of a orbit family structure, and for an existing known structure to inform the search and discovery of new orbit candidates. By following branches on a bifurcation diagram, more complicated motions can be interpreted as combinations of simple oscillations, giving a more interpretable method for analyzing periodic solutions. Promising future directions include both the continued application of machine learning methods to periodic orbits, and applying our method of systematically propagating through the orbit spectrum to a wide range of other dynamical systems.
9 ACKNOWLEDGMENTS
Z.L. and M.T. were supported by IAIFI through NSF grant PHY-2019786.
[1] M. C. Gutzwiller, Periodic orbits, in Chaos in Classical and Quantum Mechanics (Springer New York, New York, NY, 1990) pp. 75–86. [2] V. Gelfreich, Proceedings of the National Academy of Sciences 99, 13975–13979 (2002). [3] L. Sbano, Periodic orbits of hamiltonian systems, in Mathematics of Complexity and Dynamical Systems, edited by R. A. Meyers (Springer New York, New York, NY, 2011) pp. 1212–1236. [4] R. Artuso, E. Aurell, and P. Cvitanovic, Nonlinearity 3, 325–359 (1990). [5] Y. Lan, Communications in Nonlinear Science and Numerical Simulation 15, 502–526 (2010). [6] D. Lasagna, Phys. Rev. E 102, 052220 (2020). [7] Y. Saiki, Nonlinear Processes in Geophysics 14, 615 (2007). [8] S. C. Farantos, Journal of Molecular Structure: THEOCHEM 341, 91 (1995). [9] L. Drossos, O. Ragos, M. N. Vrahatis, and T. Bountis, Phys. Rev. E 53, 1206 (1996). [10] D. Pingel, P. Schmelcher, and F. K. Diakonos, Phys. Rev. E 64, 026214 (2001). [11] B. Mestel and I. Percival, Physica D: Nonlinear Phenomena 24, 172 (1987). [12] J. SÁNCHEZ and M. NET, International Journal of Bifurcation and Chaos 20, 43–61 (2010). [13] Y. Lan and P. Cvitanović, Phys. Rev. E 69, 016217 (2004). [14] D. Wang, P. Wang, and Y. Lan, Phys. Rev. E 98, 042204 (2018). [15] S. Azimi, O. Ashtari, and T. M. Schneider, Phys. Rev. E 105, 014217 (2022). [16] C. Wulff and A. Schebesch, SIAM Journal on Applied Dynamical Systems 5, 435–475 (2006). [17] M. Hénon, Generating Families in the Restricted Three-
Body Problem (Springer Berlin Heidelberg, 1997). [18] A. Wintner, Mathematische Zeitschrift 34, 321–349 (1932). [19] G. D. Birkhoff, Annali della Scuola Normale Superiore di Pisa - Classe di Scienze 5, 9 (1936). [20] M. Baranger, K. Davies, and J. Mahoney, Annals of Physics 186, 95 (1988). [21] P. Yu and Q. Bi, Journal of Sound and Vibration 217, 691–736 (1998). [22] K. Yagasaki, Nonlinear Dynamics 9, 391–417 (1996). [23] Z. Liu, X. Rao, J. Gao, and S. Ding, Chaos, Solitons & Fractals 177, 114254 (2023). [24] H. Williams, Forces in Mechanics 10, 100164 (2023). [25] A. M. Calvão and T. J. P. Penna, European Journal of Physics 36, 045018 (2015). [26] J. Jiménez-López and V. J. Garcı́a-Garrido, International Journal of Bifurcation and Chaos 34, 10.1142/s0218127424502018 (2024). [27] Z. Liu and M. Tegmark, Phys. Rev. Lett. 126, 180604 (2021). [28] K. Kaheman, J. J. Bramburger, J. N. Kutz, and S. L. Brunton, Nonlinear Dynamics 111, 7199–7233 (2023). [29] M. Net and J. Sánchez, SIAM Journal on Applied Dynamical Systems 14, 674–698 (2015). [30] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in Advances in Neural Information Processing Systems 32 , edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035. [31] M. Riedmiller and H. Braun, in IEEE international conference on neural networks (IEEE, 1993) pp. 586–591.