;;; Double Pendulum Project (1.39)
;;; James Rising
;;; A) Generate the Lagrangian
;; My technique is to split the problem up into a description of the
;; effect of gravity on two objects, convolved with a coordinate transform
;; which places the double pendulum constraints between the objects.
(define ((L-uniform-acc2 m1 m2 g) local)
(let ((q1 (ref (coordinate local) 0))
(q2 (ref (coordinate local) 1))
(v1 (ref (velocity local) 0))
(v2 (ref (velocity local) 1)))
(let ((y1 (ref q1 1))
(y2 (ref q1 1)))
(- (+ (* 1/2 m1 (square v1)) (* 1/2 m2 (square v2)))
(+ (* m1 g y1) (* m2 g y2))))))
;Value: L-uniform-acc2
(define ((dp-coordinates l1 l2) local)
(let* ((t (time local))
(theta1 (ref (coordinate local) 0))
(theta2 (+ theta1 (ref (coordinate local) 1))))
(let ((x1 (* l1 (sin theta1)))
(y1 (- (* l1 (cos theta1)))))
(let ((x2 (+ x1 (* l2 (sin theta2))))
(y2 (- y1 (* l2 (cos theta2)))))
(up (up x1 y1) (up x2 y2))))))
;Value: dp-coordinates
(define (L-double-pend m1 m2 l1 l2 g)
(compose (L-uniform-acc2 m1 m2 g)
(F->C (dp-coordinates l1 l2))))
;Value: L-double-pend
(show-expression
((L-double-pend 'm_1 'm_2 'l_1 'l_2 'g)
(->local 't (up 'theta_1 'theta_2) (up 'thetadot_1 'thetadot_2))))
(+
(* l_1 l_2 m_2 (expt thetadot_1 2) (cos (+ theta_1 theta_2)) (cos theta_1))
(* l_1 l_2 m_2 (expt thetadot_1 2) (sin theta_1) (sin (+ theta_1 theta_2)))
(* l_1 l_2 m_2 thetadot_1 thetadot_2 (cos (+ theta_1 theta_2)) (cos theta_1))
(* l_1 l_2 m_2 thetadot_1 thetadot_2 (sin theta_1) (sin (+ theta_1 theta_2)))
(* 1/2 (expt l_1 2) m_1 (expt thetadot_1 2))
(* 1/2 (expt l_1 2) m_2 (expt thetadot_1 2))
(* 1/2 (expt l_2 2) m_2 (expt thetadot_1 2))
(* (expt l_2 2) m_2 thetadot_1 thetadot_2)
(* 1/2 (expt l_2 2) m_2 (expt thetadot_2 2))
(* g l_1 m_1 (cos theta_1))
(* g l_1 m_2 (cos theta_1)))
;Unspecified return value
;; The state derivative is needed for numerical integration
(define (double-pend-state-derivative m1 m2 l1 l2 g)
(Lagrangian->state-derivative
(L-double-pend m1 m2 l1 l2 g)))
;Value: double-pend-state-derivative
;; And this is the energy of the system
(define energy-double-pend
(Lagrangian->energy (L-double-pend 1.0 3.0 1.0 .9 9.8)))
;;; B) and C) Graphing the system
;; monitor-thetas graphs both angles and the energy
(define ((monitor-thetas win1 win2 wine) state)
(let ((theta1 ((principal-value :pi) (ref (coordinate state) 0)))
(theta2 ((principal-value :pi) (ref (coordinate state) 1)))
(energy (energy-double-pend state)))
(plot-point win1 (time state) theta1)
(plot-point win2 (time state) theta2)
(plot-point wine (time state) energy)))
;Value: monitor-thetas
;; Create three windows; I will not include these lines from now on
(define plot1-win (frame 0. 50. :-pi :pi))
;Value: plot1-win
(define plot2-win (frame 0. 50. :-pi :pi))
;Value: plot2-win
(define plote-win (frame 0. 50. -1.0e-10 1.0e-10))
;Value: plote-win
;; Graph the system
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) :pi)
(up 0 0))
(monitor-thetas plot1-win plot2-win plote-win)
.001
50.
1.0e-13)
;; The energy of the system is almost perfectly conserved; the error
;; rarely goes beyond 1.0e-11
;; Now I want to save the data at .125 second intervals
;; I will save the states to a list and use them later
(define nil '())
(define states nil)
(define (save-state state)
(set! states (cons state states)))
;; Evolve the system again, this time just saving the data
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) :pi)
(up 0 0))
save-state
.125
50.
1.0e-13)
;; Reverse the list and save the data elsewhere
(define states1 (reverse states))
;; Clear the list
(define states nil)
;; D) Graph the system again, for m_2 10^-10 m higher
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) (- :pi (atan 1.0e-10 .9)))
(up 0 0))
(monitor-thetas plot1-win plot2-win plote-win)
.01
50.
1.0e-13)
;; And save the data
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) (- :pi (atan 1.0e-10 .9)))
(up 0 0))
save-state
.125
50.
1.0e-13)
(define states2 (reverse states))
(define states nil)
;; log-squared-diffs walks down the lists and generates a new list
;; of the log of the squares of the differences
(define (log-squared-diffs l1 l2)
(if (null? l1)
nil
(cons (log (square (- (ref (coordinate (car l1)) 1)
(ref (coordinate (car l2)) 1))))
(log-squared-diffs (cdr l1) (cdr l2)))))
(define states-log (log-squared-diffs states1 states2))
;; Now graph the log squared differences
(define plotl-win (frame 0. 50. -50. 50.))
;Value: plotl-win
(define (graph-points l time step win)
(if (not (null? l))
(begin
(plot-point win time (car l))
(graph-points (cdr l) (+ time step) step win))))
(graph-points states-log 0. .125 plotl-win)
;; Initially, the differences are starting very small and on the log graph
;; increase linearly (so they increase exponentially) until a certain point
;; (when they are sufficiently different) where the log difference becomes
;; constant.
;;; E) Repeat parts (B) and (D) with new initial conditions
;; Graph the system
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) 0.)
(up 0 0))
(monitor-thetas plot1-win plot2-win plote-win)
.01
50.
1.0e-13)
;; Now save the data
(define states nil)
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) 0.)
(up 0 0))
save-state
.125
50.
1.0e-13)
(define states1 (reverse states))
;; And again, slightly higher
;; graph it
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) (atan 1.0e-10 .9))
(up 0 0))
(monitor-thetas plot1-win plot2-win plote-win)
.01
50.
1.0e-13)
;; save it
(define states nil)
((evolve double-pend-state-derivative
1.0
3.0
1.0
.9
9.8)
(up 0.0
(up (/ :pi 2) (atan 1.0e-10 .9))
(up 0 0))
save-state
.125
50.
1.0e-13)
(define states2 (reverse states))
;; Generate the log plot
(define states-log (log-squared-diffs states1 states2))
(define plotl-win (frame 0. 50. -50. 50.))
;Value: plotl-win
(graph-points states-log 0. .125 plotl-win)
;; Again, the log difference increases linearly (but more slowly this time)
;; until a certain point when it becomes more or less constant.