Exercise 1.21: A dumbbell

This notebook is a solution to exercise 1.21 from Structure and Interpretation of Classical Mechanics, by Sussman and Wisdom. Visit the html version of the exercise for the surrounding context!

The uneven dumbbell. We've just made it through four exercises which embrace the idea that you can bake constraints into the coordinate transformation. But why should we believe that this is allowed?

First, load the sicmutils environment:

(require '[sicmutils.env :refer :all])

This exercise tries to do a coordinate change that is really careful about not changing the dimension of the configuration space, so that we can show that this move is allowed. Here's the setup:

The idea here is to take the distance between the particles inline_formula not implemented and treat it as a new dimension inline_formula not implemented.

Goal is to assume that Newtonian mechanics' approach to constraints, shown here, I think, is correct, and then show that the resulting equations of motion let us treat the distance coordinate inline_formula not implemented as a constant.

Multiple Particle API

Many exercises have been dealing with multiple particles so far. Let's introduce some functions that let us pluck the appropriate coordinates out of the local tuple.

If we have the velocity and mass of a particle, its kinetic energy is easy to define:

(defn KE-particle [m v]

  (* (/ 1 2) m (square v)))

This next function, extract-particle, takes a number of components – 2 for a particle with 2 components, 3 for a particle in space, etc – and returns a function of local and i, a particle index. This function can be used to extract a sub-local-tuple for that particle from a flattened list.

(defn extract-particle [pieces]

  (fn [[t q v] i]

    (let [indices (take pieces

                    (iterate

                      inc (* i pieces)))

          extract (fn [tuple]

                    (mapv (fn [i]

                            (ref tuple i))

                      indices))]

      (up t (extract q) (extract v)))))

Part A: Newton's Equations

Write Newton's equations for the balance of forces for the four rectangular coordinates of the two particles, given that the scalar tension in the rod is inline_formula not implemented.

TODO: Write these down from the notebook!

Part B: Dumbbell Lagrangian

Write the formal Lagrangian

formula not implemented

such that Lagrange's equations will yield the Newton's equations you derived in part a.

Here is how we model constraint forces. Each pair of particles has some constraint potential acting between them:

(defn U-constraint [q0 q1 F l]

  (* (/ F (* 2 l))

    (- (square (- q1 q0))

      (square l))))

And here's a Lagrangian for two free particles, subject to a constraint potential inline_formula not implemented acting between them.

(defn L-free-constrained [m0 m1 l]

  (fn [local]

    (let [extract (extract-particle 2)

          [_ q0 qdot0] (extract local 0)

          [_ q1 qdot1] (extract local 1)

          F (ref (coordinate local) 4)]

      (- (+ (KE-particle m0 qdot0)

           (KE-particle m1 qdot1))

        (U-constraint q0 q1 F l)))))

Finally, the path. This is rectangular coordinates for each piece, plus inline_formula not implemented between them.

(def q-rect

  (up

    (literal-function 'x_1)

    (literal-function 'y_0)

    (literal-function 'x_1)

    (literal-function 'y_1)

    (literal-function 'F)))

This shows the Lagrangian itself, which answers part b:

(let [L (L-free-constrained 'm_0 'm_1 'l)

      f (compose L (Gamma q-rect))]

  (f 't))

Here are the Lagrange equations, which, if you squint, are like Newton's equations from part a.

(let [L (L-free-constrained 'm_0 'm_1 'l)

      f ((Lagrange-equations L) q-rect)]

  (f 't))

Part C: Coordinate Change

Make a change of coordinates to a coordinate system with center of mass coordinates inline_formula not implemented, inline_formula not implemented, angle inline_formula not implemented, distance between the particles inline_formula not implemented, and tension force inline_formula not implemented. Write the Lagrangian in these coordinates, and write the Lagrange equations.

This is a coordinate change that is very careful not to reduce the degrees of freedom.

First, the coordinate change:

(defn cm-theta->rect [m0 m1]

  (fn [[_ [x_cm y_cm theta c F]]]

    (let [total-mass (+ m0 m1)

          m0-distance (* c (/ m1 total-mass))

          m1-distance (* c (/ m0 total-mass))]

      (up (- x_cm (* m0-distance (cos theta)))

          (- y_cm (* m0-distance (sin theta)))

          (+ x_cm (* m1-distance (cos theta)))

          (+ y_cm (* m1-distance (sin theta)))

          F))))

Then the coordinate change applied to the local tuple:

(let [local (up 't

              (up 'x_cm 'y_cm 'theta 'c 'F)

              (up 'xdot_cm 'ydot_cm 'thetadot 'cdot 'Fdot))

      C (F->C (cm-theta->rect 'm_0 'm_1))]

  (C local))

Then the Lagrangian in the new coordinates;

(defn L-free-constrained* [m0 m1 l]

  (compose (L-free-constrained m0 m1 l)

    (F->C (cm-theta->rect m0 m1))))

This shows the Lagrangian itself, after the coordinate transformation:

(let [q (up (literal-function 'x_cm)

          (literal-function 'y_cm)

          (literal-function 'theta)

          (literal-function 'c)

          (literal-function 'F))

      L (L-free-constrained* 'm_0 'm_1 'l)

      f (compose L (Gamma q))]

  (f 't))

Here are the Lagrange equations:

(let [q (up (literal-function 'x_cm)

          (literal-function 'y_cm)

          (literal-function 'theta)

          (literal-function 'c)

          (literal-function 'F))

      L (L-free-constrained* 'm_0 'm_1 'l)

      f ((Lagrange-equations L) q)]

  (f 't))

That final equation states that inline_formula not implemented. Amazing!

Part D: Substitute inline_formula not implemented

You may deduce from one of these equations that inline_formula not implemented. From this fact we get that inline_formula not implemented and inline_formula not implemented. Substitute these into the Lagrange equations you just computed to get the equation of motion for inline_formula not implemented, inline_formula not implemented, inline_formula not implemented.

We can substitute the constant value of inline_formula not implemented using a function that always returns inline_formula not implemented to get simplified equations:

(let [q (up (literal-function 'x_cm)

          (literal-function 'y_cm)

          (literal-function 'theta)

          (fn [t] 'l)

          (literal-function 'F))

      L (L-free-constrained* 'm_0 'm_1 'l)

      f ((Lagrange-equations L) q)]

  (f 't))

This is saying that the acceleration on the center of mass is 0.

The fourth equation, the equation of motion for the inline_formula not implemented, is interesting here. We need to pull in the definition of "reduced mass" from exercise 1.11:

(defn reduced-mass [m1 m2]

  (/ (* m1 m2)

     (+ m1 m2)))

If we let inline_formula not implemented be the reduced mass, this equation states:

We can verify this with Scheme by subtracting the two equations:

(let [F (literal-function 'F)

      theta (literal-function 'theta)

      q (up (literal-function 'x_cm)

          (literal-function 'y_cm)

          theta

          (fn [_] 'l)

          F)

      L (L-free-constrained* 'm_0 'm_1 'l)

      f ((Lagrange-equations L) q)

      m (reduced-mass 'm_0 'm_1)]

  (- (ref (f 't) 3)

    (- (F 't)

      (* m 'l (square ((D theta) 't))))))

Part E: New Lagrangian

Make a Lagrangian (inline_formula not implemented) for the system described with the irredundant generalized coordinates inline_formula not implemented, inline_formula not implemented, inline_formula not implemented and compute the Lagrange equations from this Lagrangian. They should be the same equations as you derived for the same coordinates in part d.

For part e, I wrote this in the notebook - it is effectively identical to the substitution that is happening on the computer, so I'm going to ignore this. You just get more cancellations.

But let's go at it, for fun.

Here's the Lagrangian of 2 free particles:

(defn L-free2 [m0 m1]

  (fn [local]

    (let [extract (extract-particle 2)

          [_ q0 qdot0] (extract local 0)

          [_ q1 qdot1] (extract local 1)]

      (+ (KE-particle m0 qdot0)

        (KE-particle m1 qdot1)))))

Then a version of cm-theta->rect where we ignore inline_formula not implemented, and sub in a constant inline_formula not implemented:

(defn cm-theta->rect* [m0 m1 l]

  (fn [[_ [x_cm y_cm theta]]]

    (let [total-mass (+ m0 m1)

          m0-distance (* l (/ m1 total-mass))

          m1-distance (* l (/ m0 total-mass))]

      (up (- x_cm (* m0-distance (cos theta)))

        (- y_cm (* m0-distance (sin theta)))

        (+ x_cm (* m1-distance (cos theta)))

        (+ y_cm (* m1-distance (sin theta)))))))

The Lagrangian:

(defn L-free-constrained2 [m0 m1 l]

  (compose (L-free2 m0 m1)

           (F->C (cm-theta->rect* m0 m1 l))))

Equations:

(let [q (up (literal-function 'x_cm)

          (literal-function 'y_cm)

          (literal-function 'theta)

          (fn [_] 'l)

          (literal-function 'F))

      L1 (L-free-constrained* 'm_0 'm_1 'l)

      L2 (L-free-constrained2 'm_0 'm_1 'l)]

  ((- ((Lagrange-equations L1) q)

     ((Lagrange-equations L2) q))

   't))

The only remaining equation is \eqref{eq:constraint-force} from above. This remains because the simplified Lagrangian ignores the inline_formula not implemented term.