%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('default')
Instructions:¶
Skim through the entire lab before beginning. Read the text carefully and complete the steps as indicated. Submit your completed lab through the D2L dropbox.
Introduction¶
We will study in this lab a strange property of the nonlinear pendulum.
Below, a modified code from Lecture 8 has been reproduced that solve either the linear or non-linear driven, damped pendulum problem.
# constant parameters
#pendulum length in metres
length = 9.8
#acceleration due to gravity
g = 9.8
def Pendulum(N=1500, dt=0.04, F_D=0.0, Ω_D=0.0, q=0.5, θ_0=0.2, nonlinear=True):
"""
Euler Cromer Solution for a damped, driven pendulum
Keyword arguments:
N number of timesteps
dt time step
F_D driving force amplitude
Ω_D driving force frequency
q damping strength
θ_0 initial displacement
Returns
t time
ω angular velocity
θ angular displacement
"""
#initializes omega, a vector of dimension N,to being all zeros
ω = np.zeros(N)
#initializes theta, a vector of dimensionN,to being all zeros
θ = np.zeros(N)
#this initializes the vector time to being all zeros
t = np.zeros(N);
#you need to have some initial displacement, otherwise the pendulum will not swing
θ[0] = θ_0
ω[0] = 0
#loop over the timesteps.
for i in range(N-1):
if nonlinear:
ω[i+1] = ω[i] + (-(g/length)*np.sin(θ[i]) - q*ω[i] + F_D*np.sin(Ω_D*t[i]))*dt
else:
ω[i+1] = ω[i] + (-(g/length)*θ[i] - q*ω[i] + F_D*np.sin(Ω_D*t[i]))*dt
temp_theta = θ[i]+ω[i+1]*dt
# We need to adjust theta after each iteration so as to keep it between +/-π
# The pendulum can now swing right around the pivot, corresponding to theta>2π.
# Theta is an angular variable so values of theta that differ by 2π
# correspond to the same position.
# For plotting purposes it is nice to keep (-π < θ < π).
# So, if theta is <-π, add 2π.If theta is > π, subtract 2π
#********************************************************************************************
if (temp_theta < -np.pi):
temp_theta = temp_theta+2*np.pi
elif (temp_theta > np.pi):
temp_theta = temp_theta-2*np.pi
#********************************************************************************************
# Update theta array
θ[i+1]=temp_theta
t[i+1] = t[i] + dt
return t, ω, θ
Linear behaviour¶
Start by carefully studying the above code. By calling this function, plot the solution of a linear pendulum with the parameters:
F_D = 1.44
Ω_D = 2/3
q = 0.5
dt = 0.01
Choose an appropriate value of N so that at least 100 seconds of the simulation is performed.
fig, axes = plt.subplots(2,1, figsize=(10,8))
t, ω, θ = Pendulum(nonlinear=False) # you'll need to use the correct keyword arguments here
plt.sca(axes[0])
plt.plot(t, θ)
plt.xlabel('t (s)')
plt.ylabel('$\\theta$ (s)')
plt.sca(axes[1])
plt.plot(t, ω)
plt.xlabel('t (s)')
plt.ylabel('$\omega$ (s)')
Recall that for the linear pendulum, the period of the solution should be the same as period of the driving force, $2\pi/\Omega_D$.
Confirm this by measuring the period from your plots, calculating the period, and comparing to the period of the driving force.
your answer here...
Non-linear behaviour¶
Now plot the motion of a non-linear pendulum with the same parameters as given above. What is the period of the non-linear pendulum?
Repeat the above analysis (run the simulation, make the plot, measure the period) with two additional simulations: F_D = 1.35 and F_D = 1.465. When measuring, double check that you have the correct period: you need to find two crests with the same height.
Relative to the the frequency of the driving force, what have you discovered about a non-linear pendulum for $F_D = $ 1.35, 1.44, and 1.465.
your answer here...
Phase plots¶
Another way of visualization the periodic nature of the pendulum is make a phase-plot as we did in Lecture 8. Below, I've already made a phase plot for the linear pendulum. Complete the other three plots to show phase plots for the non-linear pendulums.
fig, axes = plt.subplots(2,2, figsize=(10,10))
plt.sca(axes[0,0])
t, ω, θ = Pendulum(F_D= 1.35, Ω_D=2/3, q =0.5, nonlinear=False)
plt.plot(θ, ω, 'b.')
#Filter the results to exclude initial transient of 5 periods
I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]
θ[I] = float('nan')
ω[I] = float('nan')
# Remove all NaN values from the array to reduce dataset size
θ = θ[~np.isnan(θ)]
ω = ω[~np.isnan(ω)]
plt.plot(θ, ω, 'r.')
plt.xlabel('$\\theta$ (rad)')
plt.ylabel('$\omega$ (rad/s)')
plt.xlim(-np.pi, np.pi)
plt.title('F_D = 1.35 (linear)')
plt.sca(axes[0,1])
t, ω, θ = Pendulum(F_D= 1.35, Ω_D=2/3, q =0.5, nonlinear=True)
plt.plot(θ, ω, 'b.')
#Filter the results to exclude initial transient of 5 periods
I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]
θ[I] = float('nan')
ω[I] = float('nan')
# Remove all NaN values from the array to reduce dataset size
θ = θ[~np.isnan(θ)]
ω = ω[~np.isnan(ω)]
plt.plot(θ, ω, 'r.')
plt.xlabel('$\\theta$ (rad)')
plt.ylabel('$\omega$ (rad/s)')
plt.xlim(-np.pi, np.pi)
plt.title('F_D = 1.35 (linear)')
plt.sca(axes[1,0])
t, ω, θ = Pendulum()
plt.plot(θ, ω, 'b.')
#Filter the results to exclude initial transient of 5 periods
I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]
θ[I] = float('nan')
ω[I] = float('nan')
# Remove all NaN values from the array to reduce dataset size
θ = θ[~np.isnan(θ)]
ω = ω[~np.isnan(ω)]
plt.plot(θ, ω, 'r.')
plt.xlabel('$\\theta$ (rad)')
plt.ylabel('$\omega$ (rad/s)')
plt.xlim(-np.pi, np.pi)
plt.title('F_D = 1.44 (non-linear)')
plt.sca(axes[1,1])
t, ω, θ = Pendulum()
plt.plot(θ, ω, 'b.')
#Filter the results to exclude initial transient of 5 periods
I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]
θ[I] = float('nan')
ω[I] = float('nan')
# Remove all NaN values from the array to reduce dataset size
θ = θ[~np.isnan(θ)]
ω = ω[~np.isnan(ω)]
plt.plot(θ, ω, 'r.')
plt.xlabel('$\\theta$ (rad)')
plt.ylabel('$\omega$ (rad/s)')
plt.xlim(-np.pi, np.pi)
plt.title('F_D = 1.465 (non-linear)')
Better programming style¶
The above four panel plot had a lot of code-duplication. Rewrite that code by using a function to draw each subplot.
Bifurcation Diagram¶
A key observation in the this lab is that as we increase the driving force F_D the period increases. Another way to understand this is to make a plot of all the values of $\theta$ which is in-phase with the driving force. That is, consider only the value of $\theta$ that occurs when $\Omega_D t = 2 n \pi$ where $n$ is an integer. This is analogous to viewing a fast spinning object under a stroboscope.
Just run the code below; no changes are needed. Optional: Try to zoom* into the region 1.47 < F_D < 1.48 by changing the start, stop, and step of the outermost loop.*
fig, axes = plt.subplots()
Omega_D=2/3
for F_Drive_step in np.arange(1,13,0.1):
F_Drive=1.35+F_Drive_step/100;
# Calculate the plot of theta as a function of time for the current drive step
# using the function :- pendulum_function
t, omega, theta = Pendulum(F_D=F_Drive, Ω_D=Omega_D, N=10000);
#Filter the results to exclude initial transient of 10 periods, note
# that the period is 3*pi.
I = [i for i in range(len(t)) if t[i] < 3*np.pi*10]
t[I] = float('nan')
theta[I] = float('nan')
# Further filter the results so that only results in phase with the driving force
# F_Drive are displayed.
# Replace all those values NOT in phase with NaN
I = [i for i in range(len(t)) if abs(np.fmod(t[i],2*np.pi/Omega_D)) > 0.01]
t[I] = float('nan')
theta[I] = float('nan')
# Remove all NaN values from the array to reduce dataset size
t = t[~np.isnan(t)]
theta = theta[~np.isnan(theta)]
F_Drive_Array = np.zeros(len(theta))
F_Drive_Array[:] = F_Drive
#Add a column to the plot of the results
plt.plot(F_Drive_Array,theta,'ks', markersize=1)
# axes scales
plt.xlim(1.35, 1.5);
plt.ylim(1, 3);
plt.xlabel('$F_D$')
plt.ylabel('$\\theta$ (rad)')
plt.title('Bifurcation diagram for the pendulum')
The above diagram begins to show how the period changes as $F_D$ increases. As $F_D$ continues to increase on the right of the graph, there is a transition to fully chaotic behaviour. Here is much more detailed version of the bifurcation diagram (looking at $\omega$ instead of $\theta$) for the non-linear pendulum:
See this page for the details.