Some Reversing Orbits for a Rattleback Model

A physical rattleback is a toy that can exhibit counter-intuitive behavior when spun on a horizontal plate. Most notably, it can spontaneously reverse its direction of rotation. Using a standard mathematical model of the rattleback, we prove the existence of reversing motion, reversing motion combined with rolling, and orbits that exhibit such behavior repeatedly.

where x =ẋ − x × ω. Here, ω is the angular velocity, and a × b denotes the cross-product of two vectors a and b in R 3 .
In the case of the rattleback with mass m, a vertical gravitational force −mgγ acts at the center of mass G, where g is a gravitational acceleration. Suppose that the body stays in contact with a fixed horizontal plate and satisfies a no-slip condition v = r×ω. Here, v is the velocity of G and r denotes the vector from G to the point of contact C. Assuming conservation of momentum, we have mv = f − mgγ , where f is the force exerted on the body at C. Assuming conservation of angular momentum as well, we have (I ω) = r × f, where I is the inertia tensor about G.
Notice that r×f = mr×v +mgr×γ due to momentum conservation. Substituting this expression into the equation (I ω) = r×f, we end up with the equation of motion (1.1) In addition, we haveγ = γ × ω, due to the fact that γ = 0. The dynamic variables here are γ and ω. For the velocity v, we can substitute r × ω, and the vector r can be expressed in terms of γ by using the geometry of the body.
In this paper, we consider the body to be an ellipsoid in R 3 with principal semiaxes b 1 > b 2 > b 3 > 0. Consider the 3 × 3 matrix B = diag(b 1 , b 2 , b 3 ). Then, the equation for the surface of the body and the tangency condition at the point of contact C are given by (1.2) Using these equations, one easily finds that (1. 3) The inertia tensor I is assumed to be a symmetric strictly positive definite 3× 3 matrix. Then, I is invertible, and (1.1) together with the equationγ = γ × ω defines a flow on R 6 . This flow preserves the length = γ . A straightforward computation shows that another flow-invariant quantity is the total energy H = 1 2 ω I ω + 1 2 m v 2 + mgs . (1.4) The three terms on the right-hand side of this equation can be identified with the rotational kinetic energy, the translational kinetic energy, and the potential energy, respectively. We note that the no-slip condition v = r × ω is a non-holonomic constraint, so the rattleback model is not a Hamiltonian system. In what follows, we restrict to γ = 1, unless specified otherwise. Then, the phase space for our flow is S 2 × R 3 , where S 2 denotes the unit sphere in R 3 . The dimension can be reduced further from 5 to 4, if desired, by choosing an energy E > mgb 3 and restricting to the fixed-energy surface (1.5) Clearly, these invariant surfaces are all compact. So in particular, every orbit returns arbitrarily close to a point that it has visited earlier. This allows for a variety of different types of motion, including periodic, quasiperiodic, and chaotic orbits. For the parameters and energies considered in this paper, orbits that look periodic are abundant. However, finding nontrivial periodic orbits turns out to difficult, unless one focuses on reversible orbits. An important feature of the rattleback flow is reversibility. To be more precise, let be the flow for some vector field X on R n . That is, d dt t = X • t for all t ∈ R. Given an invertible map R on R n , we say that is R-reversible if R • t = −t • R for all times t. Reversible dynamical systems share many qualitative properties with Hamiltonian dynamical systems (Devaney 1976(Devaney , 1977Golubitsky et al. 1991;Lamb and Roberts 1998). But they need not preserve a volume. In fact, one of our results exploits the existence of stationary solutions that attract (or repel) nearby points with the same energy. Nontrivial attractors of the type seen in dissipative systems have been observed numerically in Borisov and Mamaev (2003), Gonchenko et al. (2005), Borisov et al. (2012).
A well-known consequence of R-reversibility is the following. Assume that is R-reversible, and that some orbit of includes two distinct points x and τ (x) that are both R-invariant. Then, the orbit is time-periodic with period 2τ . The proof is one line: This property will be used to construct symmetric periodic orbits for the rattleback flow.
It is well-known that the rattleback flow is R-reversible for the reflection (1.7) Here, and in what follows, we use the notation (x 1 , . . . , x n ) = [x 1 . . . x n ] for vectors in R n . A rattleback with ellipsoid geometry (1.2) has another symmetry: the flow commutes with the reflection S 0 (γ , ω) = (−γ , ω). Additional symmetries exist for special choices of the inertia tensor I. A standard choice in experiments is to take I 13 = I 23 = 0. Then, the system is invariant under a rotation by π about the vertical axis e 3 . In what follows, we always restrict to this situation. As a consequence, the flow commutes with the reflection S(γ , ω) = ((−γ 1 , −γ 2 , γ 3 ), (−ω 1 , −ω 2 , ω 3 )). And it commutes with the reflection S = SS 0 as well. Given that S commutes with R, our flows are RS -reversible, where and I 13 = I 23 = 0. For the gravitational acceleration we choose the value g = 40141 4096 . These parameters can be realized in a physical experiment, with the proper choice of units for length, mass, and time. (Possible units would be centimeters, decagrams, and deciseconds, respectively.) We note that the matrix I is strictly positive definite, and that the smallest possible energy of a point in M is mgb 3 = 40141 1024 = 39.2001953125. The chosen inertia tensor is of the form I = R −1 I 0 R, where I 0 is roughly the inertia tensor of a homogeneous solid ellipsoid with the given mass m and semi-axes b j , and where R is a rotation about the vertical axis e 3 by an angle δ π 20 .

Main Results
A trivial solution of the rattleback equation is the stationary solution with γ = (0, 0, 1) and ω = (0, 0, ω 3 ). If ω 3 = 0, then this corresponds to a steady rotation about the vertical axis. As mentioned earlier, one of the peculiar features of the rattleback is observed when starting with a nearby initial condition that is not a stationary point. If ω 3 is within a certain range of values, then the rotation is observed to slow down and eventually reverse direction. In order to give a precise definition of reversal, consider the column vectors α, β, and γ of the rotation matrix Q, and define the angle ψ 0 by the equation When evaluated along an orbit, this "yaw-angle" ψ 0 typically varies as a function of time. Denote by ψ : R → R a continuous lift of ψ 0 to the real line. We say that the body reverses its direction of rotation on a time interval [a, b], if there exists a time c ∈ [a, b] such that the differences ψ(c) − ψ(a) and ψ(c) − ψ(b) have the same sign and are bounded away from zero by some positive constant C. The largest such constant C will be referred to as the amplitude of the reversal, and the sign of ψ(c) − ψ(a) will be called the sign of the reversal. Our proof of this theorem is computer-assisted, in the sense that it involves estimates that have been verified (rigorously) with the aid of a computer. The same applies to the theorems stated below. The statement E = 39.683 . . . in Theorem 2.1 means that 39.683 ≤ E < 39.684. The same notation is used for other interval enclosures. We note that our actual bounds are much more accurate. The orbit mentioned in Theorem 2.1 is shown in Fig. 2. To be more precise, consider the angular velocity M = I ω−mr×v about the contact point C. This angular velocity has been used as primary variable (in place of ω) in several papers. A straightforward computation shows that (2. 2) The matrix I + m K (r) is strictly positive, so (2.2) could be used to express ω in terms of M. We note that the reflections R, S, and S commute with the change of variables (γ , ω) → (γ , M). Figure 2 shows the components of γ (left) and of M (right) as functions of time t, for the orbit described in Theorem 2.1. The R-reversibility of the orbit is equivalent to the condition that γ is an even function of t, while M is an odd function of t.
Remark 1 All of our results that refer to the parameters (1.9) hold for an open set of parameter values nearby. This is a consequence of nondegeneracy properties that are verified as part of our proofs.
Our next result concerns the existence of a reversing heteroclinic orbit between two stationary points of the form z c = (e 3 , c e 3 ). As will be shown in the next section, there exists a value c * = 1.048 . . . such that z c is repelling for c < −c * and attracting for c > c * , if the flow is restricted to the surface of fixed energy E = H(z c ). The first orbit described in this theorem is shown in Fig. 3. We note that this orbit must pass through an R-invariant point x at time t = 0. Since z −c is repelling and z c attracting (for fixed energy), every point that is sufficiently close to x and has the same energy as x lies on some heteroclinic orbit connecting z −c to z c . We expect that there exists heteroclinic orbits between z ±c for a range of values c > c * , and it is possible that such orbits exist for some range of values c < c * as well.
Numeral experiments are most often carried out for rattlebacks whose body is a cut-off elliptic paraboloid. If we replace F(r) = ρ 2 by F(r) = ρ 2 1 + ρ 2 2 − 2ρ 3 − 1, where ρ = B −1 r, then the behavior can be expected to be similar to the behavior of the ellipsoid, as long as γ stays close to e 3 . Among the features of the ellipsoid-shaped rattleback that cannot be studied in the cut-off paraboloid case is roll-over motion.
A possible definition of "rolling over e 1 " can be given in terms of the angle φ 0 defined by the equation When evaluated along an orbit, this "roll-angle" φ 0 typically varies as a function of time. Denote by φ : R → R a continuous lift of φ 0 to the real line. We say that the body rolls over e 1 on a time interval is no less than π in absolute value. The sign of φ(b) − φ(a) will be called the direction of the roll-over.

Some Simpler Solutions
After describing a periodic orbit that rolls over e 1 repeatedly in the same direction, we will discuss some stationary solutions and their stability. In fact, there exists a two-parameter family of such orbits.
The orbit described in this theorem is shown in Fig. 6. Our numerical results suggest that both the yaw-angle ψ and the roll-angle φ are monotone for this orbit, but we did not try to prove this.
We note that there exist trivial roll-over orbits as well as trivial heteroclinic orbits. Consider the manifold Fix(S ) = {(γ , ω) ∈ M : γ 3 = ω 1 = ω 2 = 0} that is invariant under the flow. At energy mgb 1 , we have a heteroclinic orbit in Fix(S ) between the points z ± = ((±1, 0, 0), 0). For energies below mgb 1 , the orbits are all closed and avoid z ± . For energies above mgb 1 , the orbits are closed and clearly roll over e 2 with the obvious definition of such a roll-over.
Next, we consider some stationary solutions. A stationary point x = (γ , ω) necessarily satisfies ω = ± ω γ , sinceγ = γ × ω has to vanish. The stability of x is best discussed in terms of the vector field X : The simplest stationary points are x j = (e j , 0), where e j is the unit vector parallel to the j-th coordinate axis. A straightforward computation shows that besides two eigenvalues zero (due to the conservation of and H), D X(x 1 ) has four real eigenvalues, D X(x 2 ) has two real and two imaginary eigenvalues, and D X(x 3 ) has four imaginary eigenvalues. This holds for any ellipsoid body with The stationary point x 3 is part of a family of stationary points z c = (e 3 , c e 3 ) parametrized by a real number c. The stability of these points has been investigated in several papers, including (Markeev 1983;Bondi 1986;Garcia and Hubbard 1988;Pascal 1994;Franti 2013). The consensus is that, for many choices of parameters, an analogue of the following holds. A pencil-and-paper proof of this lemma is possible, but tedious; so we carried out the necessary (rigorous) computations with a computer. Notice that it suffices to prove the assertions for We remark that Lemma 3.2 excludes the existence of a real analytic first integral that is independent of and H. A result on the non-existence of analytic integrals was proved also in Dullin and Tsygvintsev (2008).
Finally, let us describe two one-parameter families of RS -invariant stationary points. We have not found them discussed in the literature.
The existence of stationary points near (±e 2 , 0) or (±e 1 , 0) may be a known fact. What seems surprising is that the two critical solutions (corresponding to ω = ∞) are related via a rotation by π 2 . This part of Lemma 3.3 is not specific to the choice of parameter values (1.9).
The remaining parts of this paper are devoted to our proofs of the results stated in Sects. 2 and 3.

Integration and Poincaré Sections
The equation (1.1) determinesω as a function of x = (γ , ω). Together with the equationγ = γ × ω, this defined a vector field X = (γ ,ω) on R 6 . This vector field is considered only in a small open neighborhood of M in R 6 where it is real analytic. The resulting flow (t, x) → t (x) is then real analytic as well.
In order to construct an orbit for a non-small time-interval [0, r ], we partition this interval into m small subintervals [τ k−1 , τ k ] with τ 0 = 0 and τ m = r . On each successive subinterval, starting with k = 1, we solve the initial value problemẋ = X (x) with given initial conditions at time τ k−1 via the integral equation is a sufficiently small positive real number, then for 1 ≤ j ≤ 6, the function g defined by g(t) = x j (τ k−1 + t) is given by its Taylor series about t = 0 and has a finite norm This norm is convenient for computer-assisted proofs, since it is easy to estimate, and since the corresponding function space G ρ is a Banach algebra for the pointwise product of functions. Each function g ∈ G ρ extends analytically to the complex disk |z| < ρ and continuously to its boundary. In what follows, ρ is a fixed but arbitrary positive real number. Consider now the integral equation (4.1), with k = 1 to simplify the discussion. Let x 0 (t) = x(0). Since the vector field X defines an analytic function on some open neighborhood of x(0), the equation (4.1) can be solved order by order by iterating the transformation K given by starting with x = x 0 . That is, the Taylor polynomial x n of order n for K n (x 0 ) agrees with the Taylor polynomial of order n for x. This is of course the well-known Taylor integration method. In order to estimate the higher-order correction x − x d for some large degree d, we use a norm on G 6 ρ of the form with appropriately chosen weights w j > 0. A common approach is to apply the contraction mapping theorem on a ball centered at x d . Instead, we use Theorem 5.1 in Arioli and Koch (2015), which only requires that some closed higher-order set gets mapped into itself. In our programs, the radius ρ and the weights w j are chosen adaptively, depending on properties of X (x d ).
Computing K(x) from the Taylor series for x involves only a few basic operations like sums, products, antiderivatives, multiplicative inverses, and square roots. This is done by decomposing each function g involved into a Taylor polynomial g d of some (large) degree d and a higher-order remainder g − g d . For sums, products, and antiderivatives of functions in G ρ , it is trivial to estimate the higher-order remainder of the result.
Consider now the multiplicative inverse 1 + h of a function 1 + g, where The following is straightforward to prove.
The Taylor coefficients c n of h are given recursively by We note that the same holds if g takes values in some commutative Banach algebra X with unit. Then, the Taylor coefficients b n and c n in (4.5) are vectors in X . This fact is used when estimating the flow for initial points that depend on parameters.
Next, consider the (principal branch of the) square root of a function 1 + g.
Proposition 4.2 Let g ∈ G ρ with g < 1 2 . Then h = (1 + g) 1/2 − 1 belongs to G ρ . The Taylor coefficients c n of h are given recursively by provided that the norm on the right-hand side of this inequality does not exceed 1 4 . Proof We will use that (4.10) Verifying this identity and (4.8) is straightforward. Let now n ≥ 2. Using the power series for z → (1 + z) 1/2 − 1, and the fact that 1/2 k ≤ 1 8 for k ≥ 2, one easily finds that h n ≤ 5 8 g . This in turn yields a bound (1 + h n ) −1 ≤ 8 5 . So from (4.10), we find that δ = 4 5 h − h n satisfies (4.11) Assuming that ε ≤ 1 4 , this implies the bound (4.9). Next, we consider the problem of constructing reversible orbits for the given flow.
In what follows, we will use M as a primary variable instead of ω. The equation of motion in the variables (γ , M) is given bẏ (4.12) Here, r andṙ are obtained from (1.3), while ω is determined from r and M via (2.2). In order to simplify notation, we will now write x = (γ , M) and X = (γ ,Ṁ). Recall that R and RS commute with the change of variables (γ , ω) → (γ , M).
For the construction of periodic orbits, it is convenient to consider return maps to some codimension 1 surface . The surfaces used in our analysis are for j = 1 or j = 3. Since we are exploiting reversibility, only half-orbits or quarterorbits need to be considered. Each partial orbit starts at some symmetric (meaning R-invariant or RS -invariant) point. The goal is to determine such a point x, as well as a positive time τ = τ (x), such that τ (x) is again symmetric. To this end, we first determine a symmetric numerical approximationx for x and an approximationτ for τ . After choosing a real number τ slightly smaller thanτ , the associated Poincaré map P is then defined by setting for all (symmetric) starting points x in some neighborhood ofx.
Consider now the problem of constructing the orbit described in Theorem 2.1. The starting point at time t = 0 is R-invariant and thus of the form x = (γ , 0). Restricting to γ = 1, the possible starting points are parametrized by a vector γ = (γ 1 , γ 2 ) in R 2 of length less than 1. For the Poincaré section, we choose = 1 . Then, Our proof of this lemma is computer-assisted and will be described in Sect. 6. Notice that if γ ∈ B g is a solution of P(γ 1 , γ 2 ) = (0, 0), and if we set x = (γ , 0), then the point τ (x) is R-invariant for τ = τ (x). Thus, as described earlier, this implies that In order to construct the orbit described in Theorem 2.3, we use a Poincare map P with = 3 . The starting point x is again R-invariant, but the desired pointx = P(x) is RS -invariant, meaning that γ 3 = M 3 = 0. So the goal is to find zeros of the function P defined by P(γ ) =γ 3 .

Lemma 4.4
There exists a vectorγ ∈ S 2 such that the following holds. Letx = (γ , 0). Then, the Poincare map P with = 3 and τ = 63 is well-defined and real analytic in an open neighborhood B g × B M ofx in M. When restricted to B g , the associated function P is real analytic and takes the value 0 at some RS -invariant point. Furthermore, all orbits with starting points in B g × {0} have Poincaré time τ (x) = 63.57172 . . ., energy E = 42.0308 . . . and reverse/roll-over as described in Theorem 2.3. The same holds for a two-parameter family of RS -invariant initial points.
Our proof of this lemma will be described in Sect. 6. Notice that if γ ∈ B g is a solution of P(γ 1 , γ 2 ) = 0, and if we set x = (γ , 0), then the point τ (x) is RS -invariant for τ = τ (x). Thus, by R-reversibility, the point x = −τ (x) is RSinvariant as well. As described earlier, this implies that T (x ) = x with T = 4τ (x).

So Theorem 2.3 follows from Lemma 4.4.
Remark 2 A lemma analogous to Lemma 4.4 holds for the orbit described in Theorem 3.1, with RS -invariant starting points x. Choosing again = 3 and P =γ 3 , the equation that needs to be solved is P(γ 1 , M 1 , M 2 ) = 0. Here, the value of τ used in (4.14) is τ = 8.5.
In both cases, the goal is to prove that there exists a time τ > 0 such that τ (x) belongs to an open neighborhood ofx in M E that is attracted tox under the flow. To this end, consider the map P E : M E → R 4 given by P E (x) = Px, where P is as defined at the beginning of this section. Then, the equation of motion on M E near the origin is conjugate via P E to the equatioṅ in some open neighborhood ofx in M E . The stationary point for the associated flow isȳ = 0.
Notice that DY (0) = L(ω 3 ). Using Lemma 3.2, we have chosenx in such a way that ω 3 > c * . So we know that all eigenvalues of L(ω 3 ) have a negative real part. We expect that all eigenvalues are simple. Then, there exists an inner product . , . on R 4 such that the matrix is strictly negative definite, where L(ω 3 ) * denotes the adjoint of L(ω 3 ) with respect to the above-mentioned inner product. Assume for now that is strictly negative definite, meaning that u, u is negative for every nonzero vector u ∈ R 4 . Then, the derivative d dt y, y = 2 y, y + 2 y, N (y) , is negative, if the nonlinear part N (y) is sufficiently small compared to y. Let now y τ = P τ (x). In order to prove that y τ is attracted to zero by the flow associated with Y , it suffices to show that y τ belongs to a ball B ⊂ R 4 that is centered at the origin, with the property that | y, N (y) | < | y, y | for all nonzero y ∈ B. This property is equivalent to Lemma 5.1 Let y τ = P τ (x), with τ = 100 and x as described above (either the R-invariant or the RS-invariant choice). Then, there exists an inner product . , . on R 4 such that u, u is negative for every nonzero u ∈ R 4 . Moreover, there exists δ > 0 such that y τ belongs to the ball B = y ∈ R 4 : | y, y | 1/2 < δ , and such that the condition (5.10) holds. Furthermore, the orbit for x has the energy and reversing property described in Theorem 2.2.
Our proof of this lemma is computer-assisted and will be described in Sect. 6. Notice that t (x) →x as t → ∞, since the norm of y(t) = P t (x) tends to zero by (5.9) and (5.10). Furthermore Sx =x. So by reversibility, t (x) converges to Rx = RSx as t → −∞. In other words, we have a heteroclinic orbit connecting Rx = RSx = (e, −ω 3 e) to x = (e, ω 3 e). So Theorem 2.2 follows from Lemma 5.1.
In the remaining part of this section, we give a proof of Lemma 3.3, based in part on (trivial) estimates that have been carried out with the aid of a computer (Supplementary material). These estimates are specific to the choice of parameters (1.9), but analogous estimates should work for many other choices. The remaining arguments only use that b 1 = b 2 and I 13 = I 23 = 0.
Sketch of a proof of Lemma 3.3. We consider the equation for a stationary solution (γ , ω) with the property that γ 3 = ω 3 = 0. Then, ω = (ω 1 , ω 2 ) must be parallel to γ = (γ 1 , γ 2 ). So ω = ± ω γ . Consider also the conditionṀ = 0. From (4.12), we see that the first two components ofṀ vanish automatically. And the conditioṅ M 3 = 0 becomes This condition can be written as an equation for r = (r 1 , r 2 ) by using that γ j = −sb −1 j r j and ω = ± ω γ . To be more specific, we define two functions P and Q by the equation A straightforward computation shows that (5.11) reduces to If we stay away from the zeros of Q, then the condition is satisfied for some value of ω if and only if P(r ) and Q(r ) have opposite signs. In addition to (5.13), we also have the ellipse condition (r 1 /b 1 ) 2 +(r 2 /b 2 ) 2 = 1. So define p(θ ) = P(r ) and q(θ ) = Q(r ), using r 1 = b 1 sin(θ/2) and r 2 = b 2 cos(θ/2). Both p and q are 2π -periodic functions, since P and Q are even functions of r . In the remaining part of this paragraph, we consider just the parameters (1.9). Restricting θ to the interval [−π, π], the sign of p(θ ) is just the sign of −θ . So it suffices to determine the sign of q(θ ). This is easily done by using interval arithmetic. By estimating q and its derivative q on subintervals, one finds that q has exactly two zeros. Finally, using a (rigorous) Newton method, the zeros are located at values θ * = −0.242951 . . . and θ * = 3.13056 . . .. For details, we refer to the source code of the program RSp_Stat in Supplementary material.
Notice that γ 1 γ 2 = b 2 b 1 tan(θ/2). When computing these ratios numerically, it appears that the the vector γ for the angle θ = θ * is orthogonal to the vector γ for the angle θ * . The following argument confirms this observation.
Consider now Q as a function of γ , say Q(r ) = Q(γ ). Let γ be a solution of Q(γ ) = 0. This property of γ is equivalent to M 2 ω 1 − M 1 ω 2 = 0, meaning that So ω is an eigenvector of I + m K (r ). Equivalently, ω is an eigenvector of m −1 I − rr . But ω is parallel to γ , so γ is an eigenvector as well.
for some real number λ. This property is equivalent to the condition Q(γ ) = 0.

Computer Estimates
What remains to be done is to prove Lemmas 4.3, 4.4, and 5.1. (Our proof of the lemma referred to in Remark 2 is analogous to the proof of Lemma 4.4, so we will not discuss it separately here.) The necessary estimates are carried out with the aid of a computer. This part of the proof is written in the programming language Ada [24] and can be found in Supplementary material. The following is a rough guide for the reader who wishes to check the correctness of our programs.

Enclosures and Data Types
By an enclosure for (or bound on) an element x in a space X we mean a set X ⊂ X that includes x and is representable as data on a computer. For points in R n , this could be rectangles that contains x. Working rigorously with such enclosures is known as interval arithmetic. What we need here are enclosures for elements in Banach spaces, such as functions g(t) = n b n t n in the spaces G ρ described earlier. In addition, when considering orbits that depend on parameters (such as initial conditions), the coefficients b n can be functions themselves.
In our programs, enclosures are associated with a data type. Let X be a commutative real Banach algebra with unit 1. Our data of type Ball are pairs B = (B.C, B.R), where B.C and B.R are representable numbers, with B.R ≥ 0. The enclosure associated with a Ball B is the ball B X = {x ∈ X : x − (B.C)1 ≤ B.R}. For specific spaces X , other types of enclosures will be described below. In all cases, enclosures are closed convex subsets of X that admit a canonical finite decomposition where each x n is a representable element in X , and where each B(n) is a Ball centered at 0 or at 1.
Assume that X carries a type of enclosures named Scalar. For vectors in X 3 , we use a Scalar-type enclosure for each component. The corresponding data type SVector3 is simply an array(1..3) of Scalar. Our type Point defines enclosures for points x = (γ , M) with γ , M ∈ X 3 . But a Point P is in fact a 7-tuple P=(P.Alpha, P.Beta, P.Gamma, P.M, P.Energy, P.YawPi, P.RollPi), where the first four components are of type SVector3. The component P.Energy is a Scalar that defines an enclosure for the energy of a point, while P.YawPi and P.RollPi are integers. More specifically, P.YawPi = (ψ − ψ 0 )/π and P.RollPi = (φ − φ 0 )/π , where ψ is the lifted yaw-angle and φ the lifted roll-angle for points x in the enclosure given by P. The type Point is defined in the Ada package Rattleback.
Consider now a function g : D → X on a disk D = {z ∈ C : |z| < ρ} with representable radius ρ > 0. Denote by G the space of all such functions that admit a Taylor series representation g(z) = ∞ n=0 b n z n and have a finite norm g = ∞ n=0 b n ρ n . Here, b n ∈ X for all n. A large class of enclosures for functions in this space is determined by the type Taylor1, which is defined in the Ada package Taylors1. Since this type has been used several times before, we refer to (Arioli and Koch 2018) for a rough description and to Supplementary material for details.
Our integration method uses a much simpler type named Taylor. A Taylor P is an array(0..d) of Scalar, where d is some fixed positive integer. The associated enclosure is the set Here, P(d) G is obtained from S = P(d) X by replacing each ball B(n) X in the decomposition (6.1) of S by the corresponding ball B(n) G . The first d terms in (6.2) provide enclosures for the polynomial part g d−1 of a vector g ∈ G, as defined in (4.5), while the last term provides an enclosure for both the coefficient b d and the remainder g − g d . A precise definition of the type Taylor is given in the package ObO (an abbreviation for order-by-order). For analytic curves with values in X 3 , we use Taylor-type enclosures for each component via a type TVector3, which is simply an array(1..3) of Taylor.
Enclosures for real analytic curves t → x(t) on D are defined by the type Curve that is introduced in the package Rattleback.Flows. In our programs, a Curve C is a quadruple (C.Alpha,C.Beta,C.Gamma,C.M) whose four components are of type TVector3. These enclosures are used in our bounds on the integral operator K defined by (4.3).
We note that the types Point, Taylor, and Curve depend on choice of the Banach algebra X via the type Scalar. In the case X = R, we instantiate the package Rattleback and others with Scalar => Ball. For our analysis of the characteristic polynomial (5.6), which depends on two parameters ω 3 and λ, we use an instantiation of Rattleback with Scalar => TTay, where TTay defines enclosures for real analytic functions of two variables. (Hurwitz.TTay is a Taylor series in λ whose coefficients are Taylor series in ω 3 .) Another Banach algebra T that is very useful consists of pairs (u, u ), where u ∈ X and u ∈ X n . Addition and multiplication by scalars is as in X n+1 . The product of (u, u ) with (v, v ) is defined as (uv, uv + u v). If one thinks of u and v as being functions of n parameters, then u and v transform like gradients. Enclosures for element in T use a data type Tangent that is defined in the package Tangents. They are used to obtain bounds on the derivative of Poincaré maps (by using Scalar => Tangent) without first having to determine a formula for the derivative.

Bounds and Procedures
The next step is to implement bounds on maps between the various spaces. By a bound on a map f : X → Y, we mean a function F that assigns to a set X ⊂ X of a given type (say Xtype) a set Y ⊂ Y of a given type (say Ytype), in such a way that y = f (x) belongs to Y whenever x ∈ X . In Ada, such a bound F can be implemented by defining an appropriate procedure F(X: in Xtype; Y: out Ytype). In practice, the domain of F is restricted: if X does not belong to the domain of F, the F raises an Exception which causes the program to abort.
The The Ada package that defines a certain type also defines (usually) bounds on the basic operations that involve this type. In particular, bounds on the maps g → g −1 and g → g 1/2 on G are implemented by the procedures Inv and Sqrt, respectively, in the package Obo that defines the type Taylor. In the spirit of order-by-order computations, these procedures include an argument Deg for the order (degree) that needs to be processed. At the top degree, which corresponds to the last term in (6.2), the procedures Inv and Sqrt determine bounds on the higher order terms, using the estimate given in Propositions 4.1 and 4.2, respectively.
Bounds involving the type Point are defined mostly in Rattleback. This includes a procedure Ham that implements a bound on the energy function H, and a procedure VecField that implements a bound on the vector field (γ , M) → (γ ,Ṁ). Several other procedures deal with the construction of points (initial conditions) with prescribed properties; their role is described by short comments in our programs.
The package Rattleback.Flows implements bounds on the time-t maps t and various Poincaré maps. The first few procedures deal with the order-by-order computation of cross products and other basic operations. They maintain temporary data, so that lower order computations do not have to be repeated. And some of them can run sub-tasks in parallel, using the standard tasking facilities that are part of Ada [24]. The procedure VecField combines these computations into a bound on the vector field x →ẋ as maps between enclosures of the type Curve.
A bound on the solution of the integral equation K(x) = x is implemented by the procedure Integrate. After the polynomial part x d of the solution x has been determined, a bound on x − x d is obtained by first guessing a possible enclosure S ⊂ G 6 for this function, and then checking that x d + S is mapped into itself by the operator K. Using Theorem 5.1 in Arioli and Koch (2015), this guarantees that K has a unique fixed point in x d + S. We note that Integrate first determines a proper value of the domain parameter ρ for the space G = G ρ . This defines the time-increments τ k − τ k−1 used in (4.1).
Poincaré maps are now straightforward to implement. The type Flt_Affine specifies an affine functional F : X 6 → X whose zero defines a Poincaré section . To be more specific, F(γ , M) only depends on M. Besides an argument F that specifies F, the procedure Sign_Poincare also includes an argument TNeed for the time τ that enters the definition (4.14). Now, Sign_Poincare uses (an instantiation of) the procedure Generic_Flow to iterate Integrate, until t (x) with t ≥ τ lies on . A bound on the zero of t → F( t (x)) is determined by using the Newton-based procedure ObO.FindZero. We note that t is of type Scalar, so the stopping time τ = τ (x) can depend on parameters, if X = R.

Main Programs
Our proof of Lemma 4.3 is organized in the programs R_Der and R_Point. The initial point x = (γ , 0) is determined from data of type Point that are read from a file (Supplementary material). It suffices to control the map P described before Lemma 4.3 on a square centered at γ = (γ 1 , γ 2 ). The chosen square is 2ε × 2ε, with ε = 2 −2000 . This square also determines a domain B g ⊂ S 2 via the constraint γ = 1. After instantiating the necessary packages, the program R_Der computes an enclosure for the derivative D P on R and saves the result to a file. It also verifies that B g × B M belongs to the domain of the Poincaré map for some open neighborhood B M of the origin in R 3 . The program R_Point uses the above-mentioned enclosure for D P to verify that a quasi-Newton map associated with P maps R into its interior.
An additional program RSpR_Der can be used (optionally) to prove that D P is nonzero. This implies that the two-parameter family mentioned in Lemma 4.4 is real analytic.
The bounds referred to in Remark 2 are verified via the program Roll_Point. This program is analogous to RSpR_Point. And there is an analogue Roll_Der of RSpR_Der.
The bounds needed for Lemma 5.1 are organized by the programs Het, HetRS, and Basin. Both Het and HetRS run Plain_Flow for a time τ = 100. The initial point x is as described in Sect. 5. Enclosures for x and τ (x) are saved to data files. These files are then read by the procedure Check in Basin.
An upper bound LambdaMax on the spectrum of the (negative) linear operator defined by (5.8) is determined and shown by Basin.Show_Linear. This is done by via approximate diagonalization. The matrix that diagonalizes approximately also defines the inner product used in (5.10). Then, Basin.Show_NonLinear computes and shows an upper bound on the absolute value of the ratio y, ϑ −1 N (ϑ y) / y, y for y ∈ ∂B. By construction, this bound is non-decreasing in ϑ, so it suffices to consider ϑ = 1. At the end, (5.10) can be (and has been) checked by inspecting the output from Basic.
All of these programs were run successfully on a standard desktop machine, using a public version of the gcc/gnat compiler (https://www.gnu.org/software/gnat/). Instructions on how to compile and run these programs can be found in the file README that is included with the source code in Supplementary material. We note that the running times are rather long-days for some programs. This is due to the fact that our orbits are quite long, and that we need to use MPFR and rather high Taylor orders to get the accuracy needed.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.