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