Mathematica Activity 2:
The Swinging Atwood's Machine (SAM)

by Phillip Kent, METRIC Project, 
for M3A16/M4A16, 1999-2000.
Version 1.2, 21 November 1999.

The Mathematica notebook containing this document (Activity2.nb), plus related materials can be accessed at

http://metric.ma.ic.ac.uk/dynamics/

Note that you can save yourself a LOT of typing by downloading this notebook!

This activity forms part of your assessed coursework. As in Activity 1, the different pieces of "portfolio work" that we would like you to produce are clearly indicated by the text cells with grey background. Please show all the Mathematica code that you use and explain fully the steps in your calculations.

This coursework should be submitted by FRIDAY 14th JANUARY.

0. Loading the dynamics toolkit

On the Mathematics Department Windows NT network, you can load the toolkit into Mathematica with the following commands:

AppendTo[$Path,  "//OPHELIA/OPHELIA_D"]

<< dynapack.m

To install the toolkit file locally in NT, or on your own computer, see the instructions on the web page.

1. The Swinging Atwood's Machine

Atwood's Machine is a classic device constructed by George Atwood in the eighteenth century in order to measure the gravitational acceleration, g. Two masses are attached to either end of a (massless) string which is passed over a (frictionless) pulley. When the two masses are slightly unequal, the heavier mass falls in a slowed-down free fall motion, which allowed Atwood to measure g using a clockwork timing device. 

The Swinging Atwood's Machine (SAM) was introduced (and actually built as a lab experiment) by Tufillaro et al [American Journal of Physics 52, 10, 895-903, 1984]:

Here, one of the masses, m, is constrained to swing on a (frictionless) vertical plane, and suppose that the string passes through a hole in the planar surface, with frictionless contact, and that the string is long enough so that the counter-mass, M, never strikes the pulley.

SAM's Hamiltonian, which you should have already derived in Problem Sheet 2, Question 7, is

H =((((Subscript[p, r])^2))/(2  m(1 + μ)))+((((Subscript[p, θ])^2))/(2  m (r^2)))+m g r (μ - cosθ),       μ =M/m.

As already noted in Activity 1, subscripted forms such as (Subscript[p, r]) and (Subscript[p, θ]) are best avoided so we'll use "pr" and "pθ" instead in the Mathematica definition. And for simplicity we'll set m = 1 and g = 10.

\(samHam = 
      pr[t]\^2\/\(2 \((1 + μ)\)\) + pθ[t]\^2\/\(2\ r[t]\^2\) + 
        10\ r[t] \((μ - Cos[θ[t]])\)\)

Then the equations of motion are calculated via Hamilton's equations:

samEqns[μ_] = {r'[t] == D[samHam, pr[t]], θ'[t] == 
          D[samHam, pθ[t]], pr'[
            t] == -D[samHam, r[t]], pθ'[
            t] == -D[samHam, θ[t]]}

It is useful to define the equations as a function of μ in this way since much use will be made of that dependence. For example, here is a Trajectory for μ = 1.1185:

MakeTrajectory[
      samEqns[1.1185], {r[0] == 1.0, θ[0] == π/2, 
        pr[0] == 0.0, pθ[0] == 0.0}, {r, θ, pr, 
        pθ},  {t,  10}]

SAM's simplicity of construction belies its very complex dynamical behaviour when M > m, since a dynamical equilibrium can be achieved by the centrifugal force on the swinging mass m counterbalancing the free-fall pull of M. There are a lot of interesting behaviours that occur, and in this activity we'll look at only a few of these.

2. Periodic solutions

Portfolio Work 2.1

In Tufillaro et al's numerical study of SAM, they found a class of periodic solutions, all starting from rest at r = 1, θ = π/2, for the following (approximate) values of μ:
                     μ = 1.1185, 1.1727, 1.555, 1.665, 2.394, 2.812, 4.745.
Using MakeTrajectory and TrajPlot, produce pictures of these solutions plotted in (x, y) space*. A few of these are given below for reference. For each solution, determine the approximate period.

*You will need to use the Transform option to TrajPlot, and define a coordinate transformation function -- cf. the Kepler system plots in Activity 1.

Movies

For fun (and not to be assessed), you might like to generate some animations based on these periodic solutions, using the following code. Note that you'll need to generate a lot of frames to get a smooth-looking movie (so it's advisable to delete any movie outputs before saving your notebook).

First, a Trajectory to look at, and an (x, y) plot of it. The appropriate coordinate transformation function is left for you to insert!

traj4movie = 
      MakeTrajectory[
        samEqns[1.1185], {r[0] == 1.0, θ[0] == π/2, 
          pr[0] == 0.0, pθ[0] == 0.0}, {r, θ, pr, 
          pθ},  {t,  10}]

plot4movie = 
      TrajPlot[traj4movie, {r,  θ},  
        Transform -> <INSERT YOUR FUNCTION HERE>, 
        PlotStyle -> GrayLevel[0.7]]

This function extracts from a Trajectory the position of the mass m in the (x, y) plane at time t.

orbitPoint[traj_Trajectory, 
        t_] := TrajComponent[traj, r][
          t]*{Sin[TrajComponent[traj, θ][
              t]], -Cos[TrajComponent[traj, θ][t]]}

And this command builds the movie.

Table[
      Show[plot4movie, 
Graphics[{Line[{{0, 0}, 
                orbitPoint[traj4movie, t]}], PointSize[0.02], 
            Point[orbitPoint[traj4movie, t]]}], 
        AspectRatio -> Automatic, Axes -> True, 
        PlotRange -> {{-3, 3}, {-6.5, 0}}], {t, 0, 10, 
        0.2}]

For different Trajectories, you'll need to alter the "PlotRange", the total time (10 above) and time-step for the movie (0.2 above).

3. Zero velocity curves

In the Problem Sheet question you derived the result that the zero velocity curves (ZVCs) for the SAM system are given by

E =m g r (μ  -  
              cosθ).

With our simplifying assumptions of m = 1 and g = 10, this turns out as

r=(E/(10 (μ  -  cosθ))).

We will consider here just the case of μ > 1, for which all motions are bounded and the ZVCs are ellipses.

Portfolio Work 2.2

Use ParametricPlot to generate ZVCs for three of the periodic orbits you generated in PW2.1. Make combined plots (with the Show command) to illustrate that the orbits are indeed bounded by the ZVCs.

The ZVC should be a bounding curve for any solution of SAM of course, periodic or otherwise. Illustrate that this is the case for three different arbitrary choices of initial conditions and μ (> 1). It will probably be worth running MakeTrajectory for fairly large times (up to 50 or so), for which you will need to set the MaxSteps option to a value of 2000-3000 or higher.

4. The case of μ = 3: The invariant J

Tufillaro [American Journal of Physics, 54 (1986), 142-143] discovered (initially by numerical observation) that a second integral of the motion (besides the Hamiltonian) exists for the case of μ = 3. As you have already seen in the Problem Sheets (2, Q7, and 5, Q7) this invariant quantity is

J = 
      10 (r^2) Cos[(θ/2)]^2 Sin[(θ/2)] + 
        (1/4) pθ (pr Cos[(θ/2)] - ((2 pθ Sin[(θ/2)])/r))

A numerical check for invariance

Portfolio Work 2.3

Choose three initial points somewhere in the (r, θ) plane, make a Trajectory for μ = 3.0 starting at each point, and use CheckQuantity and ListPlot to show graphs of the variation -- or hopefully lack of! -- J along each Trajectory.

For one of your initial points, make a combined plot showing J as a function of time for three Trajectories where μ is equal to, just greater than, and just smaller than, 3.0. Distinguish the lines using different PlotStyles.

A symbolic check for invariance

Mathematica can be used to check the invariance of J(t) symbolically. First, set up an expression for J that includes the time dependencies:

Jt = 
      J /. {r -> r[t], θ -> θ[t], pr -> pr[t], 
          pθ -> pθ[t]}

The time derivative of this is:

D[Jt, t]

We want to substitute into this expression the equations of motion. However, Mathematica needs rules ("... -> ...") to do this rather than equations ("... == ..."). A little Mathematica substitution can get us there in one go:

subs = samEqns[3] /. Equal -> Rule

Then,

D[Jt, t] /. subs

Is this really all equal to zero? Generally, to get proper simplification with Mathematica you must use one of the simplification commands, and the basic one is Simplify:

Simplify[D[Jt, t] /. subs]

Terminating singular trajectories: Tears, hearts and rabbit's ears

Another interesting class of orbits discovered by Tufillaro et al. is those which start by firing the mass "out of the hole", or at least (for the purposes of numerical simulation) having r(0) very small. They called these "terminating singular trajectories", because they start from the hole (where the equations of motion are singular) and at some later time they terminate back at the hole. (Actually, the approximate numerical solutions aren't able to terminate, instead they get very close to r = 0 and then "fire" back out again.)

The singular trajectory that occurs at μ = 3 is rather special. It produces a symmetrical terminating trajectory no matter what the values of θ(0) and pr(0) are, whereas those at other μ values only exist for certain initial values.

Portfolio Work 2.4

Take μ = 3, r(0) = 0.001, pθ(0) = 0.0 and pr(0) = 4.0. Make Trajectories for values of θ(0) ranging from about 0.5 up to about π, and produce xy-plots illustrating how the Trajectory shape changes from a "tear" to a "heart" to "rabbit's ears".

Illustrate, for a few choices of θ(0) and pr(0), that this terminating trajectory persists no matter what the value of pr(0). (However, note that numerical errors accumulate rather badly when pr(0) gets small, and the numerical solutions become unreliable).