Lab.04-Non-linear-Pendulum.ipynb

{
 "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",
    "![](http://www.thphys.uni-heidelberg.de/~gasenzer/pendelplots/ppbif_q=2_g=0.99-1.51.jpeg)\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
}