{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('default')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Nonlinear Pendulum\n",
"## Lab 4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Name:\n",
"### Student ID:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Instructions:\n",
"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.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will study in this lab a strange property of the nonlinear pendulum. \n",
"\n",
"Below, a modified code from [Lecture 8](Lecture.08-Forced-Pendulums.ipynb) has been reproduced that solve either the linear or non-linear driven, damped pendulum problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# constant parameters\n",
"\n",
"#pendulum length in metres\n",
"length = 9.8\n",
"\n",
"#acceleration due to gravity\n",
"g = 9.8"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def Pendulum(N=1500, dt=0.04, F_D=0.0, Ω_D=0.0, q=0.5, θ_0=0.2, nonlinear=True):\n",
" \"\"\"\n",
" Euler Cromer Solution for a damped, driven pendulum\n",
" \n",
" Keyword arguments: \n",
" N number of timesteps\n",
" dt time step\n",
" F_D driving force amplitude\n",
" Ω_D driving force frequency\n",
" q damping strength\n",
" θ_0 initial displacement\n",
" \n",
" Returns\n",
" t time\n",
" ω angular velocity\n",
" θ angular displacement\n",
" \"\"\"\n",
" \n",
" #initializes omega, a vector of dimension N,to being all zeros\n",
" ω = np.zeros(N)\n",
" #initializes theta, a vector of dimensionN,to being all zeros\n",
" θ = np.zeros(N)\n",
" #this initializes the vector time to being all zeros\n",
" t = np.zeros(N);\n",
" \n",
" #you need to have some initial displacement, otherwise the pendulum will not swing\n",
" θ[0] = θ_0\n",
" ω[0] = 0\n",
"\n",
" #loop over the timesteps.\n",
" for i in range(N-1):\n",
" \n",
" if nonlinear:\n",
" ω[i+1] = ω[i] + (-(g/length)*np.sin(θ[i]) - q*ω[i] + F_D*np.sin(Ω_D*t[i]))*dt\n",
" else:\n",
" ω[i+1] = ω[i] + (-(g/length)*θ[i] - q*ω[i] + F_D*np.sin(Ω_D*t[i]))*dt\n",
" \n",
" temp_theta = θ[i]+ω[i+1]*dt\n",
" \n",
" # We need to adjust theta after each iteration so as to keep it between +/-π\n",
" # The pendulum can now swing right around the pivot, corresponding to theta>2π.\n",
" # Theta is an angular variable so values of theta that differ by 2π\n",
" # correspond to the same position.\n",
" # For plotting purposes it is nice to keep (-π < θ < π).\n",
" # So, if theta is <-π, add 2π.If theta is > π, subtract 2π\n",
" #********************************************************************************************\n",
" if (temp_theta < -np.pi):\n",
" temp_theta = temp_theta+2*np.pi \n",
" elif (temp_theta > np.pi):\n",
" temp_theta = temp_theta-2*np.pi\n",
" #********************************************************************************************\n",
" \n",
" # Update theta array\n",
" θ[i+1]=temp_theta\n",
" \n",
" t[i+1] = t[i] + dt\n",
" \n",
" return t, ω, θ"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linear behaviour"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Start by carefully studying the above code. By calling this function, plot the solution of a *linear* pendulum with the parameters:\n",
"\n",
" F_D = 1.44\n",
" Ω_D = 2/3\n",
" q = 0.5\n",
" dt = 0.01\n",
"\n",
"Choose an appropriate value of `N` so that at least 100 seconds of the simulation is performed."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, axes = plt.subplots(2,1, figsize=(10,8))\n",
"\n",
"t, ω, θ = Pendulum(nonlinear=False) # you'll need to use the correct keyword arguments here\n",
"\n",
"plt.sca(axes[0])\n",
"plt.plot(t, θ)\n",
"plt.xlabel('t (s)')\n",
"plt.ylabel('$\\\\theta$ (s)')\n",
"\n",
"plt.sca(axes[1])\n",
"plt.plot(t, ω)\n",
"plt.xlabel('t (s)')\n",
"plt.ylabel('$\\omega$ (s)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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$.\n",
"\n",
"Confirm this by measuring the period from your plots, calculating the period, and comparing to the period of the driving force."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*your answer here...*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Non-linear behaviour\n",
"\n",
"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?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"*your answer here...*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Phase plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another way of visualization the periodic nature of the pendulum is make a *phase-plot* as we did in [Lecture 8](lec8.ipynb). 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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(2,2, figsize=(10,10))\n",
"\n",
"plt.sca(axes[0,0])\n",
"t, ω, θ = Pendulum(F_D= 1.35, Ω_D=2/3, q =0.5, nonlinear=False)\n",
"plt.plot(θ, ω, 'b.')\n",
"#Filter the results to exclude initial transient of 5 periods\n",
"I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]\n",
"θ[I] = float('nan')\n",
"ω[I] = float('nan')\n",
"# Remove all NaN values from the array to reduce dataset size \n",
"θ = θ[~np.isnan(θ)]\n",
"ω = ω[~np.isnan(ω)]\n",
"plt.plot(θ, ω, 'r.')\n",
"plt.xlabel('$\\\\theta$ (rad)')\n",
"plt.ylabel('$\\omega$ (rad/s)')\n",
"plt.xlim(-np.pi, np.pi)\n",
"plt.title('F_D = 1.35 (linear)')\n",
"\n",
"plt.sca(axes[0,1])\n",
"t, ω, θ = Pendulum(F_D= 1.35, Ω_D=2/3, q =0.5, nonlinear=True)\n",
"plt.plot(θ, ω, 'b.')\n",
"#Filter the results to exclude initial transient of 5 periods\n",
"I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]\n",
"θ[I] = float('nan')\n",
"ω[I] = float('nan')\n",
"# Remove all NaN values from the array to reduce dataset size \n",
"θ = θ[~np.isnan(θ)]\n",
"ω = ω[~np.isnan(ω)]\n",
"plt.plot(θ, ω, 'r.')\n",
"plt.xlabel('$\\\\theta$ (rad)')\n",
"plt.ylabel('$\\omega$ (rad/s)')\n",
"plt.xlim(-np.pi, np.pi)\n",
"plt.title('F_D = 1.35 (linear)')\n",
"\n",
"plt.sca(axes[1,0])\n",
"t, ω, θ = Pendulum()\n",
"plt.plot(θ, ω, 'b.')\n",
"#Filter the results to exclude initial transient of 5 periods\n",
"I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]\n",
"θ[I] = float('nan')\n",
"ω[I] = float('nan')\n",
"# Remove all NaN values from the array to reduce dataset size \n",
"θ = θ[~np.isnan(θ)]\n",
"ω = ω[~np.isnan(ω)]\n",
"plt.plot(θ, ω, 'r.')\n",
"plt.xlabel('$\\\\theta$ (rad)')\n",
"plt.ylabel('$\\omega$ (rad/s)')\n",
"plt.xlim(-np.pi, np.pi)\n",
"plt.title('F_D = 1.44 (non-linear)')\n",
"\n",
"plt.sca(axes[1,1])\n",
"t, ω, θ = Pendulum()\n",
"plt.plot(θ, ω, 'b.')\n",
"#Filter the results to exclude initial transient of 5 periods\n",
"I = [i for i in range(len(t)) if t[i] < 3*np.pi*5]\n",
"θ[I] = float('nan')\n",
"ω[I] = float('nan')\n",
"# Remove all NaN values from the array to reduce dataset size \n",
"θ = θ[~np.isnan(θ)]\n",
"ω = ω[~np.isnan(ω)]\n",
"plt.plot(θ, ω, 'r.')\n",
"plt.xlabel('$\\\\theta$ (rad)')\n",
"plt.ylabel('$\\omega$ (rad/s)')\n",
"plt.xlim(-np.pi, np.pi)\n",
"plt.title('F_D = 1.465 (non-linear)')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Better programming style\n",
"\n",
"The above four panel plot had a lot of code-duplication. Rewrite that code by using a function to draw each subplot."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Bifurcation Diagram\n",
"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.\n",
"\n",
"*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.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, axes = plt.subplots()\n",
"\n",
"Omega_D=2/3\n",
"for F_Drive_step in np.arange(1,13,0.1):\n",
" F_Drive=1.35+F_Drive_step/100;\n",
" # Calculate the plot of theta as a function of time for the current drive step\n",
" # using the function :- pendulum_function\n",
" t, omega, theta = Pendulum(F_D=F_Drive, Ω_D=Omega_D, N=10000);\n",
" \n",
" #Filter the results to exclude initial transient of 10 periods, note\n",
" # that the period is 3*pi. \n",
" I = [i for i in range(len(t)) if t[i] < 3*np.pi*10]\n",
" t[I] = float('nan')\n",
" theta[I] = float('nan')\n",
" \n",
" # Further filter the results so that only results in phase with the driving force\n",
" # F_Drive are displayed.\n",
" # Replace all those values NOT in phase with NaN\n",
" I = [i for i in range(len(t)) if abs(np.fmod(t[i],2*np.pi/Omega_D)) > 0.01]\n",
" t[I] = float('nan')\n",
" theta[I] = float('nan')\n",
" \n",
" # Remove all NaN values from the array to reduce dataset size \n",
" t = t[~np.isnan(t)]\n",
" theta = theta[~np.isnan(theta)]\n",
" \n",
" F_Drive_Array = np.zeros(len(theta))\n",
" F_Drive_Array[:] = F_Drive\n",
"\n",
" #Add a column to the plot of the results\n",
" plt.plot(F_Drive_Array,theta,'ks', markersize=1)\n",
"\n",
"# axes scales\n",
"plt.xlim(1.35, 1.5);\n",
"plt.ylim(1, 3);\n",
"plt.xlabel('$F_D$')\n",
"plt.ylabel('$\\\\theta$ (rad)')\n",
"\n",
"plt.title('Bifurcation diagram for the pendulum')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:\n",
"\n",
"\n",
"\n",
"See [this page](http://www.thphys.uni-heidelberg.de/~gasenzer/index.php?n1=teaching&n2=chaos) for the details."
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 1
}