Week 6 Computer Lab

In this lab we will simulate the behavior of a simple pendulum (i.e. with no forces other than gravity and the pendulum arm's tension acting on it). The differential equation describing the motion of a pendulum shown below can be easily derived from Newton's Laws, but has no closed form solution. Therefore a simplifying assumption is often made that the displacement angle of the pendulum (how wide it swings up from vertical) is very small which allows a reasonable approximation to be made by simplifying the equation. With this simplifying assumption the equation becomes that of a simple harmonic oscillator and is easily solvable with a cosine function as the solution.

Where g is acceleration due to gravity, L is the length of the pendulum's arm and θ is the angular displacement from vertical.

But real pendulums have much more interesting behavior. With a computer program it becomes possible to explore this behavior by numerically solving the of the equation with Euler's method and simulating the result.

A pendulum can be started by displacing it and either dropping it or giving it a push to make it go higher. If it is pushed hard enough, it can actually loop around in a complete circle. If it has a stiff rod for an arm, and if it is pushed pushed just the right amount, it could theoretically be made to balance exactly straight over its pivot point. Our ultimate goal in this lab is to find within 3 decimal places what starting velocity we need to release our simulated pendulum to get it to balance on top. Chances are you will not be able to find an exact value for that critical velocity, but you can find two numbers, one a little too fast and one a little to slow, that only differ in the fourth decimal place.

You only need to turn in a listing of your final program or as much as you were able to get working. Please turn something in even if you got stuck so we can help you get unstuck. The more you can write about what you were having trouble with or what lame Python error message you couldn't get past the better

1. First we need a program to do an Euler approximation of the differential equation given above. This will look just like most all of the programs we have done for physics and calculus, a loop that steps along the function using the previous point and the slope at that point to compute the the next point on the curve. In this case our function will be the angular displacement theta of our pendulum as a function of time. dθ/dt, the derivative of this function we will call ω, the angular velocity. The rate of change in velocity over a given time interval is the acceleration which is the derivative of velocity or dω/dt which is also the second derivative of θ. This can be written d²θ/dt². So we can see our differential equation gives is the angular acceleration. This should work the same as our projectile motion program but we are using different variable names for angular displacement, angular velocity and angular acceleration and we have a different equation for acceleration than the uniform acceleration problem. Modify SIRPLOT or one of your other existing programs to graph the solution to the pendulum equation assuming the following conditions:
   • Acceleration due to gravity is 9.8m/s² downwards.
   • The length of the pendulum arm is 2 meters.
   • The initial position of the pendulum when you release it is straight down (θ=0).
   • Start out using an initial speed of around 3m/s to release the pendulum.
   We are assuming no friction and a Mass doesn't not appear in the motion equation so you may choose any mass you like or ignore it completely. Simulate for a long enough period of time to get two or three complete swings back and forth. Choose an appropriate δt stepsize to give a fairly smooth curve without making your program too slow. You will probably want to graph the resulting curve of θ versus time so you can verify it looks plausible. You can also make 3D model of the pendulum and see if it seems to look realistic to you (see examples from class).
2. You should now have everything you need to find the necessary starting angular velocity to balance the pendulum, however a helpful way of visualizing the behavior of a slightly complex system such as this one is to make a phase space plot which is a graph that shows all of the possible states a system could be in. In our case the important variables are θ the angular displacement and ω the angular velocity. Make a phase space plot with θ on the x axis and ω on the y axis. You should see the relationship between position and speed in phase space as a circle or ellipse centered on the origin. You might want to reduce the time period you simulate to only go about one complete swing cycle. If you push the pendulum hard enough for it to loop over the top, how will this appear in phase space? How would it look if we had friction damping the pendulum?
3. Now try zeroing in on the critical starting velocity, ω₀ we are searching for by continuously finding values that are closer and closer either just above or just below. You can do this manually by changing the starting value and rerunning the simulation You can have the computer assist you by enclosing the simulation in a loop that quickly runs your simulation testing a range of starting values for ω₀. You will need to continuously decrease the step size of this loop to increase accuracy to 3 decimal places. Again don't expect to find an exact value that makes the pendulum balance. The actual value is likely to be an irrational number or at least on that has more decimal places than the computer can represent, but it can probably get us closer than experimenting with a real pendulum.

Optional more difficult exercises

from visual.graph import *

graphics=False # turn graphing code on or off

armlength=2
g = 9.8  # acceleration due to gravity
t=0

#draw the pendulum graphics
def make_pendulum():
    global pendulum
    top=vector(0,20,0)
    basesize=vector(.2,.1,.2)
    base = box(size=basesize, pos=top+vector(0,basesize.y/2.0,0))
    pendulum = frame(pos=top,axis=(0,-1,0))
    arm = cylinder(frame=pendulum, length=armlength,  radius=0.01)
    bob = sphere(frame=pendulum, pos=(armlength,0,0), color=color.red, radius=.1)

    scene.center = (0,top.y-armlength/2.0,0)
    scene.range = armlength

#this function causes the graphics windows to appear
def enable_graphics():
    global graphics, thetagraph, frame_rate
    graphics=True
    thetagraph=gdisplay(title="Theta")
    make_pendulum()
    frame_rate=24

#theta is angle  relative to straight down
def setPendulumAngle(pendulum,theta):
    x=armlength*sin(theta)
    y=-armlength*cos(theta)
    pendulum.axis=(x,y,0)

#Run a complete simulation and return the final theta 
def simulate(theta,w,time_limit=4,dt=1.0/128):
    t=0
    theta_max=0 # maximum value of theta reached so far
    if graphics: thetacurve=gcurve(gdisplay=thetagraph)
    while t<=time_limit:

      #Euler's method
      a   = (-g/armlength)*sin(theta)
      w  += a*dt 
      theta += w*dt
      t += dt

      if theta>theta_max: theta_max=theta
      if graphics:
          if t % (1.0/frame_rate) < dt:
              setPendulumAngle(pendulum,theta)
              thetacurve.plot(pos=(t,theta))
              rate(frame_rate)

    return theta_max


dt=1.0/256
upper=5
lower = w_estimate = 3 # starting guess
count=0
last_estimate=0

while last_estimate != w_estimate:
    last_estimate = w_estimate
    theta_max = simulate(0,w_estimate,time_limit=30,dt=dt)

    #do a binary search looking to see if pendulum went over the top
    # or swung back
    if theta_max > 2*pi:
        upper = w_estimate
        print "Overshot"
    else:
        lower = w_estimate
        print "Undershot"

    w_estimate = (upper+lower)/2.0
    count += 1

    print "trial:",count,"estimated w:", w_estimate

enable_graphics()
simulate(0,w_estimate,time_limit=30,dt=dt)

import sys
sys.exit(0)