9 Double Pendulum
The most beautiful chaos
9.1 Introduction
The double pendulum is one of the most famous examples of chaos. Enjoy making your own animations!
9.2 The code
comment the lines below if you want interactive mode, i.e., if you want to see the animation in real time.
define equations of motion and other functions
def equations(t, y, args):
""" the equations for the double pendulum """
x1 = y[0] # x1 = theta1, angle
x2 = y[1] # x2 = theta2, angle
p1 = y[2] # p1 = omega1, angular velocity
p2 = y[3] # p2 = omega2, angular velocity
l1,l2,m1,m2,g = args
x1_eq = p1
x2_eq = p2
p1_eq = -((g*(2*m1+m2)*np.sin(x1)+m2*(g*np.sin(x1-2*x2)+2*(l2*p2**2+l1*p1**2*np.cos(x1-x2))*np.sin(x1-x2)))/(2*l1*(m1+m2-m2*(np.cos(x1-x2))**2)))
p2_eq = ((l1*(m1+m2)*p1**2+g*(m1+m2)*np.cos(x1)+l2*m2*p2**2*np.cos(x1-x2))*np.sin(x1-x2))/(l2*(m1+m2-m2*(np.cos(x1-x2))**2))
return [x1_eq, x2_eq, p1_eq, p2_eq]
def calculate_trajectory(args,time,y0):
""" uses scipy's ode itegrator to simulate the equations """
t0,t1,dt = time
r = ode(equations).set_integrator('dopri5')
r.set_initial_value(y0, t0).set_f_params(args)
data=[[t0, y0[0], y0[1], y0[2], y0[3] ]]
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
data.append([r.t, r.y[0], r.y[1], r.y[2], r.y[3] ])
return np.array(data)
def from_angle_to_xy(args,angles):
""" converts angles into xy positions """
l1,l2,m1,m2,g = args
time,theta1,theta2 = angles.T
x1 = l1*np.sin(theta1)
y1 = -l1*np.cos(theta1)
x2 = l2*np.sin(theta2) + x1
y2 = -l2*np.cos(theta2) + y1
return np.array([time,x1,y1,x2,y2]).Tparameters
here the magic happens
Let’s plot stuff, and make a nice movie.
Requrement: ffmpeg
make_movie=True
params = {'backend': 'ps',
'font.size': 20,
'font.family':'serif',
'font.serif':['Computer Modern Roman'], # Times, Palatino, New Century Schoolbook, Bookman, Computer Modern Roman
'ps.usedistiller': 'xpdf',
'text.usetex': True,
}
plt.rcParams.update(params)
plt.ion()
fig = plt.figure(figsize=(9.6,5.4),dpi=100) # 1920x1080
fig.subplots_adjust(left=0, right=1, top=1, bottom=0,hspace=0.02,wspace=0.02)
ax = fig.add_subplot(111)
ax.axis('off') # no frame
def plot_last_seconds(data,index):
""" Plots a line with the trajectory of the tip of pendulum 2 (x2,y2)
"""
how_long = 1.0 # seconds
n = int(how_long/time[2])
to_plot = data[:index,:]
if index < n:
prepend = np.tile(data[0],(n-index,1))
to_plot = np.vstack([prepend,to_plot])
index = n
colormap = plt.cm.Greys_r
colors = [colormap(i) for i in np.linspace(0.0, 1.0, n-1)]
plots = []
for j in np.arange(n-1):
p, = ax.plot(to_plot[index-j-1:index-j+1,3],to_plot[index-j-1:index-j+1,4],
color=colors[j], zorder=-1)
plots.append(p)
return plots
# "plot" returns a tuple of line objects, thus the comma
t,x1,y1,x2,y2 = data_TXY[0]
line1, = ax.plot([0.0,x1], [0.0,y1], 'r-')
line2, = ax.plot([x1,x2], [y1,y2], 'r-')
circ1, = ax.plot([x1], [y1], 'ro',markersize=10)
circ2, = ax.plot([x2], [y2], 'ro',markersize=10)
sizeY = 1.2
ax.axis([-sizeY*16/9,sizeY*16/9,-sizeY,sizeY])
frame_names = []
tex=ax.text(0.0,0.85,'',ha="center")
for i,v in enumerate(data_TXY):
t,x1,y1,x2,y2 = v
# print("t={:.2f}".format(t)) # you might want to know how things are going...
line1.set_data([0.0,x1],[0.0,y1])
line2.set_data([x1,x2],[y1,y2])
circ1.set_data([x1],[y1])
circ2.set_data([x2],[y2])
# plot_last_seconds considerably slows down the simulation,
# but makes it much prettier...
pls = plot_last_seconds(data_TXY,i+1)
tex.set_text(r"$t={:.3f}$ s".format(t))
fig.canvas.draw()
if make_movie:
fname = "_tmp{:05d}.png".format(i)
frame_names.append(fname)
fig.savefig(fname,bbox_inches='tight')
for k in pls:
k.remove()
if make_movie:
frames = "_tmp%5d.png"
frames = "_tmp%5d.png"
movie_command = "ffmpeg -y -r {:} -i {:} double.mp4".format(fps,frames)
os.system(movie_command)
for fname in frame_names:
os.remove(fname)