Physica D 156 (2001) 247–259
Nonlinear oscillations in the rolling motion of Euler’s disk A.A. Stanislavsky∗ , K. Weron1 Institute of Physics, Wroclaw University of Technology, Wyb. Wyspianskiego 27, Wroclaw 50-370, Poland Received 22 February 2001; received in revised form 17 April 2001; accepted 17 April 2001 Communicated by Y. Kuramoto
Abstract The fixed points of the dynamical system describing the rolling motion of a uniform disk or a uniform circular hoop on a rough horizontal plane without dissipation are analyzed. All fixed points of the system are stated to be steady motion states which are either saddle points or centers. The spectral analysis of the sound amplitude time series emitted by a toy of the rolling disk shows that the similar nonlinear oscillations are expected in the motion with dissipation. If under adiabatic approximation a finite-time singularity occurs, the bifurcation in the disk motion can appear. © 2001 Elsevier Science B.V. All rights reserved. PACS: 02.30; 05.45; 46 Keywords: Nonlinear oscillation; Euler’s disk motion; Fixed point; Stability analysis; Computer simulation
1. Introduction The rolling motion of a thin disk (called Euler’s disk [1]) on a horizontal plane is a classic problem of mechanics and has been treated by many authors [1–4]. If in the classic (nondissipative) theory the motion persists forever, the common experience shows that a circular disk (e.g., a coin) is spun upon a table, then it comes to rest. Recently, the identification of the dominant dissipative mechanism, for a given disk and a given surface, becomes the great interest [5–8]. There are a number of possible dissipative mechanisms for the rolling disk: rolling friction due to plastic deformation at the point of rolling contact, vibration of the supporting surface, viscous dissipation in the surrounding air, etc. To distinguish between these various possibilities, a careful analysis of experimental data is necessary. The point is that the dynamical system of Euler’s disk is nonlinear and its time evolution may be very surprising. In this paper we focus on the analysis of Euler’s disk motion without dissipation from the point of view of dynamical systems theory. The question “what place is taken by this dynamical system” is not only interesting itself, but may have a deeper raison d’être in the identification of the dominant dissipative mechanism for the rolling disk. As each dynamical system is characterized by its fixed (equilibrium) points, their comparison in the dissipative and nondissipative regimes of the system under consideration suggests a strategy for formulating the equations of ∗ Corresponding author. Permanent address: Institute of Radio Astronomy, 4 Chervonopraporna St., Kharkov 61002, Ukraine. E-mail addresses:
[email protected] (A.A. Stanislavsky),
[email protected] (K. Weron). 1 Fax: +48-7122-7751.
0167-2789/01/$ – see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 0 1 ) 0 0 2 8 1 - 0
248
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
motion in complex cases. We start, in Section 2, with the presentation of the equations of Euler’s disk motion without dissipation, and then show that they can be solved in separable variables which leads to the first integral of the system. In Section 3 we present the proof that all steady motion states are fixed points of the system, which form a surface in the three-dimensional space of system variables. We show that the dynamical system of Euler’s disk type has no other fixed points than steady states. The main result of Section 4 is that these fixed points can only be saddle points or centers. The properties of the fixed points do not change under linearization so that in Section 5 we can formulate the theorem characterizing the dynamical system of interest. The direct computer simulation of equations of Euler’s disk motion confirms our conclusions. We also have traced nonlinear oscillations in the recording of sound emitted by a toy of the rolling disk (with dissipation) commercially available. In Section 6 we discuss the Euler’s disk motion with dissipation. Such a motion can create a number of interesting phenomena like bifurcation.
2. The equations of motion A uniform disk rolling on a rough horizontal plane is a nonholonomic system, i.e., the dynamical system that has three degrees of freedom, but to define its configuration we need five Lagrangian coordinates. The procedure of finding the equations of motion from the Lagrangian function and equations of constraint is well known (for details, see [4]). We will not repeat the derivation here. Our starting point is these equations of motion themselves. In our consideration we will use McDonald’s notations [6]. As it is shown in Fig. 1, the axis zˆ is vertically ˆ 2, ˆ 3) ˆ are the moving axes (not fixed in the disk). Axis upwards. A right-handed coordinate triad of unit vectors (1, 1ˆ being along the symmetry axis of disk is chosen so that the component ω1 of the angular velocity vector ω of the disk about this axis is positive. Axis 3ˆ is directed from the center of the disk to the point of contact with the horizontal plane, and makes angle α to the horizontal, where 0 ≤ α ≤ π . Having the direction of the velocity of the point of contact, axis 2ˆ = 3ˆ × 1ˆ lies in both the plane of disk and the horizontal plane. The angular velocity of the moving axes about the vertical zˆ is denoted by Ω. In this case the motion of the disk (of radius a) rolling without slipping on a rough horizontal plane is described by the following equations: (2k + 1)ω˙ 1 + αΩ ˙ sin α = 0,
(1)
g kΩ 2 sin α cos α + (2k + 1)ω1 Ω sin α − (k + 1)α¨ = cos α, a
(2)
Ω˙ sin α + 2αΩ ˙ cos α + 2ω1 α˙ = 0,
(3)
Fig. 1. Sketch of the system: a disk of radius a rolls without slipping on a rough horizontal plane.
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
249
where g is the acceleration of gravity, and k is defined by the inertia tensor: k = 21 for a disk with mass m concentrated at the rim, k = 41 for a uniform disk, etc. As usual, any dotted symbol denotes the derivative with respect to time. According to Pars [4], we can derive the equations for ω1 and Ω as functions of α. Substituting dω1 /dt = (dω1 /dα)(dα/dt) and dΩ/dt = (dΩ/dα)(dα/dt), we can rewrite Eqs. (1) and (3) in the form (2k + 1)
dω1 + Ω sin α = 0, dα
(4)
dΩ sin α + 2Ω cos α + 2ω1 = 0. dα
(5)
Denoting p = cos α, we get dω1 1 = Ω, dp 2k + 1
(6)
d [(1 − p 2 )Ω] = 2ω1 . dp
(7)
Eliminating Ω, we obtain d 2 2 dω1 (1 − p ) ω1 = 0. − dp 2k + 1 dp
(8)
One can easily check that the latter equation is invariant to the change of sign of ω1 (see [4]). Eq. (8) is a differential equation of Legendre’s type; it determines ω1 as a function of p. The value of the coefficient 2/(2k+1) is 1 for a hoop √ and 43 for a disk. It should be noted that in this case, the index of the Legendre’s function is − 21 ±i (2k − 7)/(2k + 1) √ √ (for a disk it equals − 21 ±i 14/13, and a hoop: − 21 ±i( 3/2)). This means that ω1 is expressed in conical functions having the index − 21 +iλ. It follows from Eq. (6) that we can consider Ω as a derivative of ω1 , i.e., of the corresponding Legendre’s function. Eq. (2) can be written as kΩ 2 cos α + (2k + 1)ω1 Ω − V
dV (k + 1) g = cot α, dα sin α a
(9)
where V is the symbol substituting dα/dt, and d2 α/dt 2 = dV /dt = (dV /dα)(dα/dt) = V (dV /dα). In p-coordinates this equation takes the form kΩ 2 p + (2k + 1)ω1 Ω + (k + 1)V
g p dV = , dp a 1 − p2
where Ω and ω1 are the functions of p. Integrating Eq. (10), we have 1+k 2 k 2k + 1 2 g V + (1 − p 2 )Ω 2 + ω1 + 1 − p 2 = h, 2 2 2 a
(10)
(11)
where h is a constant, and the first and second terms 1 − p2 2k + 1 2 Ω 2 p dp = ω1 Ω dp = Ω − (2k + 1)ω12 , ω1 2 2 have been derived by using Eqs. (6) and (7). In α-coordinates, Eq. (11) reads 2k + 1 2 g 1 + k dα 2 k + sin2 αΩ 2 + ω1 + sin α = h. 2 dt 2 2 a
(12)
250
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
Expression (11) (or equivalently Eq. (12)) is the first integral of both Eq. (10) and the whole system of Eqs. (1)–(3). It shows the simple fact that the total energy of motion does not change with time, and the motion continues indefinitely. √ √ Since Ω and ω1 depend on cos α, Eq. (12) takes the form dα/dt = f (α), whose solution is t = dα/ f (α)+ C, where C is a constant. Therefore, the system of nonlinear differential Equations (1)–(3) may be successfully resolved in separable variables.
3. Steady motion states and fixed points As it is well known, for steady motion of this system we have α˙ = 0 and α¨ = Ω˙ = ω˙ 1 = 0. In this case Eqs. (1) and (3) become trivial, and Eq. (2) gives the following relationship between the constant, equilibrium values α0 , Ω0 and ω10 : kΩ02 sin α0 cos α0 + (2k + 1)ω10 Ω0 sin α0 =
g cos α0 . a
(13)
Let us check whether the steady states are fixed points of Eq. (9). Under the fixed point of a differential equation dv P (x, v) = , dx Q(x, v) it is understood [9] that such a point (x0 , v0 ) for which both functions P (x0 , v0 ) and Q(x0 , v0 ) take the value zero. If we write Eq. (9) in the form (k + 1)
dV kΩ 2 sin α cos α + (2k + 1)ω1 Ω sin α − (g/a) cos α = , dα V
(14)
the condition V = dα/dt = 0 leads automatically to ω˙ 1 = Ω˙ = α¨ = 0 so that the numerator of the fraction on the right-hand side of (14) also becomes zero. This means that the steady motion states are really fixed points of the system under consideration. Moreover, this dynamical system has only such fixed points.
4. Small oscillations about steady motion Now we suppose that α, Ω and ω1 undergo small oscillations about their equilibrium states: α = α0 + (α − α0 ) = α0 + α,
(15)
Ω ≈ Ω0 + a1 α,
(16)
ω1 ≈ ω10 + a2 α,
(17)
where a1 and a2 are some constants. Let α be small so that cos α ≈ 1 and sin α ≈ α. Inserting (16) and (17) in Eqs. (4) and (5) and equating terms of first order of smallness, we find the coefficients: a1 = −
2 (Ω0 cos α0 + ω10 ), sin α0
(18)
a2 = −
Ω0 sin α0 . 2k + 1
(19)
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
251
Using the equality (13) and the fact that ω10 = (b/a)Ω0 , where b is the distance from the axis of the instantaneous motion to the center of mass of the disk (see Fig. 1), we get V
dV dV =V = B(α − α0 ), dα d(α − α0 )
where g B= a
(20)
4k + 1 2k − 1 b b2 sin α0 − + cos 2α0 + 2(4k + 1) cos α0 + 2(2k + 1) 2 2 2 a a
is a constant for the given steady state. Let us observe that Eq. (20) is nothing else as a standard equation of linear oscillations. When B > 0, it has a saddle fixed point, and for B < 0, the fixed point becomes a center. As it is well known [9], the linear approximation of a nonlinear system is not sufficient for defining the global phase plane plot in the latter case. We shall show in the next section that all steady state solutions of our dynamical system are stable globally, and all states in the neighborhood of the steady state will oscillate around them, being attracted by ones and repelled from others. The nonlinear terms, which were neglected in the above stability analysis, do not change this fact. Clearly, for this dynamical system the period of small oscillations does not depend on the amplitude. When the oscillations have a large amplitude this will be wrong (see also the next section).
5. General analysis of oscillations Our analysis of the dynamical system (1)–(3) shows that, in fact, we investigate the motion being governed by the differential equation of type x¨ + f (x) = 0. In the general case, this equation describes free nonlinear oscillations [10] and may be written as dv −f (x) = , dx v
(21)
where v = dx/dt. Assume that the function f (x) is expanded in the terms of Taylor series around a fixed point x0 (in which v = 0 too): f (x) = b1 (x − x0 ) + b2 (x − x0 )2 + · · ·
(b1 = 0),
(22)
where b1 , b2 , . . . are some constants. Here, we suppose that b1 = 0. Otherwise the point x0 will be a fixed point of higher order. The integration of Eq. (21) gives 1 2 2v
+ F (x) = H,
where F (x) =
f (x) dx =
b1 b2 (x − x0 )2 + (x − x0 )3 + · · · , 2! 3!
and H is a constant. As it was noted above, the fixed points correspond to the system’s equilibrium. Owing to properties of Eqs. (21) and (22), we have dv dv −b1 (x − x0 ) + · · · . = = d(x − x0 ) v dx
(23)
Thus, for b1 < 0 the fixed point will be a saddle point, and for b1 > 0 it will be a center. Inclusion of the nonlinear terms does not change the character of fixed points because Eq. (21) has always the first integral, i.e., the total
252
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
energy does not depend on time (the kinetic energy is transformed in potential one and vice versa). In this case the oscillations cannot be damped. On the basis of the above consideration we can formulate the following assertion. Theorem. The steady motion states of the dynamical system Eqs. (1)–(3) are fixed points that may be only either saddle points or centers. Since α and V = dα/dt are functions of t, the closed trajectories on the phase plane (α, V ) will be related to periodic solutions. The period T of these oscillations may be expressed by the curvilinear integral
dα V
T =
(24)
taken along an integral curve with respect to increasing time. As V in Eq. (11) is squared, we may write T =2
αmax
αmin
√
dα , h1 − Ψ (α)
(25)
where Ψ (α) =
k 2k + 1 2 2g Ω 2 (α) sin2 α + ω (α) + sin α, k+1 k+1 1 a(k + 1)
(26)
and h1 = 2h/(k + 1) is defined by initial conditions from Eq. (12). Consequently, the period T of nonlinear oscillations of the system of interest depends on their amplitude. The examples of computer simulation for the Euler’s disk motion are presented in Figs. 2–4. They confirm the results obtained above. The disk’s radius a has been taken to be equal to 3.75 cm. We have integrated numerically the system of differential equations (1)–(3) using the Runge–Kutta formulas of fifth order with high accuracy (no problems with convergence of the numerical integration). As an example, we have taken the fixed point α0 = π/6, Ω0 = 40 and ω10 ≈ 1.7782. Fig. 2 presents the regime of small oscillations, while Fig. 3 shows the case when the nonlinearity is strong. In Fig. 4 one can see how the spectrum of the angular velocity changes for different initial conditions. The above results allow us to conclude that in the case of Euler’s disk motion with dissipation, similar oscillations are expected. To confirm our conclusion, we have analyzed the sound time series emitted during the rolling motion of a Tangent toy “Euler’s disk” [11]. We have observed the following. On the early stage, when the disk spins slowly, the rotation of the sound emission direction with respect to the microphone at rest is clearly detected; we hear sounds similar to flips. When the angular velocity of the disk motion increases, the flips unite in the uniform sound in which the nonlinear oscillations can be found. The final stage of motion is characterized by a shudder and a whirring sound of rapidly increasing frequency (Ω increases). The sound time series was long, about 12.3 s. On the interval between 9 and 11 s we have observed that the sound spectrum consists of a number of equidistant “harmonics” (Fig. 5). The only explanation of their appearance is the nonlinear oscillations of the disk motion. In Fig. 5, one can see that the angular velocity Ω increased. Since the oscillation is nonlinear, the time evolution of its “harmonics” is complex (see also Fig. 4). It should be noted that the result presented in Fig. 3 of the McDonald’s paper [6] may be interpreted as a signature of nonlinearity in the disk motion. The oscillations have the same amplitude in the interval between 7 and 7.2 s, and their form is characterized by sharp peaks and dished bottoms.
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
253
Fig. 2. Computer simulation of Euler’s disk motion: (a) time series of α, Ω, ω1 ; (b) phase portrait (α, dα/dt). Initial conditions: α(0) = π/6, α(0) ˙ = 0, Ω(0) = 40, ω1 (0) = 1.82 (α˙ = dα/dt).
254
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
Fig. 3. Computer simulation of Euler’s disk motion: (a) time series of α, Ω, ω1 ; (b) phase portrait (α, dα/dt). Initial conditions: α(0) = π/6, α(0) ˙ = 0, Ω(0) = 40, ω1 (0) = 10 (α˙ = dα/dt).
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
255
Fig. 4. Spectrum of the angular velocity Ω for different initial conditions: for top plot, α(0) = π/6, α(0) ˙ = 0, Ω(0) = 40, ω1 (0) = 2; for ˙ = 0, Ω(0) = 40, ω1 (0) = 2 (α˙ = dα/dt). central plot, α(0) = π/96, α(0) ˙ = 0, Ω(0) = 40, ω1 (0) = 2; for bottom plot, α(0) = π/384, α(0)
6. Dissipation of energy in Euler’s disk motion Eqs. (1)–(3) describe free (undamped) oscillations of the system with a nonlinear restoring force. In any real situation, the system will be dissipative. How to include a damping force in these equations? Since the system is similar to free nonlinear oscillations of type x¨ + f (x) = 0, one could consider the following form of damped oscillations x¨ + φ(x) ˙ + f (x) = 0
(xφ( ˙ x) ˙ ≥ 0).
The condition xφ( ˙ x) ˙ ≥ 0 is necessary for the resistance force φ to be really a damping force. Taking x˙ = v as a new variable, the above equation can be written in the form dv −f (x) − φ(v) = , dx v
(27)
in which the term φ(v) does not allow us to separate the variables. However, some conclusions about the behavior of the system can be formulated. Clearly, the damping force φ(x) ˙ has to be zero for x˙ = 0. This means that Eq. (27) has the same fixed points as Eq. (21) (at least when x˙ = 0), and the time evolution of such a damped system will direct to one of the fixed points. In any case the system will be in a fixed point exactly, when t → ∞. The fixed (equilibrium) points of Eq. (9) form a surface in the three-dimensional space (α, Ω, ω1 ). From common experience we know that if the Euler’s disk is
256
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
Fig. 5. Experimental data (normalized sound amplitude time series and corresponding spectra).
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
257
Fig. 5. (Continued).
a dissipative system, it has only one point of stable equilibrium (disk lies on the table at rest). Thus, the dissipative processes change the configuration of fixed points of the dynamical system (1)–(3). In this case one can suppose the following scenario of dissipation. The viscous friction (depending on α) ˙ lands the system’s evolution on the surface (13) and then another (or others) dissipative mechanism drives it on the surface to α = 0.
258
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
If one assumes that α is a slowly varying function of time (|α| ˙ Ω [5]), the energy of the motion E of the disk of mass m (the sum of the kinetic energy, 21 mΩ 2 [ 23 b2 + 41 a 2 sin2 α], and the potential energy, mga sin α) satisfies the equation dE/dt = −Φ, where Φ is the rate of dissipation. Using Eq. (13) and taking a small value of α, one can write the energy E as B 1 , (28) E ≈ mga Aα + 2 α where A = 2 + 1/(6(b/a) + 1), B = 6(b/a)2 /(6(b/a) + 1). Obviously, for b = 0 Eqs. (28) and (13) coincide with Eqs. (1) and (2) in [5]. In our case (and in [5]) the production Ω 2 α is constant so that when the angle α decreases, Ω increases. Considering the dissipation in the thin layer of air between the disk and the table, the rate of dissipation is Φ ≈ π µga2 /α 2 , where µ is the viscosity of the surrounding air (Φ → ∞ as α → 0) [5]. Then we have 1 B dα πµga2 . (29) mga A − 2 ≈− 2 dt α α2 The integration of this equation gives t0 − t A 3 , α − Bα = 2π 3 t1
(30)
where t1 = m/µa, t0 is a constant of integration determined by the initial condition: t0 = (Aα03 /3 − Bα0 )(t1 /2π ), where α = α0 at t = 0. Eq. (30) is cubic and has three roots. They can be either real or complex. From the physical point of view, the angle α takes only nonnegative values. Therefore, the complex and negative roots of (30) are unphysical. For 0 ≤ t < t0 the equation has only one positive solution. At t = t0 we obtain two nonnegative roots. One of them, α → 0, is the stable equilibrium state that is realized only at t = t0 . Let us stress that it is the only √ stable state in the investigated dynamical system. Another root, α = 3B/A, is unstable. When the disk reaches α , the cusp catastrophe, well known in physics [12], occurs. For this reason at t = t0 the disk leaps from α to zero omitting intermediate values. Then the disk stops in a trice because of dry friction. The behavior described by Eq. (30) persists until α is reached. At α = 0 the mechanism of energy dissipation for Euler’s disk changes: dry friction makes the dominant contribution, so any free rolling of the disk on the table ends. The cause of bifurcation is the finite-time singularity of Euler’s disk motion with dissipation [5]. Of course, such a singularity cannot be realized in practice. On the other hand, bifurcation may be one of physical effects which prevents its occurrence. Observe that in the final stage of Euler’s disk motion the whirring sound breaks off suddenly. The disk seems to float for a moment, then jumps down on the horizontal plane. Although this effect may be explained also by the upward contact force having vanished [6], this does not exclude the above-mentioned bifurcation. In practice, they can concur. If the rate of energy is proportional to a power of α with the exponent being less than 1, then a finite-time singularity (for which Ω becomes infinite) under the adiabatic approximation will occur [8]. Consequently, the bifurcation will be possible too. For a toy “Euler’s Disk” commercially available in Tangent Toys (Sausalito, CA) [11], m = 400 g, and a = 3.75 cm. With these values and with µ = 1.78×10−4 g cm−1 s, the characteristic constant t1 = m/µa ≈ 0.6×106 s. √ If we take α0 = 0.1(≈ 6◦ ) and b/a = 0.002, we find t0 ≈ 95 s. At α = 3B/A ≈ 0.5 × 10−2 , Ω ≈ 460 s−1 . From α0 to α , the adiabatic approximation is well satisfied. Moffatt took in his calculations the same values of a, α0 , m, µ as we did. Assuming b = 0, he has found t0 ≈ 100 s. Then, Moffatt has written in [5] the following: “This is indeed the order of magnitude (to within ± 20%) of the observed settling time in many repetitions of the spinning of the disk (with quite variable and ill-controlled initial conditions), that is, there is no doubt that dissipation associated with air friction is sufficient to account for the observed behaviour.” But the experiments performed up to now show that the process is about 10 times shorter. McDonald has mentioned that the total time of spin of a
A.A. Stanislavsky, K. Weron / Physica D 156 (2001) 247–259
259
Tangent Toy “Euler’s Disk” lies between 7.26 and 7.28 s [6]. Van den Engh et al. [7] have written that “coins in vacuo spun on average for 12.5 s; coins in air spun on average for 10.5 s (average of 10 observations each)”. The results clearly indicate that viscosity dissipation in the surrounding air is not the dominant dissipative mechanism for the rolling disk. Unfortunately, at present the Moffatt’s adiabatic approximation is the only efficient method for theoretical analysis of the disk motion with dissipation, and another approach is not known. Perhaps, if one uses another dissipative mechanism (rather than air friction) with another rate of the energy dissipation in the adiabatic approximation, the correspondent quantitative agreement will be obtained.
7. Discussion We have shown that the system of nonlinear differential equations describing Euler’s disk motion without dissipation may be transformed to the form which allow us to obtain closed differential equations for each of values ω1 , Ω, α separately. In other words, the task can be reached in separable variables. Using the dynamical systems theory (in particular, its stability analysis) we have stated that steady states of the dynamical system are fixed points. The fixed points of this system have been found to be only either saddle points or centers. When we depart from them (e.g., by choosing the initial conditions), the nonlinear free oscillations are, in general, obtained. The period of these nonlinear oscillations depends on their amplitude. It is very likely that nonlinear oscillations will appear in the case of Euler’s disk motion with dissipation. Our analysis of the sound time series recorded experimentally does not confirm that dissipation associated with air friction is sufficient to account for the observed behavior of the rolling disk. It is clear that the dynamical systems (of Euler’s disk type) without and with dissipation have different fixed (equilibrium) points. Since in the limit of very small dissipative effects the surfaces of fixed points for the two cases must coincide, the formulation of a more comprehensive mathematical model of Euler’s disk motion with dissipation and its analysis are of great interest for a future work.
Acknowledgements A.A. Stanislavsky acknowledges the NATO Science Fellowships Programme for support, he is a postdoctoral fellow at the programme. References [1] L. Euler, Theoria Motus Corporum Solidorum Seu Rigidorum, Greifswald, 1765. [2] E.J. Routh, The Advanced Part of a Treatise on the Dynamics of a System of Rigid Bodies, 6th Edition, Macmillan, London, 1905 (reprinted by Dover, New York, 1955). [3] E.A. Milne, Vectorial Mechanics, Interscience, New York, 1948. [4] L.A. Pars, Treatise on Analytical Dynamics, Heinemann, London, 1965. [5] H.K. Moffatt, Nature 404 (2000) 833. [6] K.T. McDonald, physics/0008237 (http://xxx.lanl.gov). [7] G. van den Engh, P. Nelson, J. Roach, Nature 408 (2000) 540. [8] H.K. Moffatt, Reply, Nature 408 (2000) 540. [9] D.K. Arrowsmith, C.M. Place, Ordinary Differential Equations: A Qualitative Approach with Applications, Chapman & Hall, London, 1982. [10] J.J. Stoker, Nonlinear Vibrations in Mechanical and Electrical Systems, Interscience, New York, 1950. [11] J. Bendik, The Official Euler’s Disk Website (http://www.eulerdisk.com/), Tangent Toy Co., Sausalito, CA (http://www.tangenttoy.com/). [12] T. Poston, I. Steward, Catastrophe Theory and its Applications, Pitman, London, 1978.