In the previous modules on pendulums, we let our hanging mass be attached to a rigid rod. We will now find the equations of motion for a pendulum whose mass is connected to a spring (with spring constant $k$) instead of a rigid rod. We will be ignoring friction and drag. We will also invoke Largrangian mechanics for the first time. The lagrangian of a system is defined as its kinetic energy $T$ minus its potential energy $U$.
\begin{equation*}
\mathcal{L} = T - U
\end{equation*}
Then, one might use the Euler-Lagrange equation to find the equations of motion.
\begin{equation}\label{eqn:eulerLagrange}
\frac{\partial \mathcal{L}}{\partial q} - \frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{q}}\right)=0
\end{equation}
Where $q$ is a generalized coordinate, and $\dot{q}$ is its derivative with respect to time.
We begin by finding the kinetic energy of such a system.
\begin{equation*}\label{eqn:kineticEnergy}
T = \frac{1}{2}m\dot{r}^2 + \frac{1}{2}mr^2\dot{\theta}^2
\end{equation*}
where $r(t)$ is the length of our spring, as a function of time, and $\theta$ is the angle the pendulum makes with the vertical.
Next, we define the potential energy. Our system has both graviational potential energy, and elastic potential energy.
\begin{equation*}\label{eqn:potentialEnergy}
U = \frac{1}{2}k(l-r)^2 - mgr\cos\theta
\end{equation*}
where $l$ is the equilibrium length of our spring, and $g$ is acceleration due to gravity. Thus, our Lagrangian is
\begin{equation}
\mathcal{L} = \frac{1}{2}m\dot{r}^2 + \frac{1}{2}mr^2\dot{\theta}^2 + mgr\cos\theta - \frac{1}{2}k(l-r)^2
\end{equation}
Now we may start taking derivatives to plug into equation (1). First, we consider $r$:
\begin{equation}\label{eqn:wrtR}
\frac{\partial\mathcal{L}}{\partial r} = mr\dot{\theta}^2 +k(l-r) + mg\cos\theta
\end{equation}
\begin{equation*}
\frac{\partial\mathcal{L}}{\partial\dot{r}}=m\dot{r}
\end{equation*}
\begin{equation}\label{eqn:wrtRdot}
\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{r}}\right)=m\ddot{r}
\end{equation}
Now, we consider $\theta$:
\begin{equation}\label{eqn:wrtTheta}
\frac{\partial\mathcal{L}}{\partial \theta} = -mgr\sin\theta
\end{equation}
\begin{equation*}
\frac{\partial\mathcal{L}}{\partial\dot{\theta}}=mr^2\dot{\theta}
\end{equation*}
\begin{equation}\label{eqn:wrtThetaDot}
\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{\theta}}\right)=2mr\dot{r}\dot{\theta} + mr^2\ddot{\theta}
\end{equation}
Finally, we may plug equations (3) and (4) into equation (1), take equations (5) and (6) and plug those into equation (1), do some rearranging, and we have our equations of motion.
$$
\ddot{r} = r\dot{\theta}^2 + \frac{k}{m}(l-r) +g\cos\theta
$$
$$
\ddot{\theta} = -\frac{g}{r}\sin\theta - \frac{2\dot{r}\dot{\theta}}{r}
$$
To solve these numerically, we let
$$
\dot{r} = v \hspace{1.5cm} \text{and} \hspace{1.5cm} \dot{\theta}=\omega
$$
Such that
$$
\dot{v} = r\omega^2+\frac{k}{m}(l-r) + g\cos\theta
$$
$$
\dot{\omega} = -\frac{g}{r}\sin\theta - \frac{2v\omega}{r}
$$
Just like that, we have turned two second order differential equation into four coupled first order equations. Incredible.
I decided to solve these equations using the fourth-order Runge-Kutta method. See the code block below to see how I animated the pendulums and graphs to come. I will not provide any more analysis of this system, I will just show a few animations starting with different initial conditions. The initial conditions of interest here are the intitial angular displacement, and the intitial stretch/compression of the spring. Playing around with the spring constant can also yield interesting results.
I think coding is gross!
# solve the simple pendulum using fourth order Runge-Kutta method
# 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 a function to calculate our slopes at a given position
# eta = (theta, omega, r, v)
def f(eta):
# what is our angle equal to?
theta = eta[0]
# what is our angular speed equal to?
omega = eta[1]
# what is the length of our spring/rod?
r = eta[2]
# at what speed is the length of our spring changing?
v = eta[3]
# slope of our angle
f_theta = omega
# slope of r
f_r = v
# slope of our angular speed
f_omega = -2*v*omega/r - (g/r)*sin(theta)
# slope of v
f_v = g*cos(theta) + r*omega**2 + (k/m)*(l-r)
# return an array of our slopes
return(array([f_theta,f_omega,f_r,f_v],float))
# set up the parameters of our situation
m = 1.00 # mass of bob
k = 24.0 # spring constant
l = 1.00 # length of pendulum
g = 9.81 # acceleration due to gravity
omega0 = sqrt(g/l) # resonant frequency
theta_initial_deg = 30.0 # initial angle in degrees
omega_initial_deg = 0.0 # initial angular speed
r_initial = l*2.60 # initial stretch of spring/rod (l=equilibrium)
v_initial = 0.0 # intitial tangential speed
theta_initial = theta_initial_deg*pi/180 # convert initial angle into degrees
omega_initial = omega_initial_deg*pi/180 # convert initial anglular speed into degrees
# set up a domain (time interval of interest)
a = 0.0 # interval start
b = 40.0 # interval end
dt = 0.005 # timestep
t_points = arange(a,b,dt) # array of times
# initial conditions eta = (theta, omega, r, v)
eta = array([theta_initial, omega_initial, r_initial, v_initial],float)
# create empty sets to update with values of interest, then invoke Runge-Kutta
theta_points = []
omega_points = []
r_points = []
v_points = []
for t in t_points:
# add current conditions to lists
theta_points.append(eta[0])
omega_points.append(eta[1])
r_points.append(eta[2])
v_points.append(eta[3])
# calculate where we think we are going from slopes of current position
k1 = dt*f(eta)
# calculated where we think we are going from slopes at the midpoint of
# where we are and where we just calculated we would end up
k2 = dt*f(eta + 0.5*k1)
# another mid point calculation
k3 = dt*f(eta + 0.5*k2)
# final time through we use the slope of where we think we are going
k4 = dt*f(eta + k3)
# take a weighted average of all 4 estimations to update our conditions
eta += (k1 + 2*k2 + 2*k3 + k4)/6
# convert values to degrees to plot
theta_points_deg = []
omega_points_deg = []
for i in range(len(t_points)):
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,(r_points[index])*sin(theta_points[index])],\
[0,-(r_points[index])*cos(theta_points[index])],\
color='white',lw=3)
# our mass on the end
bob, = ax_pend.plot((r_points[index])*sin(theta_points[index]),\
-(r_points[index])*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=1.6)
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=2.8)
# add pictures to ims list
ims.append([ln, bob, phase_curve, phase_dot, ang])
# only show every 8 frames
index += 8
# save animations
ani = animation.ArtistAnimation(fig, ims, interval=100)
writervideo = animation.FFMpegWriter(fps=60)
ani.save('../elasticPendulum/30-2_6-elasticPend.mp4', writer=writervideo)