Assignment 3

assignment3.scm

;; Ethan Sherbondy
;; 6.946, Fall 2012
;; Project Exercise 1.39

;; The Double Pendulum

(define ((L-rect-pend m1 m2 g) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (let* ((v1 (ref v 0))
	   (v2 (ref v 1))
	   (q1 (ref q 0))
    	   (q2 (ref q 1))
           (y1 (ref q1 1))
    	   (y2 (ref q2 1)))
      (+ (* m1 (- (* 1/2 (square v1)) (* g y1)))
	 (* m2 (- (* 1/2 (square v2)) (* g y2)))))))

(define ((p->r l1 l2) local)
  (let* ((t (time local))
         (theta (coordinate local))
	 (theta1 (ref theta 0))
	 (theta2 (ref theta 1)))
    (let* ((x1 (* l1 (sin theta1)))
	   (x2 (+ x1 (* l1 (sin (+ theta1 theta2)))))
           (y1 (- 0 (* l1 (cos theta1))))
	   (y2 (- y1 (* l2 (cos (+ theta1 theta2))))))
      (up (up x1 y1) (up x2 y2))))))

(define (L-dpend m1 m2 l1 l2 g)
  (compose (L-rect-pend m1 m2 g)
           (F->C (p->r l1 l2))))

(define gamma (up 't (up 'theta_1 'theta_2) (up 'thetadot_1 'thetadot_2)))

(se
 ((L-dpend 'm1 'm2 'l1 'l2 'g)
  (->local 't (up 'theta_1 'theta_2) (up 'thetadot_1 'thetadot_2))))

(define (pend-state-derivative m1 m2 l1 l2 g)
  (Lagrangian->state-derivative
   (L-dpend m1 m2 l1 l2 g)))

(pe
 ((pend-state-derivative 'm1 'm2 'l1 'l2 'g)
  gamma))

((pend-state-derivative 1.0 3.0 1.0 0.9 9.8) (up 0 (up :pi/2 :pi) (up 0 0)))

(define plot-win (frame 0. 50. :-pi :pi))

(define ppi (principal-value :pi))

;; This is yielding floating point overflows right now...
;; L-dpend is probably incorrectly defined
(define ((monitor-theta win) state)
  (let ((theta (coordinate state)))
    (let ((theta1 (ppi (ref theta 0)))
	  (theta2 (ppi (ref theta 1))))
      (plot-point win (time state) theta1)
      (plot-point win (time state) theta2))))

((evolve pend-state-derivative
	 1.0      ;m1 = 1.0kg
	 3.0      ;m2 = 3.0kg
	 1.0      ;l1 = 1.0m
	 0.9      ;l2 = 0.9m
	 9.8)     ;g  = 9.8m/s^2
 (up 0 ; t0 = 0
     (up :pi/2    ;theta0 = pi/2 radians
	 :pi)     ;theta1 = pi radians
     (up 0 0))    ;thetadot 0 and 1 are 0 rad/s
 (monitor-theta plot-win)
 0.01
 50
 1.0e-4)


;; Problem 1.29:
;; A particle of mass m slides off a horizontal cylinder of radius R 
;; in a uniform gravitational field with acceleration g. If the particle starts 
;; close to the top of the cylinder with zero initial speed, with what angular 
;; velocity does it leave the cylinder?

;; The particle will fall off the cylinder when the acceleration due to
;; gravity is equal to the centripetal acceleration.

;; The particle only actually exhibits constrained behavior while it's
;; on the surface of the cylinder. Luckily, we're interested in the
;; angular velocity at the very moment it falls off, so the constraints are valid.
(define ((L-cyl m g) local)
  (let ((coords (coordinate local))
     	(v (velocity local)))
    (let ((y (ref coords 1)))
      (- (* 1/2 m (square v))
	 (* m g y)))))

(define (p->r local)
  (let ((coords  (coordinate local)))
    (let ((r     (ref coords 0))
	  (theta (ref coords 1)))
      (let ((y (* r (sin theta)))
	    (x (* r (cos theta))))
	(up x y)))))


(define (L-cyl-polar m g)
  (compose (L-cyl m g) (F->C p->r)))


(se ((L-cyl-polar 'm 'g) 
     (->local 't (up 'r 'theta) (up 'rdot 'thetadot))))
#|
(+ (* 1/2 (expt r 2) m (expt thetadot 2))
   (* -1 g r m (sin theta))
   (* 1/2 m (expt rdot 2)))
|#

;; Angular momentum is conserved
(se (((partial 2) 
      (L-cyl-polar 'm 'g)) 
     ((Gamma (up (lf 'r) (lf 'theta))) 't)))
#|
(down (* m ((D r) t)) 
      (* m (expt (r t) 2) ((D theta) t)))
|#


(se
 (((Lagrange-equations
    (L-cyl-polar 'm 'g))
   (up (lf 'r)
       (lf 'theta)))
  't))
#|
(down
 (+ (* -1 m (r t) (expt ((D theta) t) 2))
    (* g m (sin (theta t)))
    (* m (((expt D 2) r) t)))
 (+ (* g m (cos (theta t)) (r t))
    (* m (expt (r t) 2) (((expt D 2) theta) t))
    (* 2 m (r t) ((D theta) t) ((D r) t))))
|#

;; The point mass will fall off when the force due to gravity
;; exceeds the centripetal force.
;; This happens when v^2/R = gsin(theta)
;; v = R*(D theta)

;; The following is equal to mgR by conservation of energy
((Lagrangian->energy 
  (L-cyl-polar 'm 'g)) 
 ((Gamma (up (lf 'r) (lf 'theta))) 't))
#|
(+ (* 1/2 m (expt (r t) 2) (expt ((D theta) t) 2))
   (* g m (r t) (sin (theta t)))
   (* 1/2 m (expt ((D r) t) 2)))
|#

;; the first expression of the sum is equivalent to
;; 1/2mgsin(theta)*R + mgRsin(theta) = mgR
;; so sin(theta) = 2/3
;; from here, we can plug theta into the Lagrange equations to
;; determine D(theta) in terms of m, g, and r at the point of separation
#|
 (+ (* -1 m (r t) (expt ((D theta) t) 2))
    (* g m (sin (theta t)))
    (* m (((expt D 2) r) t))) = 0
=
  (+ (* -1 m r (expt ((D theta) t) 2))
     (* g m 2/3)) = 0

AND HERE'S OUR SOLUTION:
(D theta) = (sqrt (* 2/3 (/ g r)))
|#


;; Problem 1B
;; Here I've arbitrarily decided to make gravity act on the x-axis

(define ((L-ellipsoid-rect m g) local)
  (let ((velo   (velocity local))
	(coords (coordinate local)))
    (let ((x (ref coords 0)))
      (- (* 1/2 m (square coords))
	 (* m g x)))))

(define ((s->r a b c) local)
  (let ((coords (coordinate local)))
    (let ((u (ref coords 0))
	  (v (ref coords 1)))
      (let ((x (* a (cos u) (cos v)))
	    (y (* b (cos u) (sin v)))
	    (z (* c (sin u))))
	(up x y z)))))

(define (L-ellipsoid-s m g a b c)
  (compose (L-ellipsoid-rect m g)
	   (F->C (s->r a b c))))

(define L-eq-bc
  (L-ellipsoid-s 'm 'g 'a 'b 'b))

(se 
 (L-eq-bc
  (->local 't (up 'u 'v) (up 'udot 'vdot))))
#|
(+ (* 1/2 (expt a 2) m (expt (cos v) 2) (expt (cos u) 2))
   (* -1/2 (expt b 2) m (expt (cos v) 2) (expt (cos u) 2))
   (* -1 a g m (cos v) (cos u))
   (* 1/2 (expt b 2) m))
|#

;; Since this ellipsoid has rotational symmetry about the x-axis,
;; I expect R_x to capture the symmetry

(define coords
  (->local 't (up 'u 'v) (up 'udot 'vdot)))

(define rotated-coords
  ((Rx 's) ((s->r 'a 'b 'b) coords)))
#|
(up (* a (cos v) (cos u))
    (+ (* b (cos u) (sin v) (cos s)) (* -1 b (sin u) (sin s)))
    (+ (* b (cos u) (sin v) (sin s)) (* b (sin u) (cos s))))
|#

(square rotated-coords)
#|
(+ (* (expt a 2) (expt (cos v) 2) (expt (cos u) 2))
   (* -1 (expt b 2) (expt (cos v) 2) (expt (cos u) 2))
   (expt b 2))
|#

(square ((s->r 'a 'b 'b) coords))
#|
(+ (* (expt a 2) (expt (cos v) 2) (expt (cos u) 2))
   (* -1 (expt b 2) (expt (cos v) 2) (expt (cos u) 2))
   (expt b 2))
|#

;; The squared coordinates, are identical across the
;; transformation. We expect velocities to be preserved also.
;; and the x-coordinate is preserved (our GPE coordinate),
;; so L' = L

(define (DRx local)
 (((D Rx) 0) 
  ((s->r 'a 'b 'b) local)))


(define Noether-integral
  (let ((L L-eq-bc))
    (* ((partial 2) L) DRx)))

(Noether-integral
 (->local 't (up 'u 'v) (up 'udot 'vdot)))

#|
(up (down 0 0) (down 0 0) (down 0 0))
|#