Integrability of Nonholonomically Coupled Oscillators

We study a family of nonholonomic mechanical systems. These systems consist of harmonic oscillators coupled through nonholonomic constraints. In particular, the family includes the so called contact oscillator, which has been used as a test problem for numerical methods for nonholonomic mechanics. Furthermore, the systems under study constitute simple models for continuously variable transmission gearboxes. The main result is that each system in the family is integrable reversible with respect to the canonical reversibility map on the cotangent bundle. This result may explain previous numerical observations, that some discretisations of the contact oscillator have favourable structure preserving properties.


Introduction
In this paper we study a family of continuous dynamical systems consisting of nonholonomically coupled oscillators. The main result is that the family is integrable and that integrability is stable under reversible perturbations. Although numerical discretisations are not discussed in the paper, our main motivation is to obtain a rigorous explanation of the numerical observations by McLachlan and Perlmutter [16], that some discretisations for nonholonomic problems have structure preserving properties.
We continue the introduction with some background in Hamiltonian, nonholonomic, and reversible dynamics. In § 2 we give a presentation of the family of problems under study, and we state the main results. These results are then proved in § 3.
The concept of Liouville-Arnold integrability is fundamental in the study of continuous Hamiltonian dynamical systems. Through Kolmogorov-Arnold-Moser (KAM) theory, Liouville-Arnold integrability implies a foliation of the phase space in invariant tori, not only for the system itself, but also for "nearby" Hamiltonian systems (see Arnold [2,Ch. 10] and references therein).
By backward error analysis a numerical discretisation method for a continuous dynamical system can, at least formally, be interpreted as the exact solution of a nearby system [10,Ch. IX]. Roughly speaking, the KAM result then states that symplectic methods applied to Liouville-Arnold integrable systems preserve perturbed invariant tori [20,21]. This explains the excellent performance of symplectic methods for Hamiltonian systems.
In nonholonomic mechanics, there is no corresponding KAM theory, so the foliation in invariant tori are in general not stable with respect to perturbations (see the discussion in the end of [1, § 5.4.2]). Nevertheless, integrability of nonholonomic systems is a well established research topic. In the current literature there are three main techniques to obtain complete integrability of nonholonomic systems. One technique uses reduction by symmetry, and examples can be found in [5]. Next, the special case of systems of Chaplygin type is treated using a nonholonomic version of the Hamilton-Jacobi theory in [11,7,17]. Finally, one can appeal to the existence of an invariant measure, which is the technique used in [14] and [1, § 5.4]. All these techniques depend on some special property of the nonholonomic system at hand, and our case is no exception. Indeed, the essential property we use is that the underlying differential equation is a non-autonomous linear system whose fundamental matrix belongs to SO(n).
For discussions on numerical discretisations of nonholonomic problems we refer to [16,8,9,12,13], to the monograph by Cortés Monforte [6, § 7] and to the references therein. Due to the lack of KAM theory for nonholonomic mechanics, it is unclear which structure to preserve when discretising nonholonomic systems in order for the modified differential equation to remain integrable.
Let Q be a smooth configuration manifold of dimension n. The corresponding phase space is given by the cotangent bundle T * Q. Recall that T * Q comes with a canonical symplectic form, which, in local canonical coordinates (q i , p i ), is given by ω = i dq i ∧ dp i . In addition, T * Q comes with a canonical reversibility map (as does any vector bundle), namely the bundle map ρ : T * Q → T * Q given in local coordinates by ( A continuous dynamical system on T * Q with flow map ϕ t is called reversible if the flow map is reversible. If X is the vector field that generates ϕ t , then the system is reversible if and only if T ρ • X = −X • ρ, or equivalently, ρ * X = −X, where ρ * X is the pushforward of the vector field (see, e.g., [15, § 4.3]). A reversible system is called integrable reversible if T * Q can be foliated into invariant tori, which are also invariant with respect to the reversibility map ρ. That is, if there exists a reversible change of coordinates into action-angle variables [10, Def. XI.1.1].
There are various ways to define nonholonomic systems. For an overview, see the monograph by Bloch [4,Ch. 5]. In our setting, a nonholonomic system is defined by a triple (Q, D, H), where Q is a configuration manifold, D is a regular distribution on Q, i.e., a vector subbundle of T Q, and H is a Hamiltonian, i.e., a smooth function on T * Q. The constraints are given byq ∈ D. Equivalently, if D * ⊂ T * Q is the annihilator of D, then the constraints are D * ,q = 0. If the one-forms τ 1 , . . . , τ k span the codistribution D * , then the Hamiltonian version of the Lagrange-d'Alembert principle yields the governing equations (see [4, § 5.2]) where the quantities τ ij are defined by τ j = i τ ij dq i , and λ j are Lagrange multipliers. We assume that the differential algebraic equation (1) defines an underlying ordinary differential equation (ODE) on the manifold Notice that if H is regular, which we assume, then M is a fibre bundle over Q. M is a vector subbundle of T * Q if and only if H is quadratic in p (which is typically the case). It is closed under the canonical reversibility map ρ on T * Q if and only if the Hamiltonian is reversible, i.e., H • ρ = H, and in this case the corresponding dynamical system on M is reversible, i.e., ρ • ϕ t = ϕ −t • ρ. We say that the nonholonomic system (1) is integrable reversible if its corresponding system on M is integrable reversible.

Main results
We consider a family of nonholonomic systems on the configuration space Q = R 3 , with a distribution of the form where f is a smooth function on R, and a Hamiltonian function on T * R 3 of the form Figure 1: Illustration of a continuously variable transmission gearbox for which the dynamics is governed by the nonholonomic system (5). The belt between the two cones is translated along the shafts in accordance with the (q 3 , p 3 ) subsystem, thus providing a varying transmission ratio. The belt is kept in a plane perpendicular to the shafts, so that the belt keeps a constant length. The variables q 1 and q 2 denote the angular deflections of the shafts.
where m i > 0, k i > 0, ε ≥ 0 are constants and F is a smooth, reversible and regular Hamiltonian function on T * R. From (1) we get the governing equations aṡ From a mechanical point of view, this system consists of two linear oscillators, coupled through a single nonholonomic constraint which depends on the independent Hamiltonian subsystem described by the Hamiltonian F (q 3 , p 3 ).
Remark 2.1. System (5) is a simple model for a continuously variable transmission (CVT) gearbox, see Figure 1 for an illustration. Indeed, the shafts are attached to spiral springs that are fixed to a chassis, and the system governing q 3 is an arbitrary Hamiltonian system with Hamiltonian F , which explains the form of the Hamiltonian (4). Now the belt imposes a constraint between the rotation velocitiesq 1 andq 2 which, assuming that both cones have the same size, is determined by If we assume that q 3 < 1 (which correspond to assuming that the gear ratio is finite), then we obtain a distribution of the form (3). As a result, the system (5) describes a CVT gearbox where the "driver" is continuously and periodically shifting according to the independent subsystem (q 3 , p 3 ), and f (q 3 ) describes the transmission ratio.
The conceptualisation of CVTs is attributed to Leonardo da Vinci in 1490. In 1879, Milton Reeves invented a CVT for saw milling, which he later used also for cars. Today, CVT gearboxes are used in a variety of vehicles, particularly small tractors, snowmobiles, and motorscooters.
The proof is given in § 3 below. Consider now a perturbed Hamiltonian H ε = H 0 + εG, where ε ≥ 0 and G is a reversible smooth function on T * R 3 . We need a Lemma on perturbations of reversible nonholonomic systems, which we specialise to our case. Lemma 2.3. For small enough ε, the underlying ODE of the nonholonomic system defined by (R 3 , D, H ε ) is isomorphic to a reversible perturbation of the underlying ODE of (R 3 , D, H 0 ).
Proof. First, notice that H ε is reversible and for small enough ε it is regular. Let M ε be the subbundle given by (2) with H = H ε , and let X ε be the corresponding reversible vector field on M ε that generates the dynamics. For small enough ε the map ∂Hε ∂p : M ε → D is a reversible bundle automorphism. Thus, X ε induces a reversible vector field X 0,ε on M 0 by This vector field is an ε perturbation of X 0 , i.e., X 0,ε = X 0 + O(ε).
An integrable system is called KAM stable if "almost every" invariant torus, in the sense made precise e.g. in [18,19], is stable under either Hamiltonian or reversible perturbations. By using Theorem 2.2 and Lemma 2.3 we can now use the reversible KAM theorem on the system (5) to show KAM stability under reversible perturbations. then the foliation into tori of the underlying ODE of (5) is KAM stable under reversible perturbations.
Proof. If the (q 3 , p 3 ) subsystem is integrable reversible it follows from Theorem 2.2 that the underlying ODE of system (5) is integrable reversible. If f, F are real analytic it follows from the proof of Theorem 2.2 (see §3 below) that ω and ξ are real analytic. From Lemma 2.3 it follows that a reversible perturbation of (5) corresponds to a reversible perturbation of the underlying ODE of (5). The KAM theorem, as given in [19,Thm. 1], then states that the system (6)

Proof of Theorem 2.2
The proof involves the following steps: 1. Derive the ordinary differential equation constituting the dynamical system on M 0 .

Derive action-angle variables.
3. Show that the transformation to action-angle variables preserves reversibility.

Ordinary differential equation
Our aim is to reduce the constrained system (5) into an ordinary differential equation. From the governing equation (5) it follows that p 2 /m 2 = −f (q 3 )p 1 /m 1 , which implies that the energy can be written By differentiating the identity p = 1 α 1 (q 3 ) p 1 we obtain dp = α 1 (q 3 ) dp 1 m 1 + α 2 (q 3 ) dp 2 m 2 , By replacingṗ 1 andṗ 2 in the equation above, and by substituting for p in the equation forq 1 andq 2 , we now obtain the following system.
We have thus obtained an ordinary differential equation describing the system (5) in the variables (q 1 , q 2 , q 3 , p, p 3 ) ∈ R 3 × R 2 , which parametrize M (which we think of as a trivial fibre bundle over Q = R 3 ). Since the Hamiltonian is reversible it follows that the canonical reversibility map on T * R 3 restricts to the map ρ : (q 1 , q 2 , q 3 , p, p 3 ) → (q 1 , q 2 , q 3 , −p, −p 3 ), and that system (11) is reversible with respect to this map.

Integrability
First, since we assumed that the Hamiltonian subsystem (q 3 , p 3 ) is integrable reversible, let a ∈ R denote the action variable and θ ∈ R/Z the angle variable, and ω(a) the corresponding angular velocity, so that the canonical reversibility map (q 3 , p 3 ) → (q 3 , −p 3 ) becomes (a, θ) → (a, −θ). The system (11), with obvious abuse notation can now be writtenq Let so(3) denote the Lie algebra of skew symmetric 3 × 3 matrices. The structure of this system becomes evident if we introduce the vector variable u = ( √ k 1 m 1 q 1 , √ k 2 m 2 q 2 , p) and the so(3) valued function Then, in the variables (u, a, θ) ∈ R 3 × R × R/Z, system (12) takes the forṁ Our aim below is to show that equation (14) is integrable for each fixed value of a, which then proves that equation (11) is integrable. If ω(a) = 0, then B(a, θ) is a constant rotation matrix, so the system has invariant tori for these values of a. If ω(a) = 0 then we can rescale the time variable so that the angular velocity of θ is 1. The system (14) can then be writtenu where A a (θ) := B(a, θ)/ω(a).
That this system is integrable is a consequence of the following more general result. where (u, θ) ∈ R n × R/Z are the dependent variables and A : R/Z → so(n) is a smooth mapping. There is a smooth change of variables such that system (17) expressed in the new variables (v, θ) takes the forṁ v = Av for a constant matrix A ∈ so(n), whose spectrum σ(A) is such that Proof. One defines the flow map Φ(τ ) as the solution operator of the differential equation defined for all τ ∈ R by dw(τ ) dτ = A(τ )w(τ ) (21) with initial condition at τ = 0 -the initial time matters because the differential equation is not autonomous. This means that if w is a solution of (21), then and vice versa.
Since A(τ ) ∈ so(n) for all τ , the flow map after one period, i.e., Φ(1), is an element of SO(n). As a result, there exists a matrix A ∈ so(n), with the restriction (20) on the spectrum, such that Φ(1) = exp(A).
Recall that, one of the pillars of Floquet's theory is that for any integer n ∈ Z and any τ ∈ R we have ([3, § 28]) Now, v(u, τ ) is periodic in τ , of period one, because, due to the Floquet property (28), and the definition (23) of A, As a result, the mapping Ψ induces a mapping Ψ : R n × R/Z → R n × R/Z which has the desired properties.

Reversibility
We equip the space R n × R/Z with a linear involution R, defined as for a given linear orthogonal involution ρ. One verifies by a straightforward calculation that the differential equation (17) is reversible with respect to R if and only if Theorem 3.2. Assume that (31) holds, and that θ → A(θ) is periodic of period one. Assume further that the spectrum σ(A) of the matrix A defined in Theorem 3.1 fulfils Then the mapping Ψ : R n × R/Z → R n × R/Z, defined in Theorem 3.1, preserves reversibility, i.e., R • Ψ = Ψ • R.
Now, recalling the definition of A in (23), we notice that (34) implies that exp(−A) = ρ exp(A)ρ = exp(ρAρ) and, using the assumption on the spectrum of A, and the fact that the eigenvalues of ρAρ are the same as those of A, we deduce that −A = ρAρ. We therefore obtain exp(−Aτ )ρ = ρ exp(Aτ ).
Remark 3.3. The condition (32) is a non-resonance assumption. We already know that the matrix A is chosen to satisfy (20), but we cannot exclude the extreme case of one or more half rotations.

Wrapping up
We finish the proof by putting all the pieces together. Notice that with the canonical mapping ρ : R 3 → R 3 defined as (q 1 , q 2 , p) → (q 1 , q 2 , −p), the property (31) is fulfilled for the matrix A a defined in (16). We thus obtain that, in the special case of R 3 , the map Ψ defined in Theorem 3.1 transforms the original system to an integrable one, and the map Ψ is reversible. In R 3 , all rotations have one axis of rotation, so for all the values of the invariant a, we obtain two more invariants, and one angle variable. In the end, we thus obtain three invariants and two angle variables.