In [15]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
from scipy.integrate import odeint
from IPython.html.widgets import interact, fixed
The equations of motion for a simple pendulum of mass $m$, length $l$ are:
$$ \frac{d^2\theta}{dt^2} = \frac{-g}{\ell}\sin\theta $$When a damping and periodic driving force are added the resulting system has much richer and interesting dynamics:
$$ \frac{d^2\theta}{dt^2} = \frac{-g}{\ell}\sin\theta - a \omega - b \sin(\omega_0 t) $$In this equation:
When $a=0$ and $b=0$, the energy/mass is conserved:
$$E/m =g\ell(1-\cos(\theta)) + \frac{1}{2}\ell^2\omega^2$$Here are the basic parameters we are going to use for this exercise:
In [16]:
g = 9.81 # m/s^2
l = 0.5 # length of pendulum, in meters
tmax = 50. # seconds
t = np.linspace(0, tmax, int(100*tmax))
Write a function derivs
for usage with scipy.integrate.odeint
that computes the derivatives for the damped, driven harmonic oscillator. The solution vector at each time will be $\vec{y}(t) = (\theta(t),\omega(t))$.
In [17]:
def derivs(y, t, a, b, omega0):
"""Compute the derivatives of the damped, driven pendulum.
Parameters
----------
y : ndarray
The solution vector at the current time t[i]: [theta[i],omega[i]].
t : float
The current time t[i].
a, b, omega0: float
The parameters in the differential equation.
Returns
-------
dy : ndarray
The vector of derviatives at t[i]: [dtheta[i],domega[i]].
"""
theta = y[0]
dtheta = y[1]
dw = -(g/l)*np.sin(theta) - a*dtheta - b*np.sin(omega0*t)
return [dtheta,dw]
In [18]:
assert np.allclose(derivs(np.array([np.pi,1.0]), 0, 1.0, 1.0, 1.0), [1.,-1.])
In [19]:
def energy(y):
"""Compute the energy for the state array y.
The state array y can have two forms:
1. It could be an ndim=1 array of np.array([theta,omega]) at a single time.
2. It could be an ndim=2 array where each row is the [theta,omega] at single
time.
Parameters
----------
y : ndarray, list, tuple
A solution vector
Returns
-------
E/m : float (ndim=1) or ndarray (ndim=2)
The energy per mass.
"""
theta = y[0]
omega = y[1]
if y.ndim == 1:
theta = y[0]
omega = y[1]
EperM = g*l*(1-np.cos(theta))+.5*(l**2)*omega**2
return EperM
if y.ndim == 2:
theta = y[:,0]
omega = y[:,1]
EperM = g*l*(1-np.cos(theta))+.5*(l**2)*omega**2
return EperM
In [20]:
assert np.allclose(energy(np.array([np.pi,0])),g)
assert np.allclose(energy(np.ones((10,2))), np.ones(10)*energy(np.array([1,1])))
Use the above functions to integrate the simple pendulum for the case where it starts at rest pointing vertically upwards. In this case, it should remain at rest with constant energy.
atol
and rtol
arguments of odeint
until $E/m$, $\theta(t)$ and $\omega(t)$ are constant.Anytime you have a differential equation with a a conserved quantity, it is critical to make sure the numerical solutions conserve that quantity as well. This also gives you an opportunity to find other bugs in your code. The default error tolerances (atol
and rtol
) used by odeint
are not sufficiently small for this problem. Start by trying atol=1e-3
, rtol=1e-2
and then decrease each by an order of magnitude until your solutions are stable.
In [21]:
y0 = [np.pi,0]
a = 0
b = 0
omega0 = 0
soln = odeint(derivs,y0,t,args=(a,b,omega0),atol=1e-5, rtol=1e-4)
In [22]:
theta = soln[:,0]
omega = soln[:,1]
In [23]:
plt.plot(t,energy(soln));
plt.title('Energy per Mass vs time');
In [31]:
plt.plot(t,theta,label='$\Theta(t)$');
plt.title('Theta and Omega vs Time');
plt.ylim((-np.pi,2*np.pi));
plt.plot(t,omega,label='$\omega(t)$');
plt.legend();
plt.xlabel('Time');
plt.ylabel('Omega,Theta');
In [25]:
assert True # leave this to grade the two plots and their tuning of atol, rtol.
Write a plot_pendulum
function that integrates the damped, driven pendulum differential equation for a particular set of parameters $[a,b,\omega_0]$.
atol
and rtol
even futher and make sure your solutions have converged.
In [42]:
def plot_pendulum(a=0.0, b=0.0, omega0=0.0):
"""Integrate the damped, driven pendulum and make a phase plot of the solution."""
y0 = [-np.pi + .1,0]
soln = odeint(derivs,y0,t,args=(a,b,omega0),atol=1e-5, rtol=1e-4)
theta = soln[:,0]
omega = soln[:,1]
plt.figure(figsize=(10,6))
plt.plot(theta,omega)
plt.title('Pendlum Motion')
plt.xlabel('Theta')
plt.ylabel('Omega')
Here is an example of the output of your plot_pendulum
function that should show a decaying spiral.
In [43]:
plot_pendulum(0.5, 0.0, 0.0)
Use interact
to explore the plot_pendulum
function with:
a
: a float slider over the interval $[0.0,1.0]$ with steps of $0.1$.b
: a float slider over the interval $[0.0,10.0]$ with steps of $0.1$.omega0
: a float slider over the interval $[0.0,10.0]$ with steps of $0.1$.
In [44]:
interact(plot_pendulum,a=(0.0,1.0,.1),b=(0.0,10.0,.1),omega0=(0.0,10.0,.1));
Use your interactive plot to explore the behavior of the damped, driven pendulum by varying the values of $a$, $b$ and $\omega_0$.
Describe the different classes of behaviors you observe below.
'a' is the damping coefficient, so as we increase 'a' with the other two parameters 0, the pendulum will spiral into the center sooner. Corresponding to the pendulum stopping sooner.
With a fixed at 0.5, if we increase b and $\omega$ together, we are increasing the amplitude of our driving force and changing the intial radial velocity. At some values b and $\omega$ we notice resonance where the pendulum motion and the driving force are in phase.
I worked with Hunter, Jessica, and Brett.