The Double Pendulum


In [1]:
(require '[clojupyter.misc.helper :as helper])
(run! helper/add-dependencies '[[org.clojure/data.json "0.2.6"]
                                [net.littleredcomputer/sicmutils "0.12.0-SNAPSHOT"]
                                [thi.ng/geom "0.0.908"]])
(ns double-pendulum
    (:refer-clojure :exclude [partial zero? + - * / ref])
    (:require [sicmutils.env :refer :all]
              [clojupyter.misc.display :as display]
              [thi.ng.geom.viz.core :as viz]
              [thi.ng.geom.svg.core :as svg]))
(defn tex [x] (-> x tex$$ display/latex))


Out[1]:
#'double-pendulum/tex

In [2]:
(defn V
  [m_1 m_2 l_1 l_2 g]
  (fn [[_ [theta phi] _]]
    (let [y_1 (- (* l_1 (cos theta)))
          y_2 (- y_1 (* l_2 (cos phi)))]
      (+ (* m_1 g y_1)
         (* m_2 g y_2)))))


Out[2]:
#'double-pendulum/V

In [3]:
(defn T
  [m_1 m_2 l_1 l_2 _]
  (fn [[_ [theta phi] [thetadot phidot]]]
    (let [v1sq (* (square l_1) (square thetadot))
          v2sq (* (square l_2) (square phidot))]
      (+ (* 1/2 m_1 v1sq)
         (* 1/2 m_2 (+ v1sq
                      v2sq
                      (* 2 l_1 l_2 thetadot phidot (cos (- theta phi)))))))))


Out[3]:
#'double-pendulum/T

In [4]:
(def L (- T V))
(tex ((L 'm_1 'm_2 'l_1 'l_2 'g) 
       (up 't 
           (up 'theta 'phi) 
           (up 'thetadot 'phidot))))


Out[4]:
$$l_1\,l_2\,m_2\,\dot {\phi}\,\dot {\theta}\,\cos\left(- \phi + \theta\right) + \frac{1}{2}\,{l_1}^{2}\,m_1\,{\dot {\theta}}^{2} + \frac{1}{2}\,{l_1}^{2}\,m_2\,{\dot {\theta}}^{2} + \frac{1}{2}\,{l_2}^{2}\,m_2\,{\dot {\phi}}^{2} + g\,l_1\,m_1\,\cos\left(\theta\right) + g\,l_1\,m_2\,\cos\left(\theta\right) + g\,l_2\,m_2\,\cos\left(\phi\right)$$

In [5]:
(defn state-derivative [m_1 m_2 l_1 l_2 g]
    (Lagrangian->state-derivative
        (L m_1 m_2 l_1 l_2 g)))


Out[5]:
#'double-pendulum/state-derivative

In [6]:
(def equations-of-motion
    ((state-derivative 'm_1 'm_2 'l_1 'l_2 'g)
               (up 't
                   (up 'theta 'phi)
                   (up 'thetadot 'phidot))))


Out[6]:
#'double-pendulum/equations-of-motion

In [7]:
(tex (nth equations-of-motion 2))


Out[7]:
$$\begin{pmatrix}\frac{- l_1\,m_2\,{\dot {\theta}}^{2}\,\cos\left(- \phi + \theta\right)\,\sin\left(- \phi + \theta\right) - l_2\,m_2\,{\dot {\phi}}^{2}\,\sin\left(- \phi + \theta\right) + g\,m_2\,\cos\left(- \phi + \theta\right)\,\sin\left(\phi\right) - g\,m_1\,\sin\left(\theta\right) - g\,m_2\,\sin\left(\theta\right)}{l_1\,m_2\,{\sin}^{2}\left(- \phi + \theta\right) + l_1\,m_1}\\\frac{l_2\,m_2\,{\dot {\phi}}^{2}\,\cos\left(- \phi + \theta\right)\,\sin\left(- \phi + \theta\right) + l_1\,m_1\,{\dot {\theta}}^{2}\,\sin\left(- \phi + \theta\right) + l_1\,m_2\,{\dot {\theta}}^{2}\,\sin\left(- \phi + \theta\right) + g\,m_1\,\cos\left(- \phi + \theta\right)\,\sin\left(\theta\right) + g\,m_2\,\cos\left(- \phi + \theta\right)\,\sin\left(\theta\right) - g\,m_1\,\sin\left(\phi\right) - g\,m_2\,\sin\left(\phi\right)}{l_2\,m_2\,{\sin}^{2}\left(- \phi + \theta\right) + l_2\,m_1}\end{pmatrix}$$

In [8]:
(def out (atom []))
(defn observe [t [_ [theta phi]]] (swap! out conj [t theta phi]))
(tex ((evolve state-derivative 1 1 1 1 9.8)
      (up 0 (up (/ pi 2) 0) (up 0 0))
      observe
      1/60
      10
      1e-6
      :compile true))


Out[8]:
$$\begin{pmatrix}9.999999999999954\\\begin{pmatrix}0.11375091645060717\\-2.223620451765166\end{pmatrix}\\\begin{pmatrix}-1.7171709187549014\\0.5300899328860362\end{pmatrix}\end{pmatrix}$$

In [13]:
(def spec
  {:x-axis (viz/linear-axis
            {:domain [(first (first @out)) (first (last @out))]
             :range  [50 (- 985 10)]
             :major 1
             :pos    550})
   :y-axis (viz/linear-axis
            {:domain      [(- Math/PI) Math/PI]
             :range       [550 20]
             :major       1
             :minor       0.2
             :pos         50
             :label-dist  15
             :label-style {:text-anchor "end"}})
   :grid   {:attribs {:stroke "#caa"}
            :minor-x true
            :minor-y false}
   :data   [{:values  (map #(vector (nth % 0) (nth % 1)) @out)
             :attribs {:stroke "#0af" :stroke-width "2pt"}
             :layout  viz/svg-line-plot}
            {:values  (map #(vector (nth % 0) (nth % 2)) @out)
             :attribs {:stroke "#f60" :stroke-width "2pt"}
             :layout  viz/svg-line-plot}]})


Out[13]:
#'double-pendulum/spec

In [14]:
(defn export-viz
  [spec]
  (->> spec
       (viz/svg-plot2d-cartesian)
       (svg/svg {:width 985 :height 600})
       (svg/serialize)))

(-> spec
    export-viz
    display/make-html)


Out[14]:
0.001.002.003.004.005.006.007.008.009.0010.00-3.00-2.00-1.000.001.002.003.00

In [ ]: