The initial boilerplate to start a sicmutils investigation:
In [1]:
(require '[clojupyter.misc.helper :as helper])
(run! helper/add-dependencies '[[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]
[sicmutils.examples.pendulum :as pendulum]
[clojupyter.misc.display :as display]
[thi.ng.geom.viz.core :as viz]
[thi.ng.geom.svg.core :as svg]
[thi.ng.color.core :as color]))
(defn tex [x] (-> x tex$$ display/latex))
Out[1]:
The pendulum library in the examples directory contains the Lagrangian for the simple pendulum with moving pivot. For this investigation, we want a pendulum with a fixed pivot, so we supply a constant function for the motion of the pivot. Let's double-check that the Lagrangian for the pendulum agrees with what we find on Wikipedia.
In [2]:
;; The coordinates of the pivot as a function of time
(defn fixed [t] (up 0 0))
;; Obtain a fixed-pivot Lagrangian from the moving-pivot Lagrangian
(defn L [m l g] (pendulum/L m l g fixed))
;; Bind the contants to symbolic values, and supply symbolic initial data
(tex ((L 'm 'l 'g) (up 't 'theta 'thetadot)))
Out[2]:
That looks correct. Now we form the state derivative:
In [4]:
(defn state-derivative [m l g]
(Lagrangian->state-derivative (L m l g)))
(tex ((state-derivative 'm 'l 'g) (up 't 'theta 'thetadot)))
Out[4]:
So far so good. The next step is to integrate the motion for various choices of initial displacement. We'll set $m = l = 1, g = 9.8$ here, and call the initial angle $\theta_0$. We will run the simulation for 4 seconds, sampling every 1/60 second. Following that is some plotting boilerplate. We resume commentary at the graph below.
In [5]:
(def g 9.8)
(defn p [theta_0]
(integrate-state-derivative state-derivative [1 1 g] (up 0 theta_0 0) 4 1/60))
Out[5]:
In [6]:
(def spec
{:x-axis (viz/linear-axis
{:domain [0 4]
: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 (for [a (range 0.15 2 0.15)]
{:values (p a)
:attribs {:stroke (color/as-css (color/rgba 0.8 (/ a 2) 0 1)) :stroke-width "2pt"}
:layout viz/svg-line-plot})
})
Out[6]:
Now let's look at the output. For a small displacement, the high-school (linearized) calculation shows we should expect a period of $2\pi\sqrt{\frac{l}g}\approx 2.0$ in our case with $l=1$ and $g=9.8$. As the initial angular displacement increases, however, the period begins to rise.
In [7]:
(->> spec
viz/svg-plot2d-cartesian
(svg/svg {:width 985 :height 600})
svg/serialize
display/html)
Out[7]:
In [8]:
(tex ((Lagrangian->energy (L 'm 'l 'g)) (up 't 'theta 'thetadot)))
Out[8]:
The energy of the system is conserved. Assuming we raise the bob to an angle $\theta_0$ and release it, the total energy must remain equal to the potential energy initially present:
In [9]:
(tex ((pendulum/V 'm 'l 'g fixed) (up 't 'theta_0 'thetadot_0)))
Out[9]:
Equating these expressions allows us to write $\dot\theta$ in terms of $\theta$ and $\theta_0$. We can then form the integral for the period $$\mathrm{period}(\theta_0)=4\int_0^{\theta_0}\frac{1}{\dot\theta}d\theta$$ We find that $$\mathrm{period}(\theta_0)=\frac{4}{\sqrt{2g}}\int_0^{\theta_0}\frac{1}{\sqrt{\cos\theta - \cos\theta_0}}d\theta$$
In [10]:
(defn integrand [theta_0]
(fn [theta] (/ (sqrt (- (cos theta) (cos theta_0))))))
(defn period-broken [theta_0]
(* (/ 4 (* 2 g)) (definite-integral (integrand theta_0) 0 (- theta_0 0.00001) )))
Out[10]:
We interrupt this broadcast to inform the viewers that the integrator we are currently using is not up to the task of computing this improper integral (the integrand blows up at the right-hand edge of the interval). I tried a few other integrators in Apache Commons Math, without success. Scmutils itself is equipped with integrators specially designed to handle integrals with singularities at the endpoints, but we have not ported them here. So we will take another tack: consider the solution of the integral in terms of $F$, the elliptic integral of the first kind: $$\mathrm{period}(\theta_0)=\frac{8}{\sqrt{2g}}\frac{F(\frac{\theta_0}{2},\csc\frac{\theta_0}{2})}{\sqrt{1-\cos(\theta_0)}}$$
In [11]:
(defn period [theta_0]
(let [t (/ theta_0 2)]
(/ (* 8 (elliptic-f t (csc t)))
(* (sqrt (* 2 g)) (sqrt (- 1 (cos theta_0)))))))
Out[11]:
In [12]:
(display/hiccup-html
[:table
[:thead [:th "θ<sub>0</sub><br/>radians"] [:th "period"]]
(for [A (range 0.15 2 0.15)]
[:tr [:td (format "%8.02f" A)] [:td (format "%8.04f" (period A))]])])
Out[12]:
In [ ]: