Last edited 11dec12 by aamurth2@illinois.edu
Find this document at http://new.math.uiuc.edu/math198/aamurth2
Double Pendulum 2012
Abstract
After reading about many different types of dynamical systems on the internet, I was intrigued by one particular dynamical system: the Double Pendulum and I decided to base my project on this system. This is the dynamical system of one pendulum connected at the end of another one but free to rotate. Its motion, which can be calculated with ordinary differential equations, can be chaotic depending on the external characteristics of the system.
Project Proposal
After deciding upon which dynamical system to focus on in my project, I attempted to design a project that would put the double pendulum to use. I have decided to model a double pendulum. I will do this by first tracing the position of the bar as it moves. Then I will plot the angles of each bar to the vertical using both torus and a wrap around screen.
The Math and Physics behind a Double Pendulum
The physics behind a single simple pendulum is fairly easy to calculate:
x = .5*l*sin θ and y = .5*l*cos θ
where l represents the length of the pendulum, θ represents the angle between the limb and the vertical, x represents the x-coordinate of the end of the pendulum, and y represents the y-coordinate of the end of the pendulum.
The physics behind the second pendulum is more difficult to model and the exact position of the end of the pendulum is given by the following equations:
x=l(sinθ+.5*l*sinφ)
and
y=−l(cosθ+.5*l*cosφ)
where l represents the length of the pendulum, θ represents the angle between the first limb and the vertical, φ represents the angle between the second limb and the vertical, x represents the x-coordinate of the end of the pendulum, and y represents the y-coordinate of the end of the pendulum.
By plotting the angles the bars make with the vertical at any time t, I will make it easier to visualize the chaotic motion of a double pendulum.
My Final Project
I successfully modified Bruce Sherwood’s double pendulum VPython program to draw the path of the end of the first pendulum, and that of the end of the second pendulum. In addition, I have learned to create a torus and a ball on VPython, and I have enlarged the VPython window and added these objects next to Sherwood's double pendulum. Then I used the ball to plot the angles of the double pendulum (angle 1 and angle 2) at any given time t. I also created a wrap around screen that just like the torus and ball model plots the angles of the double pendulum at any given time t.
My Project can be found in the following locations:
RTICA without Wrap-Around Screen
class198f12/aamurth2/finalproject/toroidalpendulum.py
RTICA with Wrap-Around Screen
class198f12/aamurth2/finalproject/boardpendulum.py
IMG="toruspic1"
IMG="toruspic2"
IMG="toruspic3"
IMG="doublepenpic"
boardpendulum.py
from __future__ import print_function, division
from visual.graph import *
# Double pendulum in a Torus
# Akshay Murthy
# An adaptation of Bruce Sherwood's program
# Draws the motion of the pendulum
# Plots the upper (theta1) and lower angles (theta2) within a torus
scene.title = 'Double Pendulum - Configuration Space'
scene.height = scene.width = 1000
scene.autoscale = true
g = 9.8 # acceleration due to gravity
M1 = 2.0 # mass of upper bar
M2 = 1.0 # mass of lower bar
L1 = 0.5 # length of upper bar
L2 = 1.0 # length of lower bar
I1 = M1*L1**2/12 # moment of inertia of upper assembly
I2 = M2*L2**2/12 # moment of inertia of lower bar
d = 0.05 # thickness of each bar
gap = 2*d # distance between two parts of upper, U-shaped assembly
L1display = L1+d # show upper assembly a bit longer than physical, to overlap axle
L2display = L2+d/2 # show lower bar a bit longer than physical, to overlap axle
hpedestal = 1.3*(L1+L2) # height of pedestal
wpedestal = 0.1 # width of pedestal
tbase = 0.05 # thickness of base
wbase = 8*gap # width of base
offset = 2*gap # from center of pedestal to center of U-shaped upper assembly
top = vector(0,0,0) # top of inner bar of U-shaped upper assembly
scene.center = top-vector(-2,(L1+L2)/2.,0) # creates the scene
theta1 = 2.1 # initial upper angle (from vertical)
theta2 = 2.7 # initial lower angle (from vertical)
theta1dot = 0 # initial rate of change of theta1
theta2dot = 0 # initial rate of change of theta2
pedestal = box(pos=top-vector(0,hpedestal/2.,offset),
height=1.1*hpedestal, length=wpedestal, width=wpedestal,
color=(0.4,0.4,0.5)) # creates the pedestal from which the pendulums hang
base = box(pos=top-vector(0,hpedestal+tbase/2.,offset),
height=tbase, length=wbase, width=wbase,
color=pedestal.color) # creates the base of the pedestal
axle = cylinder(pos=top-vector(0,0,gap/2.-d/4.), axis=(0,0,-offset),
radius=d/4., color=color.yellow) # creates the axle that connects pendulums to pedestal
frame1 = frame(pos=top) # sets the position of the frame
bar1 = box(frame=frame1, pos=(L1display/2.-d/2.,0,-(gap+d)/2.),
size=(L1display,d,d), color=color.red) # creates first pendulum
bar1b = box(frame=frame1, pos=(L1display/2.-d/2.,0,(gap+d)/2.),
size=(L1display,d,d), color=color.red) # creates another pendulum object
axle1 = cylinder(frame=frame1, pos=(L1,0,-(gap+d)/2.), axis=(0,0,gap+d),
radius=axle.radius, color=axle.color) # creates an axle coming off the first pendulum
frame1.axis = (0,-1,0) # sets the direction of frame 1
frame2 = frame(pos=frame1.axis*L1) # creates a new frame with bar 2
bar2 = box(frame=frame2, pos=(L2display/2.-d/2.,0,0),
size=(L2display,d,d), color=color.green) # creates second pendulum
frame2.axis = (0,-1,0) # sets the direction of frame 2
frame1.rotate(axis=(0,0,1), angle=theta1) # gives the first pendulum the ability to rotate
frame2.rotate(axis=(0,0,1), angle=theta2) # gives the second pendulum the ability to rotate
frame2.pos = top+frame1.axis*L1 # sets the position of frame 2
frame1.trail = curve(color=color.orange) # draws trail for end of second pendulum
frame2.trail = curve(color=color.blue) # draws trail for end of first pendulum
dt = 0.00005 # deltaT- the frequency at which the frames are updated (smaller numbers give more precise movement)
t = 0. # set time = 0 when the program is first run
C11 = (0.25*M1+M2)*L1**2+I1 # creates matrices needed for Lagrangian coordinates
C22 = 0.25*M2*L2**2+I2
ring(pos=(4,0,0), axis=(0,0,1), radius=1.0, thickness=.30,
opacity=.5, color=color.red) # creates torus
board = box(pos=(2.019,2,0), axis=(0,0,1), size=(0.05,1.25,1.25), opacity=.5, color=color.green) # creates wrap around screen
boardplotter = sphere(make_trail=True, trail_type="points",
interval=1, radius=.001, color=color.blue) # creates drawer that will draw on this screen
ball1 = sphere(pos=(4,0,0), radius=0.05, color=color.cyan) # creates a ball that will plot theta1 at any time t
ball1.trail = curve(color=color.white) # creates a method that will trace the ball's path
ball2 = sphere(pos=(4,0,0), radius=0.05, color=color.yellow) # creates a ball that will plot theta2 at any time t
ball2.trail = curve(color=color.green) # creates a method that will trace the ball's path
ball3 = sphere(pos=(4,0,0), radius=0.05, color=color.red) # creates a ball that will plot theta1 and theta2 at any time t
ball3.trail = curve(color=color.yellow) # creates a method that will trace the ball's path
while True:
rate(1/dt)
# Calculate accelerations of the Lagrangian coordinates:
C12 = C21 = 0.5*M2*L1*L2*cos(theta1-theta2)
Cdet = C11*C22-C12*C21
a = .5*M2*L1*L2*sin(theta1-theta2)
A = -(.5*M1+M2)*g*L1*sin(theta1)-a*theta2dot**2
B = -.5*M2*g*L2*sin(theta2)+a*theta1dot**2
atheta1 = (C22*A-C12*B)/Cdet
atheta2 = (-C21*A+C11*B)/Cdet
# Update velocities of the Lagrangian coordinates:
theta1dot += atheta1*dt
theta2dot += atheta2*dt
# Update Lagrangian coordinates:
dtheta1 = theta1dot*dt # rate of change of theta1
dtheta2 = theta2dot*dt # rate of change of theta2
theta1 += dtheta1 # updates theta1
theta2 += dtheta2 # updates theta2
boardplotter.pos = (1.4+(theta1 % 6.283185306)/5, 1.4+ (theta2 % 6.283185306)/5, 0)
# sets position of drawer with respect to the wrap around screen
if ((theta1 % 6.283185306)/5 < .01 or (theta2 % 6.283185306)/5 < .01):
boardplotter.make_trail=False # stops drawing when drawer hits edge of screen
else : boardplotter.make_trail=True # continues to draw if drawer does not hit edge of screen
ball1.pos = (4+cos(theta1-3.141592653/2),sin(theta1-3.141592653/2), 0) # sets position of ball 1
ball2.pos = (4,1+.3*sin(theta2), .3*cos(theta2)) # sets position of ball 2
ball2.trail.append(pos=ball2.pos) # draws trail for ball 2
ball3.pos = (4+cos(theta1-3.141592653/2),sin(theta1-3.141592653/2)+.3*sin(theta2), .3*cos(theta2)) # sets position of ball 3
ball3.trail.append(pos=ball3.pos) # draws trail for ball 3
frame1.rotate(axis=(0,0,1), angle=dtheta1) # updates position of frame 2
frame2.pos = top+frame1.axis*L1 # updates position of frame 2
frame2.rotate(axis=(0,0,1), angle=dtheta2) # rotates frame 2 to reflect the motion of the pendulums
frame1.trail.append(pos=frame2.pos+vector(cos(theta2-3.141592653/2)*L2,
sin(theta2-3.141592653/2)*L2,0)) # draws trail for end of second pendulum
frame2.trail.append(pos=frame2.pos) # draws trail for end of first pendulum
t = t+dt # updates the time