We are finally at a point in our pendulum investigations where we might tackle a truly interesting case: the so-called driven pendulum. Studying these systems was one of my introductions in seeing science as an artform (actually, this might be a poor representation of the art of science, and more of a display of the gorgeuous happenings in nature that require some digging to find). You will find that one can make beautiful pictures from these systems.
The basic idea is simple. Take the damped pendulum equation, and add a sinusoidal driving force, that periodically exerts a force onto the pendulum. When we do this, we get the following result.
$$
\ddot{\theta} + 2\beta\dot{\theta} + \omega_0^2\sin\theta = C\cos\left(\Omega t\right)
$$
Where C and Ω are constants. Our decision to use a cosine term here is fairly arbitrary. We go about this in the usual way by breaking it down into two coupled first order differential equations. We let
$$
\dot{\theta} = \omega
$$
such that
$$
\dot{\omega} = -\omega_0^2\sin\theta -2\beta\dot{\theta}+C\cos\left( \Omega t\right)
$$
I choose to use the Euler-Cromer method to solve these equations. See the code block below to see how I go about it -- note that the only difference here from my code for the damped case is the inclusion of the driving term.
# Driven pendulum
# Import packages needed
import matplotlib.animation as animation
from numpy import sqrt, sin, cos, arange, pi, append, array
from pylab import plot, xlabel, ylabel, title, show, axhline, savefig, subplots_adjust,\
figure, xlim, rcParams, rc, rc_context, subplot, tight_layout, axvline, legend
# define parameters for our problem
m = 1 # mass of bob at end of pendulum
l = 2.25 # length of pendulum
g = 9.81 # acceleration due to gravity
b = 0.25 # friction coefficient
omega0 = sqrt(g/l) # resonant frequency
beta = b/(2*m)
OMEGA = 8.0 # frequency of driving term
C = 6.0 # strength of driving term
t = 0.0 # time start
t_f = 40.0 # time end
dt = 0.005 # time step
# set up domain
t_points = arange(t,t_f+dt,dt)
# create sets with initial conditions
theta_initial_deg = 0.0
theta_initial = theta_initial_deg*pi/180
omega_initial_deg = 0.0
omega_initial = omega_initial_deg*pi/180
theta_points = [theta_initial]
omega_points = [omega_initial]
# invoke Euler-Cromer
for i in range(len(t_points)):
# find our angular velocity at the next step
omega_new = omega_points[i] + dt*(-2*beta*omega_points[i]\
-omega0**2 * sin(theta_points[i])\
+ C*cos(OMEGA*t_points[i])*cos(theta_points[i]))
# use our newly found angular velocity to find our angle
# at the next step
theta_new = theta_points[i] + dt*omega_new
# update our conditions
theta_points.append(theta_new)
omega_points.append(omega_new)
# convert values to degrees to plot
theta_points_deg = []
omega_points_deg = []
for i in range(len(t_points)+1):
theta_points_deg.append(theta_points[i]*180/pi)
omega_points_deg.append(omega_points[i]*180/pi)
# animate our results
# start with styling options
rcParams.update({'font.size': 18})
rc('axes', linewidth=2)
with rc_context({'axes.edgecolor':'white', 'xtick.color':'white', \
'ytick.color':'white', 'figure.facecolor':'darkslategrey',\
'axes.facecolor':'darkslategrey','axes.labelcolor':'white',\
'axes.titlecolor':'white'}):
# set up a figure that will be ultimately be 2x2
fig = figure(figsize=(10,9))
fig.subplots_adjust(hspace=0.3,wspace=0.5)
# create subplots for each animation
# first the pendulum itself
ax_pend = subplot(2,2,1, aspect='equal',adjustable='datalim')
# get rid of axis ticks
ax_pend.tick_params(axis='both', colors="darkslategrey")
# now our phase diagram
ax_phase = subplot(2,2,2)
axvline(color='k', lw=1)
axhline(color='k', lw=1)
title("Phase Diagram")
xlabel("Angle (deg)")
ylabel("Angular Speed (deg/s)")
# now our angle vs. time graph
# we want it to be on the bottom row, two "plots" long
ax_ang = subplot(2,2,(3,4))
axhline(color='k', lw=1)
title("Angle Vs. Time")
xlabel("Time (s)")
ylabel("Angle (deg)")
### finally we animate ###
# create a list to input images in for each time step
ims = []
index = 0
while index <= len(t_points)+1:
### pendulum ###
# the rigid rod of length l
ln, = ax_pend.plot([0,l*sin(theta_points[index])],\
[0,-l*cos(theta_points[index])], color='white',lw=3)
# our mass on the end
bob, = ax_pend.plot(l*sin(theta_points[index]),\
-l*cos(theta_points[index]),'o',\
markersize=22,color="darkturquoise",zorder=100)
### phase ###
phase_curve, = ax_phase.plot(theta_points_deg[:(index+1)],\
omega_points_deg[:(index+1)],\
color="coral",lw=2.5)
phase_dot, = ax_phase.plot(theta_points_deg[index:(index+1)],\
omega_points_deg[index:(index+1)],\
color="darkturquoise",marker="o",markersize=14)
### angle vs. time ###
ang, = ax_ang.plot(t_points[:(index+1)], theta_points_deg[:(index+1)],\
color="darkturquoise",lw=3)
# add pictures to ims list
ims.append([ln, bob, phase_curve, phase_dot, ang])
# only show every 5 frames
index += 5
# save animations
ani = animation.ArtistAnimation(fig, ims, interval=100)
writervideo = animation.FFMpegWriter(fps=60)
ani.save('../drivenPendulum/4_drivenPend.mp4', writer=writervideo)
The first system is what the code above produces, the other two are produced by playing around with the parameters $\beta$, $C$, and $\Omega$. One can spend quite a bit of time fidgeting with these and seeing how things change. Let's make it a little more interesting by making the driving term have a $\theta$ dependence as follows $$ \ddot{\theta} + 2\beta\dot{\theta} + \omega_0^2\sin\theta = C\cos\theta\cos\left(\Omega t\right) $$ The change to our code is trivial, so I will not show it. I will also not bother with boring you with more words -- the rest of this page is dedicated to showing some of the interesting and beautiful behaviors of these systems. Some are not too pretty (the next three) but are fun to watch regardless, but be sure to watch the phase diagrams until the end. I hope you enjoy.
Now we replace our cosine term with a sine term. $$ \ddot{\theta} + 2\beta\dot{\theta} + \omega_0^2\sin\theta = C\cos\theta\sin\left(\Omega t\right) $$
How wonderful. Next we will investigate the elastic pendulum, which also produces quite lovely pictures.