Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
simple pendulum
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import glob
import os
def simple_pendulum_deriv(x, t, m, g, l):
"""The simple pendulum subject to zero damping and zero control input
Based on material from the MIT OCW subject "Underactuated Robotics"
"lecture 2: the simple pendulum" given by Russell Tedrake, the
maths and some of the related diagrams are given at about 51 minutes
in.
"""
nx = np.zeros(2)
nx[0] = x[1]
nx[1] = -(m * g * l * np.sin(x[0]))
return nx
def plot_orbit( x0 ):
"""Plot the phase space of the pendulum
"""
# Pick ts to be sufficiently large (by inspection) to make sure
# that all of the cases of interest will have reached their
# homoclinic orbit
ts = np.linspace(0.0, 50.0, 1001)
ys = odeint(simple_pendulum_deriv, x0, ts, args = (1.0, 9.8, 1.0))
plt.plot(ys[:,0], ys[:,1])
def pendulum_video( q0, ts, output_file):
ts = np.linspace(0.0, 5.0, 101)
qs = odeint(simple_pendulum_deriv, q0, ts, args = (1.0, 9.8, 1.0))
pendulum_video_sequence(qs, output_file)
def plot_pendulum(q, length = 1.0):
xs = [ 0.0, length * np.sin(q[0]) ]
ys = [ 0.0, -length * np.cos(q[0]) ]
plt.plot( xs, ys, 'b-', linewidth = 3 )
plt.plot( xs[1], ys[1], 'ro' )
plt.gca().set_xlim([-1.2, 1.2])
plt.gca().set_ylim([-1.2, 1.2])
plt.gca().set_aspect('equal')
def encode_video(input_pattern, output_filename):
# taken (almost) directly from the matplotlib website:
# http://matplotlib.sourceforge.net/faq/howto_faq.html#make-a-movie
os.system("mencoder 'mf://%s' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec
def pendulum_video_sequence(qs, output_filename):
for n, q in enumerate(qs):
print n, q
plt.clf()
plot_pendulum( q )
plt.savefig( '_tmp%03d.png' % n )
encode_video('_tmp*.png', output_filename)
for fname in glob.glob('_tmp*.png'):
os.remove(fname)
def plot_interesting_orbits():
plt.clf()
# start with zero velocity at a number of differing angles
# scipy odeint does some pretty weird things if you have enough
# velocity to loop around
plot_orbit(np.array([np.pi * 0.25, 0.0]))
plot_orbit(np.array([np.pi * 0.5, 0.0]))
plot_orbit(np.array([np.pi * 0.9, 0.0]))
plot_orbit(np.array([np.pi * 0.999, 0.0]))
plt.title('homoclinic orbits for an undamped simpled pendulum')
plt.xlabel(r'$\theta$')
plt.ylabel(r'$\dot{\theta}$')
plt.show()
def interesting_videos():
ts = np.linspace(0.0, 50.0, 501)
pendulum_video( np.array([np.pi/2.0, 0.0]), ts, 'undamped_simple_pendulum01.mp4' )
pendulum_video( np.array([3*np.pi/4.0, 0.0]), ts, 'undamped_simple_pendulum02.mp4' )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment