{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Forced Pendulums\n", "## Lecture 8" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution for a damped, driven pendulum\n", "\n", "In the damped pendulum, the motion decay from an initial perturbation until it comes to rest.\n", "\n", "To keep the pendulum swinging, we need to provide a source of energy to drive the motion.\n", "\n", "One way of modeling such a system is to drive it with a driving amplitude $F_D$ and driving angular frequency $\\Omega_D$\n", "\n", "\\begin{align}\n", "\\frac{d^2\\theta}{dt^2} &= - \\frac{g}{L} \\theta -q \\frac{d \\theta}{dt} + F_D \\sin(\\Omega_D t)\\;, \\\\\n", "\\end{align}\n", "\n", "This driving force will pump energy into (or out of) the system. The external forcing frequency, $\\Omega_D$, \"competes\" with the natural frequency of the pendulum, $\\Omega = \\sqrt{g/L}$." ] }, { "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):\n", " \"\"\"\n", " Euler Cromer Solution for 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", " ω[i+1] = ω[i] + (-(g/length)*θ[i] - q*ω[i] + F_D*np.sin(Ω_D*t[i]))*dt\n", " θ[i+1] = θ[i] + ω[i+1]*dt \n", " t[i+1] = t[i] + dt\n", " \n", " return t, ω, θ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plots the numerical solution for F_Drive = 0.5" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t, ω, θ = Pendulum(N=15000, F_D=0.5, Ω_D=2/3)\n", "\n", "fig, axes = plt.subplots(figsize=(8,3))\n", "\n", "plt.plot(t, θ, 'b-')\n", "# label axes\n", "plt.xlabel('$t$ (s) ')\n", "plt.ylabel('$\\\\theta$ (rad)')\n", "# axes scales\n", "plt.xlim(0, 60);\n", "plt.ylim(-2.0, 2.0);\n", "\n", "plt.title('Results from Physical pendulum, F_D = 0.5')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plots the numerical solution for F_Drive = 1.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t, ω, θ = Pendulum(N=15000, F_D=1.2, Ω_D=2/3)\n", "\n", "fig, axes = plt.subplots(figsize=(8,3))\n", "\n", "plt.plot(t, θ, 'b-')\n", "# label axes\n", "plt.xlabel('$t$ (s) ')\n", "plt.ylabel('$\\\\theta$ (rad)')\n", "# axes scales\n", "plt.xlim(0, 60)\n", "plt.ylim(-1.5, 1.5)\n", "\n", "plt.title('Results from Physical pendulum, F_D = 1.2')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "> 1. What is the natural frequency of a pendulum?\n", "> 2. Write code using a loop to investigate how the steady-state amplitude of driven, damped pendulum depends on\n", "> the forcing frequency (while keeping $q$ and $F_D$ unchanged, see below!)\n", "> 3. What happens when the forcing frequency is near to the natural frequency?\n", "> 4. Explore the effect of changing the damping parameter, $q$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "fig, axes = plt.subplots(2,1, figsize=(8,6))\n", "\n", "for Ω_D in np.arange(0.2, 0.2, 0.1): ## <-- FIX THIS\n", " dt = 0.04\n", " tmax = 100\n", " N = int(tmax // dt)\n", " t, ω, θ = Pendulum(N=N, dt=dt, F_D=1.2, \n", " Ω_D=Ω_D, \n", " q=0.5, \n", " θ_0=0.2)\n", " \n", " # plot the solution on the upper axes\n", " axes[0].plot(t, θ,'-')\n", " \n", " ## EXERCISE: \n", " # Consider the long term (n>N//2) behaviour only\n", " max_θ = 0 ## REPLACE THIS with an expression involving θ\n", " \n", " # plot the maximum value on the lower axes\n", " axes[1].plot(Ω_D, max_θ, 'ko')\n", " \n", "# label axes\n", "plt.sca(axes[0])\n", "plt.xlabel('$t $(s) ')\n", "plt.ylabel('$\\\\theta$ (rad)')\n", "plt.ylim(-4, 4)\n", "\n", "plt.sca(axes[1])\n", "plt.xlabel('$\\Omega_D$ (s$^{-1}$) ')\n", "plt.ylabel('Amplitude (maximum angle)')\n", "plt.axvline(1, color='k', linestyle=':')\n", "plt.ylim(0, 4)\n", "\n", "plt.suptitle('Amplitude for different values of $\\Omega_D$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Discussion: Resonance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Non-linear Pendulums" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The final improvement to model a more realistic pendulum is to no longer assume the amplitude is *small*. Recall back to the model we made when deriving the equation of motion of a simple harmonic oscillator for the pendulum. We assumed\n", "\n", "$$ \\sin(\\theta) \\approx \\theta $$\n", "\n", "when $\\theta << 1$. The assumption of small oscillations is what leads to the *simple* harmonic oscillator model.\n", "\n", "For forced pendulums, we observed that the oscillations can exceed small values. Indeed, for large enough values of the driving force, especially near resonance, the values of $\\theta$ are much larger than $\\pi$! That is the pendulum is actually swinging completely around!\n", "\n", "More realistically we should consider the equation of motion:\n", "\n", "\\begin{align}\n", "\\frac{d^2\\theta}{dt^2} &= - \\frac{g}{L} \\sin(\\theta) \\;, \\\\\n", "\\end{align}\n", "\n", "\n", "And if we include both damping and forcing we have the following equation of model:\n", "\\begin{align}\n", "\\frac{d^2\\theta}{dt^2} &= - \\frac{g}{L} \\sin(\\theta) -q \\frac{d \\theta}{dt} + F_D \\sin(\\Omega_D t)\\;, \\\\\n", "\\end{align}\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solution for a non-linear, damped, driven pendulum:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#pendulum length in metres\n", "length = 9.8\n", "#acceleration due to gravity\n", "g = 9.8\n", "\n", "def Nonlinear_Pendulum(N=1500, dt=0.04, F_D=0.0, Ω_D=0.0, q=0.5, θ_0=0.2):\n", " \"\"\"\n", " Euler Cromer Solution for non-linear, 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", " #initialize vectors of dimension N to being all zeros\n", " ω = np.zeros(N)\n", " θ = np.zeros(N)\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", " ω[i+1] = ω[i] + (-(g/length)*np.sin(θ[i]) - q*ω[i] + F_D*np.sin(Ω_D*t[i]))*dt\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", " t[i+1] = t[i] + dt\n", " \n", " return t, ω, θ" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical solution for F_D = 0.5" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "t, ω, θ = Nonlinear_Pendulum(F_D=0.5, q=0.5, Ω_D=2/3)\n", "\n", "fig, axes = plt.subplots(figsize=(8,3))\n", "\n", "plt.plot(t, θ, 'b-')\n", "# label axes\n", "plt.xlabel('t (s) ')\n", "plt.ylabel('$\\\\theta$ (rad)')\n", "# axes scales \n", "plt.xlim(0, 60)\n", "plt.ylim(-1, 1)\n", "\n", "plt.title('Results from Physical pendulum, F_D = 0.5')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The pendulum oscillates with a uniform period. Let's estimate this period.\n", "\n", "\n", "Choose two peaks in the plot and estimate the period" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "t1 = 30.0 # CHANGE THIS SO THAT\n", "t2 = 40.0 # WE PICK OUT PEAKS IN THE PLOT\n", "\n", "fig, axes = plt.subplots(figsize=(8,3))\n", "plt.plot(t, θ, 'b-')\n", "plt.xlabel('t (s) ')\n", "plt.ylabel('$\\\\theta$ (rad)') \n", "plt.xlim(0, 60)\n", "plt.ylim(-1, 1)\n", "\n", "plt.axvline(t1)\n", "plt.axvline(t2)\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "T = t2 - t1\n", "print(T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the frequency is " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(2*np.pi/T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Observe that this is the driven frequency $\\Omega_D$.\n", "\n", "This is the same as we saw with the *linear* pendulum." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Numerical solution for F_Drive = 1.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t, ω, θ = Nonlinear_Pendulum(N=15000, dt=0.04, F_D=1.2, Ω_D=2/3)\n", "\n", "fig, axes = plt.subplots(figsize=(8,3))\n", "\n", "plt.plot(t, ω, 'b-');\n", "# label axes\n", "plt.xlabel('t (s) ');\n", "plt.ylabel('$\\\\theta$ (rad)');\n", "# axes scales \n", "plt.xlim(0, 60);\n", "#plt.ylim(-4, 4);\n", "\n", "plt.title('Results from Physical pendulum, F_D = 1.2')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the motion for the non-linear pendulum is completely different compared to the linear pendulum. In fact, the motion appears random and unpredictable. However the motion is completely deterministic -- run the code again and you'll get exactly the same plot.\n", "\n", "Here we have an example of a **chaotic** system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "> Explore the behaviour of the non-linear pendulum by adjusting $F_D$ and $\\Omega_D$. \n", "\n", "> Changing the parameters in the above plot from the original values of `F_D=1.2` and `Ω_D=2/3`.\n", "\n", "> What kinds of motion can you discover?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What is chaos?\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "One way of understanding chaotic systems is to study how small changes in the initial conditions lead to potentially large difference in the solution.\n", "\n", "First consider two solutions that differ only slightly in their initial conditions for a non-chaotic solution. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "θ_0 = 0.2 # intitial condition 1\n", "Δθ = 0.001 # initial condition 2: θ_0 + Δθ\n", "\n", "F_D = 0.5 # driving force\n", "\n", "t, ω1, θ1 = Nonlinear_Pendulum(N=15000, dt=0.04, F_D=F_D, Ω_D=2/3, θ_0 = θ_0)\n", "t, ω2, θ2 = Nonlinear_Pendulum(N=15000, dt=0.04, F_D=F_D, Ω_D=2/3, θ_0 = θ_0 + Δθ)\n", "\n", "fig, axes = plt.subplots()\n", "\n", "plt.plot(t, abs(θ2-θ1), 'b-')\n", "\n", "# label axes\n", "plt.xlabel('t (s) ')\n", "plt.ylabel('$\\Delta\\\\theta$ (rad)')\n", "\n", "# axes scales \n", "plt.yscale('log')\n", "plt.xlim(0, 100)\n", "\n", "plt.title('Comparison between two solution, F_D = {:.1f}'.format(F_D))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The error decreases exponentially. This solution is *stable*.\n", "\n", "Now repeat the exercise with a larger driving force that leads to a chaotic solution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "θ_0 = 0.2 # intitial condition 1\n", "Δθ = 0.001 # initial condition 2: θ_0 + Δθ\n", "\n", "F_D = 1.2 # driving force\n", "\n", "t, ω1, θ1 = Nonlinear_Pendulum(N=15000, dt=0.04, F_D=F_D, Ω_D=2/3, θ_0 = θ_0)\n", "t, ω2, θ2 = Nonlinear_Pendulum(N=15000, dt=0.04, F_D=F_D, Ω_D=2/3, θ_0 = θ_0 + Δθ)\n", "\n", "fig, axes = plt.subplots()\n", "\n", "plt.plot(t, abs(θ1-θ2), 'b-')\n", "\n", "# label axes\n", "plt.xlabel('t (s) ')\n", "plt.ylabel('$\\Delta\\\\theta$ (rad)')\n", "\n", "# axes scales \n", "plt.yscale('log')\n", "plt.xlim(0, 120)\n", "\n", "plt.title('Comparison between two solution, F_D = {:.1f}'.format(F_D))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, the difference between the two solution diverges exponentially. The divergent behaviour is key attribute of chaotic system." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Phase plots\n", "\n", "There is another way of visualizing the motion of a pendulum. We can plot the solution in *phase-space*." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t, ω, θ = Nonlinear_Pendulum(N=1500, dt=0.04, F_D=0.5, Ω_D=2/3, θ_0=0.75)\n", "\n", "fig, axes = plt.subplots()\n", "\n", "plt.plot(θ, ω, 'r.', markersize=0.5)\n", "# label axes\n", "plt.xlabel('$\\\\theta$ (rad) ')\n", "plt.ylabel('$\\omega$ (rad/s)')\n", "# axes scales\n", "plt.xlim(-1, 1)\n", "plt.ylim(-0.8, 0.8)\n", "\n", "plt.title('Results from Physical pendulum, F_D = 0.5')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plots the numerical solution for F_Drive = 1.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t, ω, θ = Nonlinear_Pendulum(N=15000, dt=0.04, F_D=1.2, Ω_D=2/3, θ_0=0.2)\n", "\n", "fig, axes = plt.subplots()\n", "\n", "plt.plot(θ, ω, 'r.', markersize=0.5)\n", "# label axes\n", "plt.xlabel('$\\\\theta$ (rad) ')\n", "plt.ylabel('$\\omega$ (rad/s)')\n", "\n", "# axes scales\n", "plt.xlim(-np.pi, np.pi)\n", "plt.ylim(-2.5, 2.5)\n", "\n", "plt.title('Results from Physical pendulum, F_D = 1.2')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want higher resolution, simply increase the resolution by changing `N`." ] } ], "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 }