Computation of minimal periods for ordinary differential equations (original) (raw)

(1Division of Mathematics, University of Dundee, Dundee, DD1 4HN, United Kingdom
jparker002@dundee.ac.uk)

Abstract

We consider the problem of finding the shortest possible period for an exactly periodic solution to some given autonomous ordinary differential equation. We show that, given a pair of Lyapunov-like observable functions defined over state space of the corresponding dynamical system and satisfying a certain pointwise inequality, we can obtain a global lower bound for such periods. We give a method valid for the case of bounding the period of only those solutions which are invariant under a symmetry transformation, as well as bounds for general periodic orbits. If the governing equations are polynomial in the state variables, we can use semidefinite programming to find such auxiliary functions computationally, and thus compute lower bounds which can be rigorously validated using rational arithmetic. We apply our method to the Lorenz and Hénon-Heiles systems. For both systems we are able to give validated bounds which are sharp to several decimal places.

1 Introduction

In a sufficiently regular vector field, it is natural to expect a lower bound on the period of closed orbits: an orbit must return to its starting point, while the local speed and curvature of the flow are limited. This intuition can be formalized by appealing to Fenchel’s theorem, which states that the total curvature of any closed curve is at least 2​π2\pi.

Yorke [25] used Fenchel’s theorem to prove an important result: for a real, finite dimensional system governed by an ordinary differential equation x˙=f​(x)\dot{x}=f(x) with a Lipschitz constant LL such that |f​(x)−f​(y)|≤L​|x−y||f(x)-f(y)|\leq L|x-y|, the period of a limit cycle is no less than 2​πL\frac{2\pi}{L}. Various different extensions and relaxations of this result, for example to general Hilbert and Banach spaces, have been found [2, 23, 27, 3, 26, 19, 15], and the factor of 2​π2\pi is optimal for general Hilbert spaces, in the sense that it is possible to find vector fields in which this bound is attained. However, for any particular system of interest, this bound is not expected to be sharp, i.e. there is not usually a limit cycle of period 2​πL\frac{2\pi}{L}, and the constant can be improved given certain restrictions on the field [15]. One notable application of these ideas is that of Kukavica [16], who showed that there is a finite lower bound on the period of solutions of the incompressible 3D Navier-Stokes equations.

The present work was inspired by the observation that the total curvature of a periodic orbit as it appears in Yorke’s application of Fenchel’s theorem can be reformulated as a time-average over a function defined over the state space of the dynamical system. Recent developments in other areas of applied dynamical systems have seen sum-of-squares optimization, a numerical method based on semidefinite programming (SDP), applied to find and prove bounds on the time averages of quantities in soltuions to ordinary differential equations (ODEs) [8, 5]. Thus applying these bounding methods to Fenchel’s theorem gives a bound on periods, which is a special case of the more general approach we describe in section 2.

We are motivated by applications of differential equations in mathematical models of real-world problems. The question we aim to answer is, given a particular ODE, what is the shortest possible period for a periodic solution? The computational methods discussed in this paper, like other related methods based on semidefinite programming, are described for ODE systems whose right-hand side is a polynomial function of the state variables. The two example systems we study are indeed of this form, but in practical applications of ODEs – from biological models, classical mechanics, planetary dynamics and so on – it is not typical to have purely polynomial ODEs. In cases where the right-hand side contains, for example, trigonometic functions, rational functions, exponentials or square roots, a straightforward change of variables can be made to transform the system into a polynomial one. Even with so-called ‘differentially transcendental’ functions where this is not possible, locally we can approximate the dynamics to arbitrary precision using polynomials. Therefore, the methods we present here are broadly applicable.

In section 2, we give a method for finding lower bounds of orbit period via the existence of state space observables satisfying a specific pointwise inequality. Unlike previous work, this does not rely on a Lipschitz constant, though the existence of a Lipschitz constant gives a special case of our result which is generally suboptimal. In sections 3 and 4 the method is applied computationally to the Lorenz and Hénon-Heiles systems respectively, in each case finding bounds which are very close to the period of known solutions. Code is provided which rigorously validates these bounds. Section 5 gives concluding remarks.

2 A Lyapunov-like method for lower bounds on periods

Firstly we recall the following form of Wirtinger’s inequality for functions (a Poincaré inequality with optimal constant on an interval) [13, section 7.7]:

Lemma 2.1.

For T>0T>0 and f∈C1​([0,T],ℝm)f\in C^{1}([0,T],\mathbb{R}^{m}), if ∫0Tf​(t)​dt=0\int_{0}^{T}f(t)\mathrm{d}t=0 and f​(T)=f​(0)f(T)=f(0) then

| ∫0T|f​(t)|2​dt≤(T2​π)2​∫0T|f˙​(t)|2​dt.\int_{0}^{T}|f(t)|^{2}\mathrm{d}t\leq\left(\frac{T}{2\pi}\right)^{2}\int_{0}^{T}\left|\dot{f}(t)\right|^{2}\mathrm{d}t. | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |

The notation C1​(A,B)C^{1}(A,B) means the set of continuously differentiable functions from AA to BB, and |⋅|\left|\cdot\right| is the usual Euclidean norm. The inequality becomes an equality if and only if all components of ff are purely sinusoidal with period TT.

Our main result, which is inspired by previous work on bounding time averages in continuous-time dynamical systems via auxiliary functions VV [8], relies on this existence of some function ϕ\phi defined over state space which vanishes on average over periodic trajectories of interest. The result is:

Proposition 2.2.

Consider a system defined over a domain Ω⊂ℝn\Omega\subset\mathbb{R}^{n} and governed by

d​xd​t=f​(x)\frac{\mathrm{d}x}{\mathrm{d}t}=f(x) (1)

for f∈C1​(Ω,ℝn)f\in C^{1}(\Omega,\mathbb{R}^{n}). Suppose there exists C>0C>0, V∈C1​(Ω,ℝ)V\in C^{1}(\Omega,\mathbb{R}) and ϕ∈C1​(Ω,ℝm)\phi\in C^{1}(\Omega,\mathbb{R}^{m}) such that

| C​|ϕ​(x)|2−|ℒf​ϕ​(x)|2+ℒf​V​(x)≥0C\left|\phi(x)\right|^{2}-\left|\mathcal{L}_{f}\phi(x)\right|^{2}+\mathcal{L}_{f}V(x)\geq 0 | (2) | | ---------------------------------------------------------------------------------------------------------------------------------------- | --- |

for all x∈Ωx\in\Omega, where ℒf​ϕ​(x)=D​ϕ​(x)​f​(x)\mathcal{L}_{f}\phi(x)=D\phi(x)\,f(x) is the Lie derivative at xx of the observable ϕ\phi with respect to the flow. Then for any solution of 1 with x​(t+T)=x​(T)x(t+T)=x(T) for some T>0T>0 and ∫0Tϕ​(x​(t))​dt=0\int_{0}^{T}\phi(x(t))\mathrm{d}t=0, either ϕ​(x​(t))=0\phi(x(t))=0 for all tt or T≥2​π/CT\geq 2\pi/\sqrt{C}.

Proof.

Suppose we have CC and VV satisfying the above criteria. Consider a periodic orbit x​(t)x(t) with period TT on which ϕ\phi vanishes on average, and write ϕ​(t)=ϕ​(x​(t))\phi(t)=\phi(x(t)), V​(t)=V​(x​(t))V(t)=V(x(t)). Then for all tt,

| C​|ϕ​(t)|2−|ϕ˙​(t)|2+V˙​(t)≥0,C\left|\phi(t)\right|^{2}-\left|\dot{\phi}(t)\right|^{2}+\dot{V}(t)\geq 0, | (3) | | ------------------------------------------------------------------------------------------------------------------ | --- |

and so integrating over one period,

| C​∫0T|ϕ​(t)|2​dt−∫0T|ϕ˙​(t)|2​dt≥0.C\int_{0}^{T}\left|\phi(t)\right|^{2}\mathrm{d}t-\int_{0}^{T}\left|\dot{\phi}(t)\right|^{2}\mathrm{d}t\geq 0. | (4) | | --------------------------------------------------------------------------------------------------------------------------------------------------------------- | --- |

Assuming that ϕ​(t)\phi(t) does not vanish everywhere,

| ∫0T|ϕ˙​(t)|2​dt∫0T|ϕ​(t)|2​dt≤C.\frac{\int_{0}^{T}\left|\dot{\phi}(t)\right|^{2}\mathrm{d}t}{\int_{0}^{T}\left|\phi(t)\right|^{2}\mathrm{d}t}\leq C. | (5) | | -------------------------------------------------------------------------------------------------------------------------------------------------------------------- | --- |

Lemma 2.1 tells us that

| (2​πT)2≤∫0T|ϕ˙​(t)|2​dt∫0T|ϕ​(t)|2​dt,\left(\frac{2\pi}{T}\right)^{2}\leq\frac{\int_{0}^{T}\left|\dot{\phi}(t)\right|^{2}\mathrm{d}t}{\int_{0}^{T}\left|\phi(t)\right|^{2}\mathrm{d}t}, | (6) | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | --- |

and the result follows. ∎

To derive a bound on the period of periodic solutions using this result, it is necessary to find a ϕ\phi with the property that ϕ\phi does not vanish identically over any periodic orbit, but whose average does vanish for each orbit. The obvious choice is the state space velocity itself, ϕ=f\phi=f. In this case, if ff is globally Lipschitz with constant LL, then it follows that |ℒf​f|≤‖D​f‖​|f|≤L​|f||\mathcal{L}_{f}f|\leq\left\|Df\right\||f|\leq L|f|, and so we can take V=0V=0 and C=L2C=L^{2} to satisfy 2, giving T≥2​π/LT\geq 2\pi/L for any periodic orbit. This is exactly the classical result of Yorke [25]. The additional freedom to choose V≠0V\neq 0 potentially gives sharper bounds and allows for the case where ff is not globally Lipschitz. In fact, given any polynomial ϕ\phi for a system with polynomial ff, we can directly optimize the bound with the following sum-of-squares program

| max(x,T)∈𝒜(2​πT)2\displaystyle\max_{(x,T)\in\mathcal{A}}\left(\frac{2\pi}{T}\right)^{2} | ≤infV∈C1​(Ω,ℝ)sup(x,T)∈𝒜∫0T(|ℒf​ϕ​(x​(t))|2−ℒf​V​(x​(t)))​𝑑t∫0T|ϕ​(x​(t))|2​𝑑t\displaystyle\leq{\inf_{V\in C^{1}(\Omega,\mathbb{R})}\sup_{(x,T)\in\mathcal{A}}\frac{\int_{0}^{T}\left(|\mathcal{L}_{f}\phi(x(t))|^{2}-\mathcal{L}_{f}V(x(t))\right)dt}{\int_{0}^{T}|\phi(x(t))|^{2}dt}} | (7) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | --- | | ≤infV∈C1​(Ω,ℝ)C>0C​ s.t. ​C​|ϕ​(x)|2−|ℒf​ϕ​(x)|2+ℒf​V​(x)≥0​∀x∈Ω\displaystyle\leq\inf_{\begin{subarray}{c}V\in C^{1}(\Omega,\mathbb{R})\\ C>0\end{subarray}}C\text{ s.t. }C\left|\phi(x)\right|^{2}-\left|\mathcal{L}_{f}\phi(x)\right|^{2}+\mathcal{L}_{f}V(x)\geq 0\,\forall x\in\Omega | | | | ≤infV​ polynomialC>0C​ s.t. ​C​|ϕ|2−|ℒf​ϕ|2+ℒf​V∈S​O​S\displaystyle\leq\inf_{\begin{subarray}{c}V\text{ polynomial}\\ C>0\end{subarray}}C\text{ s.t. }C\left|\phi\right|^{2}-\left|\mathcal{L}_{f}\phi\right|^{2}+\mathcal{L}_{f}V\in SOS | | |

where the notation g∈S​O​Sg\in SOS means that the polynomial function gg can be expressed as a sum of squares of a polynomials. Here

| 𝒜={(x,T)∈C1​(ℝ,Ω)×(0,∞):x˙​(t)=f​(x​(t))∀t,x​(t+T)=x​(t)∀t,∫0Tϕ​(x​(t))​dt=0,∫0T|ϕ​(x​(t))|2​dt≠0}\mathcal{A}=\left\{(x,T)\in C^{1}(\mathbb{R},\Omega)\times(0,\infty):\begin{aligned} &\dot{x}(t)=f(x(t))\quad\forall t,\\ &x(t+T)=x(t)\quad\forall t,\\ &\int_{0}^{T}\phi(x(t))\,\mathrm{d}t=0,\\ &\int_{0}^{T}\lvert\phi(x(t))\rvert^{2}\,\mathrm{d}t\neq 0\end{aligned}\right\} | (8) | | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | --- |

is the set of all periodic orbits of interest – which is all nonstationary periodic orbits of the system in the case that ϕ=f\phi=f.

However, this is not guaranteed to give a sharp bound on the period for any given ϕ\phi. Instead we would like to simultaneously optimize the choice of ϕ\phi and VV. Suppose that we have a ‘library’ of mm linearly independent functions gi∈C0​(Ω,ℝ)g_{i}\in C^{0}(\Omega,\mathbb{R}) which each have the desired property of integrating to zero over periodic orbits of interest, and furthermore that at least one of the gig_{i} is non-zero over each periodic orbit of interest. One option is gi=ℒf​pig_{i}=\mathcal{L}_{f}p_{i} for non-constant polynomials pip_{i}, which generalizes the above choice ϕ=f\phi=f corresponding to gi=ℒf​xig_{i}=\mathcal{L}_{f}x_{i}. Other possibilities may arise from symmetries of the system, as discussed in section 3.

Then we let ϕ=U​g\phi=Ug where UU is an m×mm\times m invertible matrix and gg is the vector of gig_{i}. This, together with the conditions given for gg, implies that ∫0T|ϕ|2​𝑑t≠0\int_{0}^{T}|\phi|^{2}dt\neq 0 for periodic orbits of interest. Furthermore, ϕ​(x)=U​g​(x)\phi(x)=Ug(x) for some invertible UU if and only if there exists positive definite QQ such that |ϕ​(x)|2=g​(x)𝖳​Q​g​(x)|\phi(x)|^{2}=g(x)^{\mathsf{T}}Qg(x), via Cholesky decomposition. This implies that |ℒf​ϕ​(x)|2=(ℒf​g​(x))𝖳​Q​(ℒf​g​(x))|\mathcal{L}_{f}\phi(x)|^{2}=\left(\mathcal{L}_{f}g(x)\right)^{\mathsf{T}}Q\left(\mathcal{L}_{f}g(x)\right), which leads directly to:

Proposition 2.3.

For a system d​xd​t=f​(x)\frac{\mathrm{d}x}{\mathrm{d}t}=f(x) with ff polynomial, and given g∈C0​(Ω,ℝm)g\in C^{0}(\Omega,\mathbb{R}^{m}) with linearly independent polynomial components, consider the set of periodic orbits

| 𝒜={(x,T)∈C1​(ℝ,Ω)×(0,∞):x˙​(t)=f​(x​(t))∀t,x​(t+T)=x​(t)∀t,∫0Tg​(x​(t))​dt=0,∫0T|g​(x​(t))|2​dt≠0}.\mathcal{A}=\left\{(x,T)\in C^{1}(\mathbb{R},\Omega)\times(0,\infty):\begin{aligned} &\dot{x}(t)=f(x(t))\quad\forall t,\\ &x(t+T)=x(t)\quad\forall t,\\ &\int_{0}^{T}g(x(t))\,\mathrm{d}t=0,\\ &\int_{0}^{T}\lvert g(x(t))\rvert^{2}\,\mathrm{d}t\neq 0\end{aligned}\right\}. | (9) | | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ | --- |

The periods of these orbits are bounded below via

max(x,T)∈𝒜(2​πT)2≤infQ​ positive definiteV​ polynomialC>0C s.t. Cg𝖳Qg−(ℒfg)𝖳Q(ℒfg)+ℒfV∈SOS.\max_{(x,T)\in\mathcal{A}}\left(\frac{2\pi}{T}\right)^{2}\leq\inf_{\begin{subarray}{c}Q\text{ positive definite}\\ V\text{ polynomial}\\ C>0\end{subarray}}C\text{ s.t. }Cg^{\mathsf{T}}Qg-\left(\mathcal{L}_{f}g\right)^{\mathsf{T}}Q\left(\mathcal{L}_{f}g\right)+\mathcal{L}_{f}V\in SOS. (10)

In practice the SOS condition is enforced via another positive (semi-)definite matrix, and the polynomial VV is chosen from some finite basis – typically a basis of monomials. Therefore the whole constraint for the infimum is in written terms of unknown semidefinite matrices, the coefficients of VV and the value CC. By fixing a potential value of CC, this constraint becomes a tractable feasibility problem (equivalently, a linear matrix inequality problem) which can be solved directly via an SDP solver. The value of CC is then reduced until the problem becomes infeasible. Note that an additional constraint must be imposed to fix the scaling; we enforce that tr​Q=1\mathrm{tr}\,Q=1.

3 Minimal periods for symmetric trajectories in the Lorenz system

The celebrated Lorenz [17] system is governed by the ordinary differential equations

d​x1d​t=σ​(x2−x1),d​x2d​t=x1​(ρ−x3)−x2,d​x3d​t=x1​x2−β​x3.\displaystyle\begin{split}\frac{\mathrm{d}x_{1}}{\mathrm{d}t}&=\sigma(x_{2}-x_{1}),\\ \frac{\mathrm{d}x_{2}}{\mathrm{d}t}&=x_{1}(\rho-x_{3})-x_{2},\\ \frac{\mathrm{d}x_{3}}{\mathrm{d}t}&=x_{1}x_{2}-\beta x_{3}.\end{split} (11)

The right-hand side of these equations is equivariant under the transformation x↦(−x1,−x2,x3)x\mapsto(-x_{1},-x_{2},x_{3}), with the result that every periodic orbit is either (a) invariant under this transformation, or (b) has an equivalent orbit found by applying the transformation. We focus on the former case, attempting to find lower bounds for the period of symmetric orbits. At the standard parameters of σ=10\sigma=10, β=8/3\beta=8/3 and ρ=28\rho=28, the shortest known orbit is a symmetric orbit, called ‘LR’ because of its symbolic dynamics, and has a rigorously validated period T≈1.5587T\approx 1.5587 [1]. LR is depicted in figure 1.

Refer to caption

Figure 1: The LR orbit, the shortest known periodic solution in the Lorenz system at the standard parameters, with period T≈1.5587T\approx 1.5587.

3.1 Analytical result

The right-hand side of 11 is not globally Lipschitz so we cannot easily obtain a lower bound for periods using the method of Yorke [25]. Compact globally attracting sets are known [7], so it would be possible to compute a Lipschitz constant on these, but our new method proposition 2.3 can easily and analytically give us a lower bound, by considering a very simple choice of gg.

Observe that for any orbit which is invariant under the symmetry, ∫0Tx1​(t)​dt=0\int_{0}^{T}x_{1}(t)\,\mathrm{d}t=0. Furthermore, if x1​(t)=0x_{1}(t)=0 for all tt then it is straightforward to see that the only possible solution is the fixed point x=0x=0. That is, with g=x1g=x_{1}, proposition 2.3 concerns

| 𝒜={(x,T)∈C1​(ℝ,Ω)×(0,∞):x˙​(t)=f​(x​(t))∀t,x​(t+T)=x​(t)∀t,∫0Tx1​(t)​dt=0,∫0T|x1​(t)|2​dt≠0},\mathcal{A}=\left\{(x,T)\in C^{1}(\mathbb{R},\Omega)\times(0,\infty):\begin{aligned} &\dot{x}(t)=f(x(t))\quad\forall t,\\ &x(t+T)=x(t)\quad\forall t,\\ &\int_{0}^{T}x_{1}(t)\,\mathrm{d}t=0,\\ &\int_{0}^{T}\lvert x_{1}(t)\rvert^{2}\,\mathrm{d}t\neq 0\end{aligned}\right\}, | (12) | | ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

which certainly contains any symmetric orbit.

Take C=σ2​(ρ−1)C=\sigma^{2}(\rho-1), Q=1Q=1 and

V​(x)=12​(σ​(ρ−2)​x12−σ2​x22−σ2​x32).V(x)=\frac{1}{2}\left(\sigma\left(\rho-2\right)x_{1}^{2}-\sigma^{2}x_{2}^{2}-\sigma^{2}x_{3}^{2}\right). (13)

Then

C​g2−ℒf​g2+ℒf​V​(x)=σ2​(ρ−1)​x12−σ2​(x2−x1)2+σ2​(ρ−2)​x1​(x2−x1)−σ2​x2​(x1​(ρ−x3)−x2)−σ2​x3​(x1​x2−β​x3)=β​σ2​x32≥0,\displaystyle\begin{split}Cg^{2}-\mathcal{L}_{f}g^{2}+\mathcal{L}_{f}V(x)&=\sigma^{2}(\rho-1)x_{1}^{2}-\sigma^{2}(x_{2}-x_{1})^{2}+\sigma^{2}(\rho-2)x_{1}(x_{2}-x_{1})\\ &\qquad-\sigma^{2}x_{2}\left(x_{1}\left(\rho-x_{3}\right)-x_{2}\right)-\sigma^{2}x_{3}(x_{1}x_{2}-\beta x_{3})\\ &=\beta\sigma^{2}x_{3}^{2}\geq 0,\end{split} (14)

and therefore proposition 2.3 implies that, for any symmetric orbit, T≥2​π/σ2​(ρ−1)T\geq 2\pi/\sqrt{\sigma^{2}(\rho-1)}. At the standard parameters, this gives a lower bound of approximately 0.1210.121.

3.2 Validated numerical results using SDP

In order to find a sharper lower bound for the period of symmetric orbits, we turn to numerics. For reasons of numerical conditioning, we rescale time by a factor of 66 and space by a factor of 2525 to instead study the system

d​x1d​t=60​(x2−x1),d​x2d​t=168​x1−150​x3−6​x2,d​x3d​t=150​x1​x2−16​x3.\displaystyle\begin{split}\frac{\mathrm{d}x_{1}}{\mathrm{d}t}&=60(x_{2}-x_{1}),\\ \frac{\mathrm{d}x_{2}}{\mathrm{d}t}&=168x_{1}-150x_{3}-6x_{2},\\ \frac{\mathrm{d}x_{3}}{\mathrm{d}t}&=150x_{1}x_{2}-16x_{3}.\end{split} (15)

This is a polynomial system with fully integer coefficients, giving neater validation, and also (empirically) better conditioned SDPs than the standard Lorenz system. Bounds for orbit periods for this system are equivalent to bounds for 11 at the standard parameters after multiplication by 6.

We apply proposition 2.3 with gg being the vector of monomials up to degree dgd_{g} which change sign under the symmetry transformation (x1,x2,x3)↦(−x1,−x2,x3)(x_{1},x_{2},x_{3})\mapsto(-x_{1},-x_{2},x_{3}). For example, with dg=2d_{g}=2,

g={x1,x2,x1​x3,x2​x3}.g=\{x_{1},x_{2},x_{1}x_{3},x_{2}x_{3}\}. (16)

As these are monomials, this collection is automatically linearly independent. The auxiliary function V=v𝖳​aV=v^{\mathsf{T}}a, with v∈ℝkv\in\mathbb{R}^{k}, is a polynomial constructed from a basis aa mostly of monomials, with a particular form of polynomial for the highest order terms, see the accompanying code and Goluskin [10] for details. It is sufficient to consider only VV which is invariant under the symmetry transformation. The maximum degree of the basis aa is dad_{a}. The SOS constraint in proposition 2.3 is enforced by requiring that

C​g𝖳​Q​g−(ℒf​g)𝖳​Q​(ℒf​g)+ℒf​(v𝖳​a)−be𝖳​Pe​be−bo𝖳​Po​bo=0,Cg^{\mathsf{T}}Qg-\left(\mathcal{L}_{f}g\right)^{\mathsf{T}}Q\left(\mathcal{L}_{f}g\right)+\mathcal{L}_{f}(v^{\mathsf{T}}a)-b_{e}^{\mathsf{T}}P_{e}b_{e}-b_{o}^{\mathsf{T}}P_{o}b_{o}=0, (17)

where PeP_{e} and PoP_{o} are positive semidefinite matrices and beb_{e} and bob_{o} are vectors of monomials up to degree dbd_{b} which are respectively invariant and change sign with respect to the symmetry transformation. This is a very simple form of block-diagonalisation with respect to the symmetry, see Gatermann and Parrilo [9] for an in-depth discussion of this topic.

For each choice of degrees dgd_{g}, dad_{a} and dbd_{b}, we reduce the value of CC until the SDP representing the SOS condition in proposition 2.3 becomes infeasible. Close to the limits of feasibility, the conditioning of the SDP becomes very poor, with the consequence that the metrics returned by the SDP solver are not a reliable indicator of feasibility. For this reason, we turn to rigorous numerics to validate the solution of the SDP. The approach is similar to that of Parker [21], except that we use rational arithmetic rather than interval arithmetic:

    1. For fixed CC, gg, aa, beb_{e} and bob_{o}, we attempt to solve 17 to find positive definite matrices QQ, PeP_{e} and PoP_{o} and real vector vv. To fix the scaling, we impose the additional constraint that QQ has unit trace. The problem is solved using JuMP [18] on top of the optimizer Hypatia [6], with the Double64 datatype from DoubleFloats.jl.
    1. If Hypatia indicates a likely feasible solution, we then attempt to find a rational solution to 17. The equation is converted to a large linear system A​y=cAy=c by comparing coefficients of monomials. Here yy is a vector containing the entries of QQ, PeP_{e}, PoP_{o} and vv. Note that by construction, AA and cc have only rational entries. An initial guess solution y0y_{0} is constructed by approximating the floating point result from Hypatia with rational numbers, with Julia’s arbitrary precision type BigInt used for numerator and denominator. An exact solution is constructed as
      y=y0+A+​(c−A​y0)y=y_{0}+A^{+}(c-Ay_{0}) (18)
      where A+A^{+} is the pseudoinverse constructed though QR decomposition, which has rational entries.
    1. We validate that we have an exact rational solution to 17 – this is necessary as the algorithm used to construct the pseudoinverse may be inexact.
    1. We validate that the reconstructed QQ, PeP_{e} and PoP_{o} from this yy are strictly positive definite. This is done with Sylvester’s criterion, which is exact with exact rational arithmetic.

The key step here is 18: we project the approximate, floating point solution returned by the SDP solver to an exact rational solution to 17, and then ensure that the resulting symmetric matrices QQ, PeP_{e} and PoP_{o} are indeed positive definite. If all of these validation steps pass then we have a rigorously validated bound. The full results are give in table 1. We successfully validated a lower bound which is sharp to four decimal places.

Table 1: The smallest possible CC for a solution to the SDP constructed from 17, for given polynomial degrees. In each case decreasing the final quoted digit of CC by 1 results in an SDP which fails to validate. The final column gives (to four decimal places) the corresponding lower bound for periods of symmetric orbits in the Lorenz system at the standard parameters. The underline indicates which digits are a sharp lower bound for the shortest known orbit with T≈1.55865T\approx 1.55865.

4 Minimal periods in the Hénon-Heiles system at a fixed energy level

The Hénon-Heiles [14] system is defined by the Hamiltonian

H​(x)=12​(x12+x22+x32+x42)+x12​x2−13​x23,H(x)=\frac{1}{2}\left(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}\right)+x_{1}^{2}x_{2}-\frac{1}{3}x_{2}^{3}, (19)

which gives a polynomial ODE system

d​xd​t=[x3,x4,−x1−2​x1​x2,−x2−x12+x22].\frac{\mathrm{d}x}{\mathrm{d}t}=\left[x_{3},\,x_{4},\,-x_{1}-2x_{1}x_{2},\,-x_{2}-x_{1}^{2}+x_{2}^{2}\right]. (20)

We focus on the fixed energy level H=1/7H=1/7, which is close to the loss of bounded trajectories at H=1/6H=1/6, has the advantage of having no fixed points, and has previously been studied using SOS methods [20]. The shortest period we were aware of at this energy level was T=6.9665T=6.9665 [20], but there is in fact a shorter one as we will discuss below. This system has a rotational symmetry which we do not exploit for technical reasons, but we will make use of the sign symmetry (x1,x2,x3,x4)↦(−x1,x2,−x3,x4)(x_{1},x_{2},x_{3},x_{4})\mapsto(-x_{1},x_{2},-x_{3},x_{4}) to accelerate the SDP.

To find lower bounds for the period, we use proposition 2.3 with

f​(x)=20​[x3,x4,−x1−2​x1​x2,−x2−x12+x22]f(x)=20\left[x_{3},\,x_{4},\,-x_{1}-2x_{1}x_{2},\,-x_{2}-x_{1}^{2}+x_{2}^{2}\right] (21)

which corresponds to 19 with time rescaled by a factor of 20. Let {wi}\{w_{i}\} be the set of nonconstant monomials in xx up to degree dwd_{w}, and then let{gi}\{g_{i}\} be a linearly independent basis for the space of polynomials spanned by {ℒf​wi}\{\mathcal{L}_{f}w_{i}\}, found by discarding any which are found to be dependent by QR factorisation – see the code for details. Note that each of these polynomials integrates to 0 over one period of any periodic orbit. Since this space includes ℒf​xi=fi\mathcal{L}_{f}x_{i}=f_{i}, the set 𝒜\mathcal{A} defined by 9 is exactly the set of all non-constant periodic orbits. One further optimization is to consider separately the wiw_{i} which are invariant or change sign under the sign symmetry mentioned above, which allows us to block-diagonalize QQ into QeQ_{e} and QoQ_{o} [9].

Table 2: The smallest possible CC for a solution to the SDP 22, for given polynomial degrees. In each case decreasing the final digit of CC by 1 results in an SDP which fails to validate. The final column gives (to four decimal places) the corresponding lower bound for periods of symmetric orbits in the Hénon-Heiles system with H=1/7H=1/7. Given dwd_{w}, we use db=dw+1d_{b}=d_{w}+1, da=2​dw+1d_{a}=2d_{w}+1 and dρ=2​dw−1d_{\rho}=2d_{w}-1 to ensure consistent maximum degrees. The underline indicates which digits are a sharp lower bound for the newly discovered orbit with T=6.0521T=6.0521.

To enforce H​(x)=1/7H(x)=1/7, we introduce a new unknown polynomial ρ𝖳​c\rho^{\mathsf{T}}c where cc is another basis of monomials invariant under the sign symmetry, with maximum degree dρd_{\rho}. By multiplying this unknown polynomial by (H−1/7)(H-1/7) and adding this to our SDP, we ensure that the solution is valid whenever H=1/7H=1/7 while disregarding any other value of HH. This is sometimes known as the S-procedure. Therefore the final equation we are trying to solve is

C​(ge𝖳​Qe​ge+go𝖳​Qo​go)−((ℒf​ge)𝖳​Qe​(ℒf​ge)+(ℒf​go)𝖳​Qo​(ℒf​go))+ℒf​(v𝖳​a)−be𝖳​Pe​be−bo𝖳​Po​bo+(6−42​H)​ρ𝖳​c=0,C\left(g_{e}^{\mathsf{T}}Q_{e}g_{e}+g_{o}^{\mathsf{T}}Q_{o}g_{o}\right)-\left(\left(\mathcal{L}_{f}g_{e}\right)^{\mathsf{T}}Q_{e}\left(\mathcal{L}_{f}g_{e}\right)+\left(\mathcal{L}_{f}g_{o}\right)^{\mathsf{T}}Q_{o}\left(\mathcal{L}_{f}g_{o}\right)\right)+\mathcal{L}_{f}(v^{\mathsf{T}}a)\\ -b_{e}^{\mathsf{T}}P_{e}b_{e}-b_{o}^{\mathsf{T}}P_{o}b_{o}+(6-42H)\rho^{\mathsf{T}}c=0, (22)

for positive definite QeQ_{e}, QoQ_{o}, PeP_{e} and PoP_{o}, and unknown vectors vv and ρ\rho.

Once the bases geg_{e}, gog_{o}, aa, beb_{e} and bob_{o} have been constructed, the validation proceeds exactly as per section 3.2, with one additional step: After the SDP is approximately solved, entries in the bases with close-to-zero coefficients are removed, and the SDP is repeated. This is done 5 times, which results in a much cleaner numerical result before the validation steps. We also used the Clarabel [12] optimizer instead of Hypatia as the latter reported infeasibility in all cases. With this approximate numerical solution, we project onto an exact rational solution for 22 and then check that the resulting symmetric matrices QeQ_{e}, QoQ_{o}, PeP_{e} and PoP_{o} are strictly positive definite. The full validated results are given in table 2. Observe that the bounds appear to saturate at T≈6.05T\approx 6.05, which is well below the shortest known period of 6.976.97.

Refer to caption

Figure 2: Two symmetry-related periodic orbits with period T≈6.0521T\approx 6.0521 in the Hénon-Heiles system 19 at energy level E=1/7E=1/7.

Consider the expression

C​g𝖳​Q​g−(ℒf​g)𝖳​Q​(ℒf​g)+ℒf​(v𝖳​a).Cg^{\mathsf{T}}Qg-\left(\mathcal{L}_{f}g\right)^{\mathsf{T}}Q\left(\mathcal{L}_{f}g\right)+\mathcal{L}_{f}(v^{\mathsf{T}}a). (23)

At a solution of the SDP, this expression must be everywhere nonnegative, and furthermore if CC is sharp then it must integrate to zero over the shortest periodic orbit. This means that with a good approximate solution, the value of this expression must be very close to zero around the shortest orbit. By minimizing this expression over the Poincare section x4=0x_{4}=0, we were able to find a new periodic orbit of period T≈6.0521T\approx 6.0521, starting from initial condition

x≈(0,−0.3979,0.2922,0).x\approx(0,-0.3979,0.2922,0). (24)

In fact there are two symmetry-related orbits with the same period, as shown in figure figure 2. Both trajectories are invariant under the 2​π/32\pi/3 rotation symmetry and map to each other under the reflection symmetry. These are neutrally stable (elliptical) orbits, so it is expected that the brute-force recurrency search of Oeri and Goluskin [20] did not detect them. With noticeably lower degrees than for the Lorenz system, we were able to validate bounds which agree with the period of the orbits to 4 decimal places.

5 Discussion and conclusion

In this paper we have presented a new method for computing rigorous lower bounds for periods in differential equation systems. We applied this to symmetric orbits in the Lorenz system and general orbits at a fixed energy level in the Hénon-Heiles system, in both cases proving bounds which agree to four decimal places with known periodic orbits.

Though our method is described for polynomial systems, it can be applied to non-polynomial systems through changes of variables [22] or by polynomial approximation. Similar approaches are also possible for partial differential equation systems using either a truncation to an ODE paired with rigorous bounds [11] or through weak integral formulations [4]. For discrete- rather than continuous-time dynamical systems, the analogous problem is to consider what integer periods are possible, though in that case it makes more sense to computer upper bounds for the period. This is covered in the appendix to Parker [21].

Two possibly related open questions arise directly from this work. The first is whether the inequality

| maxperiodic orbits(2​πT)2≤infϕ,V∈C1​(Ω,ℝ)C>0C s.t. C|ϕ(x)|2−|ℒfϕ(x)|2+ℒfV(x)≥0∀x∈Ω\max_{\text{periodic orbits}}\left(\frac{2\pi}{T}\right)^{2}\leq\inf_{\begin{subarray}{c}\phi,V\in C^{1}(\Omega,\mathbb{R})\\ C>0\end{subarray}}C\text{ s.t. }C\left|\phi(x)\right|^{2}-\left|\mathcal{L}_{f}\phi(x)\right|^{2}+\mathcal{L}_{f}V(x)\geq 0\;\forall x\in\Omega | (25) | | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- | ---- |

can be written as an equality. A sharpness result for related SOS problems is available [24] but this does not immediately apply in this case. If this were equality for C1C^{1} functions then we could get arbitrarily sharp bounds computationally by increasing the degrees of our polynomials. Our conjecture is that the sharp equality is possible at least for C1C^{1} Axiom A systems on a compact domain. The second question is whether we can find rigorous lower bounds close to T=1.557T=1.557 in the Lorenz system for general rather than symmetric orbits. Following an approach identical to that of section 4, we have been able to validate bounds T≥0.98T\geq 0.98, which already required significantly more effort than the symmetric case.

Data availability statement

Acknowledgements

The contributions of David Goluskin were instrumental in the preparation of this paper. JP also wishes to thank James Robinson, Warwick Tucker and Ian Tobasco for helpful discussions and references.

References