Off[General::spell1]

A Mechanical System with Constraints

The solution of a well-known problem in Classical Mechanics involves some unexpected subtleties, illustrating the often surprising complexity of constrained systems.

by Tim Pritchett

The marriage between Mathematica and physics is a very natural one, and recent years have seen publication of a number of books describing applications of Mathematica to a range of problems in theoretical physics.  Classical Mechanics, arguably the oldest branch of physics, is particularly well served by the means available in Mathematica for both symbolic and numerical calculation as well as for ready visualization of results.  This article describes how Mathematica was used to investigate a well-known problem from Classical Mechanics. The problem, described by Littlewood [Bollabás 1986] and recently revisited by Tokeida [1997], involves an asymmetrically weighted, rigid, circular hoop which initially rolls without slipping along a flat surface. For certain initial conditions, such a hoop is observed at a certain point to lose contact with the supporting surface and to hop into the air.  A stroboscopic photograph of the hop appears in [Pritchett 1999].

The Problem

Figure  depicts the mechanical system which is the subject of this article:  a perfectly rigid circular hoop of mass (1 - λ) M and radius a with an object of mass λM rigidly attached to its rim.  Here, M is the total mass of the combined system consisting of the hoop and the object.  The parameter λ is the ratio of the mass of the attached object to the total; of particular interest will be the limit λ->1, in which the mass of the hoop itself is negligible relative to that of the attached object.  The center of mass of the combined hoop-object system is indicated by the small dot in Figure ; it lies along the line connecting the center of the  hoop  and the attached object at a distance λ a from the center of the hoop.

Code for Figure 1  (DO NOT PRINT!)

\(\(\(Module[{a = 1, λ = 0.7, θ = 50\ Degree, rCM, \ 
          rCtrHoop, \ rMass, \ hoop, \ floor, \ diagramLabels}, 
	

        rCM = a\ {θ + λ\ Sin[θ], 
              1 + λ\ Cos[θ]}; 


        rCtrHoop = rCM - λ\ a\ {Sin[θ], Cos[θ]}; 


        rMass = rCtrHoop + \((a - \ hoopThickness/2)\)\ {Sin[θ], 
                Cos[θ]}; 
hoopThickness = 0.05; 


        hoop = Graphics[
		{{GrayLevel[0.7], 
                Disk[rCtrHoop, a]}, 
				Circle[rCtrHoop, 
                a], {GrayLevel[1], Disk[rCtrHoop, a - hoopThickness]}, 
              Circle[rCtrHoop, a - hoopThickness], {AbsolutePointSize[10], 
                Point[rMass]}, {AbsolutePointSize[5], 
                Point[rCM]}, {Dashing[{0.02, 0.01}], 
                Line[{rCtrHoop + {0, a}, rCtrHoop}]}, {Dashing[{0.02, 0.01}], 
                Line[{rCtrHoop, rMass}]}, 
			{{GrayLevel[1], \ 
                  Disk[rCtrHoop, \ rr = 0.025]}, \ 
                Circle[rCtrHoop, \ rr]}}]; 


        floor = Graphics[Line[{{0, 2.2\ a}, {0, 0}, {2\ a, 0}}]]; 


        diagramLabels = 
          Graphics[{{Dashing[{0.02, 0.01}], 
                Line[{rCM + {0, a/4}, rCM - {0, a/2}}]}, {Circle[rCM, 
                  a/4, {\(-\((π/2 + θ)\)\), \(-π\)/2}], 
                Circle[rCtrHoop, a/4, {π/2 - θ, π/2}]}, 
              Text[StyleForm["λ M", FontSize -> 12, 
                  FontFamily -> "Times", FontSlant -> Italic], 
                rMass + \ a/10\ {Sin[θ], Cos[θ]}, {\(-1\), 0}], 
              Text[StyleForm["θ", FontSize -> 12, 
                  FontFamily -> "Times", 
                  FontSlant -> "Italic"], 
                rCM - 0.35\ a\ {Sin[θ/2], Cos[θ/2]}]}]; 	

        Show[hoop, floor, diagramLabels, Axes -> True, 
          Ticks -> False, AxesLabel -> {"x", "y"}, 
          AspectRatio -> Automatic]];\)\(\ \)\)\)

Both the hoop and the surface upon which it rolls are assumed to be perfectly rigid; in other words, neither the hoop nor the surface suffers deformation of any kind as the former rolls along the latter.  At least initially, the hoop rolls entirely without slipping.  The problem is then to determine the subsequent motion of the hoop as a function of time t.

Mathematical Formulation

The most general motion of a rigid body consists of a translation of the center of mass, combined with a  rotation  of the body about an axis containing the center of mass. Accordingly, we  shall specify the motion of the hoop by the Cartesian coordinates {X[t],Y[t] } of the center of mass of the system and by the angle θ[t] through which the hoop has turned.

Of course, X[t], Y[t], and θ[t] are not the only variables we could choose to specify the configuration of the weighted hoop. Indeed, instead of the horizontal and vertical  coordinates of the center of mass of the system, it might seem more  natural to choose the horizontal and vertical coordinates of the center of the hoop, quantities that we will denote by (Subscript[hoop, x])[t] and (Subscript[hoop, y])[t], respectively.  Since the hoop is assumed to be perfectly rigid, the two sets of coordinates are related by

(Subscript[hoop, x])[t_] := X[t] - λ a Sin[θ[t]]\(hoop\_y[t_] := Y[t] - λ\ a\ Cos[θ[t]]
\)hoop[t_] := {(Subscript[hoop, x])[t], (Subscript[hoop, y])[t]}

We will see presently that dynamics of the system make X[t], Y[t], and θ[t] a particularly convenient choice of coordinates.

We will later have need of symbolic expressions for the velocity and acceleration of the center of the hoop:

hoopV[t_] := Evaluate[D[hoop[t], t]](Subscript[hoopA, y])[t_] := Evaluate[D[(Subscript[hoop, y])[t], {t, 2}]]

The coordinate system shown in Figure  is chosen so that at t = 0 the center of the hoop is situated at {0, a}. In general, three  parameters are required to specify the initial state of the weighted hoop: one to specify the initial orientation of the hoop, and two more to  completely specify its initial motion.  One very natural choice would be to specify the angle (Subscript[θ, 0]) = θ[0], the initial translational speed (Subscript[hoopV, x])[0] of the center of the hoop, and the initial angular speed of rotation \(TraditionalForm\`ω\_0 = θ\&.[0]\). In the present case, however, where we assume that the hoop initially rolls without slipping, the latter two parameters are not independent, but are related by  (Subscript[hoopV, x])[0] = a (Subscript[ω, 0]). We could thus specify the initial state of the weighted hoop by giving (Subscript[θ, 0]) and (Subscript[ω, 0]); in fact, it will prove to be convenient to use (Subscript[θ, 0]) and the dimensionless quantity c=a ((Subscript[ω, 0])^2)/g, where g=9.8 m/(s^2) is the constant acceleration imparted by the earth’s  gravity to every object near the earth’s surface.

As a practical matter, we must ensure that the Mathematica kernel recognizes (Subscript[θ, 0]) not as a subscript, but as a single symbol [Zeller 1997].

Needs["Utilities`Notation`"]

Symbolize[(Subscript[θ, 0])]

The Physics

As we noted above, the most general motion of a rigid body consists of a translation of the body's center of mass plus a rotation of the body about
an axis through its center of mass.  Newton's Second Law delivers the following
pair of second-order differential equations describing the translational motion.

(X^″)[t] = f[t]

(Y^″)[t] = n[t] - g

Here, g, n[t], and f[t] represent, respectively, the downward force per unit mass (acceleration) due to gravity and the horizontal and vertical components of the contact force per unit mass due to the surface on which the hoop rolls.   The vertical component of the contact force, which in this case acts normal to the supporting surface, is accordingly often referred to as the "normal force"; the horizontal component, which acts parallel to the surface, is friction.  We observe that if the hoop loses contact with the supporting surface, then  n[t] = f[t] = 0, and the above equations describe the characteristic parabolic trajectory of an object falling freely under the influence of gravity.

The contact forces n[t] and f[t] which appear in equations () and () are not known a priori, but must be determined by the effects they produce on the motion of the system; they are forces of constraint.  The friction force f[t]  is responsible for maintaining the constraint that the hoop rolls without slipping, while the normal force n[t] ensures that the vertical height of the center of the hoop, (Subscript[hoop, y])[t], is always greater than or equal to a, the radius of the hoop.  (Equality holds as long as the hoop and the supporting surface are in contact.)  Indeed, () and ()   may be thought of as defining  n[t] and f[t]  given the constrained trajectory of the system center of mass.

In order to describe the rotational motion of the hoop-object system, we must know the system's moment of inertia about its center of mass.  The entire mass of the hoop is generally taken to be concentrated in the rim.  For the case in which the attached object is a point mass — a mathematical idealization — the moment of inertia of the system is M (a^2) (1 - (λ^2)).  Of course, in the real world the attached object will not be a mathematical point, but will be of nonvanishing physical dimension,  and thus will possess a nonvanishing moment of inertia about its centroid. We will take

inertia = M (a^2) (1 - (λ^2) +  λ ϵ);

where ϵ is a dimensionless parameter (generally quite  small)  related to the physical dimensions of the attached object.

The time evolution of the angular variable θ[t] is governed by the following equation, obtained by considering the torques produced by the contact forces n[t] and f[t].  (The force per unit mass g due to gravity acts at the center of mass of the system, and thus can produce no torque about the center of mass.)

(θ^″)[t] = 
      (M/inertia) 
        (n[t] a λ Sin[θ[t]] - 
            f[t] a (1 + λ Cos[θ[t]]))

Phase 1:  Rolling Without Slipping

As long as the hoop and the supporting surface are both perfectly rigid and the former rolls over the latter entirely without slipping, the trajectory of the center of mass of the hoop-object system is a curtate cycloid:

(Subscript[cycloid, x])[θ_] := 
      a (θ - (Subscript[θ, 0]) + λ Sin[θ])
(Subscript[cycloid, y])[θ_] := a (1 + λ Cos[θ])
cycloid[θ_] := {(Subscript[cycloid, x])[θ], 
        (Subscript[cycloid, y])[θ]}
{X[t_], Y[t_]} = cycloid[θ[t]];

and the mechanical energy of the system,

energy=(1/2) M (((((∂_t) X[t]))^2) + ((((∂_t) Y[t]))^2))+(1/2) inertia ((θ'[t])^2)+M g Y[t]//Simplify

is a conserved quantity — a constant — whose value may be determined from the initial conditions.

initialEnergy=energy/.{θ[t] -> (Subscript[θ, 0]),((θ'[t])^2)->g c/a}//Simplify

The principle of energy conservation thus provides us with the following expression for θ'[t], valid as long as the cycloid constraint remains in force i.e., as long as the (perfectly rigid) hoop remains in contact with the (perfectly rigid) supporting surface and experiences no slipping.

\(θ\& . \)=FullSimplify[MapAll[Expand,Simplify[θ'[t]/.Solve[
			energy == initialEnergy,θ'[t]][[2]]]]//.\({\@a_\/\@b_ :> \@\(a\/b\), \(\@a_\) \@b_ :> \@\(a\ b\)}\)]

To solve for the motion of the system, we have only to integrate the single first-order differential equation θ'[t] = \(θ\&.\).  In the special case λ = 1, ϵ = 0, an analytic solution is possible (see below); in the general case, the equation must be integrated numerically.

For later use, we differentiate the above expression for θ'[t] to obtain an expression for the angular acceleration (θ^″)[t]:

(Overscript[θ, ¨])=Simplify[\(D[θ\&., t]\)/.θ'[t]->\(θ\&.\)];

The resulting Mathematica output is then:

((Sec[(θ/2)] (Sqrt[g (c + Sin[(θ/2)]^2)]))/(Sqrt[a]))

Forces of Constraint

As we noted above, the forces of constraint f[t] and n[t] appearing in (), (), and () are not known a priori, but must be determined by the effects they produce on the motion of the system.  We may obtain expressions for the forces of constraint required to keep the center of mass moving along a cycloid by substituting the parametric equations for the cycloid into the equations of motion  () and () of the center of mass of the system.

fCyc[ξ_]:=Evaluate[Simplify[D[(Subscript[cycloid, x])[θ[t]],  {t, 2}]/.θ''[t]->(Overscript[θ, ¨])/.θ'[t]->\(θ\&.\) /.
			θ[t] -> ArcCos[ξ]]]
fCyc[ξ]

nCyc[ξ_]:=Evaluate[Together[g + D[(Subscript[cycloid, y])[θ[t]], {t, 2}]/.θ''[t]->(Overscript[θ, ¨])/.θ'[t]->\(θ\&.\)/. θ[t] -> ArcCos[ξ]]]
		nCyc[ξ]

Onset of Slipping

The hoop rolls without slipping until |  f[t] |, the magnitude of the force required to maintain the kinematic constraint that the center of mass of the system move along a curtate cycloid, exceeds the maximum value possible for static friction, (Subscript[μ, s])  | n[t]  | .  The quantity (Subscript[μ, s]), the coefficient of static friction, depends on the nature of the surfaces involved; for most substances, (Subscript[μ, s]) lies in the range 0.1 to 1.5.

\(angleOnsetOfSlipping[initialAngle_, \ initialAnglrSpeed_] := 
      Module[{slipEqn, angles0ToPi}, 
	
        slipEqn = 
          Simplify[
            Numerator[fCyc[ξ]]\^2\  == \ 
                \(μ\_s\^2\) Numerator[nCyc[ξ]]\^2 /. {
                θ\_0 -> initialAngle, \ 
                c -> \((a/g)\)*\((initialAnglrSpeed)\)\^2}]; 
		
        angles0ToPi = ArcCos[Cases[ξ /. NSolve[slipEqn, ξ], _Real]]; 
        
Min[Select[
            Union[angles0ToPi, 2\ π\  - \ angles0ToPi, 
              2\ π\  + \ angles0ToPi, 
              4\ π\  - \ angles0ToPi], # ≥ initialAngle
                \ &]]]\)

Analytic Solution for a Special Case

Up until the onset of slipping (i.e., for θ < (Subscript[θ, Slip]) ), the motion of the system is given by a single first-order differential equation: \(TraditionalForm\`\(θ\^′\)[t] = θ\&.\), with \(TraditionalForm\`θ\&.\) as determined above.  In the case of a massless hoop (λ = 1) with attached point mass (ε = 0), this differential equation may be integrated analytically.  In this section, we use Mathematica to obtain a symbolic solution for this special case.

Clear[g, a, λ, 
        ϵ, c];

With λ = 1, ε = 0, and (Subscript[θ, 0]) = 0, the expression for \(TraditionalForm\`θ\&.\) on the cycloid reduces to

\(dThetaDt = 
      Simplify[Simplify[
              θ\&. /. {λ -> 1, \ ϵ -> 0, 
                  θ\_0 -> 0, \ θ[t] -> θ}] /. {
              \(-Cos[θ_]\) -> 2\ Sin[θ/2]\^2 - 1, 
              1 + Cos[θ_] -> 2\ \ Cos[θ/2]\^2}, \ Trig -> False]
           // PowerExpand\)

where we have made use of the trigonometric identities 1 + Cos[θ] = 2  (Cos^2)[θ/2] and 1 - Cos[θ] = 2  (Sin^2)[θ/2]. The above expression is separable, so we can straightforwardly  integrate.

integratedEqn = 
      PowerExpand /@ ((∫_0^(θ[t]))\[DifferentialD](θ/dThetaDt) == t)

For (Subscript[θ, 0])  ≠  0, the left-hand side involves elliptic functions.

eqnSimpler = 
      Exp /@(((Sqrt[g])/(2 (Sqrt[ a]))) #& /@ integratedEqn) // Simplify

sinThetaOver2 = 
      Sin[θ[t]/2] /. 
        Solve[eqnSimpler, Sin[θ[t]/2]][[1
          ]]

simplerSinThetaOver2 = 
      PowerExpand[Simplify[ExpToTrig[sinThetaOver2]]]

θ[t] /. 
      Solve[Sin[θ[t]/2] == simplerSinThetaOver2, θ[t]]
          [[1]]

The above expression for θ[t] is valid as long as the hoop remains in contact with the ground and  rolls  without slipping; in other words, as long as the attached object is constrained to move along a cycloid.

Numerical Solution

In the general case, the differential equation \(TraditionalForm\`\(θ\^′\)[t] = θ\&.\) governing the motion of the system for θ ≤ (Subscript[θ, Slip]) must be integrated numerically.  That task is accomplished by the following module.  (The parameters a, g, λ, ϵ, and (Subscript[μ, s]) are assumed to have numerical values.)

rollingWithoutSlipping[initialAngle_, initialAnglrSpeed_]:=Module[{thetaDot, theta},
		
  

UnSet[(Subscript[θ, Slip])];
(Subscript[θ, Slip]) := 
                (Subscript[θ, Slip]) = 
                  angleOnsetOfSlipping[initialAngle,  initialAnglrSpeed];
			
  
	
\(thetaDot\  = \ 
                θ\&. /. {θ\_0 -> initialAngle, \ 
                    c -> \(a\/g\) \((initialAnglrSpeed)\)\^2}\);
theta=θ/.NDSolve[{θ'[t]==thetaDot,θ[0] == initialAngle},θ,{t, 0, 500},StoppingTest -> 
                          (θ[t] > (Subscript[θ, Slip]))][[1]];
		
(*   Find the time (Subscript[t, Slip]) at which the hoop reaches the critical angle for the onset of slipping.   *)
				
	(Subscript[t, Slip]) = 
                t /. FindRoot[
                      theta[t] == (Subscript[θ, Slip]), {t, 
                        theta[[1, 1, 2
                          ]]}][[1
                    ]];
			
θSoln[tt_ /; 0 ≤ tt ≤ (Subscript[t, Slip])]
                   := theta[tt];
CMSoln[tt_ /; (0 ≤ tt ≤ (Subscript[t, Slip]))]
                   := cycloid[theta[tt]] /. (Subscript[θ, 0]) -> initialAngle;
	
  

{(Subscript[X, Slip]), (Subscript[Y, Slip])} = 
                cycloid[(Subscript[θ, Slip])] /. (Subscript[θ, 0]) -> initialAngle;
\({X\&.\_Slip, Y\&.\_Slip} = 
                D[cycloid[theta[t]], t] /. t -> t\_Slip\);
\(θ\&.\_Slip = 
                thetaDot /. θ[t] -> θ\_Slip\);			
	]

Phase 2:  Rolling With Slipping

Once the point of contact between the hoop and the supporting surface begins to slip, the magnitude of the friction force is known experimentally to be proportional to the magnitude of the normal force:   | f[t]|= (Subscript[μ, k]) |  n[t]  | .  For a given surface, the constant of proportionality, the coefficient of kinetic friction (Subscript[μ, k]), is always less than the coefficient of static friction (Subscript[μ, s]).  This is because the "cold-welding" which serves to bind the surfaces together in the static case cannot occur when there is relative motion between the two surfaces.

The direction of f[t]  is such that it opposes the relative motion of the two surfaces

fSign[hoopVx_, ω_] := Sign[a ω  -  hoopVx]

The normal force  n[t]  ≥ 0  must be sufficiently large to ensure that the height of the center of the hoop above the supporting surface is never less than the radius of the hoop, i.e., that the point of contact never lies below the x-axis.  As we explain below, this constraint is used to dynamically determine the value of n.

To obtain the motion of the system during this intermediate phase in which the hoop rolls with slipping but has not yet lost contact with the supporting surface, one must simultaneously solve equations () and (), subject to the constraint that the point of contact remain at or above the x-axis.  Since f[t] = fSign (Subscript[μ, k]) n[t], both of these differential equations depend on the normal force n[t], and n[t] must be dynamically determined from the constraint.  It is for this reason that one cannot simply use NDSolve to numerically integrate the pair of differential equations governing this phase of the motion; rather, one is forced to adopt a more complicated numerical integration scheme such as the one which we describe below.

Since our scheme employs a fixed time step Δt, it will be convenient in this second phase of the solution to parameterize the  evolution  of the system in terms not of the time t, but rather in terms of the integer variable i, which counts the number of time steps completed. We define new  variables  relevant to this slipping phase: x[i] and y[i] for the horizontal and vertical coordinates of the center of mass, vx[i] and vy[i] for the horizontal and vertical components of the center of mass  velocity, and th[i] and thDot[i] for the rotation angle θ and its derivative.

Clear[X, Y];

slipCoords = {X -> x, Y -> y, θ -> th };
slipVelocities={X'->vx,Y'->vy, θ'->thDot};
CMEqns={X''[i]->(Subscript[μ, k]) n fSign[xctr[i], thDot[i]], Y''[i]->n - g};
RotationEqn= {θ''[i]->(M/inertia) 
              (n a λ Sin[θ[i]] - 
                  (Subscript[μ, k]) n fSign[xctr[i], thDot[i]] a 
                    (1 + λ Cos[θ[i]]))};

We express the position, velocity, and vertical acceleration of the center
of the hoop in terms of these new variables:

\({xctr[i_], yctr[i_]} = hoop[i] /. slipCoords; 
    
{vxctr[i_], vyctr[i_]} = \(hoopV[i] /. slipVelocities\) /. slipCoords; 
    
ayctr[i_] = 
      \(\(\(hoopA\_y[i] /. RotationEqn\) /. CMEqns\) /. slipVelocities\) /. 
          slipCoords // Simplify; \)

At each time step, we compute a provisional value for the height of the center of the hoop, assuming that the relevant acceleration is constant over the step:  \(TraditionalForm
      \`yctr[i\  + \ 1] = \ 
        yctr[i]\  + \ vyctr[i] Δt\  + \ 
          ayctr[i] Δt\^2/2\).  Using this provisional height, we then check explicitly whether the constraint is violated: yctr[i + 1]<a.  If it is, the normal force n to increased to give yctr[i  +  1]  =  a for the time step, and, using the revised value of n,

nRevised = 
      Simplify[n /. 
          Solve[a == 
                yctr[i]  +  vyctr[i] Δt  +  
                  ayctr[i] Δ(t^2)/2, n][[1]]];

we compute the new values x[i + 1], y[i + 1], vx[i + 1], vy[i + 1], th[i + 1], and thDot[i + 1], assuming that the corresponding accelerations are constant over a time step.  For x, vx, y, and vy, whose accelerations depend only on n (which we update only once every time step), this is strictly true, and it is a very good approximation for th. The procedure continues until yctr[i + 1]  is strictly greater than a. Of course, since all of these comparisons involve floating point  values,  they can be done only to a given tolerance.  Thus, the termination condition actually takes the form yctr[i]-a ≤ tolerance.

To be able to start the integration procedure for phase two of the motion, we require the following values at the onset of slipping (i = 0):  x[0], y[0], th[0], vx[0], vy[0], and thDot[0] — all of which are computed in the module rollingWithoutSlipping.  We also require the value of the normal force n at the onset of slipping, as well as the value  of  the vertical acceleration of the center of mass one time step prior to the onset of slipping.

slipping[X0_,  Y0_, theta0_, Vx0_, Vy0_, thetaDot0_, AyPrior_, n0_]:=(
		
  

tolerance = 0.01;

  
  
	x[0] = X0;
  y[0] = Y0;
  th[0] = theta0;
	vx[0] = Vx0; 
	vy[0] = Vy0;
  thDot[0] = thetaDot0;

  

n = n0;
ay[-1] = AyPrior;

(*  This For loop performs the integration procedure.  *)

For[i = 0,yctr[i]  - a  ≤  tolerance,i++,
	
	  
	 
	If[yctr[i] + vyctr[i] Δt  +  
                      ayctr[i] Δ(t^2)/2  <  a, 
                  n = Max[0, nRevised]];
	
	(* Recalculate current accelerations to reflect the change in n. *) 
	
	ax[i] = (Subscript[μ, k]) n fSign[vxctr[i], thDot[i]];
	ay[i] = n - g;
	thDblDot[i] = 
                  (n/(a 
                          (1 - (λ^2) + λ ϵ))) 
                    (λ Sin[th[i]] + 
                        (Subscript[μ, k]) (1 + λ Cos[th[i]]));
	
    
	
	vx[i + 1] = vx[i] + ax[i]*Δt;
	vy[i + 1] = vy[i] + ay[i]*Δt;
	thDot[i + 1] = thDot[i] + thDblDot[i]*Δt;
	x[i + 1] = 
                  x[i] + (vx[i] + vx[i + 1])*Δt/2;
	y[i + 1] = 
                  y[i] + (vy[i] + vy[i + 1])*Δt/2;
	th[i + 1] = 
                  th[i] + (thDot[i] + thDot[i + 1])*Δt/2;
	
(* Estimate new value of n. *)

n = Max[0, n + (ay[i] - ay[i - 1])] ];     (* End For loop *)

  

nsteps = i;
(Subscript[t, Hop]) = (Subscript[t, Slip]) + nsteps Δt;
{(Subscript[X, Hop]), (Subscript[Y, Hop])} = {x[nsteps], y[nsteps]};
\({X\&.\_Hop, Y\&.\_Hop} = {vx[nsteps], vy[nsteps]}\);
\({θ\_Hop, θ\&.\_Hop} = {th[nsteps], thDot[nsteps]}\);
i =. ;

  

\(CMSoln[tt_ /; t\_Slip ≤ tt < t\_Hop] = {
		
              \(Interpolation[
                  Table[{j\ Δt, {x[j], vx[j]}}, {j, 0, nsteps}]]
                  \)[tt - t\_Slip], 
              \(Interpolation[
                  Table[{j\ Δt, {y[j], vy[j]}}, {j, 0, nsteps}]]
                  \)[tt - t\_Slip]}\);
θSoln[tt_ /; (Subscript[t, Slip]) ≤ tt < (Subscript[t, Hop])] = 
            Interpolation[
                Table[{j Δt, {th[j], thDot[j]}}, {j, 0, 
                    nsteps}]][tt - (Subscript[t, Slip])];nsteps
)

Phase 3:  The Hop Itself (Free Fall)

When the hoop loses contact with the supporting surface, n[t] = f[t] = 0, so equations () and () describe a parabolic trajectory, while () reduces to \(TraditionalForm\`\(θ\&..\)[t] = 0\).  Thus, the hoop simply rotates about its center of mass at a constant angular speed equal to its angular speed at the instant the contact forces went to zero while the center of mass itself falls freely under the influence of gravity.  This combined free rotation/free fall continues until the height of the center of the hoop decreases once again to a, at which time n[t] acquires the positive value required to keep the point of contact between hoop and ground from moving below ground level.

The impact of the hoop with the supporting surface is extremely complicated to model, so our simulation terminates at the instant (Subscript[t, Impact]) at which the hoop reestablishes contact with the supporting surface, completing the hop.  Unfortunately, FindRoot is not guaranteed to deliver the root which actually corresponds to (Subscript[t, Impact]); if the initial guess passed to FindRoot is wrong, the root returned may be a time prior to the initial loss of contact between the hoop and the supporting surface!  For this reason, we adopt instead the following procedure, which is guaranteed to work in all cases:  We first find a time sufficiently large that the function ctrHoop[t] [[2]] - a is negative, and then find the root by bisection.  (Alternately, one may simply plot the vertical height of the center of the hoop, ctrHoop[t] [[2]], as a function of time, and so determine by inspection values for the initial guesses passed to FindRoot which will yield the correct root for (Subscript[t, Impact]).)

\(freeFall := 
      \((
CMHop = {X\_Hop, Y\_Hop}; 
vCMHop = {X\&.\_Hop, Y\&.\_Hop}; 

        CMSoln[tt_ /; tt ≥ t\_Hop] := 
          CMHop + vCMHop\ \((tt - t\_Hop)\) + {0, \(-g\)}/2 
              \((tt - t\_Hop)\)^2; 

        θSoln[tt_ /; tt >= t\_Hop] := 
          θ\_Hop + θ\&.\_Hop\ \((tt - t\_Hop)\); 


        ctrHoop[tt_] := 
          CMSoln[tt] - 
            λ\ a\ {Sin[θSoln[tt]], \ Cos[θSoln[tt]]}; 

        		
rootInterval = 1.0; 
accuracy = 1.0*10\^\(-6\); 

        bottomHoop[tt_] := \(ctrHoop[tt]\)[\([2]\)] - a; 

        For[t1 = t\_Hop, Positive[bottomHoop[t2 = t1 + rootInterval]], 
          t1 += rootInterval]; 

        While[t2 - t1 > accuracy, tmidpt = \((t1 + t2)\)/2; 
          If[Positive[bottomHoop[tmidpt]], t1 = tmidpt, t2 = tmidpt]]; 

        t\_Impact = tmidpt; 
)\)\)

Combined Solution (Phase 1, 2, and 3)

The figures below depict schematically the results of typical numerical simulations. The values used for the parameters correspond to a system that one could potentially realize experimentally:  the hoop is not massless, and the attached object is not an idealized point mass.

a = 1;  g = 9.8; λ = .95; ϵ = 0.01; (Subscript[μ, s]) = 1; 
    (Subscript[μ, k]) =  0.8;  Δt =  0.001;

To aid one in visualizing the motion of the hoop, we generate a figure which
depicts schematically the results of the numerical simulations described above. Shown are the trajectory of the center of mass (black line) and the positions of the attached object (large dots) and the center of the hoop (small dots) at forty equally spaced times.  "Spokes" (light gray lines) connecting simultaneous positions of the object and the center of the hoop have been added as an aid in visualization.

Code for Simulation Figures  (DO NOT PRINT!)

Options[plotSoln] := {ShowHoops -> False, NumberOfStrobeFlashes -> 40}
          
plotSoln[opts___]:=(
{plotWithHoops, ntSteps}  =  
            {ShowHoops, NumberOfStrobeFlashes}  /.  {opts}  /.  
              Options[plotSoln];		

rMass[t_] := 
            ctrHoop[t] + a {Sin[θSoln[t]], Cos[θSoln[t]]};
		
dt = N[(Subscript[t, Impact])/ntSteps];
ptsCM  =  CMSoln /@ Range[0, (Subscript[t, Impact]),  dt];
ptsCtrHoop  =  ctrHoop /@ Range[0, (Subscript[t, Impact]),  dt];
ptsMass  =  rMass /@  Range[0, (Subscript[t, Impact]),  dt];
		
hoops  =  Circle[#, a]& /@ ptsCtrHoop;

		
If[plotWithHoops,
			
		  
			
			lines={{GrayLevel[0.8], RGBColor[1, .4, .4], 
                        Map[Line, Transpose[{ptsCtrHoop, ptsMass}]]},{AbsolutePointSize[4], RGBColor[1, 0, 0], 
                        Map[Point, ptsMass]},\({AbsolutePointSize[2], RGBColor[1, 0, 0], \ 

                          \ \ \ \ \ \ \ 
                          Map[If[#[\([2]\)]\  > \ a + tolerance, {
                                  RGBColor[1, 1, 0], PointSize[0.01], 
                                  Point[#]}, Point[#]]&, \ ptsCtrHoop]}\)};\(ParametricPlot[CMSoln[tt], {tt, 0, t\_Impact}, 
                  Compiled -> False, 
	
                  Prolog -> {
                      hoops, {AbsoluteThickness[2], RGBColor[1, 1, 0], 

                        			\(Circle[#, a]&\)\ /@\ 
                          Select[ptsCtrHoop, 
                            \ #[\([2]\)]\  > \ 1 + tolerance\ &]}}, 
	
                  Epilog -> {lines, AbsoluteThickness[1.5], 
			
                      RGBColor[1, 1, 0], 
                      Line[{{\(-a\), aa = a + Max[Last/@ptsCtrHoop] + .018}, 
                          \ 
					{ptsCtrHoop[\([\(-1\), 1]\)] + a, 
                            aa}}], 
		AbsoluteThickness[.5], 
                      Line\ [{{\(-a\), 0}, 
                          \ {\ ptsCtrHoop[\([\(-1\), 1]\)] + a, 0}}]}, 
                  AspectRatio -> Automatic, 
                  Background -> RGBColor[0.3, 0.4, 1], 
			
                  DefaultColor -> GrayLevel[0], 
	
                  PlotRange 
                    -> {{\(-1.2\), 
                        ptsCtrHoop[\([\(-1\), 1]\)] + a + .3}, {\(-.1\), 
                        \((a + Max[Last/@ptsCtrHoop])\) 1.1}}, \ 
                  Axes -> None, \ Frame -> True]\),

			  
					
		lines=Graphics[{{GrayLevel[0.8], 
                          Map[Line, Transpose[{ptsCtrHoop, ptsMass}]]},{AbsolutePointSize[4], Map[Point, ptsMass]},\({AbsolutePointSize[2], \ 
\ \ \ \ \ \ \ 
                            Map[Point, \ ptsCtrHoop]}\)}];
		Show[lines, 
                  ParametricPlot[CMSoln[tt], {tt, 0, (Subscript[t, Impact])}, 
                    Compiled -> False, Axes -> False, 
                    AspectRatio -> Automatic, PlotRange -> All, 
                    DisplayFunction -> Identity]]
			]     (*  End  If  *)
		)

Small Hop:  Attached Weight Initially at the Top of the Hoop

The character of the motion can change radically, depending on the initial conditions.  For instance, from an initial configuration in which the attached weight is at the top of the hoop ((Subscript[θ, 0]) = 0), the hoop executes a modest hop, the size of which depends on the initial speed of the hoop.  Here we take  (Subscript[ω, 0]) = 1.

(Subscript[θ, 0]) = 0;  (Subscript[ω, 0]) = 1;  
    c = (a/g) ((Subscript[ω, 0])^2);
rollingWithoutSlipping[(Subscript[θ, 0]), (Subscript[ω, 0])];
\(\(slipping[X\_Slip, Y\_Slip, θ\_Slip, X\& . \_Slip, 
        Y\& . \_Slip, θ\& . \_Slip, 
        nCyc[Cos[θSoln[t\_Slip - Δt]]] - g, 
        nCyc[Cos[θ\_Slip]]];\)\)
freeFall;
plotSoln[]

Figure 2.  Parametric plot showing the motion of the hopping hoop for  (Subscript[θ, 0]) = 0 and (Subscript[ω, 0]) = 1.

Large Hop:  Attached Weight Initially at the Bottom of the Hoop

If the attached weight is initially at the bottom of the hoop  ((Subscript[θ, 0]) = π), the subsequent motion of the hoop + weight system is radically different.  In the simulation which follows, the mechanical energy of the system is exactly the same as in the previous case.

(Subscript[ω, 0])=x/.Solve[(energy/.{θ[t] -> 0,θ'[t]->1})==(energy/.{θ[t] -> Pi,θ'[t]->x}),x][[2]];

(Subscript[θ, 0]) = π;   c = (a/g) ((Subscript[ω, 0])^2); 
rollingWithoutSlipping[(Subscript[θ, 0]), (Subscript[ω, 0])];
\(\(slipping[X\_Slip, Y\_Slip, θ\_Slip, X\& . \_Slip, 
        Y\& . \_Slip, θ\& . \_Slip, 
        nCyc[Cos[θSoln[t\_Slip - Δt]]] - g, 
        nCyc[Cos[θ\_Slip]]];\)\)
freeFall;
plotSoln[];  plotSoln[ShowHoops -> True];

Figure 3.  Motion of the hoop for  (Subscript[θ, 0]) = π.  The mechanical energy of the system is identical to that for the case illustrated in Figure 2.  The hoop is airborne for most of the motion depicted here; observe that the center of mass (black line) follows the parabolic trajectory characteristic of an object falling freely under the influence of gravity.

Conclusions

In the decade since its initial release, Mathematica has had a tremendous impact on physics teaching and research.  This article describes how Mathematica was used to investigate an old chestnut from Classical Mechanics, the well-known hopping hoop problem.  In exploring the motion of the hopping hoop, we were well served the the means available in Mathematica to perform numerical, algebraic, and graphical tasks in a single, unified environment.

Acknowledgement

The author wishes to thank Stan Wagon for his help and encouragement with this project, and for providing some of the code to improve the figures.

References

Bollabás, B., editor, Littlewood's  Miscellany, Cambridge University Press, London, UK, 1986, p. 37.

Pritchett, T., "The Hopping Hoop Revisited," The American Mathematical Monthly 106:7, 609-617 (1999).

Tokieda, T.,  "The Hopping Hoop," The American Mathematical Monthly 104:2, 152-154 (1997).

Zeller, T., "Frequently Asked Questions about Version 3.0," The Mathematica Journal 7:1, 35 (1997).

About the Author

Tim Pritchett is an Associate Professor of Physics at the United States Military Academy, West Point.  He makes extensive use of Mathematica in his research in nonlinear optics.

Tim Pritchett
Department of Physics
United States Military Academy
West Point, NY  10996
ht1187@usma.edu