Driven, damped, pendulum¶
- Period doubling route to chaos
- Poincare sections
Note: The notation and graphs in this notebook parallel those in Chaotic Dynamics by Baker and Gollub. (There's a copy in the department office.)
For the driven, damped, pendulum, Newton's second law gives
$$ \frac{d^2\theta}{dt^2} = -\omega_0^2\sin \theta -\gamma \frac{d\theta}{dt} + A\cos(\omega_d t), $$where $\omega_0$ is the natural frequency of the low-amplitude, undamped, and undriven pendulum, and $\omega_d$ is the frequency of the drive. Using dimensionless variables in which time is measured in units of $\omega_0^{-1}$, i.e., $t^\prime = \omega_0 t$, this equation reduces to
$$ \frac{d^2\theta}{d{t^\prime}^2} = -\sin \theta - b \frac{d\theta}{dt^\prime} + g\cos(\omega_d^\prime t^\prime), $$where $b\equiv \gamma/\omega_0 $, $\omega_d^\prime \equiv \omega_d/\omega_0$, and $g \equiv A/\omega_0^2$. In these units the period is $2\pi$. (This is essentially the parametrization of Baker and Gollub.) From this point on I will drop the primes, but retain the understanding that these are nondimensional quantities.
Before doing the numerical integration I will break the second-order equation of motion up into two coupled first-order equations:
\begin{eqnarray*} \frac{d\theta}{dt} &=& \omega \\ \frac{d\omega}{dt} &=& - b\omega -\sin\theta+ g\cos(\omega_d t) \end{eqnarray*}These equations will be entered in the function eqs() below.
Marty Ligare, August 2020
import numpy as np
from scipy import integrate
import matplotlib as mpl
import matplotlib.pyplot as plt
# Following is an Ipython magic command that puts figures in the notebook.
%matplotlib notebook
# Following sets up LateX fonts
#mpl.rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
#mpl.rc('text', usetex=True)
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size
Differential Equations¶
# Function returning derivatives of the dependent quantities u[0]
# and u[1], or more physically in this case, theta and omega.
def eqs(u,t,b,c,om_d):
th = u[0]
om = u[1]
return (om,-np.sin(th) - b*om + c*np.cos(om_d*t))
Example of period doubling¶
See Golub and Baker Fig. 3.4b
Parameters¶
These parameters result in period doubling; see graphs below.
th0 = 0.1 # Initial theta
om0 = 0 # Initial angular velocity
u0 = np.array([th0,om0])# Combine initial conditions in array
b = 0.5 # Damping parameter
g = 1.07 # Driving amplitude
om_d = 2./3 # Driving frequency
Pick times for output, and solve DEs:¶
Note: odeint returns an array:
[[th_0 om_0],
[th_1 om_1],
[th_2 om_2], ... ]
To get single list for th and single list for om we need the
transpose of the returned array. (Could also keep return as a
single array if that's more useful down the road.)
# NOTE: The time points t selected for the output are not the
# points used for the numerical evalution.
t = np.linspace(0,50*2.*np.pi,50*201)
th, om = integrate.odeint(eqs,u0,t,args=(b,g,om_d)).T
plt.figure(1)
plt.plot(t,om)
plt.plot(t,th)
plt.axhline(0)
plt.title("Example of period doubling",fontsize=14)
plt.xlabel("$t$")
plt.ylabel("$\\theta(t)$ and $\omega(t)$"); # No idea why I need \\theta
plt.figure(2)
plt.xlabel("$\\theta$")
plt.ylabel("$\omega$")
plt.title("phase space",fontsize=14)
tplot = 9000
plt.plot(th[tplot:],om[tplot:]);
Example of chaos¶
Same parameters as above, except for drive amplitude $g$.
See Golub and Baker Fig. 3.4c
th0 = 0.1 # Initial theta
om0 = 0 # Initial angular velocity
u0 = np.array([th0,om0])# Combine initial conditions in array
b = 0.5 # Damping parameter
g = 1.15 # Driving amplitude
om_d = 2./3 # Driving frequency
t = np.linspace(0,50*2.*np.pi,50*201) # NOTE: The points selected
# for plotting are not the
# points used for the numerical evalution.
th, om = integrate.odeint(eqs,u0,t,args=(b,g,om_d)).T
plt.figure(3)
plt.plot(t,om)
plt.plot(t,th)
plt.axhline(0)
plt.title("chaotic motion",fontsize=14)
plt.xlabel("$t$")
plt.ylabel("$\\theta(t)$ and $\omega(t)$"); # No idea why I need \\theta
plt.figure(4)
plt.xlabel("$\\theta$")
plt.ylabel("$\omega$")
plt.title("phase space",fontsize=14)
tplot = 1000
plt.plot(th[tplot:],om[tplot:]);
Poincare Section¶
Simply choose points for output of numerical integration (in the 'linspace()' function) that are integer multiples of the drive period:
$$ t_n = n\ \frac{2\pi}{\omega_d}. $$th0 = 0.1 # Initial theta
om0 = 0 # Initial angular velocity
u0 = np.array([th0,om0])# Combine initial conditions in array
b = 0.5 # Damping parameter
g = 1.15 # Driving amplitude
om_d = 2./3 # Driving frequency
t_d = 2*np.pi/om_d
t = np.linspace(0,1000*t_d, 1001) # NOTE: The points selected
# for plotting are not the
# points used for the numerical evalution.
th, om = integrate.odeint(eqs,u0,t,args=(b,g,om_d)).T
plt.figure(5)
plt.xlabel("$\\theta$")
plt.ylabel("$\omega$")
plt.title("Poincare section",fontsize=14)
plt.scatter(th%(2*np.pi),om,s=0.1);
Version Information¶
version_information is from J.R. Johansson (jrjohansson at gmail.com); see Introduction to scientific computing with Python for more information and instructions for package installation.
version_information is installed on the linux network at Bucknell
%load_ext version_information
version_information numpy, scipy, matplotlib