{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Damped-Driven Pendulums" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We start by reviewing how to solve the simple pendulum. In equations, the problem is\n", "\n", "\\begin{align}\n", "\\frac{d\\omega}{dt} &= - \\frac{g}{L} \\theta\\;, \\\\\n", "\\frac{d\\theta}{dt} &= \\omega\\;,\n", "\\end{align}\n", "\n", "Let's solve this first with Euler's method." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Euler calculation of motion of simple pendulum\n", "# by Kevin Berwick,\n", "# based on 'Computational Physics' book by N Giordano and H Nakanishi,\n", "# section 3.1\n", "\n", "#pendulum length in metres\n", "length = 1\n", "#acceleration due to gravity\n", "g = 9.8\n", "#Discretize time into 250 intervals\n", "N = 250\n", "#Time step in seconds\n", "dt = 0.04\n", "\n", "#initializes omega, a vector of dimension N,to being all zeros\n", "omega = np.zeros(N)\n", "#initializes theta, a vector of dimensionN,to being all zeros\n", "theta = np.zeros(N)\n", "#this initializes the vector time to being all zeros\n", "t = np.zeros(N);\n", "#you need to have some initial displacement, otherwise the pendulum will not swing\n", "theta[0]=0.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#loop over the timesteps\n", "for i in range(N-1):\n", " omega[i+1] = omega[i] - (g/length)*theta[i]*dt\n", " theta[i+1] = theta[i] + omega[i]*dt\n", " t[i+1] = t[i] + dt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "fig, axes = plt.subplots()\n", "#plots the numerical solution in red\n", "plt.plot(t, theta, 'r-');\n", "\n", "# label axes\n", "plt.xlabel('time (seconds) ');\n", "plt.ylabel('theta (radians)');\n", "plt.title('Simple Pendulum - Euler Method')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Solution using the Euler-Cromer method.\n", "This problem with growing oscillations is addressed by performing the solution using the Euler - Cromer method. The code is below\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Euler_cromer calculation of motion of simple pendulum\n", "# by Kevin Berwick,\n", "# based on 'Computational Physics' book by N Giordano and H Nakanishi,\n", "# section 3.1\n", "\n", "#pendulum length in metres\n", "length = 1\n", "#acceleration due to gravity\n", "g = 9.8\n", "#Discretize time into 250 intervals\n", "N = 250\n", "#Time step in seconds\n", "dt = 0.04" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#initializes omega, a vector of dimension N,to being all zeros\n", "omega = np.zeros(N)\n", "#initializes theta, a vector of dimensionN,to being all zeros\n", "theta = np.zeros(N)\n", "#this initializes the vector time to being all zeros\n", "t = np.zeros(N);\n", "#you need to have some initial displacement, otherwise the pendulum will not swing\n", "theta[0]=0.2" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "#loop over the timesteps\n", "for i in range(N-1):\n", " omega[i+1] = omega[i] - (g/length)*theta[i]*dt\n", " theta[i+1] = theta[i] + omega[i+1]*dt #This line is the only change between this program and standard Euler\n", " t[i+1] = t[i] + dt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "fig, axes = plt.subplots()\n", "\n", "#plots the numerical solution in red\n", "plt.plot(t, theta, 'r-');\n", "# label axes\n", "plt.xlabel('time (seconds) ');\n", "plt.ylabel('theta (radians)');\n", "# axes scales\n", "plt.xlim(0, 10);\n", "plt.ylim(-0.2, 0.2);\n", "# Figure 3. Simple Pendulum: Euler - Cromer method\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Simple Harmonic motion example using a variety of numerical approaches\n", " \n", "In this example I use a variety of approaches in order to solve the following, very simple,\n", "equation of motion. This is the general equation for all simple harmonic oscillators which all have the same mathematical structure.\n", " \n", " $$ \\frac{\\mathrm{d}^2y}{\\mathrm{d}t^2} = -y $$\n", " \n", "Here, four approaches are used to solving the equation, illustrating the use of the Euler, Euler Cromer, Second order Runge-Kutta and finally the python ODE solver integrate.odeint from SciPy. We will need to import this function." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ " \n", " The first step is to take the second order ODE equation and split it into 2 first order ODE equations that can be solved using the mentioned methods. These are\n", " \n", "$$ \\frac{\\mathrm{d}y}{\\mathrm{d}t} = v $$\n", " \n", "and\n", " \n", "$$ \\frac{\\mathrm{d}v}{\\mathrm{d}t} = -y $$\n", " \n", "These can then be solved using the preferred method, by dividing them into individual python functions that return the values of interest, in this case the time and y. Here are the functions to do the individual calculations and the needed code to plot them." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "# Simple harmonic motion - Euler, Euler Cromer and Runge Kutta methods\n", "# by Kevin Berwick,\n", "# based on 'Computational Physics' book by N Giordano and H Nakanishi,\n", "# section 3.1\n", "# Equation is d²y/dt² = -y " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def Euler(initial_displacement):\n", " N = 2500 #Discretize time into 2500 intervals\n", " dt = 0.04 #time step in seconds\n", " \n", " #Initialize v,y and time (t) to all zeros\n", " v = np.zeros(N)\n", " y = np.zeros(N)\n", " t = np.zeros(N)\n", " y[0] = initial_displacement #Add initial displacement\n", " \n", " #Euler Solution\n", " for i in range(N-1):\n", " v[i+1] = v[i] - y[i]*dt\n", " y[i+1] = y[i] + v[i]*dt\n", " t[i+1] = t[i] + dt\n", " \n", " return t,y" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def Euler_Cromer(initial_displacement):\n", " N = 2500 #Discretize time into 2500 intervals\n", " dt = 0.04 #time step in seconds\n", " \n", " #Initialize v,y and time (t) to all zeros\n", " v = np.zeros(N)\n", " y = np.zeros(N)\n", " t = np.zeros(N)\n", " y[0] = initial_displacement #Add initial displacement\n", " \n", " #Euler Cromer Solution\n", " for i in range(N-1):\n", " v[i+1] = v[i] - y[i]*dt\n", " y[i+1] = y[i] + v[i+1]*dt\n", " t[i+1] = t[i] + dt\n", " \n", " return t,y" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "def Runge_Kutta(initial_displacement):\n", " N = 2500 #Discretize time into 2500 intervals\n", " dt = 0.04 #time step in seconds\n", " \n", " #Initialize v,y along with v',y' and time (t) to all zeros\n", " v = np.zeros(N)\n", " vp = 0.0\n", " y = np.zeros(N)\n", " yp = 0.0\n", " t = np.zeros(N)\n", " y[0] = initial_displacement #Add initial displacement\n", " \n", " #Runger Kutta Solution\n", " for i in range(N-1):\n", " vp = v[i]-0.5*y[i]*dt\n", " yp = y[i]+0.5*v[i]*dt\n", " \n", " v[i+1] = v[i] - yp*dt\n", " y[i+1] = y[i] + vp*dt\n", " \n", " t[i+1] = t[i] + dt\n", " \n", " return t,y" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#With the functions written, they need to be called to obtain the results\n", "#Set the initial displacement\n", "initial_displacement = 10\n", "\n", "#Obtain the results for the three first methods\n", "[t_euler,y_euler] = Euler(initial_displacement)\n", "[t_euler_cromer,y_euler_cromer] = Euler_Cromer(initial_displacement)\n", "[t_runge_kutta,y_runge_kutta] = Runge_Kutta(initial_displacement)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#Python Integration method\n", "#To use the the python method, we need to import the SciPy integrate\n", "# module\n", "from scipy import integrate\n", "\n", "#A function is created that has information on the ODEs\n", "def Python_ODEint(variables, t):\n", " y, v = variables #Seperate the two variables\n", " dydt = [v, -y] #Create an array that holds the factors of the ODEs\n", " return dydt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#For the python method, first create an array for the time variable\n", "t_integrate_odeint = np.arange(0, 100, 0.04)\n", "\n", "#Call the integrate.odeint method with the function that was written, \n", "#giving initial conditions and the time array\n", "sol = integrate.odeint(Python_ODEint, [initial_displacement, 0], \n", " t_integrate_odeint)\n", "\n", "#Save the arrays that can be used to plot\n", "[t_integrate_odeint, y_integrate_odeint] = [t_integrate_odeint, sol[:,0]]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false, "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#Plot the four figures next to each other in a 2x2 grid.\n", "fig, plt_array = plt.subplots(2, 2, figsize=(9, 9))\n", "\n", "#Plot the Euler values\n", "plt_array[0, 0].plot(t_euler, y_euler,'r-') \n", "plt_array[0, 0].set_xlabel('Time')\n", "plt_array[0, 0].set_ylabel('Displacement')\n", "plt_array[0, 0].set_ylim(-100, 100) #Set the y range to be higher for Euler\n", "plt_array[0, 0].set_title('Euler')\n", "\n", "#Plot the Euler Cromer values\n", "plt_array[0, 1].plot(t_euler_cromer, y_euler_cromer,'b-') \n", "plt_array[0, 1].set_xlabel('Time')\n", "plt_array[0, 1].set_ylabel('Displacement')\n", "plt_array[0, 1].set_ylim(-15, 15)\n", "plt_array[0, 1].set_title('Euler Cromer')\n", "\n", "#Plot the Runge Kutta values\n", "plt_array[1, 0].plot(t_runge_kutta, y_runge_kutta,'g-') \n", "plt_array[1, 0].set_xlabel('Time')\n", "plt_array[1, 0].set_ylabel('Displacement')\n", "plt_array[1, 0].set_ylim(-15, 15)\n", "plt_array[1, 0].set_title('Runge Kutta')\n", "\n", "#Plot the integrate.odeint values\n", "plt_array[1, 1].plot(t_integrate_odeint,y_integrate_odeint,'k-')\n", "plt_array[1, 1].set_xlabel('Time')\n", "plt_array[1, 1].set_ylabel('Displacement')\n", "plt_array[1, 1].set_ylim(-15, 15)\n", "plt_array[1, 1].set_title('integrate.odeint')\n", "\n", "fig.suptitle('Simple pendulum solution using Euler, Euler Cromer, Runge Kutta and Python ODE solver.')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "fig" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Damped pendulum\n", "\n", "If we add some damping to our model, we can begin to include the effect of friction into the pendulum.\n", "\n", "The friction is assumed to be linearly proportional to the velocity:\n", "\n", "$$ F_{fric} \\propto - \\frac{d \\theta}{dt} $$\n", "\n", "We can model this with a *damping coefficient*, $q$, that represents the amount of friction in the system. The equation of motion then becomes:\n", "\n", "\n", "\\begin{align}\n", "\\frac{d^2\\theta}{dt^2} &= - \\frac{g}{L} \\theta -q \\frac{d \\theta}{dt} \\;, \\\\\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "or written as a pair of first-order equations\n", "\n", "\\begin{align}\n", "\\frac{d\\omega}{dt} &= - \\frac{g}{L} \\theta -q \\frac{d \\theta}{dt} \\;, \\\\\n", "\\frac{d\\theta}{dt} &= \\omega\\;,\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "The analytical/theoretical solution to such differential equations are typically studied in a course on Ordinary Differential Equations. \n", "\n", "\\begin{align}\n", "\\frac{d^2\\theta}{dt^2} &= - \\frac{g}{L} \\theta -q \\frac{d \\theta}{dt} \\;, \\\\\n", "\\end{align}\n", "\n", "Formally, we call this a linear, second-order differential equation with constant coefficients.\n", "\n", "Without assuming knowledge about how to solve such an equation with analytical methods, let's try a numerical approach." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Solution for a damped pendulum using the Euler-Cromer method.\n", "This solution uses q=1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "# Euler_cromer calculation of motion of simple pendulum with damping\n", "# by Kevin Berwick,\n", "# based on 'Computational Physics' book by N Giordano and H Nakanishi,\n", "# section 3.2\n", "\n", "#pendulum length in metres\n", "length = 1\n", "#acceleration due to gravity\n", "g = 9.8\n", "#damping strength\n", "q = 1\n", "#Discretize time into 250 intervals\n", "N = 250\n", "#Time step in seconds\n", "dt = 0.04" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "#initializes omega, a vector of dimension N,to being all zeros\n", "omega = np.zeros(N)\n", "#initializes theta, a vector of dimensionN,to being all zeros\n", "theta = np.zeros(N)\n", "#this initializes the vector time to being all zeros\n", "t = np.zeros(N);\n", "#you need to have some initial displacement, otherwise the pendulum will not swing\n", "theta[0]=0.2\n", "\n", "#loop over the timesteps\n", "for i in range(N-1):\n", " omega[i+1] = omega[i] - (g/length)*theta[i]*dt - q*omega[i]*dt\n", " theta[i+1] = theta[i] + omega[i+1]*dt #This line is the only change between this program and standard Euler\n", " t[i+1] = t[i] + dt" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "fig, axes = plt.subplots()\n", "\n", "#plots the numerical solution in red\n", "plt.plot(t, theta, 'r-')\n", "# label axes\n", "plt.xlabel('time (seconds) ')\n", "plt.ylabel('theta (radians)')\n", "# axes scales\n", "plt.xlim(0, 10)\n", "plt.ylim(-0.2, 0.2)\n", "plt.title('The damped pendulum using the Euler-Cromer method')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Exercise\n", "> 1. Write code using a loop to calculate and plot the motion of damped pendulum for several different values of `q`\n", "> 2. Damped pendulums exhibit three regimes that depend on `q`. What do you think this means?\n", "> * underdamped\n", "> * critically damped\n", "> * overdamped\n", "> 3. Try to estimate the value `q` when the system transitions from underdamped to overdamped; this is what critically damped means." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# in class exercise" ] }, { "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": "markdown", "metadata": {}, "source": [ "### Exercise\n", "> 1. Modify the code below so that it includes a driving force term in the equation of motion." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "# Euler Cromer Solution for damped, driven pendulum\n", "\n", "#pendulum length in metres\n", "length = 9.8\n", "#acceleration due to gravity\n", "g = 9.8\n", "#damping strength\n", "q = 0.5\n", "\n", "#The following is created as a function to be able to easily change values.\n", "def Euler_Cromer_Pendulum(N, dt, F_D, Omega_D):\n", " #initialize array storage\n", " omega = np.zeros(N)\n", " theta = np.zeros(N)\n", " t = np.zeros(N)\n", " \n", " #need to have some initial displacement, otherwise the pendulum will not swing\n", " theta[0]=0.2\n", " omega[0]=0\n", "\n", " #loop over the timesteps.\n", " for i in range(N-1):\n", " omega[i+1] = omega[i] + (-(g/length)*theta[i] - q*omega[i])*dt ## CHANGE THIS LINE\n", " theta[i+1] = theta[i] + omega[i+1]*dt \n", " t[i+1] = t[i] + dt\n", " \n", " return t, omega, theta" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "F_D = 1.2 # Driving force ampltiude\n", "Omega_D = 0.6 # Driving force frequency\n", "N = 15000 # Discretize time\n", "dt = 0.04 # Time step in seconds\n", "\n", "t, omega, theta = Euler_Cromer_Pendulum(N, dt, F_D, Omega_D)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "slideshow": { "slide_type": "subslide" } }, "outputs": [], "source": [ "fig, axes = plt.subplots()\n", "\n", "plt.plot(t, theta, 'b-');\n", "# label axes\n", "plt.xlabel('time (seconds) ');\n", "plt.ylabel('theta (radians)');\n", "# axes scales\n", "plt.xlim(0, 60);\n", "plt.ylim(-2.0, 2.0);\n", "plt.title('Results from Physical pendulum,\\nusing the Euler-Cromer method\\n' + \n", " '$F_D = {:.1f}, \\Omega_D = {:.1f}$'.format(F_D, Omega_D))\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Exercise\n", "> 1. Explore the damped, driven pendulum for different forcing amplitudes and frequencies." ] } ], "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 }