Structure and Interpretation of Classical Mechanics in Clojure

1.
Introduction

This Article is based on the sicmutils project by C. Smith.

The indicated chapters and subsections correspond to the book "Structure and Interpretation of Classical Mechanics" 2nd ed. by G.J. Sussman and J. Wisdom (freely available at MIT and github).

First, we provide the necessary Libraries (remove the # in the next cell if you open this notebook for the first time): 

# #remove the hashes at the beginning of the lines
#echo '{:deps {net.littleredcomputer/sicmutils {:mvn/version "0.12.0-SNAPSHOT"}}}' > deps.edn
#clj -Sforce -e ":ok"

0.1s

Language:Bash

(require 'sicmutils.env)
(in-ns 'sicmutils.env)
(require '[clojure.spec.alpha :as s])
(require '[clojure.spec.test.alpha :as ts])
(require '[sicmutils.numerical.minimize :as m])
(require '[sicmutils.mechanics.lagrange :refer :all])

0.0s

Language:Clojure

2.
Computing Actions (Chapter 1.4)

Calculate the action for the free particle along a path. Consider a particle moving at uniform speed along a straight line.

(defn test-path
  "See p. 20"
  [t]
  (up (+ (* 4 t) 7)
      (+ (* 3 t) 5)
      (+ (* 2 t) 1)))

(Lagrangian-action (L-free-particle 3.0) test-path 0.0 10.0)

0.0s

Language:Clojure

435

Paths of minimum Action

Show that the action is smaller along a straight-line test path than along nearby paths

(defn make-η
  "See p. 21"
  [ν t1 t2]
  (fn [t]
    (* (- t t1) (- t t2) (ν t))))

(defn varied-free-particle-action
  "See p. 21"
  [mass q ν t1 t2]
  (fn [ε]
    (let [η (make-η ν t1 t2)]
      (Lagrangian-action (L-free-particle mass)
                         (+ q (* ε η)) t1 t2))))
						 
((varied-free-particle-action 3.0 test-path (up sin cos square) 0.0 10.0) 0.01)

0.1s

Language:Clojure

564.1214285996772

Simulate lots of the paths in this manner. Proof that the minimum value of the action is the action along the straight path. For this proof it suffices that some optimization parameter is close to zero (need not be exactly zero).

(minimize (varied-free-particle-action 3.0 test-path (up sin cos square) 0.0 10.0) -2 1)

0.5s

Language:Clojure

3[-9.163058228844889e-10,435.0000000594215,6]

We provide some helper functions for visualization:

0.0s

Language:Clojure

(defn alt-range [t0 t1 nofsteps]
    (map float (flatten [t0 (linear-interpolants t0 t1 (- nofsteps 2)) t1])))

(defn make-path-txyz [fn_t t0 t1 nofsteps]
  (mapv #(cons % (.v (fn_t %)))
    (alt-range t0 t1 nofsteps)))

(defn points->plotly2 [paths [idx1 label1 idx2 label2]]
  (letfn [(f1x [u]
            (map (fn [y] (map (fn [x] (nth x y)) u)) (range 4)))
          (f4x [name w]
            [{:x (nth w idx1)
              :y (nth w idx2)
              :type "scatter"
              :name name}])
          (f5x [w]
            (map #(f4x (first %) (f1x (second %))) w))
          (f6x [data]
           {:data data
            :layout
            {:autosize false
             :width 600
             :height 500
             :xaxis {:title label1}
             :yaxis {:title label2}}})]
   (f6x (flatten (f5x paths)))))

(defn points-xy->plotly [paths] 
  (points->plotly2 paths [1 "x" 2 "y"]))

(defn points-xz->plotly [paths] 
  (points->plotly2 paths [1 "x" 3 "z"]))

(defn points-tz->plotly [paths] 
  (points->plotly2 paths [0 "t" 3 "z"]))

(comment "a test example"
 (def werte {"line1" [[0 7 5 1] [1 11 8 10]]
            "line2" [[2 9 2 4] [3 3  9 7]]})
	(points-xy->plotly werte))

0.0s

Language:Clojure

Create data to plot two straight paths in the xy plane. One path is along the x axis (name: path-along-x), the second path leads in all directions.

(defn path-along-x
  [t]
  (up (+ (* 5 t) 1)
      (* 0 t)
      (* 0 t)))

(def path-points-1 {"path-along-x" (make-path-txyz path-along-x 0 10 3)
                 "all directions" (make-path-txyz test-path 0 10 3)})

(points-xy->plotly path-points-1)

0.0s

Language:Clojure

2{:data(2):layout{5}}

In a ClojureScript cell, plot the two paths.

(.plot js/Nextjournal plotly-1↩)

0.0s

Language:ClojureScript

Create two variations of the path-along-x. Calculate the action. Show once again that the Lagrangian-action is indeed smallest for the straight path.

(defn make-varied-path [ε t0 t1]
 (+ path-along-x (* ε (make-η (up #(* 0 %) identity #(* 5.0 (sin %))) t0 t1))))

(def small-varied-path (make-varied-path 0.01 0 10))
(def large-varied-path (make-varied-path 0.02 0 10))


[(Lagrangian-action (L-free-particle 3.0) path-along-x 0.0 10.0)
 (Lagrangian-action (L-free-particle 3.0) small-varied-path 0.0 10.0)
 (Lagrangian-action (L-free-particle 3.0) large-varied-path 0.0 10.0)]

0.2s

Language:Clojure

3[375,383.8260332594603,410.3041330378416]

Create data to plot the three paths in the xz plane along with their actions.

(def path-points-2
    {"path-along-x" (make-path-txyz path-along-x 0 10 8)
     "small-varied-path" (make-path-txyz small-varied-path 0 10 24)
     "large-varied-path" (make-path-txyz large-varied-path 0 10 32)})
(points-xz->plotly path-points-2)

0.0s

Language:Clojure

2{:data(3):layout{5}}

Plot the three paths.

(.plot js/Nextjournal plotly-2↩)

0.0s

Language:ClojureScript

Finding trajectories that minimize the action

The SICM library provides a procedure that constructs a one dimensional path (along, say, the z axis) using an interpolation polynomial: (make-path t0 q0 t1 q1 qs), where q0 and q1 are the endpoints, t0 and t1 are the corresponding times, and qs is a list of intermediate points. We give an example (note that the result, initial-path, is itself a function):

(def pi-half (* 0.5 Math/PI))
(def initial-qs [0.1 0.2 0.2])
(def initial-path (make-path 0 1.0  pi-half 0.0 initial-qs))

0.0s

Language:Clojure

sicmutils.env/initial-path

Construct a parametric action that is just the action computed along that parametric path. Find approximate solution paths of a free particle and the harmonic oszillator respectively (hint: use the SICM procedure multidimensional-minimize).

(defn parametric-path-actn
  [Lagrangian t0 q0 t1 q1]
  (fn [qs]
    (let [path (make-path t0 q0 t1 q1 qs)]
      (Lagrangian-action Lagrangian path t0 t1))))

(defn fnd-path
  [Lagrangian t0 q0 t1 q1 initial-qs]
  (let [minimizing-qs
        (m/multidimensional-minimize
          (parametric-path-actn Lagrangian t0 q0 t1 q1)
          initial-qs)]
    (make-path t0 q0 t1 q1 minimizing-qs)))

(def free-path 
  (fnd-path (L-free-particle 3.0) 0.0 1.0 pi-half 0.0 initial-qs))

(def harmonic-path 
  (fnd-path (L-harmonic 1.0 1.0) 0.0 1.0 pi-half 0.0 initial-qs))

5.6s

Language:Clojure

sicmutils.env/harmonic-path

Make a plot of these one dimensional paths, this time not in configuration space (i.e. the x-z plane) but in phase space (i.e. the t-z plane). This shows that, upon optimization, the initial-path turns into a streight line and a sinusoidal curve respectively.

(defn make-path-tz [fn_t t0 t1 nofsteps]
  (map #(vector % 0 0 (fn_t %)) (alt-range t0 t1 nofsteps)))

(let [i (make-path-tz initial-path 0 pi-half 25)
      f (make-path-tz free-path 0 pi-half 25)
      h (make-path-tz harmonic-path 0 pi-half 25)]
  (points-tz->plotly {"initial-path" i "free-path" f "harmonic-path" h}))

0.0s

Language:Clojure

2{:data(3):layout{5}}

(.plot js/Nextjournal plotly-3↩)

0.0s

Language:ClojureScript

Show that your numerically attained harmonic-path approximates the well known analytic solution x(t) = cos(t). 

(points-tz->plotly 
 {"diff-cos-num" 
  (make-path-tz #(- (cos %) (harmonic-path %)) 0 pi-half 50)})

0.0s

Language:Clojure

2{:data(1):layout{5}}

(.plot js/Nextjournal plotly-4↩)

0.0s

Language:ClojureScript

Calculate the Lagrange equation of the harmonic oszillator.

;;we provide a function simplify-spec (defined at the bottom)
;;which should be called if the specs are activated 
(defmacro simplify-alt [x]
    (if (resolve 'simplify-spec)
        `(simplify-spec ~x)
        `(simplify ~x)))

0.0s

Language:Clojure

sicmutils.env/simplify-alt

(-> (((Lagrange-equations (L-harmonic 'm 'k)) (literal-function 'q)) 't)
  simplify-alt
  ->TeX	
  print)

0.6s

Language:Clojure

k\,q\left(t\right) + m\,{D}^{2}q\left(t\right)

3.
The Euler-Lagrange Equations (Chapter 1.5)

3.1.
Computing Lagrange's Equations (Section 1.5.2)

The free particle

State the dynamic equation of motion (i.e. the Lagrange equation a.k.a Newton's second law) of the free particle.

(-> (((Lagrange-equations (L-free-particle 'm)) (literal-function 'q)) 't)
  simplify-alt
  ->TeX
  print)

0.6s

Language:Clojure

m\,{D}^{2}q\left(t\right)

Check that an arbitrary straight-line path satisfies this equation, i.e. that inserting a straight line for q(t)q(t)q(t)gives identically zero (strictly speaking the zero covector of three dimensions).

(letfn [(straight-line [t]
  (up (+ (* 'a t) 'a0)
    (+ (* 'b t) 'b0)
    (+ (* 'c t) 'c0)))]
  (print (simplify-alt (((Lagrange-equations (L-free-particle 'm)) straight-line) 't))))

0.7s

Language:Clojure

The harmonic oscillator

State the dynamic equation of motion for the harmonic oszillator with arbitrary mass  and spring constant .

(-> (((Lagrange-equations (L-harmonic 'm 'k)) (literal-function 'q)) 't)
  simplify-alt
  ->TeX
  print)

0.6s

Language:Clojure

k\,q\left(t\right) + m\,{D}^{2}q\left(t\right)

Plug in a sinusoid with arbitrary amplitude  AAA, frequency ω\omegaωand phase  ϕ\phiϕand show that the only solutions allowed are ones where  ω=k/m\omega = \sqrt{k/m} ω=√​k/m​​​ 

(letfn [(proposed-solution [t]
  (* 'A (cos (+ (* 'omega t) 'φ))))]
  (-> (((Lagrange-equations (L-harmonic 'm 'k)) proposed-solution) 't)
    simplify-alt
    ->TeX
    print))

0.8s

Language:Clojure

- A\,m\,{\omega}^{2}\,\cos\left(\omega\,t + \varphi\right) + A\,k\,\cos\left(\omega\,t + \varphi\right)

Exercise 1.11: Kepler's third law

Show that a planet in circular orbit satisfies Kepler's third law n2a3=G(M1+m2)n^2a^3=G(M_1+m_2)n​2​​a​3​​=G(M​1​​+m​2​​), where nnnis the angular frequency of the orbit and aaais the distance between sun and planet. (Hint: use the reduced mass  to construct the Lagrangian)

(defn gravitational-energy [G M_1 m_2]
  (fn [r]
   (- (/ (* G M_1 m_2) r))))

(defn circle [t]
  (up 'a (* 'n t)))

(let [lagrangian (L-central-polar     
                  (/ (* 'M_1 'm_2) (+ 'M_1 'm_2))     
                  (gravitational-energy 'G 'M_1 'm_2))]
      (-> (((Lagrange-equations lagrangian) circle) 't)
        simplify-alt
        ->TeX
        print))

0.6s

Language:Clojure

\begin{bmatrix}\frac{- M_1\,{a}^{3}\,m_2\,{n}^{2} + G\,{M_1}^{2}\,m_2 + G\,M_1\,{m_2}^{2}}{M_1\,{a}^{2} + {a}^{2}\,m_2}&0\end{bmatrix}

4.
How to find Lagrangians (Chapter 1.6)

Central force field

State the dynamic equation of motion for the uniform acceleration and the central potential, the latter in rectangular as well as in polar coordinates.

(defn make-tex-matrix [lst]
  (str "\\begin{matrix} "
    (->> lst
     (map #(->TeX (ts/with-instrument-disabled (simplify %))))
     (map #(str % " \\\\ "))
     (apply str))
    "\\end{matrix}"))

(def three-equations
 (let [f1 (((Lagrange-equations
              (L-uniform-acceleration 'm 'g))
            (up (literal-function 'x)
              (literal-function 'y)))
           't)
       f2 (((Lagrange-equations
              (L-central-rectangular 'm (literal-function 'U)))
            (up (literal-function 'x)
              (literal-function 'y)))
           't)
       f3 (((Lagrange-equations
              (L-central-polar 'm (literal-function 'U)))
            (up (literal-function 'r)
              (literal-function 'phi)))
           't)]
   (make-tex-matrix [f1 f2 f3])))

(print three-equations)

1.0s

Language:Clojure

\begin{matrix} \begin{bmatrix}m\,{D}^{2}x\left(t\right)&g\,m + m\,{D}^{2}y\left(t\right)\end{bmatrix} \\ \begin{bmatrix}\frac{m\,{D}^{2}x\left(t\right)\,\sqrt {{\left(x\left(t\right)\right)}^{2} + {\left(y\left(t\right)\right)}^{2}} + x\left(t\right)\,DU\left(\sqrt {{\left(x\left(t\right)\right)}^{2} + {\left(y\left(t\right)\right)}^{2}}\right)}{\sqrt {{\left(x\left(t\right)\right)}^{2} + {\left(y\left(t\right)\right)}^{2}}}&\frac{m\,{D}^{2}y\left(t\right)\,\sqrt {{\left(x\left(t\right)\right)}^{2} + {\left(y\left(t\right)\right)}^{2}} + y\left(t\right)\,DU\left(\sqrt {{\left(x\left(t\right)\right)}^{2} + {\left(y\left(t\right)\right)}^{2}}\right)}{\sqrt {{\left(x\left(t\right)\right)}^{2} + {\left(y\left(t\right)\right)}^{2}}}\end{bmatrix} \\ \begin{bmatrix}- m\,{\left(D\phi\left(t\right)\right)}^{2}\,r\left(t\right) + m\,{D}^{2}r\left(t\right) + DU\left(r\left(t\right)\right)&2\,m\,D\phi\left(t\right)\,r\left(t\right)\,Dr\left(t\right) + m\,{\left(r\left(t\right)\right)}^{2}\,{D}^{2}\phi\left(t\right)\end{bmatrix} \\ \end{matrix}

Coordinate transformations (Section 1.6.1)

Calculate the  [x˙ y˙][\dot x \space \dot y][​x​˙​​ ​y​˙​​]velocity vector in polar coordinates.

(def vel (-> (velocity ((F->C p->r)
                   (->local 't (up 'r 'φ) (up 'rdot 'φdot))))
    simplify-alt ->TeX))
(print vel)

0.8s

Language:Clojure

\begin{pmatrix}- r\,\dot {\varphi}\,\sin\left(\varphi\right) + \dot r\,\cos\left(\varphi\right)\\r\,\dot {\varphi}\,\cos\left(\varphi\right) + \dot r\,\sin\left(\varphi\right)\end{pmatrix}

Calculate the lagrangian in polar coordinates twice. Once directly and once via the  lagrangian in rectangular coordinates.

(defn L-alternate-central-polar
  [m U]
  (compose (L-central-rectangular m U) (F->C p->r)))

(def two-lagrangians
 (let [lcp ((L-central-polar 'm (literal-function 'U))
            (->local 't (up 'r 'φ) (up 'rdot 'φdot)))
       lacp ((L-alternate-central-polar 'm (literal-function 'U))
             (->local 't (up 'r 'φ) (up 'rdot 'φdot)))]
   (make-tex-matrix [lcp lacp])))

(print two-lagrangians)

0.6s

Language:Clojure

\begin{matrix} \frac{1}{2}\,m\,{r}^{2}\,{\dot {\varphi}}^{2} + \frac{1}{2}\,m\,{\dot r}^{2} - U\left(r\right) \\ \frac{1}{2}\,m\,{r}^{2}\,{\dot {\varphi}}^{2} + \frac{1}{2}\,m\,{\dot r}^{2} - U\left(r\right) \\ \end{matrix}

Coriolis and centrifugal forces

State, in cartesian coordinates, the Lagrangian for the two dimensional free particle in a rotating coordinate system. 

(def L-free-rectangular L-free-particle)

(defn L-free-polar [m]
 (compose (L-free-rectangular m) (F->C p->r)))

(defn F [Omega]
 (fn [[t [r theta]]]
  (up r (+ theta (* Omega t)))))

(defn L-rotating-polar [m Omega]
 (compose (L-free-polar m) (F->C (F Omega))))

(defn r->p
  "rectangular to polar coordinates of state."
  [[_ [x y :as q]]]
  (up (sqrt (square q)) (atan (/ y x))))

(defn L-rotating-rectangular [m Omega]
  (compose (L-rotating-polar m Omega) (F->C r->p)))

(def L-simplify
  (simplify-alt ((L-rotating-rectangular 'm 'Omega)
       	(up 't (up 'x 'y) (up 'xdot 'ydot)))))

(def lrr (->TeX L-simplify))
(print  lrr)

1.1s

Language:Clojure

\frac{1}{2}\,{\Omega}^{2}\,m\,{x_r}^{2} + \frac{1}{2}\,{\Omega}^{2}\,m\,{y_r}^{2} + \Omega\,m\,x_r\,{\dot y}_r - \Omega\,m\,{\dot x}_r\,y_r + \frac{1}{2}\,m\,{{\dot x}_r}^{2} + \frac{1}{2}\,m\,{{\dot y}_r}^{2}

Derive the equations of motion, in which the centrifugal and the coriolis force appear. (Hint: take the expression just obtained and use this data as code via Clojure's macro feature).

(defmacro L-fn [args1 args2]
  `(fn ~args1 (fn ~args2 ~L-simplify)))

(def L-rotating-rectangular-simp (L-fn [m Omega] [[t [x y] [xdot ydot]]]))

(def leq 
  (->TeX
      (simplify-alt
        (((Lagrange-equations (L-rotating-rectangular-simp 'm 'Omega))
          (up (literal-function 'x) (literal-function 'y)))
         't))))

(print leq)

0.7s

Language:Clojure

\begin{bmatrix}- {\Omega}^{2}\,m\,x_r\left(t\right) -2\,\Omega\,m\,Dy_r\left(t\right) + m\,{D}^{2}x_r\left(t\right)&- {\Omega}^{2}\,m\,y_r\left(t\right) + 2\,\Omega\,m\,Dx_r\left(t\right) + m\,{D}^{2}y_r\left(t\right)\end{bmatrix}

Plot a clockwise rotating path. (Hints: a) Use the SICM function "Gamma" to create the triple  (t(xy)(x˙y˙))(t \: (x \: y) \: (\dot{x} \: \dot{y}))(t(xy)(​x​˙​​​y​˙​​))which can be transformed, b) the angular frequency must be negative) 

(defn test-path-2d
  [t]
  (up
   (+ (* 3 t) 7)
   (+ (* 5 t) 11)))

(defn rotate-path [path-fn]
(simplify-alt
 ((F->C p->r)
  ((F->C (F 'Omega))
   ((F->C r->p)
    ((Gamma path-fn) 't))))))

(def xy (nth (rotate-path test-path-2d) 2))
(def x (apply list (nth xy 1)))
(def y (apply list (nth xy 2)))

(defmacro P-fn [args1 args2]
  `(fn ~args1 (up (fn ~args2 ~x) (fn ~args2 ~y))))

(def rotating-path-2d (P-fn [Omega] [t]))

(let [NegativeOm -3]
(points-xy->plotly {"rotate" (make-path-txyz (rotating-path-2d NegativeOm) 0 3 64)}))

0.5s

Language:Clojure

2{:data(1):layout{5}}

(.plot js/Nextjournal plotly-5↩)

0.0s

Language:ClojureScript

Show that this path indeed satiesfies the equations of motion in a rotating coordinate system.

(simplify-alt
   (let [Om 'Omega
         NegativeOm (* -1 Om)]
     (((Lagrange-equations (L-rotating-rectangular-simp 'm Om))
       (rotating-path-2d NegativeOm))
      't)))

7.3s

Language:Clojure

3(down,0,0)

Systems with rigid constraints (Section 1.6.2)

A pendulum driven at the pivot

State Lagrange’s equation for this system.

(defn T-pend
  [m l _ ys]
  (fn [local]
    (let [[t theta thetadot] local
          vys (D ys)]
      (* 1/2 m
         (+ (square (* l thetadot))
            (square (vys t))
            (* 2 l (vys t) thetadot (sin theta)))))))

(defn V-pend
  [m l g ys]
  (fn [local]
    (let [[t theta _] local]
      (* m g (- (ys t) (* l (cos theta)))))))

(def L-pend (- T-pend V-pend))

(def theta (literal-function 'theta))
(def y (literal-function 'y))

(print (->TeX
      (simplify-alt
        (((Lagrange-equations (L-pend 'm 'l 'g y)) theta) 't))))

0.8s

Language:Clojure

g\,l\,m\,\sin\left(\theta\left(t\right)\right) + {l}^{2}\,m\,{D}^{2}\theta\left(t\right) + l\,m\,{D}^{2}y_s\left(t\right)\,\sin\left(\theta\left(t\right)\right)

State the Lagrangian

(print (->TeX
      (simplify-alt
        ((L-pend 'm 'l 'g y) (->local 't 'θ 'θdot)))))

0.6s

Language:Clojure

\frac{1}{2}\,{l}^{2}\,m\,{\dot {\theta}}^{2} + l\,m\,\dot {\theta}\,Dy_s\left(t\right)\,\sin\left(\theta\right) + g\,l\,m\,\cos\left(\theta\right) - g\,m\,y_s\left(t\right) + \frac{1}{2}\,m\,{\left(Dy_s\left(t\right)\right)}^{2}

Constraints as Coordinate Transformations (Section 1.6.3)

Derive the previous result by using coordinate transformations.

(defn L-uniform-acceleration-a [m g]
  (fn [[_ [_ y] v]]
    (- (* 1/2 m (square v)) (* m g y))))

(defn dp-coordinates [l y_s]
  (fn [[t θ]]
    (let [x (* l (sin θ))
          y (- (y_s t) (* l (cos θ)))]
      (up x y))))

(defn L-pend2 [m l g y_s]
  (comp (L-uniform-acceleration-a m g)
        (F->C (dp-coordinates l y_s))))
(print (->TeX
      (simplify-alt
        ((L-pend2 'm 'l 'g y) (->local 't 'θ 'θdot)))))

0.5s

Language:Clojure

\frac{1}{2}\,{l}^{2}\,m\,{\dot {\theta}}^{2} + l\,m\,\dot {\theta}\,Dy_s\left(t\right)\,\sin\left(\theta\right) + g\,l\,m\,\cos\left(\theta\right) - g\,m\,y_s\left(t\right) + \frac{1}{2}\,m\,{\left(Dy_s\left(t\right)\right)}^{2}

Central Forces in Three Dimensions (Section 1.8.3)

Calculate the z-component of the angular momentum of an arbitrary path in rectangular and spherical coordinates.

(def rectangular-state (up 't
                           (up 'x 'y 'z)
                           (up 'xdot 'ydot 'zdot)))

(def spherical-state (up 't
                         (up 'r 'θ 'φ)
                         (up 'rdot 'θdot 'φdot)))

(defn ang-mom-z [m]
  (fn [[_ xyz v]]
    (get (cross-product xyz (* m v)) 2)))

(print
  (make-tex-matrix
    [((ang-mom-z 'm) rectangular-state)
     ((compose (ang-mom-z 'm) (F->C s->r)) spherical-state)]))

0.7s

Language:Clojure

\begin{matrix} m\,x\,\dot y - m\,\dot x\,y \\ m\,{r}^{2}\,\dot {\varphi}\,{\sin}^{2}\left(\theta\right) \\ \end{matrix}

Using sherical coordinates, calculate the generalized forces and the generalized momenta of a planet moving in a central potential. Thus show that the momentum conjugate to the third coordinate ϕ\phiϕis (1) conserved (because the respective force is zero) and (2) identical the z-component of the angular momentum.

(def U (literal-function 'U))

(defn T3-spherical [m]
 (compose (L-free-rectangular m) (F->C s->r)))

(defn L3-central [m Vr]
  (let [Vs (fn [[_ [r]]] (Vr r))]
    (- (T3-spherical m) Vs)))

(print
  (make-tex-matrix
    [(((∂ 1) (L3-central 'm U)) spherical-state)
     (((∂ 2) (L3-central 'm U)) spherical-state)]))

1.5s

Language:Clojure

\begin{matrix} \begin{bmatrix}m\,r\,{\dot {\varphi}}^{2}\,{\sin}^{2}\left(\theta\right) + m\,r\,{\dot {\theta}}^{2} - DV\left(r\right)&m\,{r}^{2}\,{\dot {\varphi}}^{2}\,\sin\left(\theta\right)\,\cos\left(\theta\right)&0\end{bmatrix} \\ \begin{bmatrix}m\,\dot r&m\,{r}^{2}\,\dot {\theta}&m\,{r}^{2}\,\dot {\varphi}\,{\sin}^{2}\left(\theta\right)\end{bmatrix} \\ \end{matrix}

Show that the energy state function computed from the Lagrangian for a central field is in fact T + V.

(print
  (make-tex-matrix
   [(simplify-alt
     ((T3-spherical 'm) (->local 't (up 'r 'θ 'φ) (up 'rdot 'θdot 'φdot))))
    (simplify-alt
        ((Lagrangian->energy (L3-central 'm U)) spherical-state))]))

1.2s

Language:Clojure

\begin{matrix} \frac{1}{2}\,m\,{r}^{2}\,{\dot {\varphi}}^{2}\,{\sin}^{2}\left(\theta\right) + \frac{1}{2}\,m\,{r}^{2}\,{\dot {\theta}}^{2} + \frac{1}{2}\,m\,{\dot r}^{2} \\ \frac{1}{2}\,m\,{r}^{2}\,{\dot {\varphi}}^{2}\,{\sin}^{2}\left(\theta\right) + \frac{1}{2}\,m\,{r}^{2}\,{\dot {\theta}}^{2} + \frac{1}{2}\,m\,{\dot r}^{2} + V\left(r\right) \\ \end{matrix}

5.
Clojure Spec Definitions for this Notebook

The definitions are minimal and closely knit to the code above. To activate them, uncomment the last (ts/instrument) expression

(require '[sicmutils.value :as vl]
         '[sicmutils.numsymb :as ny]
         '[sicmutils.generic :as gn]
         '[sicmutils.structure :as st]
         '[sicmutils.simplify :as sp] ;;necessary to load simplify multifunction
         '[sicmutils.value :as vl]
         '[sicmutils.function :as fu] ;;necessary for ((gn/* 5 (fn[x] x)) 4)
         '[sicmutils.expression :as ex]
         '[sicmutils.calculus.derivative :as dr]
         '[sicmutils.numerical.minimize :as mn])

0.0s

Language:Clojure

0.0s

Language:Clojure

(s/def ::nu number?)
(s/def ::sy symbol?)
(s/def ::nu-sy (s/or :nu ::nu :sy ::sy))
(s/def ::expression
      (s/or :exp2 ::expression-2 :exp-q ::expression-q :exp-U ::expression-U
            :exp3 ::expression-3 :exp2-7 ::expression-gt2 :ns ::nu-sy))

(s/def ::expression-gt2
      (s/and seq?
             #(> (count %) 2)
             (s/cat :a-gt2 #{'* '+} :rest (s/+ ::expression))))
(s/def ::expression-3
      (s/and seq?
             #(= (count %) 3)
             (s/cat :a-3 #{'- '/ 'expt} :rest (s/+ ::expression))))
(s/def ::expression-U
      (s/and seq?
             #(= (count %) 2)
             (s/cat :a-U #{'U '(D U) '((expt D 2) U)}
                    :rest ::expression)))

(s/def ::expression-2
      (s/and seq?
             #(= (count %) 2)
             (s/cat :a-2 #{'sin 'cos '- 'sqrt 'atan} :rest ::expression)))

(s/def ::expression-q
      (s/and seq?
             #(= (count %) 2)
             (s/cat :a-q #{'q '(D q) '((expt D 2) q)
                           'x '(D x) '((expt D 2) x)
                           'y '(D y) '((expt D 2) y)
                           'z '(D z) '((expt D 2) z)
                           'r '(D r) '((expt D 2) r)
                           'phi '(D phi) '((expt D 2) phi)
                           'theta '(D theta) '((expt D 2) theta)}
                    :rest ::nu-sy)))

(s/def ::type #{:sicmutils.expression/numerical-expression})
(s/def ::eexpression (s/keys :req-un [::type ::expression]))
(s/def ::nu-sy-eex (s/or :ns ::nu-sy :eex ::eexpression))

(s/def ::differential (s/and #(instance?
                              sicmutils.calculus.derivative.Differential %)
                            #(s/valid? ::nu-sy-eex (second (first (.terms %))))))

(s/def ::nu-sy-eex-dr (s/or :ns ::nu-sy :dr ::differential :eex ::eexpression))

(s/def ::up-args (s/cat :a-up ::nu-sy-eex-dr-fn-up
                        :b-up ::nu-sy-eex-dr-fn-up
                        :c-up (s/? ::nu-sy-eex-dr-fn-up)))

(s/def ::up (s/and #(instance? sicmutils.structure.Structure %)
                   #(> (count %) 1)
                   #(< (count %) 4)
                   #(#{:sicmutils.structure/up} (vl/kind %))
                   #(s/valid? ::up-args (seq %))))

(s/def ::nu-sy-eex-dr-fn-up (s/or :nse ::nu-sy-eex :fn fn? :up ::up :dr ::differential))
(s/def ::dow-args (s/cat :a-dow ::nu-sy-eex-dr-fn-up
                         :b-dow ::nu-sy-eex-dr-fn-up
                         :c-dow (s/? ::nu-sy-eex-dr-fn-up)))

(s/def ::dow (s/and #(instance? sicmutils.structure.Structure %)
                    #(> (count %) 1)
                    #(< (count %) 4)
                    #(#{:sicmutils.structure/down} (vl/kind %))
                    #(s/valid? ::dow-args (seq %))))

(s/def ::nu-sy-eex-dr-fn-up-dow (s/or :nsedfu ::nu-sy-eex-dr-fn-up :dow ::dow))

(s/def ::sfunction #(instance? sicmutils.function.Function %))
(s/def ::nu-sy-eex-dr-fn-up-sfn (s/or :nsedfu ::nu-sy-eex-dr-fn-up :sf ::sfunction))
(s/def ::up-fnargs (s/cat :a-uf ::nu-sy-eex-dr-fn-up-sfn
                          :opt (s/? (s/cat
                                     :b-uf ::nu-sy-eex-dr-fn-up-sfn
                                     :c-uf (s/? ::nu-sy-eex-dr-fn-up-sfn)))))

0.0s

Language:Clojure

:sicmutils.env/up-fnargs

0.0s

Language:Clojure

(s/fdef ny/mul :args (s/cat :a-mul ::expression :b-mul ::expression))
(s/fdef gn/bin* :args (s/cat :a-bmul ::nu-sy-eex-dr-fn-up-dow :b-bmul ::nu-sy-eex-dr-fn-up-dow))
(s/fdef ny/add :args (s/cat :a-add ::expression :b-add ::expression))
(s/fdef gn/bin+ :args (s/cat :a-badd ::nu-sy-eex-dr-fn-up :b-badd ::nu-sy-eex-dr-fn-up))
(s/fdef ny/sub :args (s/cat :a-sub ::expression :b-sub ::expression))
(s/fdef gn/bin- :args (s/cat :a-sub ::nu-sy-eex-dr-fn-up-dow :b-sub ::nu-sy-eex-dr-fn-up-dow))
(s/fdef ny/div :args (s/cat :a-div ::expression :b-div ::expression))
(s/fdef gn/bin-div :args (s/cat :a-bdiv ::nu-sy-eex-dr-fn-up :b-bdiv ::nu-sy-eex-dr-fn-up))
(s/fdef gn/square :args (s/cat :a-squre (s/or :nse ::nu-sy-eex :up ::up :dr ::differential)))
(s/fdef ny/expt :args (s/cat :a-expt (s/or :dr ::differential :up ::up :exp ::expression) :expt2 int?))
(s/fdef ny/sine :args (s/cat :a (s/or :dr ::differential :ex ::expression)))
(s/fdef gn/sin :args (s/cat :a ::nu-sy-eex-dr))
(s/fdef ny/cosine :args (s/cat :a (s/or :dr ::differential :ex ::expression)))
(s/fdef gn/cos :args (s/cat :a ::nu-sy-eex-dr))
(s/fdef st/up :args ::up-fnargs)

(s/fdef lg/Lagrangian-action
  :args (s/cat :lagrangian fn? :test-path fn? :start-time number? :end-time number?))

(s/fdef mn/minimize
  :args (s/alt
         :three (s/cat :a3-min fn? :b3-min number? :c3-min number?)
         :four (s/cat :a4-min fn? :b4-min number? :c4-min number? :d4-min any?)))

(s/fdef lg/make-path :args
        (s/cat :t0 number? :q0 number?
               :t1 number? :q1 number?
               :qs (s/spec #(s/valid? (s/+ number?) (into [] %)))))

(s/fdef lg/parametric-path-action :args
        (s/cat :a1 fn? :a2 number? :a3 number? :a4 number? :a5 number?))

(s/fdef mn/multidimensional-minimize :args
        (s/alt :two (s/cat :a1 fn? :a2 (s/spec (s/+ number?)))
               :three (s/cat :a1 fn? :a2 (s/spec (s/+ number?)) :a3 any?)))

(s/fdef lg/find-path :args (s/cat :a1 fn? :a2 number?
                                  :a3 number? :a4 number? :a5 number? :a6 int?))

(s/fdef lg/linear-interpolants :args (s/cat :a1 number? :a2 number? :a3 int?))
(s/fdef lg/L-harmonic :args (s/cat :a1 ::nu-sy :a2 ::nu-sy))
(s/fdef lg/Lagrange-equations :args (s/cat :a1 fn?))
(s/fdef lg/L-uniform-acceleration :args (s/cat :a1 ::nu-sy :a2 ::nu-sy))
(s/fdef lg/L-central-rectangular :args
        (s/cat :a1 ::nu-sy-eex :a2 (s/or :o1 ::sfunction :o2 fn?)))

(s/fdef lg/L-central-polar :args
        (s/cat :a1 ::nu-sy-eex :a2 (s/or :o1 ::sfunction :o2 fn?)))

(s/fdef lg/F->C :args (s/cat :a fn?))
(s/fdef lg/L-free-particle :args (s/cat :a ::nu-sy))
(s/fdef lg/Gamma :args (s/cat :a fn?))
(s/fdef lg/->local :args (s/cat :a symbol?
                                :b (s/or :s symbol? :u ::up)
                                :c (s/or :s symbol? :u ::up)))

0.0s

Language:Clojure

lg/->local

(s/def ::up-down-expression
  (s/and seq?
         #(> (count %) 2)
         #(< (count %) 5)
         (s/cat :a-udex #{'up 'down}
                :rest (s/+ (s/or :e ::expression :udex ::up-down-expression)))))

;;simplify-spec-full is not used, it checks the result of a simplify call
;;but this check would not be deactiveted with (ts/unstrument)
(defn simplify-spec-full [x]
      (let [s (ts/with-instrument-disabled (gn/simplify x))]
          (do
              (if-not  (s/valid? (s/or :e ::expression :udex ::up-down-expression) s)
                (throw (Exception. "simplify result does not conform to spec")))
              s)))

(defn simplify-spec [x]
 (ts/with-instrument-disabled (gn/simplify x)))

(s/fdef simplify-spec :args (s/cat :a (s/or :e ::eexpression :u ::up :d ::dow)))

0.0s

Language:Clojure

sicmutils.env/simplify-spec

;;(ts/instrument)

0.0s

Language:Clojure

18[sicmutils.numsymb/mul,sicmutils.numsymb/add,sicmutils.env/simplify-spec,sicmutils.numsymb/cosine,sicmutils.generic/bin*,sicmutils.numsymb/sub,sicmutils.structure/up,sicmutils.numerical.minimize/multidimensional-minimize,sicmutils.generic/cos,sicmutils.numerical.minimize/minimize,sicmutils.numsymb/div,sicmutils.numsymb/expt,sicmutils.generic/bin+,sicmutils.generic/square,sicmutils.numsymb/sine,sicmutils.generic/sin,sicmutils.generic/bin-div,sicmutils.generic/bin-]

Shift+Enter to run

Language:Clojure

Structure and Interpretation of Classical Mechanics in Clojure

1. Introduction

2. Computing Actions (Chapter 1.4)

3. The Euler-Lagrange Equations (Chapter 1.5)

3.1. Computing Lagrange's Equations (Section 1.5.2)

4. How to find Lagrangians (Chapter 1.6)

5. Clojure Spec Definitions for this Notebook

1.
Introduction

2.
Computing Actions (Chapter 1.4)

3.
The Euler-Lagrange Equations (Chapter 1.5)

3.1.
Computing Lagrange's Equations (Section 1.5.2)

4.
How to find Lagrangians (Chapter 1.6)

5.
Clojure Spec Definitions for this Notebook